ccc
%% Define search path.
parentDir = fileparts(fileparts(mfilename('fullpath')));
saveFolder = [parentDir filesep 'data'];
load([saveFolder filesep 'Fig3-data_Tm9GCaMP6f-Tm4jRGECO1a_ONOFFPaperNew']);

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

%% Genereta figure 3b
onOffPerCell{1} = cell2mat(cellfun(@(x) cell2mat(cellfun(@(y) y, x(:, 1), 'uni', 0)), onoffData, 'Uni', 0)')';
onOffPerCell{2} = cell2mat(cellfun(@(x) cell2mat(cellfun(@(y) y, x(:, 2), 'uni', 0)), onoffData, 'Uni', 0)')';
timePoints = (0: size(onOffPerCell{1}, 2)-1) /10;

plotONOFFTracesFromCell(onOffPerCell, {'Tm4', 'Tm9'}, [figDir filesep], brewermap(2, 'PrGn'));

%% Save data from plot.
dataFileName = [figDir filesep 'Fig3.xls'];
DataTable = createRFsTable(onOffPerCell{1}, 'Time_seconds_Tm4');
DataTable(:, 1).Time_seconds_Tm4 = timePoints(:);
writetable(DataTable, dataFileName, 'Sheet', 'Fig3b', 'Range', 'A1');

DataTable = createRFsTable(onOffPerCell{2}, 'Time_seconds_Tm9');
DataTable(:, 1).Time_seconds_Tm9 = timePoints(:);
writetable(DataTable, dataFileName, 'Sheet', 'Fig3b', 'Range', ['A' num2str(size(onOffPerCell{1}, 2) + 3)])

%% Reanalyze data of dual imaging.
ccc

parentDir = fileparts(fileparts(mfilename('fullpath')));
saveFolder = [parentDir filesep 'data'];
load([saveFolder filesep 'Fig3-data_Tm9GCaMP6f-Tm4jRGECO1a_All']);

%%
figDir = [parentDir filesep 'figures' filesep 'Fig3' filesep];
if ~exist(figDir, 'dir'); mkdir(figDir); end
%%
[nFlies, nLayers, nStims] = size(dataEdges);
% New functions use cell arrays.
dataEdges = mat2cell(dataEdges, ones(1, nFlies), ones(1, nLayers), ones(1, nStims));
% Get tuning curves.
tuningStat = 'maxAbs';
% tuningStat = 'median';
barRFs = collectRFs(dataEdges, tuningStat);

% Remove blanck epoch.
barRFs = cellfun(@(x) x(:, 2: end), barRFs, 'Uni', 0);
barRFs([2, 4], :) = cellfun(@(x) x(:, 3: end - 2), barRFs([2, 4],  :), 'Uni', 0);

%%

fitStructOff = cellfun(@(x) fitGaussian1(1: size(x, 2), x, 1), barRFs(1:2, :), 'Uni', 0);
fitStructOn = cellfun(@(x) fitGaussian1(1: size(x, 2), x, -1), barRFs(3:4, :), 'Uni', 0);
fitStructBars = [fitStructOff; fitStructOn];

%% Extract the parameters for the analysis of spatial RFs.
for iChannel = 1: 2
    amp{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.a1, x, 'uni', 1), ...
                                        fitStructBars(:, iChannel), 'uni', 0));
    pos{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.b1, x, 'uni', 1), ...
                                        fitStructBars(:, iChannel), 'uni', 0));
    width{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.fit.c1, x, 'uni', 1), ...
                                          fitStructBars(:, iChannel), 'uni', 0));
    rsquare{iChannel} = cell2mat(cellfun(@(x) arrayfun(@(y) y.gof.rsquare, x, 'uni', 1), ...
                                            fitStructBars(:, iChannel), 'uni', 0));
end
pos = cellfun(@(x) 2 * (x - 1), pos, 'uni', 0);
width = cellfun(@(x) 2 * (2 * sqrt(log(2))) * x, width, 'uni', 0);

%%
% barRFs = cellfun(@(x) alignRFsToAbsMax(x, 0), barRFs, 'uni', 0);

% The screen is about 48-60 deg, we center the RFs at about the center 24
% deg, using the fit positions for higher presicion.
for iChannel = 1: 2
    lagsBars = round((pos{iChannel}(1: 4, :) - 25) / 2);
    for iBar = 1: size(barRFs, 1)
        barRFs{iBar, iChannel} = alignRFsToAbsMax(barRFs{iBar, iChannel}, 0, lagsBars(iBar, :)');
    end
end
%% Use all data for fig3c-d
splitSign = 1;
doSmoothing = 1;
rSquareCutoff = -20;
plotImage = 1;
hue = [];
useSubplots = 1;

for iChannel = 1: 2
    rSquareArray = rsquare{iChannel}';
    % Check good fit for all stim.
    [rfInds{1: 2, iChannel}] =  deal(all(rSquareArray(:, [1: 2]) >= rSquareCutoff, 2));
    [rfInds{3: 4, iChannel}] =  deal(all([rSquareArray(:, [1:2]) >= rSquareCutoff ... 
                                          rSquareArray(:, [3:4]) >= -20], 2));
end
%% Figure 3c, examples extracted from here.
hFig = createPrintFig(10 * [ 1 1]);
for iChannel = 1: 2
    subplot(2, 2, 2 * iChannel - 1);
    plotRFEnsemble(barRFs{2, iChannel}, barRFs{1, iChannel}, ...
                   splitSign, doSmoothing, rfInds{1, iChannel}, plotImage, hue, useSubplots)
    subplot(2, 2, 2 * iChannel);
    plotRFEnsemble(barRFs{4, iChannel}, barRFs{3, iChannel}, ...
                   splitSign, doSmoothing, rfInds{3, iChannel}, plotImage, hue, useSubplots)
end
figFileName = ['RFIm-ON-OFF-Noise-raw-ON-norm'];    
print(hFig, [figDir filesep figFileName '.pdf'], '-dpdf')
%% For the mean images in fig. 3d
filteredRFs = cellfun(@(x, y) x(y, :), barRFs, rfInds, 'uni', 0);

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)'), ...
%                   barRFs, 'uni', 0);
normRFs = cellfun(@(x) cell2mat(arrayfun(@(y) smooth(getRow(x, y))', 1: size(x, 1), 'uni', 0)'), ...
                  filteredRFs, 'uni', 0);
normRFs = cellfun(@(x) maxNorm(x,2), normRFs, 'uni', 0);

% Maximum of the mean tuning curve to normalize image accross stimuli.
maxPerStim = max(cellfun(@(x) max(abs(mean(x))), filteredRFs, 'uni', 1), [], 2);
maxPerPolarity = max(reshape(maxPerStim, [2 2]));
hFig = createPrintFig(15 * [ 1 1 / 3]);
for iRFType = 1:4
    subplot_tight(1, 4, iRFType);
%     imshow(getRFIm(mean(normRFs{2*iRFType}), mean(normRFs{2*iRFType - 1})))
    imshow(getRFIm(mean(filteredRFs{2*iRFType}), mean(filteredRFs{2*iRFType - 1}), ...
           maxPerPolarity(2 - mod(iRFType, 2))))
end
%%
figFileName = ['RFIm-ON-OFF-Noise-mean'];    
print(hFig, [figDir filesep figFileName '.pdf'], '-dpdf')

%% For the mean traces in fig 3d.
hFig = createPrintFig(15 * [ 1 2/3]);
labelsCell = {'Tm4 off x' 'Tm4 off y' 'Tm4 on x' 'Tm4 on y'...
             'Tm9 off x' 'Tm9 off y' 'Tm9 on x' 'Tm9 on y'};
rfsToPlot = filteredRFs; % or normRFs.
plotONOFFTracesFromCell(rfsToPlot, labelsCell, [figDir filesep], brewermap(8, 'Paired'))
hSubAx = plotLasagnaTracesFromCell(gcf, rfsToPlot, brewermap(8, 'Paired'));
ylim(hSubAx(:, 1), [-1 1]);
arrayfun(@(x) set(get(x, 'YAxis'), 'Visible', 'on'), hSubAx(1, :))
cellfun(@(x) caxis(ZeroCenteredBounds(x, [0 100])), rfsToPlot, 'uni', 0)
arrayfun(@prettifyAxes, hSubAx)
arrayfun(@offsetAxes, hSubAx)
set(hSubAx, 'Visible', 'off')
% figFileName = ['RFcurves-ON-OFF-Noise-XY-lasagna'];    
% print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%%
dataFileName = [figDir filesep 'Fig3.xls'];

DataTable = createRFsTable(rfsToPlot{1, 1}, 'Tm4_OFF_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A1')
DataTable = createRFsTable(rfsToPlot{2, 1}, 'Tm4_OFF_vertical_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A29')
DataTable = createRFsTable(rfsToPlot{3, 1}, 'Tm4_ON_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A56')
DataTable = createRFsTable(rfsToPlot{4, 1}, 'Tm4_ON_vertical_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A83')

DataTable = createRFsTable(rfsToPlot{1, 2}, 'Tm9_OFF_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A110')
DataTable = createRFsTable(rfsToPlot{2, 2}, 'Tm9_OFF_vertical_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A137')
DataTable = createRFsTable(rfsToPlot{3, 2}, 'Tm9_ON_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A164')
DataTable = createRFsTable(rfsToPlot{4, 2}, 'Tm9_ON_vertical_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A191')

%%
splitSign = 1;
doSmoothing = 1;
rSquareCutoff = 0.5;
plotImage = 1;
hue = [];
useSubplots = 1;
% for iChannel = 1: 2
%     rSquareArray = rsquare{iChannel}';
%     % Check good fit for all stim.
%     indsToPlot = all([rSquareArray(:, [1:2]) >= rSquareCutoff ... 
%                       rSquareArray(:, [3:4]) >= 0.1], 2);
%     rfInds{:, iChannel} = indsToPlot;
% end


for iChannel = 1: 2
    rSquareArray = rsquare{iChannel}';
    % Check good fit for all stim.
    [rfInds{1: 2, iChannel}] =  deal(all(rSquareArray(:, [1: 2]) >= rSquareCutoff, 2));
    [rfInds{3: 4, iChannel}] =  deal(all([rSquareArray(:, [1:2]) >= rSquareCutoff ... 
                                          rSquareArray(:, [3:4]) >= 0.2], 2));
end


%% Figure 3e
hFig = createPrintFig(15 * [1/2 1/2]);

meanWidth = {mean(width{1}(1: 2, rfInds{1, 1}))
            mean(width{2}(1: 2, rfInds{1, 2}))
            mean(width{2}(3: 4, rfInds{3, 2}))};

beanPlot(meanWidth, {'Tm4 off', 'Tm9 off', ' Tm9 on'}, flipud(brewermap(4, 'Set1')), 1, gca,'unbounded', 'vertical', false, true)
figFileName = ['RF-ON-OFF-Noise-Width-Bean'];    
print(hFig, [figDir filesep figFileName '.pdf'], '-dpdf')

%% Save data for Fig 3e
sheetName = 'Fig3e';
xlswrite(dataFileName, {'Tm4-FHWM-OFF-bars'}, sheetName, 'A1')
xlswrite(dataFileName, meanWidth{1}, sheetName, 'B1')
xlswrite(dataFileName, {'Tm9-FHWM-OFF-bars'}, sheetName, 'A2')
xlswrite(dataFileName, meanWidth{2}, sheetName, 'B2')
xlswrite(dataFileName, {'Tm9-FHWM-ON-bars'}, sheetName, 'A3')
xlswrite(dataFileName, meanWidth{3}, sheetName, 'B3')

%% For revision, reviewer 3 plot some raw traces.
% 
receptiveFieldTraces = collectRFTraces(dataEdges);

filteredRFTraces = cellfun(@(x, y) x(y, :), receptiveFieldTraces, rfInds, 'uni', 0);
rng(1001)
subsampleInds = randsample(min(cellfun(@(x) size(x, 1), receptiveFieldTraces(1, :))) - 1, 10);
filteredRFTraces = receptiveFieldTraces;

%% Attempt to plot traces correctly, this is working now!

deleteLastROI = 1;

% Define the array to compare to randomized recordings.
xOffStimFileName = 'StandingStripe_1s_XAxis_5degWide_2degSep_m1.0Con_rand_USEFRUSTUM.txt';
yOffStimFileName = 'StandingStripe_1s_YAxis_5degWide_2degSep_m1.0Con_rand_USEFRUSTUM.txt';
xOnStimFileName = 'StandingStripe_1s_XAxis_5degWide_2degSep_p1.0Con_rand_USEFRUSTUM.txt';
yOnStimFileName = 'StandingStripe_1s_YAxis_5degWide_2degSep_p1.0Con_rand_USEFRUSTUM.txt';
barStimFileNames = {xOffStimFileName, yOffStimFileName, xOnStimFileName, yOnStimFileName};

%%
% Take rows as ROIs.
hFunYOffsetTraces = @(traces, yOffset) ...
                     bsxfun(@plus, traces, yOffset * (1: size(traces, 1))');
%%
nChannels = numel(dataEdges{1}.roiTCsInterp);
% receptiveFieldsAll = cell(4, nChannels);
cellCounter = zeros(2, 4);
hFig = createPrintFig(30 * [1 1 / 2]);
axInds(1, :) = [1, 2, 5, 6];
axInds(2, :) = [3, 4, 7, 8];
for iFly = 1: size(dataEdges, 1)
%     flyPathInd = iFly;
    for jStim = 2:5
        % Find index of matching stimulus.
        rfInd = find(strcmp(dataEdges{iFly, 1, jStim}.stimulusName, barStimFileNames));
        if ~isempty(rfInd)
            for kChannel = 1: nChannels
                dataTmp = dataEdges{iFly, 1, jStim};
                nEpochs = numel(dataTmp.stimFrameInds); 
                nCells = size(dataTmp.roiTCs{kChannel}, 2);
                roiTCs = dataTmp.roiTCs{kChannel};
                % Add a fake epoch at the end.
                dataTmp.stimFrameInds{end + 1} = dataTmp.stimFrameInds{end}(end) + ...
                                                 (1: numel(dataTmp.stimFrameInds{end}));
                roiTCs(:, :, dataTmp.stimFrameInds{end}) = nan;                             
                indsCellArray = [dataTmp.stimFrameInds(1: end - 1); 
                                 repmat(dataTmp.stimFrameInds(end), 1, nEpochs)];
                interleavedCellArray = reshape(indsCellArray, 1, []);
                plotIndsVec = cat(2, interleavedCellArray{:});
                roiTraces = squeeze(roiTCs(:, 1: end - 1 * deleteLastROI, plotIndsVec));
                times = dataTmp.frameDuration * (1: numel(plotIndsVec));
                hAx(kChannel, rfInd) = subplot(2, 4, axInds(kChannel, rfInd));
                plot(times, hFunYOffsetTraces(roiTraces + cellCounter(kChannel, rfInd), 1));
                hold on
                cellCounter(kChannel, rfInd) = cellCounter(kChannel, rfInd) + nCells;
                rawTracesStruct(iFly, rfInd, kChannel).times = times;
                rawTracesStruct(iFly, rfInd, kChannel).traces = roiTraces;
            end
        end
    end
end

%% This is figure S2
close all

hFig = createPrintFig(30 * [1 1 / 2]);
axInds(1, :) = [1, 2, 5, 6];
axInds(2, :) = [3, 4, 7, 8];
flyInd = 10; % Good examples 4, 6, 7, and 10
genotypeColors = [166,118,29; 27,158,119] / 255;
for iStim = 1: 4
    for jChannel = 1: 2
        hTracesAx(jChannel, iStim) = subplot(2, 4, axInds(jChannel, iStim));
        tmpTimes = rawTracesStruct(flyInd, iStim, jChannel).times;
        tmpTraces = rawTracesStruct(flyInd, iStim, jChannel).traces;
        plot(tmpTimes, hFunYOffsetTraces(tmpTraces, 0.7), ...
             'Color', genotypeColors(jChannel, :), 'LineWidth', 1);
        axis tight
        hLines = arrayfun(@(x) refline(0, x), 0.7 * (1: size(tmpTraces, 1)));
        [hLines.Color] = deal([1 1 1] * 0.5);
        [hLines.LineWidth] = deal(0.5);
    end
end


xlabel(hTracesAx(5), 'Time (s)')
ylabel(hTracesAx(5), '\DeltaF/F_0')

setFontForThesis(hTracesAx, gcf);
arrayfun(@prettifyAxes, hTracesAx)
arrayfun(@offsetAxes, hTracesAx)
figFileName = ['Fig3-suppl-meanTraces-sortedPositions'];    
print(hFig, [figDir figFileName '.pdf'], '-dpdf')
%%
DataTable = createRFsTable(rfsToPlot{1, 1}, 'Tm4_OFF_horizontal_position_deg');
writetable(DataTable, dataFileName, 'Sheet', 'Fig3cd', 'Range', 'A1')
%%
flyInd = 10; % Good examples 4, 6, 7, and 10
sheetName = 'SupFig3ac';
startRow = 1;
stimNameStrCell = {'Tm4_OFF_horizontal_bar_time_seconds', 'Tm4_OFF_vertical_bar_time_seconds', ...
                   'Tm4_ON_horizontal_bar_time_seconds', 'Tm4_ON_vertical_bar_time_seconds'};
for iStim = 1: 4
    DataTable = createRFsTable(rawTracesStruct(flyInd, iStim, 1).traces, stimNameStrCell{iStim});
    DataTable(:, 1).(stimNameStrCell{iStim}) = rawTracesStruct(flyInd, iStim, 1).times(:);
    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
    startRow = startRow + size(DataTable, 1) + 2;
end

sheetName = 'SupFig3bd';
startRow = 1;
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'};
for iStim = 1: 4
    DataTable = createRFsTable(rawTracesStruct(flyInd, iStim, 2).traces, stimNameStrCell{iStim});
    DataTable(:, 1).(stimNameStrCell{iStim}) = rawTracesStruct(flyInd, iStim, 2).times(:);
    writetable(DataTable, dataFileName, 'Sheet', sheetName, 'Range', ['A' num2str(startRow)]);
    startRow = startRow + size(DataTable, 1) + 2;
end