close all;
clear all;
clc;

%% Define search path.
parentDir = fileparts(fileparts(mfilename('fullpath')));
saveFolder = [parentDir filesep 'data'];
load([saveFolder filesep 'Fig7-Tm9CDM-processedTableWithoutBackgroundOnOff'], 'ProcessedTable');

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

%%
% Stimuli of interest.
stimSetCellStr = getStimSetCellStr(2);
% Check only bars.
nStims = 5;
stimSetCellStr = stimSetCellStr(1: nStims);
%% 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))
    flyInds = ProcessedTable.flyInd == iFly;
    flyStims = ProcessedTable(flyInds, :).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: nStims, [flyStimInds{iFly}{:}]);
    else
        setDiff{iFly} = -1;
    end
    % For flies with all chosen stimuli order the table accordingly.
    if isempty(setDiff{iFly})
        sortFlyIndsTemp = find(flyInds);
        [~, sortInd] = sort([flyStimInds{iFly}{:}]);
        sortFlyInds{iFly} = sortFlyIndsTemp(sortInd);
    end
end

StimOrderedTable = ProcessedTable(vec(cat(1, sortFlyInds{:})), :);
%%
%% Split table into cell types.

tokens = arrayfun(@(x)  regexp(x, '.*Tm(\d).*plus(.*CDM)', 'tokens'), StimOrderedTable.timeSeriesPath);
tokens = cellfun(@(x) [x{1}{:}], tokens, 'uni', 0);
cellTypeNum = unique(tokens);
for iNeuron = 1: numel(cellTypeNum)
    StimTableNeuronCell{iNeuron} = StimOrderedTable(strcmp(tokens, 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

%%
% Maybe this should get its own file
doFit = 0;
if doFit
    flyFits = fitDifferenceOfGaussiansFromNeuronStimTables(StimTableNeuronStim);
    save([dataSaveFolder 'DoGFitsTm9TestFminCon'], 'flyFits', '-v7.3');
end

%%
load([saveFolder filesep 'Fig7-Tm9CDM-DoGFitsTm9TestFminCon'], 'flyFits');

%% From DoG

if size(flyFits, 2) == 5
    flyFits(:, 1, :) = [];
end
cdm = squeeze(flyFits(1, :, :));
noCdm = squeeze(flyFits(2, :, :));
cdm(:, all(cellfun(@isempty, cdm))) = [];
noCdm(:, all(cellfun(@isempty, noCdm))) = [];

%% Bad chunk of code to try to align data structure of main analysis with the 
% new Diff of Gaussians fit structure.
rSqThresh = 0;
respIndThresh = 0.5;
iInd = 1;
clear roisToRemoveInd
for iNeuron = [2, 1] % first control then cdm
    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
    [tcArrayCell, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
      flyIndsCell, tcExtCellArray, tcArgExtCellArray, roisToRemoveIndTmp] = ...
                                    getStimArraysForSameRois(StimTableCell, 1); 
    roisToRemoveInd{iInd} = roisToRemoveIndTmp;
    validIndsCell{iInd} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
                            rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
    iInd = iInd + 1;
end

[flyDelInd, stimDelInd] = cellfun(@(x) ind2sub(size(x), ...
                                               find(~cellfun(@isempty, x))), ...
                                               roisToRemoveInd, 'uni', 0);
%% for control
genotypeInd = 1;
roiDelInds = roisToRemoveInd{genotypeInd};
delIndsPerFly = {};
for iRoi = 1: numel(flyDelInd{genotypeInd})
    delInds = roiDelInds{flyDelInd{genotypeInd}(iRoi), ...
                         stimDelInd{genotypeInd}(iRoi)}{1};
    delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [];
    delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [delIndsPerFly{flyDelInd{genotypeInd}(iRoi)}, delInds];
end
                               
for iFly = 1: numel(delIndsPerFly)
    if ~isempty(delIndsPerFly{iFly})
        for jStim = 1: size(noCdm, 1)
           tmp = noCdm{jStim, iFly};
           tmp(unique(delIndsPerFly{iFly})) = [];
            noCdm{jStim, iFly} = tmp;
        end
    end
end
%% for cdm
clear delIndsPerFly
genotypeInd = 2;
roiDelInds = roisToRemoveInd{genotypeInd};
delIndsPerFly = {};
for iRoi = 1: numel(flyDelInd{genotypeInd})
    delInds = roiDelInds{flyDelInd{genotypeInd}(iRoi), ...
                         stimDelInd{genotypeInd}(iRoi)}{1};
    delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [];
    delIndsPerFly{flyDelInd{genotypeInd}(iRoi)} = [delIndsPerFly{flyDelInd{genotypeInd}(iRoi)}, delInds];
end
                               
for iFly = 1: numel(delIndsPerFly)
    if ~isempty(delIndsPerFly{iFly})
        for jStim = 1: size(cdm, 1)
           tmp = cdm{jStim, iFly};
           tmp(unique(delIndsPerFly{iFly})) = [];
            cdm{jStim, iFly} = tmp;
        end
    end
end

%%
% reviewer #1 analysis by fly

% this method is not aligned with the new fit structure, but the structure
% is separated by flies already so we can recreate the indices.
% the minux is to exclude the background roi.
flyIndsCdm = repelem(1: size(cdm, 2), cellfun(@numel, cdm(1, :)) - 1);
flyIndsNoCdm = repelem(1: size(noCdm, 2), cellfun(@numel, noCdm(1, :)) - 1);
              
%%
[Amp, Pos, Sigma, RSq, baseline] = deal(cell(1, 2));

[Amp{1}, Pos{1}, Sigma{1}, RSq{1}, baseline{1}] = parseFitParams(noCdm);
[Amp{2}, Pos{2}, Sigma{2}, RSq{2}, baseline{2}] = parseFitParams(cdm);

centerSurround = cellfun(@(x, y) x ./ y, Amp{1}.center, Amp{1}.surround, 'uni', 0);
%%
allAmpCenter   = cellfun(@(x) cell2mat(x.center), Amp, 'uni', 0);
allAmpSurround   = cellfun(@(x) cell2mat(x.surround), Amp, 'uni', 0);
allAmpCenterSurround = cellfun(@(x, y) x ./ y, allAmpCenter, allAmpSurround, 'uni', 0);
if any(~cellfun(@isnan, baseline))
    allAmpBaseline = cellfun(@(x) cell2mat(x), baseline, 'uni', 0);
end

allRSq = cellfun(@cell2mat, RSq, 'uni', 0);
%%
minOffRSq = cellfun(@(x) min(x(1:2, :)), allRSq, 'uni', 0);
minOnRSq  = cellfun(@(x) min(x(3:4, :)), allRSq, 'uni', 0);
minAllRSq = cellfun(@min, allRSq, 'uni', 0);
%%
respIndThresh = 0.5;
% thresholdRSq = 0.2;
thresholdRSq = -20;
goodCells = cellfun(@(x) x > thresholdRSq, minAllRSq, 'uni', 0);

%% This is to plot matching RFs to dog fit
% plotRFsAlignedWithFitGenotypes(StimTableNeuronStim([2, 1], 2: end), 0, 0.5, figDir, [cellTypeStr{[2, 1]}]);

[offX, offY, onX, onY] = deal(cell(size(StimTableNeuronStim, 1), 1)); 
validIndsCellTmp = goodCells([2, 1]); % flip from (noCDM, CDM) to conserve original order in table (CDM, noCDM).
for iNeuron = [2, 1] % 2: noCDM, 1: CDM, thus offX = [cdm, noCDM]
    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
    % Overwrite the gauss1 rsquare with gauss2 rsquare.
    [tcArrayCell, responseIndArrayCell, ~, fitParamsCell, ...
      flyIndsCell, tcExtCellArray, tcArgExtCellArray] = getStimArraysForSameRois(StimTableCell);                                          
    % Align to fit
    posCell = cellfun(@(x) arrayfun(@(y) y.b1, x, 'uni', 0), fitParamsCell, 'uni', 1);
    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);

    %% 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) = [];
    %%
    % Include the response quality index.
    validRespInds = all(cat(2, responseIndArrayCell{:}) > respIndThresh, 2)';
    offX{iNeuron} = barRFs{1}(validIndsCellTmp{iNeuron} & validRespInds, :);
    offY{iNeuron} = barRFs{2}(validIndsCellTmp{iNeuron} & validRespInds, :);
    onX{iNeuron} = barRFs{3}(validIndsCellTmp{iNeuron} & validRespInds, :);
    onY{iNeuron} = barRFs{4}(validIndsCellTmp{iNeuron} & validRespInds, :);
    validRespIndsCell{iNeuron} = validRespInds; 
end

% Include the response quality index
goodCells = cellfun(@(x, y) x & y, goodCells, validRespIndsCell([2, 1]), 'uni', 0);
%%
[offX{cellfun(@isempty, offX)}] = deal(zeros(1, size(offX{find(cellfun(@isempty, offX)==0, 1)}, 2)));
[offY{cellfun(@isempty, offY)}] = deal(zeros(1, size(offY{find(cellfun(@isempty, offY)==0, 1)}, 2)));
[onX{cellfun(@isempty, onX)}] = deal(zeros(1, size(onX{find(cellfun(@isempty, onX)==0, 1)}, 2)));
[onY{cellfun(@isempty, onY)}] = deal(zeros(1, size(onY{find(cellfun(@isempty, onY)==0, 1)}, 2)));
%% Plot selected rf traces for figure 7a-b
cellTypeStr = {'CDM', 'noCDM'};
plotLasagnaRFs(offX([2, 1]), offY([2, 1]), 'OFF', figDir,  [cellTypeStr{[2, 1]}]);
plotLasagnaRFs(onX([2, 1]), onY([2, 1]), 'ON', figDir,  [cellTypeStr{[2, 1]}]);
%% Export data from figure 7a-b

dataFileName = [figDir filesep 'Fig7.xls'];

stimNameStrCell = {'_OFF_horizontal_position_deg', '_OFF_vertical_position_deg', ...
                   '_ON_horizontal_position_deg', '_ON_vertical_position_deg'};

cellTypeStr = {'CDM', 'Control'};
receptiveFields = [offX(1) offY(1) onX(1) onY(1); offX(2) offY(2) onX(2) onY(2)];     
startRow = 1;
sheetName = ['Fig7ab'];

for iCellType = 1: numel(cellTypeStr)
    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

%% For reviewer #1 comment 19. Plotting traces with different normalizations.
% this normalizes each tuning curve, so width differences are highlighted
% regardless of amplitude differences between conditions.

% normalize to peak and smooth and normalize again.
maxAbsNorm = @(x, dim) x ./ max(abs(x), [], dim);
maxNorm = @(x, dim) x ./ max(x, [], dim);
minNorm = @(x, dim) x ./ -min(x, [], dim);
rangeNorm = @(x, dim) x ./ range(x, dim);

maxNormRFs = cellfun(@(x) maxNorm(x, 2), [offX, offY], 'uni', 0);
minNormRFs = cellfun(@(x) minNorm(x, 2), [onX, onY], 'uni', 0);
[normOffX, normOffY, normOnX, normOnY] = deal(maxNormRFs([2, 1], 1), maxNormRFs([2, 1], 2), ...
                                              minNormRFs([2, 1], 1), minNormRFs([2, 1], 2));
plotLasagnaRFs(normOffX, normOffY, 'OFF', figDir,  [cellTypeStr{[2, 1]} '-maxMinNorm']);
plotLasagnaRFs(normOnX, normOnY, 'ON', figDir,  [cellTypeStr{[2, 1]} '-maxMinNorm']);


rangeNormRFs = cellfun(@(x) rangeNorm(x, 2), [offX, offY, onX, onY], 'uni', 0);
[normOffX, normOffY, normOnX, normOnY] = deal(rangeNormRFs([2, 1], 1), rangeNormRFs([2, 1], 2), ...
                                              rangeNormRFs([2, 1], 3), rangeNormRFs([2, 1], 4));

plotLasagnaRFs(normOffX, normOffY, 'OFF', figDir,  [cellTypeStr{[2, 1]} '-rangeNorm']);
plotLasagnaRFs(normOnX, normOnY, 'ON', figDir,  [cellTypeStr{[2, 1]} '-rangeNorm']);
%% Actually quantification is averaging over orientations, so let's plot those
normOff = cellfun(@(x, y) (x + y) / 2, normOffX, normOffY, 'uni', 0);
normOn = cellfun(@(x, y) (x + y) / 2, normOnX, normOnY, 'uni', 0);

plotLasagnaRFs(normOff, normOn, 'ONOFF', figDir,  [cellTypeStr{[2, 1]} '-rangeNorm-meanXY']);
%%
offAmpCenter = cellfun(@(x, y) mean(x(1:2, y)), allAmpCenter, goodCells, 'uni', 0);
offAmpSurround = cellfun(@(x, y) mean(x(1:2, y)), allAmpSurround, goodCells, 'uni', 0);
onAmpCenter = cellfun(@(x, y) mean(x(3:4, y)), allAmpCenter, goodCells, 'uni', 0);
onAmpSurround = cellfun(@(x, y) mean(x(3:4, y)), allAmpSurround, goodCells, 'uni', 0);
%%
sigma2fwhmFactor = 2 * 2 * sqrt(log(2));
allWidthCenter   = cellfun(@(x) cell2mat(x.center) * sigma2fwhmFactor, Sigma, 'uni', 0);
allWidthSurround   = cellfun(@(x) cell2mat(x.surround) * sigma2fwhmFactor, Sigma, 'uni', 0);
% allWidthCenterSurround = cellfun(@(x, y) x ./ y, allWidthCenter, allWidthSurround, 'uni', 0);
%%
offWidthCenter = cellfun(@(x, y) mean(x(1:2, y)), allWidthCenter, goodCells, 'uni', 0);
offWidthSurround = cellfun(@(x, y) mean(x(1:2, y)), allWidthSurround, goodCells, 'uni', 0);
onWidthCenter = cellfun(@(x, y) mean(x(3:4, y)), allWidthCenter, goodCells, 'uni', 0);
onWidthSurround = cellfun(@(x, y) mean(x(3:4, y)), allWidthSurround, goodCells, 'uni', 0);

%% Start plotting
flyIndsCell = {flyIndsNoCdm, flyIndsCdm};
goodFlyIndsCell = cellfun(@(x, y) x(y), flyIndsCell, goodCells, 'uni', 0);
offWidthCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, offWidthCenter, goodCells, ...
                                 'uni', 0);
offWidthSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, offWidthSurround, goodCells, ...
                                 'uni', 0);
onWidthCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, onWidthCenter, goodCells, ...
                                 'uni', 0);
onWidthSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, onWidthSurround, goodCells, ...
                                 'uni', 0);

offAmpCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, offAmpCenter, goodCells, ...
                                 'uni', 0);
offAmpSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, offAmpSurround, goodCells, ...
                                 'uni', 0);
onAmpCenterFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, onAmpCenter, goodCells, ...
                                 'uni', 0);
onAmpSurroundFlyMeans = cellfun(@(x, y, z) accumarray(x(z)', y, [], @mean), ...
                                 flyIndsCell, onAmpSurround, goodCells, ...
                                 'uni', 0);

emptyFlyInds = cellfun(@(x, y) find(accumarray(x(y)', 1) == 0), flyIndsCell, goodCells, 'uni', 0);
% Remove empty flies form array for plotting and quantification.
for iNeuron = 1: 2
    onAmpCenterFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
    offAmpCenterFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
    onAmpSurroundFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
    offAmpSurroundFlyMeans{iNeuron}(emptyFlyInds{iNeuron}) = [];
end


%% Plot amplitudes of center vs surround. Figure 7e
clear hSubAx
hFig = createPrintFig(10 * [3/2 1/2]);
% suptitle('Receptive field amplitude of center vs surround')
hSubAx(1) = subplot(1, 3, 1);
scatter(offAmpCenter{1}, offAmpSurround{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
hold on
scatter(offAmpCenter{2}, offAmpSurround{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
lsline
% refline(-1, 0)
axis equal square
title('OFF RF Amplitude')
xlabel('Center (\Delta F / F_0)')
ylabel('Surround (\Delta F / F_0)')
legend({'Control', 'CDM'}, 'Location', 'southwest')

hSubAx(2) = subplot(1, 3, 2);
scatter(onAmpCenter{1}, onAmpSurround{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
hold on
scatter(onAmpCenter{2}, onAmpSurround{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)

lsline
% refline(-1, 0)
axis equal square
title('ON RF Amplitude')

hSubAx(3) = subplot(1, 3, 3);
scatter(onAmpCenter{1}, onAmpSurround{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
hold on
scatter(onAmpCenter{2}, onAmpSurround{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
lsline
% refline(-1, 0)
axis equal square
xlim([-0.7 -0.2]);
ylim([0 0.5]);
title('ON RF Amplitude')
arrayfun(@prettifyAxes, hSubAx);
setFontForThesis(hSubAx, hFig);


% figDir = 'P:\Documents\Figures\Paper\Tm9\CDM\';
figFileName = ['CenterSurround-Amp-Scatter'];
set(gcf,'renderer','painters')
%     set(gcf, 'PaperPositionMode', 'auto');
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% Correlations for figure 7e
[offAmpCorr, offAmpCorrP] = cellfun(@(x, y) corr(x(:), y(:)), ...
                                    offAmpCenter, offAmpSurround, 'uni', 1);
[onAmpCorr, onAmpCorrP] = cellfun(@(x, y) corr(x(:), y(:)), ...
                                    onAmpCenter, onAmpSurround, 'uni', 1);
% Bootstrap CI
[offAmpCorrCI, offAmpCorrs] = cellfun(@(x, y) bootci(1000, @corr, x(:), y(:)), ...
                                      offAmpCenter, offAmpSurround, 'uni', 0);
[onAmpCorrCI, onAmpCorrs] = cellfun(@(x, y) bootci(1000, @corr, x(:), y(:)), ...
                                      onAmpCenter, onAmpSurround, 'uni', 0);                                  

offAmpCorrCIArray = cat(2, offAmpCorrCI{:}) % each column is one condition
offAmpCorrsArray = cellfun(@mean, offAmpCorrs)
onAmpCorrCIArray = cat(2, onAmpCorrCI{:}) % each column is one condition
onAmpCorrsArray = cellfun(@mean, onAmpCorrs)

%% Export data for figure 7e.
sheetName = ['Fig7e'];
startRow = 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_OFF_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_OFF_RF_Surround_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_RF_Surround_Amplitude_dF_F0']};
        
DataTable = table(offAmpCenter{1}(:), offAmpSurround{1}(:), ...
                  'VariableNames', varNames(1: 2));
writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);

DataTable = table(offAmpCenter{2}(:), offAmpSurround{2}(:), ...
                  'VariableNames', varNames(3: 4));
    
writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['D' num2str(startRow)]);

startRow = max(cellfun(@numel, offAmpCenter)) + 3;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_ON_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_ON_RF_Surround_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_RF_Surround_Amplitude_dF_F0']};
        
DataTable = table(onAmpCenter{1}(:), onAmpSurround{1}(:), ...
                  'VariableNames', varNames(1: 2));
writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);

DataTable = table(onAmpCenter{2}(:), onAmpSurround{2}(:), ...
                  'VariableNames', varNames(3: 4));
    
writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['D' num2str(startRow)]);

 
%% Same but per fly for revision
hFig = createPrintFig(10 * [2/2 1/2]);
% suptitle('Receptive field amplitude of center vs surround')
hSubAx(1) = subplot(1, 2, 1);
scatter(offAmpCenterFlyMeans{1}, offAmpSurroundFlyMeans{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
hold on
scatter(offAmpCenterFlyMeans{2}, offAmpSurroundFlyMeans{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
lsline
% refline(-1, 0)
axis equal square
title('OFF RF Amplitude')
xlabel('Center (\Delta F / F_0)')
ylabel('Surround (\Delta F / F_0)')
legend({'Control', 'CDM'}, 'Location', 'southwest')

hSubAx(2) = subplot(1, 2, 2);
scatter(onAmpCenterFlyMeans{1}, onAmpSurroundFlyMeans{1}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
hold on
scatter(onAmpCenterFlyMeans{2}, onAmpSurroundFlyMeans{2}, 50, 'filled', 'MarkerFaceAlpha', 0.5)
lsline
% refline(-1, 0)
axis equal square
title('ON RF Amplitude')

arrayfun(@prettifyAxes, hSubAx);
setFontForThesis(hSubAx, hFig);


% figDir = 'P:\Documents\Figures\Paper\Tm9\CDM\';
figFileName = ['CenterSurround-Amp-Scatter-FlyMeans'];
set(gcf,'renderer','painters')
%     set(gcf, 'PaperPositionMode', 'auto');
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% RF params stats from DoG fits

paramsToTest = {offAmpCenter, offAmpSurround, onAmpCenter, onAmpSurround, ...
                offWidthCenter, offWidthSurround, onWidthCenter, onWidthSurround};

%%
nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
% Delete flies with no rois to analyze.
nROIsPerFlyCell{1}(nROIsPerFlyCell{1} == 0) = [];
nROIsPerFlyCell{2}(nROIsPerFlyCell{2} == 0) = [];
% 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.
sample1 = 1: numel(nROIsPerFlyCell{1});
sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{2}));
%%
permutations = 10000;
sidedness = 'both';
plotresult = 1;

pVals = nan(1, numel(paramsToTest));
obsDiffs = nan(1, numel(paramsToTest));
for iParam = 1: numel(paramsToTest) 
    % Permutation test needs row vectors
    paramsQuality = cellfun(@(x) x(:)', paramsToTest{iParam}, 'uni', 0);
    [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
                                             paramsQuality, reorderedFlyIndsCell, ...
                                             permutations, sidedness, plotresult);

end
%% Plotting RF params for figure 7c-d

positiveSupport = [0-eps 10000];
negativeSupport = [-1000 0+eps];

hFig = createPrintFig(15 * [1/2 1/3]);
colors = getONOFFXYColors;
dark2 = brewermap(8, 'Dark2');
colors = dark2([3, 2], :);
% stimCellStr = {'XOFF', 'YOFF', 'XON', 'YON'};
stimCellStr = {'Control', 'CDM'};
titleStrCell = {'OFF Center', 'Surround', 'ON Center', 'Surround'};
dataCell = {offAmpCenter, offAmpSurround, onAmpCenter, onAmpSurround};
supportCell = {positiveSupport, negativeSupport, negativeSupport, positiveSupport};
for iAx = 1: 4
    hSubAx(iAx) = subplot(1, 4, iAx);
    title(titleStrCell{iAx})
    beanPlot(dataCell{iAx}, stimCellStr, colors, 1, gca, supportCell{iAx}, 'vertical', 0);
    text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
         sprintf('p = %1.1e', pVals(iAx)), 'HorizontalAlignment', 'center');
    axis tight
end
% linkaxes(hSubAx(1:2), 'xy');
% linkaxes(hSubAx(3:4), 'xy');
ylabel(hSubAx(1), 'RF Amplitude (\Delta F/F_0)')
hFig.Color = 'none';
set(hSubAx, 'Color', 'none');
arrayfun(@prettifyAxes, hSubAx)
setFontForThesis(hSubAx, hFig);

figFileName = ['CenterSurround-Amp-Bean'];
set(gcf,'renderer','painters')
print(hFig, [figDir figFileName '.pdf'], '-dpdf')

hFig = createPrintFig(15 * [1/2 1/3]);
dataCell = {offWidthCenter, offWidthSurround, onWidthCenter, onWidthSurround};
for iAx = 1: 4
    hSubAx(iAx) = subplot(1, 4, iAx);
    title(titleStrCell{iAx})
    beanPlot(dataCell{iAx}, stimCellStr, colors, 1, gca, positiveSupport, 'vertical', 0);
     text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
         sprintf('p = %1.1e', pVals(iAx + 4)), 'HorizontalAlignment', 'center');
    axis tight
end
% linkaxes(hSubAx(1:2), 'xy');
% linkaxes(hSubAx(3:4), 'xy');
ylabel(hSubAx(1), 'RF FWHM (deg)')
hFig.Color = 'none';
set(hSubAx, 'Color', 'none');
arrayfun(@prettifyAxes, hSubAx)
setFontForThesis(hSubAx, hFig);

figFileName = ['CenterSurround-Width-Bean'];
set(gcf,'renderer','painters')
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%% Export data for figure 7c-d
dataCell = [offAmpCenter(:)', offAmpSurround(:)', onAmpCenter(:)', onAmpSurround(:)'];

sheetName = ['Fig7cd'];
startRow = 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_OFF_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_OFF_RF_Surround_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_RF_Surround_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_ON_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_RF_Center_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_ON_RF_Surround_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_RF_Surround_Amplitude_dF_F0']};

capitalStr = 'ABCDEFGHIJKLMNOPQ';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end


% Export FWHM
dataCell = [offWidthCenter(:)', offWidthSurround(:)', onWidthCenter(:)', onWidthSurround(:)'];

startRow = max(cellfun(@numel, offAmpCenter)) + 3;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_OFF_RF_Center_FWHM_deg'], ...
            [conditionStrCell{2} '_OFF_RF_Center_FWHM_deg'], ...
            [conditionStrCell{1} '_OFF_RF_Surround_FWHM_deg'], ...
            [conditionStrCell{2} '_OFF_RF_Surround_FWHM_deg'], ...
            [conditionStrCell{1} '_ON_RF_Center_FWHM_deg'], ...
            [conditionStrCell{2} '_ON_RF_Center_FWHM_deg'], ...
            [conditionStrCell{1} '_ON_RF_Surround_FWHM_deg'], ...
            [conditionStrCell{2} '_ON_RF_Surround_FWHM_deg']};

capitalStr = 'ABCDEFGHIJKLMNOPQ';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end

%% Export data for figure 7c-e stats
dataCell = goodFlyIndsCell;

sheetName = ['Fig7cd'];
startRow = 2 * (max(cellfun(@numel, goodFlyIndsCell)) + 2) + 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_Fly_Index'], ...
            [conditionStrCell{2} '_Fly_Index']};
capitalStr = 'ABCD';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end

%% This is the single Gaussian fit version used for Supplementary figure 6.
rSqThresh = -20;
respIndThresh = 0.5;
alignAll = true;
[offX, offY, onX, onY] = deal(cell(size(StimTableNeuronStim([2, 1], 2: end), 1), 1));
jTmp = 1;
for iNeuron = [2, 1]
    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
    [offX{jTmp}, offY{jTmp}, ...
     onX{jTmp}, onY{jTmp}] = getRFsAlignedWithFit(StimTableCell, ...
                                                  rSqThresh, ...
                                                  respIndThresh, alignAll);
    jTmp = jTmp + 1;
end


%% Inline from plotRFsAlignedWithFitGenotypes
jTmp = 1;
for iNeuron = [2, 1]
    StimTableCell = StimTableNeuronStim(iNeuron, 2: end);
    [~, responseIndArrayCell, rSquareArrayCell, fitParamsCell, ...
      flyIndsCell, ~, ~] = getStimArraysForSameRois(StimTableCell);                                          
    % Align to fit
    widthCell{jTmp} = 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{jTmp} = cellfun(@(x, y) x(:) > rSqThresh & y(:) > respIndThresh, ...
                                     rSquareArrayCell, responseIndArrayCell, 'uni', 0)';
    jTmp = jTmp + 1;
end
%%
validOffInds = cellfun(@(x) x{1} & x{2}, validIndsCell, 'uni', 0);
validOnInds = cellfun(@(x)  x{3} & x{4}, validIndsCell, 'uni', 0);
validAllInds = cellfun(@(x)  x{1} & x{2} & x{3} & x{4}, validIndsCell, 'uni', 0);

%%
widthCell = cellfun(@(x) cat(1, x{:}), widthCell, 'uni', 0);
%%
xOFFWidth = cellfun(@(x, y) x(1, y), widthCell, validAllInds, 'uni', 0);
yOFFWidth = cellfun(@(x, y) x(2, y), widthCell, validAllInds, 'uni', 0);
xONWidth = cellfun(@(x, y) x(3, y), widthCell, validAllInds, 'uni', 0);
yONWidth = cellfun(@(x, y) x(4, y), widthCell, validAllInds, 'uni', 0);
allWidthsCell = {xOFFWidth, yOFFWidth, xONWidth, yONWidth};

%%
for iW = 1: numel(allWidthsCell)
    allWidthsCell{iW}(cellfun(@isempty, allWidthsCell{iW}, 'uni', 1)) = {pi};
end

%%
offWidthCell = cellfun(@(x, y) (x + y) / 2, allWidthsCell{1}, allWidthsCell{2}, 'uni', 0);
onWidthCell = cellfun(@(x, y) (x + y) / 2, allWidthsCell{2}, allWidthsCell{3}, 'uni', 0);

%% Statistics
% reviewer #1 analysis by fly
% Number of ROIs differs across stimuli, but this has been accounted at
% least for bars in getStimArraysForSameRois

goodFlyIndsCell = cellfun(@(x, y) x(y), flyIndsCell(1: 2), validAllInds, 'uni', 0);

nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
% Delete flies with no rois to analyze.
nROIsPerFlyCell{1}(nROIsPerFlyCell{1} == 0) = [];
nROIsPerFlyCell{2}(nROIsPerFlyCell{2} == 0) = [];
% 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.
sample1 = 1: numel(nROIsPerFlyCell{1});
sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{2}));
%% Statistics for supplementary figure 6b
paramsCell = {offWidthCell, onWidthCell};
permutations = 10000;
sidedness = 'both';
plotresult = 1;

pVals = nan(1, numel(paramsCell));
obsDiffs = nan(1, numel(paramsCell));
for iParam = 1: numel(paramsCell) 
    % Permutation test needs row vectors
    paramsQuality = cellfun(@(x) x(:)', paramsCell{iParam}, 'uni', 0);
    [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
                                             paramsQuality, reorderedFlyIndsCell, ...
                                             permutations, sidedness, plotresult);

end
%% Plot supplementary figure 6b
hFig = createPrintFig(15 * [1/4 1/3]);
colors = brewermap(4, 'Set1');
colors = colors([1, 2, 4, 3], :);
hSubAx(1) = subplot(1, 2, 1);
hWidthAx = beanPlot(offWidthCell, fliplr(cellTypeStr), colors, ...
                    1, gca, [0-eps 100], 'vertical', 0, 1, 0);
text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
     sprintf('p = %1.1e', pVals(1)), 'HorizontalAlignment', 'center');
hSubAx(2) = subplot(1, 2, 2);                
hWidthAx = beanPlot(onWidthCell, fliplr(cellTypeStr), colors, ...
                    1, gca, [0-eps 100], 'vertical', 0, 1, 0);
text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
     sprintf('p = %1.1e', pVals(2)), 'HorizontalAlignment', 'center');                
axis(hSubAx(1), 'tight')
axis(hSubAx(2), 'tight')

% linkaxes(hSubAx, 'y');
setFontForThesis(hWidthAx, gcf);
arrayfun(@prettifyAxes, hWidthAx)
% arrayfun(@offsetAxes, hWidthAx)
figFileName = ['RF-FWHM-XYMean-bean'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')

%% Export data for supplementary figure 6b
dataCell = [offWidthCell(:)', onWidthCell(:)'];

sheetName = ['SuppFig6b'];
startRow = 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_OFF_RF_Center_FWHM_deg'], ...
            [conditionStrCell{2} '_OFF_RF_Center_FWHM_deg'], ...
            [conditionStrCell{1} '_ON_RF_Center_FWHM_deg'], ...
            [conditionStrCell{2} '_ON_RF_Center_FWHM_deg']};

capitalStr = 'ABCDEFGHIJKLMNOPQ';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end

%% Min max orientation mean
maxAmpCell = cellfun(@(x) max(x, [], 2), [offX offY onX onY], 'uni', 0);
minAmpCell = cellfun(@(x) min(x, [], 2), [offX offY onX onY], 'uni', 0);
offMaxCell = cellfun(@(x, y) (x + y) / 2, maxAmpCell(:, 1), maxAmpCell(:, 2), 'uni', 0);
onMaxCell = cellfun(@(x, y) (x + y) / 2, maxAmpCell(:, 3), maxAmpCell(:, 4), 'uni', 0);
offMinCell = cellfun(@(x, y) (x + y) / 2, minAmpCell(:, 1), minAmpCell(:, 2), 'uni', 0);
onMinCell = cellfun(@(x, y) (x + y) / 2, minAmpCell(:, 3), minAmpCell(:, 4), 'uni', 0);

%% Statistics for supplementary figure 6a

paramsCell = {offMaxCell, offMinCell, onMinCell, onMaxCell};
permutations = 10000;
sidedness = 'both';
plotresult = 1;

pVals = nan(1, numel(paramsCell));
obsDiffs = nan(1, numel(paramsCell));
for iParam = 1: numel(paramsCell) 
    % Permutation test needs row vectors
    paramsQuality = cellfun(@(x) x(:)', paramsCell{iParam}, 'uni', 0);
    [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
                                             paramsQuality, reorderedFlyIndsCell, ...
                                             permutations, sidedness, plotresult);

end

%% Plot supplementary figure 6a
hFig = createPrintFig(15 * [1/2 1/3]);
colors = getONOFFXYColors;
colors = lines(2);
stimCellStr = {'Control', 'CDM'};
titleStrCell = {'OFF Max (center)', 'OFF Min (surround)', ...
                'ON Min (center)', 'ON Max (surround)'};
positiveSupport = [0-eps 1000];
negativeSupport = [-1000 0-eps];
supportCell = {negativeSupport, positiveSupport};
supportInds = [2 1 1 2];
for iParam = 1: numel(paramsCell)
    hSubAx(iParam) = subplot(1, numel(paramsCell), iParam);
    title(titleStrCell{iParam});
    beanPlot(paramsCell{iParam}, stimCellStr, colors, 1, gca, ...
             supportCell{supportInds(iParam)}, 'vertical', 0);
    text(mean(xlim), 0.9 * getArrayElem(ylim, 2), ...
         sprintf('p = %1.1e', pVals(iParam)), 'HorizontalAlignment', 'center');
     axis tight
end

% linkaxes(hSubAx(1:2), 'xy');
% linkaxes(hSubAx(3:4), 'xy');
ylabel(hSubAx(1), 'RF max (\DeltaF/F_0)')
% hFig.Color = 'none';
set(hSubAx, 'Color', 'none');
arrayfun(@prettifyAxes, hSubAx)
setFontForThesis(hSubAx, hFig);

figFileName = ['MaxMinAmp-XYMean-Bean'];
set(gcf,'renderer','painters')
print(hFig, [figDir figFileName '.pdf'], '-dpdf')

%% Export data for supplementary figure 6b
dataCell = [offMaxCell(:)', offMinCell(:)', onMinCell(:)', onMaxCell(:)'];

sheetName = ['SuppFig6a'];
startRow = 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_OFF_RF_Maximum_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_RF_Maximum_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_OFF_RF_Minimum_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_RF_Minimum_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_ON_RF_Minimum_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_RF_Minimum_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_ON_RF_Maximum_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_RF_Maximum_Amplitude_dF_F0']};

capitalStr = 'ABCDEFGHIJKLMNOPQ';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end

%% Data for figure 7f-g
for iNeuron = 1: size(StimTableNeuronStim, 1)
    allOnOffTracesCell = cat(1, StimTableNeuronStim{iNeuron, 1}.paddedEpochStacks{:});
    allOnOffTracesCell = cellfun(@(x) x', allOnOffTracesCell, 'uni', 0);
    allOffTraces = cat(1, allOnOffTracesCell{:, 1});
    allOnTraces = cat(1, allOnOffTracesCell{:, 2});
    allOnOffTracesArray = [allOffTraces allOnTraces];
    allNeuronsOnOffTraces{iNeuron} = allOnOffTracesArray;
end
%%

% reviewer #1 analysis by fly
flyIndsCdm = repelem(1:size(StimTableNeuronStim{1, 1}), ... 
                  cellfun(@(x) size(x, 2), StimTableNeuronStim{1, 1}.responseIndex));
flyIndsNoCdm = repelem(1:size(StimTableNeuronStim{2, 1}), ... 
                  cellfun(@(x) size(x, 2), StimTableNeuronStim{2, 1}.responseIndex));

%% Plot figure 7g
cellTypeStr = {'CDM', 'noCDM'};
desiredCellTypesStr = {'noCDM', 'CDM'};
cellTypeInds = cellfun(@(x) find(strcmp(cellTypeStr, x)), desiredCellTypesStr);
% cellTypeInds = 1:2;
out = getONOFFParamsFromTable(StimTableNeuronStim);

paramsCell = {out.offMean, out.onMean, out.responseInd};
supportCell = {'unbounded', 'unbounded', [-1.01 1.01], [-1 1], [0 2], [-eps 2.01], [0 1]};
paramsNames = {'Mean OFF', 'Mean ON', 'Polarity index', 'Sustenance index', ...
               'Half-rise time (s)', 'Time to extreme (s)', 'Response quality'};
qualityThreshold = 0.5;
paramsCell = cellfun(@(x) x(cellTypeInds), paramsCell, 'uni', 0);
colors = repmat([1 0.75 1] * 0.25, numel(cellTypeInds), 1);
colors = lines(2);
hFig = createPrintFig(15 * [1/3 1]);
for iParam = 1: numel(paramsCell) - 1
    hSubAx(iParam) =  subplot(2, 1, iParam);
    title(paramsNames{iParam});
    paramsQuality = cellfun(@(x, y) x(y > qualityThreshold), paramsCell{iParam}, paramsCell{end}, 'uni', 0);
    paramsToPlot = paramsQuality;
%     paramsToPlot = paramsCell{iParam};
    if 1
        beanPlot(paramsToPlot, desiredCellTypesStr, colors, 1, gca, ...
                 supportCell{iParam}, 'vertical', 0, 0, 1);
        arrayfun(@prettifyAxes, gca);
        arrayfun(@offsetAxes, gca);
        setFontForThesis(gca, hFig)
    else
        boxInd = repelem(1:numel(paramsToPlot), cellfun(@numel, paramsToPlot));
        boxplot(cat(2, paramsToPlot{:}), boxInd, 'Notch', 'on')
    end
    if numel(paramsToPlot) > 1
        for jNeuron = 2: numel(paramsToPlot) 
            [pVal(iParam, jNeuron - 1), ...
             obsDiff(iParam, jNeuron - 1)] = permutationTest(paramsToPlot{1}, ...
                                                             paramsToPlot{jNeuron}, ...
                                                             50000, 'plotresult', 0);
        end
    end
end

setFontForThesis(gca, hFig)
figFileName = ['OnOff-Params-'];
set(gcf,'renderer','painters')
%     set(gcf, 'PaperPositionMode', 'auto');
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
print(hFig, [figDir figFileName '.png'], '-dpng', '-r300')

%% Export data for figure 7g
paramsQuality = arrayfun(@(z) cellfun(@(x, y) x(y > qualityThreshold), paramsCell{z}, paramsCell{end}, 'uni', 0), 1:2,'uni', 0);

dataCell = [paramsQuality{1}(:)', paramsQuality{2}(:)'];

sheetName = ['Fig7g'];
startRow = 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_OFF_Fullfield_Mean_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_OFF_Fullfield_Mean_Amplitude_dF_F0'], ...
            [conditionStrCell{1} '_ON_Fullfield_Mean_Amplitude_dF_F0'], ...
            [conditionStrCell{2} '_ON_Fullfield_Mean_Amplitude_dF_F0']};

capitalStr = 'ABCDEFGHIJKLMNOPQ';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end

%% plot figure 7f
onOffTraces = cellfun(@(x, y) x(y > qualityThreshold, :), ...
                        allNeuronsOnOffTraces(cellTypeInds), out.responseInd(cellTypeInds), 'uni', 0);

colors = repmat([1 0.75 1] * 0.25, numel(cellTypeStr), 1);
% colors = brewermap(numel(desiredCellTypesStr), 'Dark2');
colors = brewermap([], 'Dark2');
colors = colors([8,  (1: numel(cellTypeStr) - 1) + 3*0], :);
% For cdm
colors = lines(2);
plotONOFFTracesFromCell(onOffTraces, desiredCellTypesStr, figDir, colors)
%% Export data for figure 7f
dataFileName = [figDir filesep 'Fig7.xls'];

timePoints = (1: numel(onOffTraces{1}(1, :))) / 10; % Data interpolated at 10Hz

sheetName = 'Fig7f';
stimNameStrCell = {'Tm9_Control_OFFON_fullfield_time_seconds', ...
                   'Tm9_CDM_OFFON_fullfield_time_seconds'};
startRow = 1;

for iStim = 1: 2
    DataTable = createRFsTable(onOffTraces{iStim}, stimNameStrCell{iStim});
    DataTable(:, 1).(stimNameStrCell{iStim}) = timePoints(:);
    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
    startRow = startRow + size(DataTable, 1) + 2;
end
%% More sophisticated test for reviewer #1
flyIndsCell = {flyIndsNoCdm flyIndsCdm};
goodFlyIndsCell = cellfun(@(x, y) x(y > qualityThreshold), flyIndsCell, paramsCell{end}, 'uni', 0);

flyGroups = [goodFlyIndsCell{1} max(goodFlyIndsCell{1}) + goodFlyIndsCell{2}];
genotypeGroups = [repmat({'control'}, 1, numel(goodFlyIndsCell{1})) ...
                  repmat({'cdm'}, 1, numel(goodFlyIndsCell{2}))];
groups = {genotypeGroups, flyGroups};
% I think the fly level is a random effect, right?
pValsCell = cell(1, numel(paramsCell));
for iParam = 1: numel(paramsCell) 
    paramsQuality = cellfun(@(x, y) x(y > qualityThreshold), paramsCell{iParam}, paramsCell{end}, 'uni', 0);

    [pValsCell{iParam}, ~, statsCell{iParam}] = ...
            anovan(cat(2, paramsQuality{:}), ...
                   groups, ...
                   'nested',[0 0; 1 0], ...
                   'varnames', {'Genotype', 'Fly'}, ...
                   'random', [2], 'display', 'on');
end
%%
hFig = createPrintFig(10 * [1 1/2]);
hIm = imagesc([], [1, 2], cat(2, pValsCell{:}));
colormap([plasma; 0.5 0.5 0.5]);
caxis([0 0.051]);
set(gca, 'XTick', 1: numel(paramsNames), 'XTickLabels', paramsNames, 'XTickLabelRotation', 90);
set(gca, 'YTick', [1 2], 'YTickLabels', {'Pooling across cells', 'Pooling across flies'}, ...
   'XTickLabelRotation', 30)
title('p-values from nested anova FFF');
colorbar;
arrayfun(@prettifyAxes, gca)
setFontForThesis(gca, hFig);

figFileName = ['NestedAnova-pValues-ONOFF-FFF-'];
set(gcf,'renderer','painters')
print(hFig, [figDir figFileName '.pdf'], '-dpdf');

%% Still do the nested permutation
nROIsPerFlyCell = cellfun(@(x) accumarray(x(:), 1), goodFlyIndsCell, 'uni', 0);
% Delete flies with no rois to analyze.
nROIsPerFlyCell{1}(nROIsPerFlyCell{1} == 0) = [];
nROIsPerFlyCell{2}(nROIsPerFlyCell{2} == 0) = [];
% 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.
sample1 = 1: numel(nROIsPerFlyCell{1});
sample2 = numel(sample1) + (1: numel(nROIsPerFlyCell{2}));
%% Statistics for figure 7g
permutations = 10000;
sidedness = 'both';
plotresult = 1;

pVals = nan(1, numel(paramsCell));
obsDiffs = nan(1, numel(paramsCell));
for iParam = 1: numel(paramsCell) 
    paramsQuality = cellfun(@(x, y) x(y > qualityThreshold), paramsCell{iParam}, paramsCell{end}, 'uni', 0);
               
    [obsDiffs(iParam), pVals(iParam)] = nestetPermutationTest(sample1, sample2, ...
                                             paramsQuality, reorderedFlyIndsCell, ...
                                             permutations, sidedness, plotresult);

end

%% Plot the pvalues as bars
hFig = createPrintFig(15 * [1 3/4]);

bar(pVals, 'FaceColor', [1 1 1] * 0.15, 'EdgeColor', 'none');
arrayfun(@(x) text(x, 0.02 + pVals(x), num2str(pVals(x), '%1.2e'), ...
                   'HorizontalAlignment', 'center'), 1: numel(pVals))
ylabel('p-value')
hRefLine = refline(0, 0.05);
hRefLine.Color = [1 1 1] * 0.5;
hRefLine.LineStyle = ':';
set(gca, 'XTick', 1: numel(paramsNames), ...
    'XTickLabels', paramsNames, 'XTickLabelRotation', 30);
set(gca, 'YTick', [0, 0.01, 0.05, 0.1, 0.2]);
title('p-values from nested permutation test');
arrayfun(@prettifyAxes, gca)
setFontForThesis(gca, hFig);

figFileName = ['ONOFF-NestedPermutations-pValues'];
set(gcf,'renderer','painters')
print(hFig, [figDir figFileName '.pdf'], '-dpdf');

%% Export data for figure 7g stats
dataCell = goodFlyIndsCell;

sheetName = ['Fig7g'];
startRow = 1;

conditionStrCell = {'Tm9_Control', 'Tm9_CDM'};
varNames = {[conditionStrCell{1} '_Fly_Index'], ...
            [conditionStrCell{2} '_Fly_Index']};
capitalStr = 'EFGHIJKLMNOPQ';
% Export response amplitudes
for iVar = 1: numel(dataCell)
    DataTable = table(dataCell{iVar}(:), 'VariableNames', varNames(iVar));
    writetable(DataTable, dataFileName, 'Sheet', sheetName, ...
               'Range', [capitalStr(iVar) num2str(startRow)]);
end

%%
function [Amp, Pos, Sigma, RSq, varargout] = parseFitParams(fitStructCellArray)

TmpGofCell = cellfun(@(y) arrayfun(@(x) x.gof.rsquare, y(1: end-1)), fitStructCellArray, 'uni', 0);

if isstruct(fitStructCellArray{1, 1}(1).fit)
    fitParams = fieldnames(fitStructCellArray{1, 1}(1).fit);
    ParamCell = cell(1, numel(fitParams));
    for iField = 1: numel(fitParams)
        % Delete the last ROI from background
        fieldCell = cellfun(@(y) arrayfun(@(x) x.fit.(fitParams{iField}), y(1: end-1)), fitStructCellArray, 'uni', 0);
        ParamCell{iField} = cell(1, size(fieldCell, 2));
        for jFly = 1: size(ParamCell{iField}, 2)
            ParamCell{iField}{jFly} = cat(1, fieldCell{:, jFly});
        end
    end
else
    % Params are given as an array.
    ParamCell = cell(1, numel(fitStructCellArray{1, 1}(1).fit));
    for iField = 1: numel(fitStructCellArray{1, 1}(1).fit)
        % Delete the last ROI from background
        fieldCell = cellfun(@(y) arrayfun(@(x) x.fit(iField), y(1: end-1))', fitStructCellArray, 'uni', 0);
        fieldCell = squeeze(fieldCell);
        ParamCell{iField} = cell(1, size(fieldCell, 2));
        for jFly = 1: size(ParamCell{iField}, 2)
            ParamCell{iField}{jFly} = cat(1, fieldCell{:, jFly});
        end
    end
    TmpGofCell = cellfun(@(x) x', TmpGofCell, 'uni', 0);
end

for jFly = 1: size(TmpGofCell, 2)
    RSq{jFly} = cat(1, TmpGofCell{:, jFly});
end

if numel(ParamCell) > 5
    [Amp.center, Pos, Sigma.center, Amp.surround, Sigma.surround, baseline] = deal(ParamCell{:});
else
    [Amp.center, Pos, Sigma.center, Amp.surround, Sigma.surround] = deal(ParamCell{:});
end
if nargout > 4
    if numel(ParamCell) > 5
        varargout{1} = baseline;
    else
        varargout{1} = nan;
    end
end



end