ccc
%% Define search path.
parentDir = fileparts(fileparts(mfilename('fullpath')));
saveFolder = [parentDir filesep 'data'];
load([saveFolder filesep 'Fig6-T5Tm9-processedTable'], 'ProcessedTable');

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

%%
stimSetCellStr = getStimSetCellStr(4);

StimTable = getStimSubsetTable(ProcessedTable, stimSetCellStr, 7);

tuningMethods = StimTable.nonParamTuning{1}(1).tuningMethods;

fieldNamesToUse = {'spacing', 'temporalFrequency'};
[ stimSize, sortedFieldNameInds, fieldNames, stimParamsFlat ] = ...
    StimulusDomain_ND( StimTable.stimParams{1}{1}, fieldNamesToUse );

%%
spacing = [StimTable.stimParams{1}{1}(2:end).spacing];
stimtrans = [StimTable.stimParams{1}{1}(2:end).stimtrans];
temporalFreq = unique([stimtrans.mean] ./ spacing);
spatialFreq = 1 ./ unique(spacing);

%%
selectedTuningMethod = '1F_amplitude';

tcRoiCell = cellfun(@(x) ...
    getTuningCurvePerROI(x(1).tc, tuningMethods, '1F_amplitude'), ...
    StimTable.nonParamTuning, 'uni', 0);
tm9TCs = [tcRoiCell{:}];
%%
tcRoiCell = cellfun(@(x) ...
    getTuningCurvePerROI(x(2).tc, tuningMethods, '1F_amplitude'), ...
    StimTable.nonParamTuning, 'uni', 0);
t5TCs = [tcRoiCell{:}];
%%

[xGrid, yGrid] = meshgrid(spatialFreq, temporalFreq);

%% These are the maps which average is shown in figure 6g
respIndex = cell2mat(StimTable.responseIndex')';

validInds = all(respIndex > 0.5, 2);
figure
hAx = createSquarishSubplotGrid(sum(validInds));
iAx = 1;
for iRoi = find(validInds)'
    imagesc(hAx(iAx),  t5TCs{iRoi})
    iAx = iAx + 1;
end
figure
hAx = createSquarishSubplotGrid(sum(validInds));
iAx = 1;
for iRoi = find(validInds)'
    imagesc(hAx(iAx),  tm9TCs{iRoi})
    iAx = iAx + 1;
end
%% Export maps for figure 6g
% Columns are increasing temporal frequency, rows are increasing spatial
% frequency. But transpose the array so it matches the spatial and temporal frequency
% axes. Then columns are increasing spatial frequency and rows increasing
% temporal frequency.
[spatialFreqGrid, temporalFreqGrid] = meshgrid(spatialFreq, temporalFreq); % xGrid-space, yGrid-time
% Create table to append columns with ROI tunint to it.
DataTable = table(spatialFreqGrid(:), temporalFreqGrid(:), 'VariableNames', {'spatial_frequency_cycles_per_deg', 'temporal_frequency_Hz'});

for iRoi = find(validInds)'
    tm9Map = tm9TCs{iRoi}';
    t5Map = t5TCs{iRoi}';
    DataTable.(sprintf(['ROI%d_Tm9_F1_amplitude'], iRoi)) = tm9Map(:);
    DataTable.(sprintf(['ROI%d_T5_F1_amplitude'], iRoi)) = t5Map(:);
end

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

sheetName = 'Fig6g';
writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', 'A1');

%% Figure 6h
hFig = createPrintFig(10*[1/3 1]);        
support = [-1-eps 1+eps];
plotOrientation = 'vertical'; % Can change to plotHorizontalBool.
halveBean = false;
useAlpha = true;
bwOption = true;
mapCorrs = cellfun(@(x,y) corr2(x, y), t5TCs(validInds), tm9TCs(validInds));
beanPlot({mapCorrs},{'Tm9-T5'},'k', 1, gca, support, plotOrientation, halveBean);
% beanPlot({mapCorrs},{'Tm4-T5'},'k', 1, gca, support, plotOrientation, halveBean);
print(gcf, [figDir 'tuningMap-corr-bean' '.pdf'], '-dpdf', '-r300');
%% Export pairwise correlation values ffrom figure 6h

labelStrCell = sprintfc(['ROI%d_T5vsTm9_F1_Amplitude_CorrelationCoefficient'], find(validInds));
for iRow = 1: numel(labelStrCell)
    xlswrite(dataFileName, labelStrCell(iRow), 'Fig6h', ['A' num2str(iRow)])
    xlswrite(dataFileName, {mapCorrs(iRow)}, 'Fig6h', ['B' num2str(iRow)])
end
%% Shuffle control for correlation of maps for confidence interval in Fig 6h
validT5TCs = t5TCs(validInds);
validTm9TCs = tm9TCs(validInds);
mapCorrs = cellfun(@(x,y) corr2(x, y), validT5TCs, validTm9TCs);
nShuffles = 1000000;
% exactShuffles = nchoosek(2 * sum(validInds), sum(validInds)); %even 14
% gives 40116600 shuffles, so 1e6 is ok approximation
shuffleMapCorrs = zeros(nShuffles, 1);

for iShuffle = 1: nShuffles
    shuffleInds = randperm(numel(validTm9TCs));
    shuffleMapCorrs(iShuffle) = mean(cellfun(@(x,y) corr2(x, y), ...
                                             validT5TCs, ...
                                             validTm9TCs(shuffleInds)));
end
pValueCorrMaps = sum(shuffleMapCorrs > mean(mapCorrs)) / nShuffles;
corrConfInt = arrayfun(@(x) quantile(shuffleMapCorrs, x), [0.05, 0.95]);
meanShuffleCorr = mean(shuffleMapCorrs);

%%
% Plot spatial tuning by collapsing over temporal frequencies.
spatialTuning{1} = cell2mat(cellfun(@(x) mean(x, 2)', t5TCs, 'uni', 0)');
spatialTuning{2} = cell2mat(cellfun(@(x) mean(x, 2)', tm9TCs, 'uni', 0)');
% Plot temporal tuning by collapsing over spatial frequencies.
temporalTuning{1} = cell2mat(cellfun(@(x) mean(x, 1), t5TCs, 'uni', 0)');
temporalTuning{2} = cell2mat(cellfun(@(x) mean(x, 1), tm9TCs, 'uni', 0)');
%%
spatialTuning = cellfun(@(x) x./max(x, [], 2), spatialTuning, 'uni', 0);
temporalTuning = cellfun(@(x) x./max(x, [], 2), temporalTuning, 'uni', 0);
validInds = all(respIndex > 0.5, 2); % Response quality threshold

%% From here we get mean maps in figure 6g
selectedTuningMethod = '1F_amplitude';
createPrintFig(10*[1 1]);
plotROISpatioTemporalTuning(spatialFreq, temporalFreq, mean(cat(3, t5TCs{validInds}), 3))
print(gcf, [figDir 'tuningMap-roiMean-0p5respInd' selectedTuningMethod 'T5' '.pdf'], '-dpdf', '-r300');

createPrintFig(10*[1 1]);
plotROISpatioTemporalTuning(spatialFreq, temporalFreq, mean(cat(3, tm9TCs{validInds}), 3))
print(gcf, [figDir 'tuningMap-roiMean-0p5respInd' selectedTuningMethod 'Tm9' '.pdf'], '-dpdf', '-r300');

%%
ProcessedTable = ProcessedTable(ProcessedTable.typeOfProblem == 0, :);

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: 7, [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(cat(1, sortFlyInds{:}), :);
StimOrderedTable = StimOrderedTable(StimOrderedTable.flyInd~=7,:);
%% 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
%%
stimInd = 1;
nStims = numel(stimSetCellStr);
for iNeuron = 1: numel(cellTypeNum)
    for jStim = 1: nStims
        StimTableNeuronStim{iNeuron, jStim} = getStimSubsetTable(StimTableNeuronCell{iNeuron}, stimSetCellStr, jStim);
    end
end

%% Extract the parameters for the analysis of spatial RFs.
clear amp pos width rsquare responseIndex
jBar = 1;
for iStim = 2:5
    dummyTable = StimTableNeuronStim{1, iStim};
    for jChannel = 1: 2
        responseIndex{jChannel}(jBar, :) = cell2mat(cellfun(@(x, y) x(jChannel, 1: end - 1*0), dummyTable.responseIndex, 'uni', 0)');        
        paramsCell = cellfun(@(x) x(jChannel, 1: end - 1*0), dummyTable.paramTuning, 'uni', 0);
        amp{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.a1, x, 'uni', 1), paramsCell, 'uni', 0)');
        pos{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.b1, x, 'uni', 1), paramsCell, 'uni', 0)');
        width{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.c1, x, 'uni', 1), paramsCell, 'uni', 0)');
        rsquare{jChannel}(jBar, :) = cell2mat(cellfun(@(x) arrayfun(@(y) y.gof.rsquare, x, 'uni', 1), paramsCell, 'uni', 0)');      
    end
    jBar = jBar + 1;
end
pos = cellfun(@(x) 2 * (x - 1), pos, 'uni', 0);
width = cellfun(@(x) 2 * (2 * sqrt(log(2))) * x, width, 'uni', 0);
%%
flyInds = repelem(1:size(StimTableNeuronStim{1}), ... 
                  cellfun(@(x) size(x, 2), StimTableNeuronStim{1}.responseIndex));
responseIndexOff = cat(1, responseIndex{1}(1:2, :), responseIndex{2}(1:2, :));
rSquareOff = cat(1, rsquare{1}(1:2, :), rsquare{2}(1:2, :));

minOffRSquare = cell2mat(cellfun(@(x) min(x(1:2, :)), rsquare, 'uni', 0)');
minOnRSquare = cell2mat(cellfun(@(x) min(x(3:4, :)), rsquare, 'uni', 0)');
thresholdRSquare = 0.25;
minRSquare = min(minOffRSquare);
minRespInd = min(responseIndexOff);
isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;

%% Figure 6f and supplementary figure 5a
responseIndexOff = cat(1, responseIndex{1}(1:2, :), responseIndex{2}(1:2, :));
responseIndexOn = cat(1, responseIndex{1}(3:4, :), responseIndex{2}(3:4, :));

rSquareOff = cat(1, rsquare{1}(1:2, :), rsquare{2}(1:2, :));
rSquareOn = cat(1, rsquare{1}(3:4, :), rsquare{2}(3:4, :));

minOffRSquare = cell2mat(cellfun(@(x) min(x(1:2, :)), rsquare, 'uni', 0)');
minOnRSquare = cell2mat(cellfun(@(x) min(x(3:4, :)), rsquare, 'uni', 0)');

plotOFF = true;
if plotOFF
    minRespInd = min(responseIndexOff);
    thresholdRSquare = 0.25;
    minRSquare = min(minOffRSquare);
    isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
else
    if 0
        minRespInd = max(responseIndexOn);
        thresholdRSquare = 0;
        minRSquare = max(minOnRSquare);
        isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
    else
        minRespInd = min(responseIndexOff);
        thresholdRSquare = 0.25;
        minRSquare = min(minOffRSquare);
        isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
    end
end

clear hSubAx
hFig = createPrintFig(20 * [1 2/3]);
set(hFig,'color','w');

colorLims = [-1 1] + 0*thresholdRSquare;
fitParamsCell = {width , pos};
for jOrientation = 1: 2
    for iParam = 1: numel(fitParamsCell)
        iAx = numel(fitParamsCell) * (jOrientation - 1) + iParam;
        hSubAx(iAx) = subplot(2, numel(fitParamsCell), iAx);
        if plotOFF
            x = fitParamsCell{iParam}{1}(jOrientation, isAboveThreshold);
            y = fitParamsCell{iParam}{2}(jOrientation, isAboveThreshold);
        else
            x = fitParamsCell{iParam}{1}(2 + jOrientation, isAboveThreshold);
            y = fitParamsCell{iParam}{2}(2 + jOrientation, isAboveThreshold);
        end
        hScatter = scatter(x, y, [], [1 1 1]*0.25, 'filled', 'MarkerFaceAlpha', 1);

        minVal = min([hSubAx(iAx).XLim(1) hSubAx(iAx).YLim(1)]);
        maxVal = min([hSubAx(iAx).XLim(2) hSubAx(iAx).YLim(2)]);
        lsline
        line([minVal maxVal], [minVal maxVal], 'Color', 0.5 * [1 1 1]);
        hChildren = get(hSubAx(iAx), 'Children');
        set(gca, 'Children', [hChildren(end); hChildren(1: end - 1)]);
        axis square;
        [corrCoef, p] = corr(x(:), y(:));
        legStr = ['n = ' num2str(numel(x)) ', corr = ' num2str(corrCoef, 2) ', p = ' num2str(p, 2)];
        hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
        hLeg.Box = 'off';
    end
end
linkaxes(hSubAx([1,3]), 'xy')
linkaxes(hSubAx([2,4]), 'xy')
hYLabel(1) = ylabel(hSubAx(1), 'T5 OFF RF FWHM Horizontal Bar (deg)');
hXLabel(1) = xlabel(hSubAx(1), 'Tm9 OFF RF FWHM Horizontal Bar (deg)');
hYLabel(2) = ylabel(hSubAx(2), 'T5 OFF RF Center Horizontal Bar (deg)');
hXLabel(2) = xlabel(hSubAx(2), 'Tm9 OFF RF Center Horizontal Bar (deg)');
hYLabel(3) = ylabel(hSubAx(3), 'T5 OFF RF FWHM Vertical Bar (deg)');
hXLabel(3) = xlabel(hSubAx(3), 'Tm9 OFF RF FWHM Vertical Bar (deg)');
hYLabel(4) = ylabel(hSubAx(4), 'T5 OFF RF Center Vertical Bar (deg)');
hXLabel(4) = xlabel(hSubAx(4), 'Tm9 OFF RF Center Vertical Bar (deg)');

[hXLabel.FontSize] = deal(10);
[hYLabel.FontSize] = deal(10);
hLeg.FontSize = 8;
[hSubAx.FontSize] = deal(8);
[hSubAx.Color] = deal('none');

arrayfun(@prettifyAxes, [hSubAx]); 
if plotOFF
    print(gcf, [figDir 'rfFitParamsScatter-Off-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
else
    print(gcf, [figDir 'rfFitParamsScatter-On-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
end
%% Export data for figure 6f ans supplementary figure 6a
tm9XFWHM = width{1}(1, isAboveThreshold);
t5XFWHM = width{2}(1, isAboveThreshold);
tm9YFWHM = width{1}(2, isAboveThreshold);
t5YFWHM = width{2}(2, isAboveThreshold);

tm9XPos = pos{1}(1, isAboveThreshold);
t5XPos = pos{2}(1, isAboveThreshold);
tm9YPos = pos{1}(2, isAboveThreshold);
t5YPos = pos{2}(2, isAboveThreshold);


DataTable = table(tm9XFWHM(:), t5XFWHM(:), tm9YFWHM(:), t5YFWHM(:), ...
                  tm9XPos(:), t5XPos(:), tm9YPos(:), t5YPos(:), ...
                  'VariableNames', ...
                   {'Tm9_OFF_RF_FWHM_Horizontal_Bar_deg', ...
                    'T5_OFF_RF_FWHM_Horizontal_Bar_deg', ...
                    'Tm9_OFF_RF_FWHM_Vertical_Bar_deg', ...
                    'T5_OFF_RF_FWHM_Vertical_Bar_deg', ...
                    'Tm9_OFF_RF_Position_Horizontal_Bar_deg', ...
                    'T5_OFF_RF_Position_Horizontal_Bar_deg', ...
                    'Tm9_OFF_RF_Position_Vertical_Bar_deg', ...
                    'T5_OFF_RF_Position_Vertical_Bar_deg', ...
                   });
   

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

sheetName = 'Fig6f';
writetable(DataTable(:, 1: 4), dataFileName, 'Sheet', sheetName, 'Range', 'A1');

sheetName = 'SuppFig6ab';
writetable(DataTable(:, 5: 8), dataFileName, 'Sheet', sheetName, 'Range', 'A1');

%% Reply to reviwer#1 comment 5
doBootstrap = 0;

xH = vec(width{1}(1, isAboveThreshold));
yH = vec(width{2}(1, isAboveThreshold));
xV = vec(width{1}(2, isAboveThreshold));
yV = vec(width{2}(2, isAboveThreshold));
tmpInds = flyInds(isAboveThreshold);
nValidFlies = numel(unique(tmpInds));
iAx = 1;


clear hSubAx
mSize = 15;
hFig = createPrintFig(25 * [1 1/3]);
for iFly = unique(tmpInds)   
    hSubAx(iAx, 1) = subplot(2, nValidFlies, iAx);
    x = xH(iFly == tmpInds);
    y = yH(iFly == tmpInds);
    scatter(xH, yH, mSize, 'filled', 'MarkerFaceAlpha', 0.1, 'MarkerFaceColor', 'k');
    hold on;
    hScatter = scatter(x, y, mSize);
    hLines = lsline;
    delete(hLines(2));
%     refline(1,0)
    if doBootstrap && numel(x) > 1
        [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
        legStr = sprintf('n = %d\ncorr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
                          numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
    else
        [corrCoef, p] = corr(x(:), y(:));
        legStr = ['n = ' num2str(numel(x)) ...
                  ', \n corr = ' num2str(corrCoef, 2) ', \n p = ' num2str(p, 2)];
    %     hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
    end
%     hLeg.Box = 'off';    
    title(sprintf(legStr));
    corrArray(iAx, 1) = corrCoef;

    hSubAx(iAx, 2) = subplot(2, nValidFlies, nValidFlies + iAx);
    x = xV(iFly == tmpInds);
    y = yV(iFly == tmpInds);
    scatter(xV, yV, mSize, 'filled', 'MarkerFaceAlpha', 0.1, 'MarkerFaceColor', 'k');
    hold on;
    hScatter = scatter(x, y, mSize);
    hLines = lsline;
    delete(hLines(2));
    if doBootstrap && numel(x) > 1
        [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
        legStr = sprintf('n = %d\ncorr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
                          numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
    else
        [corrCoef, p] = corr(x(:), y(:));
        legStr = ['n = ' num2str(numel(x)) ...
                  '\n corr = ' num2str(corrCoef, 2) '\n p = ' num2str(p, 2)];
    %     hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
    end
    title(sprintf(legStr))
    corrArray(iAx, 2) = corrCoef;
    iAx = iAx + 1;
end

axis(hSubAx, 'square');
linkaxes(hSubAx, 'xy')
clear hXLabel hYLabel
hXLabel(1) = xlabel(hSubAx(1, 1), 'Tm9 fwhm hor. (deg)');
hYLabel(1) = ylabel(hSubAx(1, 1), 'T5 fwhm hor. (deg)');
hXLabel(2) = xlabel(hSubAx(1, 2), 'Tm9 fwhm vert. (deg)');
hYLabel(2) = ylabel(hSubAx(1, 2), 'T5 fwhm vert. (deg)');

[hXLabel(1).FontSize] = deal(10);
[hYLabel(1).FontSize] = deal(10);
% hLeg.FontSize = 8;
[hSubAx.FontSize] = deal(8);
[hSubAx.Color] = deal('none');
arrayfun(@prettifyAxes, [hSubAx]); 

print(gcf, [figDir 'rfFwhmScatterPerFly-Off-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');

%% average values within fly and plot
% fly means
xHFlyMeans = accumarray(tmpInds', xH, [], @mean);
yHFlyMeans = accumarray(tmpInds', yH, [], @mean);
xHFlyStds = accumarray(tmpInds', xH, [], @std);
yHFlyStds = accumarray(tmpInds', yH, [], @std);
xHFlyMins = accumarray(tmpInds', xH, [], @min);
yHFlyMins = accumarray(tmpInds', yH, [], @min);
xHFlyMaxs = accumarray(tmpInds', xH, [], @max);
yHFlyMaxs = accumarray(tmpInds', yH, [], @max);

[x, y] = deal(xHFlyMeans(xHFlyMeans~=0), yHFlyMeans(xHFlyMeans~=0));
[xStd, yStd] = deal(xHFlyStds(xHFlyMeans~=0), yHFlyStds(xHFlyMeans~=0));
[xMin, yMin] = deal(xHFlyMins(xHFlyMeans~=0), yHFlyMins(xHFlyMeans~=0));
[xMax, yMax] = deal(xHFlyMaxs(xHFlyMeans~=0), yHFlyMaxs(xHFlyMeans~=0));

hFig = createPrintFig(7 * [1 1]);
hScatter = errorbar(x, y, xMin - x, xMax - x, yMin - y, yMax - y, ...
                    'LineStyle', 'none', 'CapSize', 0, 'LineWidth', 1, ...
                    'Color', [1 1 1] * 0.4);
hold on
hScatter = scatter(x, y, 25, 'filled', 'MarkerFaceColor', [1 1 1] * 0.20);
lsline
if doBootstrap && numel(x) > 1
    [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
    legStr = sprintf('n = %d, corr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
                      numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
else
    [corrCoef, p] = corr(x(:), y(:));
    legStr = ['n = ' num2str(numel(x)) ...
              ', \n corr = ' num2str(corrCoef, 2) ', \n p = ' num2str(p, 2)];
end
title(sprintf(legStr));
xlabel('Tm9 fwhm hor. (deg)');
ylabel('T5 fwhm hor. (deg)');
% hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
% hLeg.Box = 'off';   
prettifyAxes(gca);

print(gcf, [figDir 'rfFwhmScatterFlyMeans-XOff-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');
%% Same fo vertical
xHFlyMeans = accumarray(tmpInds', xV, [], @mean);
yHFlyMeans = accumarray(tmpInds', yV, [], @mean);
xHFlyStds = accumarray(tmpInds', xV, [], @std);
yHFlyStds = accumarray(tmpInds', yV, [], @std);
xHFlyMins = accumarray(tmpInds', xV, [], @min);
yHFlyMins = accumarray(tmpInds', yV, [], @min);
xHFlyMaxs = accumarray(tmpInds', xV, [], @max);
yHFlyMaxs = accumarray(tmpInds', yV, [], @max);

[x, y] = deal(xHFlyMeans(xHFlyMeans~=0), yHFlyMeans(xHFlyMeans~=0));
[xStd, yStd] = deal(xHFlyStds(xHFlyMeans~=0), yHFlyStds(xHFlyMeans~=0));
[xMin, yMin] = deal(xHFlyMins(xHFlyMeans~=0), yHFlyMins(xHFlyMeans~=0));
[xMax, yMax] = deal(xHFlyMaxs(xHFlyMeans~=0), yHFlyMaxs(xHFlyMeans~=0));

hFig = createPrintFig(7 * [1 1]);
hScatter = errorbar(x, y, xMin - x, xMax - x, yMin - y, yMax - y, ...
                    'LineStyle', 'none', 'CapSize', 0, 'LineWidth', 1, ...
                    'Color', [1 1 1] * 0.4);
% hScatter = errorbar(x, y, -xStd, xStd, -yStd, yStd, ...
%                     'LineStyle', 'none', 'CapSize', 0, 'LineWidth', 1, ...
%                     'Color', [1 1 1] * 0.2);
hold on
hScatter = scatter(x, y, 25, 'filled', 'MarkerFaceColor', [1 1 1] * 0.20);
lsline
if doBootstrap && numel(x) > 1
    [bootstrapCorrCI, bootstrapCorrs] = bootci(10000, @corr, x(:), y(:));
    legStr = sprintf('n = %d, corr = %.2f\n0.95-CI = [%.2f, %.2f]', ...
                      numel(x), mean(bootstrapCorrs), bootstrapCorrCI);
else
    [corrCoef, p] = corr(x(:), y(:));
    legStr = ['n = ' num2str(numel(x)) ...
              ', \n corr = ' num2str(corrCoef, 2) ', \n p = ' num2str(p, 2)];
end
title(sprintf(legStr));
xlabel('Tm9 fwhm vert. (deg)');
ylabel('T5 fwhm vert. (deg)');
% hLeg = legend(hScatter, legStr, 'Location', 'northoutside');
% hLeg.Box = 'off';   
prettifyAxes(gca);

print(gcf, [figDir 'rfFwhmScatterFlyMeans-YOff-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');

%% Check RF position distribution

xH = vec(pos{1}(1, isAboveThreshold));
yH = vec(pos{2}(1, isAboveThreshold));
xV = vec(pos{1}(2, isAboveThreshold));
yV = vec(pos{2}(2, isAboveThreshold));
tmpInds = flyInds(isAboveThreshold);
nValidFlies = numel(unique(tmpInds));
iAx = 1;

hFig = createPrintFig(7 * [ 1 2/3]);

hHAx(1) = subplot(1, 2, 1);
hHist(1) = histogram(xH - yH);
hHAx(2) = subplot(1, 2, 2);
hHist(2) = histogram(xV - yV);

[hHist.EdgeColor] = deal('w');
[hHist.FaceColor] = deal([1 1 1] * 0.0);

title(hHAx(1), {'Hor. RF Pos.', '(Tm9 - T5)'});
title(hHAx(2), {'Hor. RF Pos.', '(Tm9 - T5)'});

xlabel(hHAx(1), 'RF Pos. Diff. (deg)');

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

print(gcf, [figDir 'rfPosDiffScatter-Off-rSquare-' num2str(thresholdRSquare) '.pdf'], '-dpdf', '-r300');

%% Extract the parameters for the analysis of spatial RFs.
% clear amp pos width rsquare responseIndex
tuningMethod = 'median';
StimTableTmp = StimTableNeuronStim{2};
tuningInd = strcmp(tuningMethod, StimTableTmp.nonParamTuning{1}(1).tuningMethods);
jBar = 1;
for iStim = 2:5
    dummyTable = StimTableNeuronStim{1, iStim};
    for jChannel = 1: 2
        responseIndex{jChannel}(jBar, :) = cell2mat(cellfun(@(x, y) x(jChannel, 1: end - 1*0), dummyTable.responseIndex, 'uni', 0)');        
        % Hard coded getting rid of background ROI, but this should also be
        % imporved to handle the roiMeta field to discard invalid ROIs
        % (rois with NaNs)
        tuningCurves{jChannel}(jBar, :) = cellfun(@(x) squeeze(x(jChannel).tc{tuningInd}(:, 1: end - 1, :)), ...
                                                  dummyTable.nonParamTuning, 'uni', 0);
    end
    jBar = jBar + 1;
end
%% Plot tuning curves.
for iStim = 1: 4
    tm9AllTCs{iStim} = cat(1, tuningCurves{1}{iStim, :});
    t5AllTCs{iStim} = cat(1, tuningCurves{2}{iStim, :});
end
% Align to fit
% pos is already original parameter b1 as 2*(b1 - 1)

lagsCell = cellfun(@(x, y) round(x/2+1 - size(y, 2) / 2), ...
                   mat2cell(pos{1}, [1 1 1 1])', tm9AllTCs, 'uni', 0);
Tm9BarRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), tm9AllTCs, lagsCell, 'uni', 0);

lagsCell = cellfun(@(x, y) round(x/2+1 - size(y, 2) / 2), ...
                   mat2cell(pos{2}, [1 1 1 1])', t5AllTCs, 'uni', 0);
T5BarRFs = cellfun(@(x, y) alignRFsToAbsMax(x, 0, y'), t5AllTCs, lagsCell, 'uni', 0);


%% Blank epoch was already deleted, now delete vertical bars on edges of screen.
Tm9BarRFs{2}(:, end-1: end) = [];
Tm9BarRFs{4}(:, end-1: end) = [];
Tm9BarRFs{2}(:, 1: 2) = [];
Tm9BarRFs{4}(:, 1: 2) = [];

T5BarRFs{2}(:, end-1: end) = [];
T5BarRFs{4}(:, end-1: end) = [];
T5BarRFs{2}(:, 1: 2) = [];
T5BarRFs{4}(:, 1: 2) = [];
%% Plot only analyzed OFF RFs.
minRespInd = min(responseIndexOff);
thresholdRSquare = 0.25;
minRSquare = min(minOffRSquare);
isAboveThreshold = minRSquare > thresholdRSquare & minRespInd > 0.3;
Tm9BarRFs = cellfun(@(x) x(isAboveThreshold, :), Tm9BarRFs, 'uni', 0);
T5BarRFs = cellfun(@(x) x(isAboveThreshold, :), T5BarRFs, '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 1/4]);
hSubAx = plotLasagnaTracesFromCell(gcf, [Tm9BarRFs, Tm9BarRFs(1)], brewermap(6, 'Paired'));
colorbar
arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx([1, 3], 1));
arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'off'), hSubAx(2: end, 2));
arrayfun(@prettifyAxes, hSubAx)
arrayfun(@offsetAxes, hSubAx)
arrayfun(@(x) set(x, 'Color', 'none'), hSubAx);

% set(hSubAx, 'Visible', 'off')
linkaxes(hSubAx(1: 2, 1), 'y')
linkaxes(hSubAx(3: 4, 1), 'y')
hFig.Color = 'none';
figFileName = ['Tm9-RFcurves-OFF-XY-lasagna'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%%
hFig = createPrintFig(15 * [ 1 1/4]);
hSubAx = plotLasagnaTracesFromCell(gcf, [T5BarRFs, T5BarRFs(1)], brewermap(6, 'Paired'));
colorbar
arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx([1, 3], 1));
arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'off'), hSubAx(2: end, 2));
arrayfun(@prettifyAxes, hSubAx)
arrayfun(@offsetAxes, hSubAx)
arrayfun(@(x) set(x, 'Color', 'none'), hSubAx);

% set(hSubAx, 'Visible', 'off')
linkaxes(hSubAx(1: 2, 1), 'y')
linkaxes(hSubAx(3: 4, 1), 'y')
hFig.Color = 'none';
figFileName = ['T5-RFcurves-OFF-XY-lasagna'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')

%% Export data for figure 6b-e

receptiveFields = [Tm9BarRFs; T5BarRFs];
dataFileName = [figDir filesep 'Fig6.xls'];
cellTypeStr = {'Tm9', 'T5'};
sheetNamesCell = {'bc', 'de'};

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 = ['Fig6' 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
%% Tests if T5 ON response is 0 for reviewer #2

xH = vec(amp{1}(3, isAboveThreshold));
yH = vec(amp{2}(3, isAboveThreshold));
xV = vec(amp{1}(4, isAboveThreshold));
yV = vec(amp{2}(4, isAboveThreshold));
% use minimum for negative ON response.
minPeaks = cellfun(@(x) min(x, [], 2), T5BarRFs(3: 4), 'uni', 0);
tcStd = cellfun(@(x) std(x, [], 2), T5BarRFs(3: 4), 'uni', 0);
peakStdRatio = cellfun(@(x, y) x ./ y, minPeaks, tcStd, 'uni', 0);

% figure
% for iAx = 1: 2
%     subplot(1, 2, iAx);
%     histogram(peakStdRatio{iAx})
% end
 % check if the fitted amplitude is higher than the tuning curve values
 % (noise), signrank uses x - y, so tc - amp shuould be positive for a
 % smaller amplitude (larger negative peak), so we need right tail test.
[pValCell, isSignificant] = cellfun(@(x, y) arrayfun(@(z) signrank(x(z, :), y(z), 'tail', 'right'), 1: numel(y)), ...
                   T5BarRFs(3: 4), {yH, yV}, 'uni', 0);

proportionSignificant = cellfun(@(x) sum(x) / numel(x), isSignificant);

clear hAx
hFig = createPrintFig(15 * [3/4 1]);
hAx(1) = subplot(3, 2, 1);
histogram(yH);
hAx(2) = subplot(3, 2, 2);
histogram(yV);


for iAx = 1: 2
    hAx(2 + iAx) = subplot(3, 2, 2  * iAx + 1);
    histogram(pValCell{iAx}, 'NumBins', 100);
end
hAx(5) = subplot(3, 2, [4, 6]);
hBar = bar([proportionSignificant; 1 - proportionSignificant]', 'stacked');
hLeg = legend(hAx(5), {'p < 0.05', 'p >= 0.05'}, 'Location', 'best');
hAx(5).XTickLabels = {'Hor.',  'Vert.'};
ylim([0 1])

suptitle('T5 Responses to ON bars')

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

xLabelCell = {'Amplitude Hor. (\Delta F / F_0)', ...
              'Amplitude Vert. (\Delta F / F_0)', ...
              '', 'p-Value', 'Bar orientation'};
yLabelCell = {'No. ROIs', '', 'No. ROIs', 'No. ROIs', 'Proportion of ROIs'};

arrayfun(@(x, y) xlabel(x, y), hAx, xLabelCell);
arrayfun(@(x, y) ylabel(x, y), hAx, yLabelCell);

figFileName = ['T5-ON-Amp-XY-stats'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')

%%

function tcRoiCell = getTuningCurvePerROI(tc, tuningMethods, selectedTuningMethod)
%% Expose single Roi tuning maps.
% The epochs are varying first in the spacing then in the temporal
% frequency, so tc consists of 1 x nTFreqs and each element is of shape
% 1 x nRois x nSFreqs.
% Here we first extract per Roi the nSFreqs tuning, as column vector, and
% concatenate the nTFreqs column vectors arising this way to get a tuning
% matrix of nSFreqs x nTFreqs.
% tc1F = tc(:, strcmp(tuningMethods, '1F_amplitude'));
tc = tc(:, strcmp(tuningMethods, selectedTuningMethod));
tcRoiCell = cell(1, size(tc{1}, 2));

for iRoi = 1: size(tc{1}, 2)
    tcRoiCell{iRoi} = cell2mat(cellfun(@(x) squeeze(x(:, iRoi, :)), tc, 'uni', 0)');
end

if  0
% Try to reshape the tc into one tc array for each roi. is equivalent but
% maybe less clear. First concatenate along stim dimension to get an array
% of dims 1 x nRois x (nSFreqs x nTFreqs), squeeze it to get an array of
% dims nRois x (nSFreqs x nTFreqs), and reshape it into an array of dims
% nRois x nSFreqs x nTFreqs. Now we slice it along the first dimension into
% a cell array of nRois x 1, each cell being nSFreqs x nTFreqs.
tcRoiCell = mat2cell(reshape(squeeze(cat(3, tc{:})), sizeTC(2), sizeTC(3), []), ...
                    ones(1, sizeTC(2)), stimSizeTC(1), stimSizeTC(2));
end                

end