%% Calculate and plot summary data for hemolymph transfusion experiments

clear;

% Directory to look in when selecting data file
defaultDir = uigetdir();

dataSet = 'healthy'; % Specify either 'sens' for sensitized hemo experiments or 'healthy' for experiments with healthy recipients

[expName, expFile, expDir, irange, irangeby2, cutoff] = deal(cell(1,3));

switch dataSet
    case 'healthy'       
        [supName,expFolder,expName, expFile, expDir, irange, irangeby2, cutoff] = loadHealthyHemoTrxData;
    case 'sens'
        [supName,expFolder,expName, expFile, expDir, irange, irangeby2, cutoff] = loadSensHemoTrxData;
end

%% Settings

% Select which variable to plot
var = "Speed"; %% options are "Centroid" (y position) or "Speed" (total speed).

% Options
plotPairs = 0;
doBlind = 0; % option to plot pairs blind (0) or not (1)
drop_NS_XI = 0; % option to differentiate (plotAllInd) or exclude (plotMean) non-summiting flies whose donors died from infection on the day experiment was performed
metric = 'mean'; % metric for summary var plot (choices = mean or median)
numHours = 2; % upper x bound for summary plots
minWindow = 30; % window width (minutes) for smoothing speed plots
speedConvert = pixel2mm; % pull conversion from pixels to mm 

% Decide method for determining time of first movement
useCalledFirst = 1; % (1) if use manually called time, (0) for auto detection (first movement = first detected movement)

%% Figure out the length of longest experiment so all can be aligned to it
maxLength = 0;

for i=1:length(expName)
    % Load in .mat file with expmt struct and grab corresponding survival data
    load([expDir{i},expFile{i}]);
    if length(expmt.Centroid.data) > maxLength
        maxLength = length(expmt.Centroid.data);
    end
end

%% Compile all data
% Initialize empty variables
[allSummit, allNonSummit, allSumwalk, allNonSumWalk] = deal(NaN(96,maxLength));
[allSiRunning, allNsiRunning] = deal(NaN(length(expDir),maxLength));

baseSI = 1;
baseNSI = 1;

for i=1:length(expName)
    % Load in .mat file with expmt struct and grab corresponding survival data
    load([expDir{i},expFile{i}]);
    disp(expFile{i});
    expDate = returnExpDate(string(expFile{i}));

    % Automatically determine filename for survival Excel sheet
    survivalFile = strcat(expDir{i},expDate,"-survival_data.xlsx");

    % Load Excel sheet
    survival = readtable(survivalFile);
    
    %% Fill all NaN values so that areaBar will plot SE
    expmt.Centroid.data = fillmissing(expmt.Centroid.data,'previous','EndValues','nearest');
    expmt.Speed.data = fillmissing(expmt.Speed.data,'constant',0);

    % Grab first mvmts from survival wkst
    firstMoves = grabFirstMoves(expmt, survival, irange{i}, useCalledFirst);
    
    % Load in summit and nonsummit data for experiment i
    [summit, nonsummit, si, nsi] = getHemoTraces (survival, expmt, var, irange{i}, firstMoves, drop_NS_XI);
    [sumWalk, siRunning, nonSumWalk, nsiRunning] = findFracWalking(cutoff{i}, summit, nonsummit, si, nsi);
    
    % Compile all summit traces aligned by time of first movement, either
    % centroid or speed depending on option selected above
    allSummit(baseSI:baseSI+si-1,:) = horzcat(summit, NaN(si,maxLength-length(summit)));
    allNonSummit(baseNSI:baseNSI+nsi-1,:) = horzcat(nonsummit, NaN(nsi,maxLength-length(nonsummit)));
    
    % Compile walking status at each point across experiment (walking >
    % moving above chosen threshold)
    allSumWalk(baseSI:baseSI+si-1,:) = horzcat(sumWalk, NaN(si,maxLength-length(sumWalk)));
    allNonSumWalk(baseNSI:baseNSI+nsi-1,:) = horzcat(nonSumWalk, NaN(nsi,maxLength-length(nonSumWalk)));
    
    % Compile how many flies were being tracked at each frame across
    % experiment
    allSiRunning(i,:) = horzcat(siRunning,NaN(1,maxLength-length(siRunning)));
    allNsiRunning(i,:) = horzcat(nsiRunning,NaN(1,maxLength-length(nsiRunning)));
    
    baseSI = baseSI + si;
    baseNSI = baseNSI + nsi;
end

if var == "Speed"
    allSummit = allSummit * speedConvert * 3;
    allNonSummit = allNonSummit * speedConvert * 3;    
end

% Correct count # for summiting and nonsummiting
baseSI = baseSI - 1;
baseNSI = baseNSI - 1;

% Sum how many flies were being tracked at each point in time across
% experiment
siRunning = nansum(allSiRunning); 
nsiRunning = nansum(allNsiRunning); 

%% Plot summary metric (mean/median etc) +/- SE of each group using areaBar
% Only makes sense if you set var == speed
if var == "Speed"
    plotHemoSpeedHist(allSummit, allNonSummit, expFolder, var, supName);
    plotDistTraveled(allSummit, allNonSummit, siRunning, nsiRunning, metric, numHours, drop_NS_XI, expFolder, supName, minWindow*3*60); 
end

plotFracWalking(allSumWalk, allNonSumWalk, siRunning, nsiRunning, numHours, expFolder, var, supName);

plotMeanHemoVar(allSummit, allNonSummit, siRunning, nsiRunning, var, metric, numHours, drop_NS_XI, expFolder, supName, minWindow*3*60); 

%% P sweep 
% Run t-test between summiting and non-summiting recipients, integrating
% over desired interval of experiment
minFrame = 1;
nHours = 5/60:5/60:2;
pSweep = NaN(1,length(nHours));

for i=1:numel(nHours)    
    maxFrame = 10800 * nHours(i);
    [sumDist] = integrateHemoSpeed(allSummit, minFrame, maxFrame);
    [nonSumDist] = integrateHemoSpeed(allNonSummit, minFrame, maxFrame);
    [h,p,CI,stats] = ttest2(sumDist,nonSumDist);
    pSweep(i) = p;
end

logpSweep = -1*log10(pSweep);
cutoff = -1*log10(0.05);

figure;
hold on;
scatter(nHours(logpSweep<cutoff),logpSweep(logpSweep<cutoff),15,'MarkerFaceColor',[0.5 0.5 0.5],'MarkerEdgeColor','none');
scatter(nHours(logpSweep>=cutoff),logpSweep(logpSweep>=cutoff),20,'k*');
plot(nHours,repmat(cutoff,1,length(nHours)));
ylabel('p-value');
xlabel('Integration window (hours)');
title(supName);
set(gca,'xtick',0:0.5:2,'TickLength',[0 0]);
ylim([0 1.8]);
set(gcf,'PaperUnits','inches','PaperPosition',[0 0 7.5 6]);
outName = strcat(expFolder,'pSweep');
print(strcat(outName,'.pdf'),'-dpdf'); 
print(strcat(outName,'.png'),'-dpng'); 
close();