close all;
clear all;
clc;

%% Define search path.
parentDir = fileparts(fileparts(mfilename('fullpath')));
saveFolder = [parentDir filesep 'data'];
load([saveFolder filesep 'Fig2-Tm9-processedTableWithBackgroundOnOff'], 'ProcessedTable');

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

%%
% Stimuli of interest.
stimSetCellStr = getStimSetCellStr(2);
%%
StimTable = ProcessedTable;
StimSubTablesCell = arrayfun(@(x) getStimSubsetTable(ProcessedTable, ...
                                                     stimSetCellStr, x), ...
                             1: numel(stimSetCellStr), 'uni', 0);

%% Plot tuning curves.
StimTableCell = StimSubTablesCell(2: end);
% This extracts the unaligned tuning curves (TCs).
% Retrieve table data in arrays.
[tcArrayCell, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
 tcExtCellArray, tcArgExtCellArray] = cellfun(@getTuningCurveArrayFromTable, ...
                                              StimTableCell, 'uni', 0);
% Align to fit
posCell = cellfun(@(x) arrayfun(@(y) y.b1, x, 'uni', 0), fitParamsCell, 'uni', 1);
widthCell = 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.

lagsCell = cellfun(@(x, y) round(x - size(y, 2) / 2), posCell, tcArrayCell, 'uni', 0);
barRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), tcArrayCell, lagsCell, 'uni', 0);
validIndsCell = cellfun(@(x, y) x(:) > 0.2 & y(:) > 0.5, ...
                        rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
%% Blank epoch was already deleted, now delete vertical bars on edges of screen.
barRFs{2}(:, end-1: end) = [];
barRFs{4}(:, end-1: end) = [];
barRFs{2}(:, 1: 2) = [];
barRFs{4}(:, 1: 2) = [];
%%
offX = barRFs{1}(validIndsCell{1}, :);
offY = barRFs{2}(validIndsCell{2}, :);
onX = barRFs{3}(validIndsCell{3}, :);
onY = barRFs{4}(validIndsCell{4}, :);

% normalize to peak and smooth and normalize again.
maxNorm = @(x, dim) x ./ max(abs(x), [], dim);
getRow = @(x, nRow) x(nRow, :);
normRFs = cellfun(@(x) cell2mat(arrayfun(@(y) smooth(getRow(maxNorm(x,2), y))', 1: size(x, 1), 'uni', 0)'), ...
                  {offX, offY, onX, onY}, 'uni', 0);
normRFs = cellfun(@(x) maxNorm(x,2), normRFs, 'uni', 0);
%%
% Plot lasagna plots, separate figures for off and on, and reapeated panel
% is to have a colorbar that does not distort existing axis of interest.
hFig = createPrintFig(15 * [ 1 2/3]);
hSubAx = plotLasagnaTracesFromCell(gcf, {offX, offY, offY}, brewermap(6, 'Paired'));
colorbar
% arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
arrayfun(@prettifyAxes, hSubAx)
arrayfun(@offsetAxes, hSubAx)
% set(hSubAx, 'Visible', 'off')
figFileName = ['RFcurves-OFF-XY-lasagna'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%%
colorbar
hFig = createPrintFig(15 * [ 1 2/3]);
hSubAx = plotLasagnaTracesFromCell(gcf, {onX, onY, onY}, brewermap(6, 'Paired'));
colorbar
% arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
arrayfun(@prettifyAxes, hSubAx)
arrayfun(@offsetAxes, hSubAx)
% set(hSubAx, 'Visible', 'off')
figFileName = ['RFcurves-ON-XY-lasagna'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% For nat comms editorial policies source data
dataFileName = [figDir filesep 'Fig2.xls'];
DataTable = createRFsTable(offX, 'OFF_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig2b')
DataTable = createRFsTable(offY, 'OFF_vertical_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig2b', 'Range', 'A28:IV100')
DataTable = createRFsTable(onX, 'ON_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig2b', 'Range', 'A56:IV100')
DataTable = createRFsTable(onY, 'ON_vertical_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig2b', 'Range', 'A84:IV100')

%% This produces panel 2

validWidthCell = getCellArrayFromIndexThreshold(widthCell, validIndsCell', 0);
hFig = createPrintFig(15 * [ 1/3 1/3]);
hWidthAx = beanPlot(validWidthCell, {'xOFF' 'yOFF' 'xON' 'yON'}, getONOFFXYColors, ...
                    1, gca, [0-eps 100], 'vertical', 0, 1, 0);
setFontForThesis(hWidthAx, gcf);
arrayfun(@prettifyAxes, hWidthAx)
arrayfun(@offsetAxes, hWidthAx)
figFileName = ['RF-FWHM-bean'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% Write data for fig2c
xlswrite(dataFileName, {'FHWM-OFF-bars-horizontal'}, 'Fig2c', 'A1')
xlswrite(dataFileName, validWidthCell{1}, 'Fig2c', 'B1')
xlswrite(dataFileName, {'FHWM-OFF-bars-vertical'}, 'Fig2c', 'A2')
xlswrite(dataFileName, validWidthCell{2}, 'Fig2c', 'B2')
xlswrite(dataFileName, {'FHWM-ON-bars-horizontal'}, 'Fig2c', 'A3')
xlswrite(dataFileName, validWidthCell{3}, 'Fig2c', 'B3')
xlswrite(dataFileName, {'FHWM-ON-bars-vertical'}, 'Fig2c', 'A4')
xlswrite(dataFileName, validWidthCell{4}, 'Fig2c', 'B4')

%% Range and mean of ON RF fwhm
min([validWidthCell{3}(:); validWidthCell{4}(:)])
max([validWidthCell{3}(:); validWidthCell{4}(:)])
mean([validWidthCell{3}(:); validWidthCell{4}(:)])
std([validWidthCell{3}(:); validWidthCell{4}(:)])

%%
threshold = 0.5;
rSquareThreshold = 0.2;
qualityMeasure = 'responseQuality';
widthCell = plotParamsThresholdedDistribution(StimTableCell, threshold, qualityMeasure, figDir, rSquareThreshold);
%% Do stats for panel 2c
% Not really normal lly distributed
cellfun(@lillietest, widthCell);
% Test offX vs onX and offY vs onY
[pX, meanXDiff, ~] = permutationTest(widthCell{1}, widthCell{3}, 100000);
[pY, meanYDiff, ~] = permutationTest(widthCell{2}, widthCell{4}, 100000);
%% Nested permutations
flyIndsCell = cellfun(@(y) repelem(y.flyInd, cellfun(@numel, y.responseIndex)), ...
                      StimTableCell, 'uni', 0);
for iStim = 1: size(StimTableCell, 2)
    flyIndsCellArray{iStim} = mat2cell(flyIndsCell{iStim}, ...
                cellfun(@numel, StimTableCell{iStim}.responseIndex));
    
    inputCellArray = flyIndsCellArray{iStim};
    StimTable = StimTableCell{iStim};
    
    nonParamTunings = getRoisFeaturePerChannel(cellfun(@(x) x', StimTable.nonParamTuning, 'uni', 0));

    domainCellArray = arrayfun(@(x) x.domain, nonParamTunings, 'uni', 0);
    domainSize = cellfun(@numel, domainCellArray);
    outputCellArray{iStim} = cell2mat(inputCellArray(domainSize == mode(domainSize)));

end
correctedFlyInds = outputCellArray;
%%
goodFlyIndsCell = getCellArrayFromIndexThreshold(correctedFlyInds, validIndsCell', 0);

nFlies = cellfun(@(x) numel(unique(x)), goodFlyIndsCell);
%% Store the fly indices
xlswrite(dataFileName, {'FlyID-OFF-bars-horizontal'}, 'Fig2c', 'A6')
xlswrite(dataFileName, goodFlyIndsCell{1}', 'Fig2c', 'B6')
xlswrite(dataFileName, {'FlyID-OFF-bars-vertical'}, 'Fig2c', 'A7')
xlswrite(dataFileName, goodFlyIndsCell{2}', 'Fig2c', 'B7')
xlswrite(dataFileName, {'FlyID-ON-bars-horizontal'}, 'Fig2c', 'A8')
xlswrite(dataFileName, goodFlyIndsCell{3}', 'Fig2c', 'B8')
xlswrite(dataFileName, {'FlyID-ON-bars-vertical'}, 'Fig2c', 'A9')
xlswrite(dataFileName, goodFlyIndsCell{4}', 'Fig2c', 'B9')

%% Still do the nested permutation
nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
% Delete flies with no rois to analyze.
for iStim = 1: numel(nROIsPerFlyCell)
    nROIsPerFlyCell{iStim}(nROIsPerFlyCell{iStim} == 0) = [];
end% generate continuos numbering.
reorderedFlyIndsCell = cellfun(@(x) repelem(1: numel(x), x), nROIsPerFlyCell, 'uni', 0);
%%
% samples are cells, but we will permute all cells of a fly together.
% Test only within orientation differences
sample1 = 1: numel(nROIsPerFlyCell{1});
sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{3}));

permutations = 10000;
sidedness = 'both';
plotresult = 1;

[obsDiffX, pValX] = nestetPermutationTest(sample1, sample2, ...
                                        widthCell([1, 3]), ...
                                        reorderedFlyIndsCell([1, 3]), ...
                                        permutations, sidedness, plotresult);

sample1 = 1: numel(nROIsPerFlyCell{2});
sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{4}));

[obsDiffY, pValY] = nestetPermutationTest(sample1, sample2, ...
                                        widthCell([2, 4]), ...
                                        reorderedFlyIndsCell([2, 4]), ...
                                        permutations, sidedness, plotresult);
                                    
%% Plot raw traces for supplementary figure.
hFig = createPrintFig(15 * [1 1]);

drawRandomSubset = 1;
for iStim = 1: size(StimTableCell, 2)
    
    inputCellArray = StimTableCell{iStim}.paddedEpochStacks;
    StimTable = StimTableCell{iStim};
    % To remove flies with different number of epochs than the majority.
    nonParamTunings = getRoisFeaturePerChannel(cellfun(@(x) x', StimTable.nonParamTuning, 'uni', 0));
    
    domainCellArray = arrayfun(@(x) x.domain, nonParamTunings, 'uni', 0);
    domainSize = cellfun(@numel, domainCellArray);
    outputCellArray{iStim} = inputCellArray(domainSize == mode(domainSize));

    allRFTraces = cell2mat(cellfun(@(x) cat(1, x{:})', ...
                                   outputCellArray{iStim}, 'uni', 0));

    size(allRFTraces)
    rng('default')
    % imagesc(allRFTraces(validIndsCell{1}, :))
    hTracesAx(iStim) = subplot(2, 2, iStim);
    cla
    plotInds = find(validIndsCell{iStim});
    if drawRandomSubset
        plotInds = plotInds(randperm(10));
    end
    barResponseTraces{iStim} = allRFTraces(plotInds, :);
    timePoints{iStim} = (1: numel(allRFTraces(plotInds(1), :))) / 10; % Data interpolated at 10Hz
    plot(timePoints{iStim}, bsxfun(@plus, allRFTraces(plotInds, :)', ...
                0.9 * (1: numel(plotInds))), 'LineWidth', 1, 'Color', 0.25 * [1 1 1]);
    hold on;
    axis tight
    hLines = arrayfun(@(x) refline(0, x), 0.9 * (1: numel(plotInds)));
    [hLines.Color] = deal([1 1 1] * 0.5);
    [hLines.LineWidth] = deal(0.5);
    timePointsPerEpoch = size(allRFTraces, 2) / mode(domainSize);
    hLines = arrayfun(@(x) line([x, x], ylim), ...
                      timePointsPerEpoch / 10 * (1/3: 1: mode(domainSize)) + 1/3);
    [hLines.Color] = deal([1 1 1] * 0.75);
    [hLines.LineWidth] = deal(0.5);
    [hLines.LineStyle] = deal(':');

    % Send the reference lines to the back.
    chH = get(gca,'Children');
    set(gca,'Children',flipud(chH));
end

xlabel(hTracesAx(3), 'Time (s)')
xlabel(hTracesAx(4), 'Time (s)')
ylabel(hTracesAx(1), '\DeltaF/F_0')
ylabel(hTracesAx(3), '\DeltaF/F_0')

setFontForThesis(hTracesAx, gcf);
arrayfun(@prettifyAxes, hTracesAx)
arrayfun(@offsetAxes, hTracesAx)
figFileName = ['Fig2-suppl-meanTraces-sortedPositions'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% Export data for Supplementary Fig 2
dataFileName = [figDir filesep 'Fig2.xls'];

sheetName = 'SupFig2';
stimNameStrCell = {'Tm9_OFF_horizontal_bar_time_seconds', 'Tm9_OFF_vertical_bar_time_seconds', ...
                   'Tm9_ON_horizontal_bar_time_seconds', 'Tm9_ON_vertical_bar_time_seconds'};
startRow = 1;

for iStim = 1: 4
    DataTable = createRFsTable(barResponseTraces{iStim}, stimNameStrCell{iStim});
    DataTable(:, 1).(stimNameStrCell{iStim}) = timePoints{iStim}(:);
    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
    startRow = startRow + size(DataTable, 1) + 2;
end

%% Now pick a value an compare stimuli. And do some statistics.
plotBootstrapCoeffOfVariationTest(widthCell, figDir, threshold, rSquareThreshold)
%% Variance plot
plotBarVariance(widthCell, figDir, threshold, rSquareThreshold)