% 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 = collectRFTraces(dataEdges, 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};

%% Get the shortest time points of the interpolated traces, to have a common
% reference frame.
barStimInds = cellfun(@(x) find(strcmp(x.stimulusName, barStimFileNames)), ...
                      squeeze(dataEdges), 'uni', 0);
barStimInds(:, 1) = [];                  
invalidFlies = any(cellfun(@isempty, barStimInds), 2);
tStart = cellfun(@(x) x.times(1), squeeze(dataEdges));
tEnd = cellfun(@(x) x.times(end), squeeze(dataEdges));
tStart(:, 1) = [];
tEnd(:, 1) = [];
tStart(invalidFlies, :) = [];
tEnd(invalidFlies, :) = [];
barStimInds(invalidFlies, :) = [];
barStimInds = cell2mat(barStimInds);
for iFly = 1: size(barStimInds, 1)
    [~, sortInds] = sort(barStimInds(iFly, :));
    tStartSorted(iFly, :) = tStart(iFly, sortInds);
    tEndSorted(iFly, :) = tEnd(iFly, sortInds);
end
maxTStart = max(tStartSorted, [], 1);
minTEnd = min(tEndSorted, [], 1);

%%
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}.roiTCsInterp{kChannel}, 2);
                roiTCs = dataEdges{iFly, 1, jStim}.roiTCsInterp{kChannel};
                startInd = find(dataEdges{iFly, 1, jStim}.times == maxTStart(rfInd));
                % For better numerical tolerance.
                endInd = find(abs(dataEdges{iFly, 1, jStim}.times - minTEnd(rfInd)) < 1e-5);
                receptiveFields = squeeze(roiTCs);
                receptiveFields = receptiveFields(1: end - 1 * deleteLastROI, startInd : endInd);
                receptiveFieldsAll{rfInd, kChannel}  = [receptiveFieldsAll{rfInd, kChannel}; receptiveFields];                
            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