close all;
clear all;
clc;

%% Define search path.
parentDir = fileparts(fileparts(mfilename('fullpath')));
saveFolder = [parentDir filesep 'data'];
load([saveFolder filesep 'Fig4-Tm1Tm2Tm4-processedTableWithBackground'], 'ProcessedTable');

%%
figDir = [parentDir filesep 'figures' filesep 'Fig4' filesep];
if ~exist(figDir, 'dir'); mkdir(figDir); end

%%
% Stimuli of interest.
stimSetCellStr = getStimSetCellStr(2);
%% Split table into cell types.
tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), ProcessedTable.timeSeriesPath);
cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
cellTypeNum = unique(cellTypeNumArray);
for iNeuron = 1: numel(cellTypeNum)
    StimTableNeuronCell{iNeuron} = ProcessedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
end
%%
stimInd = 1;
nStims = numel(stimSetCellStr);
for iNeuron = 1: numel(cellTypeNum)
    for jStim = 1: nStims
        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
    end
end

for iFly = 1: numel(unique(ProcessedTable.flyInd))
    flyIndsCell = ProcessedTable.flyInd == iFly;
    flyStims = ProcessedTable(flyIndsCell, :).stimParamFileName;
    % Find if stim is within the chosen subset.
    flyStimInds{iFly} = (cellfun(@(x) find(strcmp(x, stimSetCellStr)), flyStims, 'uni', 0));
    % If it contains stimuli, check that it includes the 5 chosen.
    if ~isempty(flyStimInds{iFly})
        setDiff{iFly} = setdiff(1: 5, [flyStimInds{iFly}{:}]);
    else
        setDiff{iFly} = -1;
    end
    % For flies with all chosen stimuli order the table accordingly.
    if isempty(setDiff{iFly})
        sortFlyIndsTemp = find(flyIndsCell);
        [~, sortInd] = sort([flyStimInds{iFly}{:}]);
        sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
    end
end

StimOrderedTable = ProcessedTable(vec([sortFlyInds{:}]), :);

%% Split table into cell types.

tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*', 'tokens'), StimOrderedTable.timeSeriesPath);
cellTypeNumArray = cellfun(@(x) str2num(x{1}{1}), tokens);
cellTypeNum = unique(cellTypeNumArray);
for iNeuron = 1: numel(cellTypeNum)
    StimTableNeuronCell{iNeuron} = StimOrderedTable(cellTypeNumArray == cellTypeNum(iNeuron), :);
end
%%
nStims = numel(stimSetCellStr);
for iNeuron = 1: numel(cellTypeNum)
    for jStim = 1: nStims
        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
    end
end

%% Now can use this table cell array into the other functions.

StimTable = ProcessedTable;
StimSubTablesCell = arrayfun(@(x) getStimSubsetTable(ProcessedTable, ...
                                                     stimSetCellStr, x), ...
                             1: numel(stimSetCellStr), 'uni', 0);  

%% Hand stitched data for Tmneurons and Tm9 saved.
load([saveFolder filesep 'Fig4-TmNeurons-allNeuronsNeuronStimTableOnOff'], 'tmAll');
StimTableNeuronStim = tmAll;
%% Plot tuning curves. Do not discard the on.

alignToOff = true;
alignToAll = false;
respQualityThreshold = 0.5;
rSquareThreshold = 0.2;
padWithNaNs = 0;
qualityMeasure = 'responseQuality';

%%
cellTypeStr = {'Tm1', 'Tm2', 'Tm4', 'Tm9'};
% plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), 0.0, 0.5, figDir, [cellTypeStr{:}]);
plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), 0.2, 0.5, figDir, [cellTypeStr{:}]);
% plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), -20.0, 0, figDir, [cellTypeStr{:}]);
%% Export excel data figure 4b-i
receptiveFields = plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), 0.2, 0.5, figDir, [cellTypeStr{:}]);
%% Write data per cell type
dataFileName = [figDir filesep 'Fig4.xls'];
sheetNamesCell = {'bf', 'cg', 'dh', 'ei'};

stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
               
for iCellType = 1: numel(cellTypeStr)
    startRow = 1;
    sheetName = ['Fig4' sheetNamesCell{iCellType}];
    for iStim = 1: 4
        parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
        DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
        writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
        startRow = startRow + size(DataTable, 1) + 2;
    end
end

%% Export data from Supplementary Figure 4
receptiveFields = plotRFsAlignedWithFitGenotypes(StimTableNeuronStim(:, 2: end), -20, 0, figDir, [cellTypeStr{:}]);
%%
dataFileName = [figDir filesep 'Fig4.xls'];
sheetNamesCell = {'ae', 'bf', 'cg', 'dh'};

stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};
               
for iCellType = 1: numel(cellTypeStr)
    startRow = 1;
    sheetName = ['SuppFig4' sheetNamesCell{iCellType}];
    for iStim = 1: 4
        parameterLabel = [cellTypeStr{iCellType} stimNameStrCell{iStim}];
        DataTable = createRFsTable(receptiveFields{iCellType, iStim}, parameterLabel);
        writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
        startRow = startRow + size(DataTable, 1) + 2;
    end
end

%% Inline from plotRFsAlignedWithFitGenotypes
respIndThresh = respQualityThreshold;
rSqThresh = rSquareThreshold;
clear flyIndsCell
for iNeuron = 1: size(StimTableNeuronStim(:, 2: end), 2)
    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
    [~, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
      flyIndsCell{iNeuron}, ~, ~] = getStimArraysForSameRois(StimTableCell);                                          
    % Align to fit
    widthCell{iNeuron} = cellfun(@(x) arrayfun(@(y) 2 * (2 * sqrt(log(2))) * y.c1, x, 'uni', 0), ...
                        fitParamsCell, 'uni', 1); % The 2 factor is because of bar separation.
    validIndsCell{iNeuron} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
                                     rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
end
%%
validOffInds = cellfun(@(x) x{1} & x{2}, validIndsCell, 'uni', 0);
validOnInds = cellfun(@(x)  x{3} & x{4}, validIndsCell, 'uni', 0);
%%
widthCell = cellfun(@(x) cat(1, x{:}), widthCell, 'uni', 0);
%%
xOFFWidth = cellfun(@(x, y) x(1, y), widthCell, validOffInds, 'uni', 0);
yOFFWidth = cellfun(@(x, y) x(2, y), widthCell, validOffInds, 'uni', 0);
xONWidth = cellfun(@(x, y) x(3, y), widthCell, validOnInds, 'uni', 0);
yONWidth = cellfun(@(x, y) x(4, y), widthCell, validOnInds, 'uni', 0);
allWidthsCell = {xOFFWidth, yOFFWidth, xONWidth, yONWidth};
%%
for iW = 1: numel(allWidthsCell)
    allWidthsCell{iW}(cellfun(@isempty, allWidthsCell{iW}, 'uni', 1)) = {pi};
end
%%
hFig = createPrintFig(15 * [ 1 1/4]);
colors = brewermap(4, 'Set1');
colors = colors([1, 2, 4, 3], :);
for iAx = 1: numel(allWidthsCell)
    hSubAx(iAx) = subplot(1, 4, iAx);
    hWidthAx = beanPlot(allWidthsCell{iAx}, cellTypeStr, colors, ...
                        1, gca, [0-eps 100], 'vertical', 0, 1, 0);
end
linkaxes(hSubAx, 'y');
setFontForThesis(hWidthAx, gcf);
arrayfun(@prettifyAxes, hWidthAx)
arrayfun(@offsetAxes, hWidthAx)
figFileName = ['RF-FWHM-bean'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% Export data from figure 4 j-m
%%
dataFileName = [figDir filesep 'Fig4.xls'];

stimNameStrCell = {'_OFF_bars_horizontal_deg', '_OFF_bars_vertical_deg', ...
                   '_ON_bars_horizontal_deg', '_ON_bars_vertical_deg'};
startRow = 1;
sheetName = ['Fig4j-m'];
for iCellType = 1: numel(cellTypeStr)
    for iStim = 1: 4
        parameterLabel = [cellTypeStr{iCellType} '_FWHM' stimNameStrCell{iStim}];
        xlswrite(dataFileName, {parameterLabel}, sheetName, ['A' num2str(startRow)])
        if allWidthsCell{iStim}{iCellType} == pi % pi was placeholder for no data
            xlswrite(dataFileName, {''}, sheetName, ['B' num2str(startRow)])
        else
            xlswrite(dataFileName, allWidthsCell{iStim}{iCellType}, sheetName, ['B' num2str(startRow)])
        end
        startRow = startRow + 1;
    end
    startRow = startRow + 1;
end
%%
meanONWidth = cellfun(@mean, cellfun(@(x, y) (x + y) / 2, xONWidth, yONWidth, 'uni', 0));