% plotMultiDataStruct.m
%
% Purpose-plot several experiments (from data structs) on the same axes for
% easy comparison. Use data aligned by local time.

clear
tic;

%% Load data structs

%defaultDir = "D:\Harvard\Manuscripts\Summit paper\Figures\Figure 2\";
defaultDir = "D:\Harvard\Manuscripts\Summit paper\Figures\Fig 2 Supp\";

expChoice = 'DN1p'; % see switch logic below for options

% Temperature file to use to establish *common* color scale for plot
tempRef = readtable('D:\Harvard\Summit behavior\Behavior data\SUM6\2022\03-28-2022-15-53-55__Circadian_x907F1s-BoardC15_MF_None_1-128_Day5\CNE behavior room_export_202203311027_manually_adjusted.csv');

switch expChoice
    case 'PICA'
        structFiles = ["D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 6 - TrpA1 R19G10 in TM6B bg (2022-4-3)\R19G10_TrpA1(Tb+)_MixedSex_last12 hours.mat",...
            "D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 6 - TrpA1 R19G10 in TM6B bg (2022-4-3)\TrpA1_+;TM6B_+_MixedSex_last12 hours.mat"];
        tempLog = readtable('D:\Harvard\Summit behavior\Behavior data\SUM6\2022\03-28-2022-15-53-55__Circadian_x907F1s-BoardC15_MF_None_1-128_Day5\CNE behavior room_export_202203311027_manually_adjusted.csv');
        h=1; % "Rx" that flies are denoted as in data struct
        %peakSpeed = [0.918,0.1644];
        %preSpeed = [0.1583,0.09342]
    case 'DN1p'
        structFiles = ["D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 4 - R19G10, tutl, Clk4.1, R57C10\Clk4_1_TrpA1 (Sna[sco])_MixedSex_last12 hours.mat",...
            "D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 4 - R19G10, tutl, Clk4.1, R57C10\Sna[sco]_UAS-TrpA1; TM6B_+_MixedSex_last12 hours.mat"];
        tempLog = readtable('D:\Harvard\Summit behavior\Behavior data\SUM5\2021\07-25-2021-18-39-47__Circadian_x747F1s-BoardC17_MF_TempShift_1-128_Day1\CNE behavior room_manuallyadjusted.csv');
        h=3; % "Rx" that flies are denoted as in data struct
        %peakSpeed = [1.129,0.1402]; % determined by selecting points on plot in
        %debugging mode
        %preSpeed = [0.3406,0.1137]
    case 'PICA-F'
        structFiles = ["D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 6 - TrpA1 R19G10 in TM6B bg (2022-4-3)\R19G10_TrpA1(Tb+)_Female_last12 hours.mat",...
            "D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 6 - TrpA1 R19G10 in TM6B bg (2022-4-3)\TrpA1_+;TM6B_+_Female_last12 hours.mat"];
        tempLog = readtable('D:\Harvard\Summit behavior\Behavior data\SUM6\2022\03-28-2022-15-53-55__Circadian_x907F1s-BoardC15_MF_None_1-128_Day5\CNE behavior room_export_202203311027_manually_adjusted.csv');
        h=1; % "Rx" that flies are denoted as in data struct
    case 'PICA-M'
        structFiles = ["D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 6 - TrpA1 R19G10 in TM6B bg (2022-4-3)\R19G10_TrpA1(Tb+)_Male_last12 hours.mat",...
            "D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 6 - TrpA1 R19G10 in TM6B bg (2022-4-3)\TrpA1_+;TM6B_+_Male_last12 hours.mat"];
        tempLog = readtable('D:\Harvard\Summit behavior\Behavior data\SUM6\2022\03-28-2022-15-53-55__Circadian_x907F1s-BoardC15_MF_None_1-128_Day5\CNE behavior room_export_202203311027_manually_adjusted.csv');
        h=1; % "Rx" that flies are denoted as in data struct
    case 'DN1p-F'
        structFiles = ["D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 4 - R19G10, tutl, Clk4.1, R57C10\By sex\Clk4_1_TrpA1 (Sna[sco])_Female_last12 hours.mat",...
            "D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 4 - R19G10, tutl, Clk4.1, R57C10\By sex\Sna[sco]_UAS-TrpA1; TM6B_+_Female_last12 hours.mat"];
        tempLog = readtable('D:\Harvard\Summit behavior\Behavior data\SUM6\2022\03-28-2022-15-53-55__Circadian_x907F1s-BoardC15_MF_None_1-128_Day5\CNE behavior room_export_202203311027_manually_adjusted.csv');
        h=3; % "Rx" that flies are denoted as in data struct
    case 'DN1p-M'
        structFiles = ["D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 4 - R19G10, tutl, Clk4.1, R57C10\By sex\Clk4_1_TrpA1 (Sna[sco])_Male_last12 hours.mat",...
            "D:\Harvard\Summit behavior\Analysis\Screen\Boxes v3\TrpA1\Exp 4 - R19G10, tutl, Clk4.1, R57C10\By sex\Sna[sco]_UAS-TrpA1; TM6B_+_Male_last12 hours.mat"];
        tempLog = readtable('D:\Harvard\Summit behavior\Behavior data\SUM6\2022\03-28-2022-15-53-55__Circadian_x907F1s-BoardC15_MF_None_1-128_Day5\CNE behavior room_export_202203311027_manually_adjusted.csv');
        h=3; % "Rx" that flies are denoted as in data struct        
end
            
meanWindow = 10800; % # of frames over which to take average for smoothing

% Graphing options
C = [0 0 0; 0 0 0]; % Struct colors
lineStyle = {'-','--'};
xInt = 1; 
legFont = 12;                    % size of legend font, delta AUC shown on graph
tickFont = 12;
labelFont = 14;
titleFont = 10;
alpha = 0.1;
hrInc = 2; % number of hours to increment on x-axis

%% Determine where to begin aligning data, and which ZT0 to align to
minStart = 13;
refZT0 = 25; % by default, align to DST (EDT)

%% Pull temp data from log
tempHr = tempLog{:,'Hour'};
tempC = tempLog{:,'Temperature_Celsius'};
tempStart = 19;
startRow = find(tempHr == tempStart); % find temp at minStart in temp log data
tempRamp = repelem(tempC(startRow:length(tempC)),180);
tempRamp = vertcat(NaN(180*(tempStart-minStart)*60,1),tempRamp);
tempRamp = downsample(tempRamp, 100);

% Establish color scale to represent temperature range
tempRefC = tempRef{:,'Temperature_Celsius'};
tempRange = min(tempRefC):0.1:max(tempRefC);
cScale = jet(numel(tempRange+10));
tempAlpha = 0.8; %transparency of temperature color gradient

%%
[timeStamp,dateStamp] = Stamps;

speedConv = speedToMmPerS;

figure;
hold on;

keys = [];

for a=1:numel(structFiles)
    try
        load(structFiles(a), 'data');
    catch
        disp(strcat('No data struct in file ',num2str(a),': ',structFiles(a)));
    end
    
    color = C(a,:);
    
    if data.params.zt0 ~= refZT0
        missingFrames = zeros(1,round(3600*3*(data.startTime - minStart + 1)));
    else
        missingFrames = zeros(1,round(3600*3*(data.startTime - minStart)));
    end
    
    frmsPerHr = (3600*data.params.fr)/data.params.downSamp;
    dataSex = split(data.filePre,'-');
    dataSex = dataSex(length(dataSex));
        
    % If you're on first data struct, plot temperature as color gradient 
    if a == 1
        for j=1:length(tempRamp)
            [n,idx] = min(abs(tempRamp(j)-tempRange));
            area([j j+1],[5 5],'FaceColor',cScale(idx,:),'EdgeColor','none','FaceAlpha',tempAlpha,'HandleVisibility','off');
        end
    end       
    
    speedData = horzcat(missingFrames,nanmean(transpose(data.(data.rXes{h}).speed.allMatrix)));
    speedData = downsample(smoothdata(speedData,'movmean',meanWindow),data.params.downSamp)*speedConv*data.params.fr;
    speedSe = horzcat(missingFrames,(nanstd(transpose(data.(data.rXes{h}).speed.allMatrix))/sqrt(size(data.(data.rXes{h}).speed.allMatrix,2))));
    speedSe = downsample(smoothdata(speedSe,'movmean',meanWindow),data.params.downSamp)*speedConv*data.params.fr;

    % Replace NaN values with 0 to ensure error bars are plotted with
    % areaBar
    speedData(isnan(speedData)) = 0;
    speedSe(isnan(speedSe)) = 0;        

    try  
        areaBar(1:length(speedData),speedData,speedSe,color,color,lineStyle{a},alpha);
        keys = [keys,strcat(data.geno,'-',dataSex,' (N=',string(size(data.(data.rXes{h}).ypos.toDmatrix,2)),')')];
    catch
        keys = [keys,strcat(data.geno,'-',dataSex,' (N=0)')];
    end
end

ylabel('Speed (mm/s)','fontsize',labelFont);
ylim([0 2]);
xlabel('Time (ZT)','fontsize',labelFont);
xtickangle(0);        
set(gca,'xtick',1:hrInc*(3600*data.params.fr)/data.params.downSamp:ceil((2)*24*3600*data.params.fr/data.params.downSamp));
set(gca,'xticklabel',minStart-24-1:hrInc:(2)*24,'fontsize',tickFont,'TickLength',[0 0]);

zt0Frame = (refZT0-minStart)*3600*data.params.fr/data.params.downSamp;
zt12Frame = zt0Frame + (12*3600*data.params.fr/data.params.downSamp);

xBuffer = 2; % # hours to show on either side of x-axis
xlim([zt0Frame-(xBuffer*3600*data.params.fr/data.params.downSamp) zt12Frame+(xBuffer*3600*data.params.fr/data.params.downSamp)]);

legend(keys,'Location','best');

set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 3]);
%set(gcf,'PaperUnits','inches','PaperPosition',[0 0 6 4.5]);
print(strcat(defaultDir,expChoice,'-TrpA1.png'),'-dpng');
print(strcat(defaultDir,expChoice,'-TrpA1.pdf'),'-dpdf');
close();

% figure;
% hold on;
% for i=1:numel(tempRange)
%     scatter(1,tempRange(i),'s','MarkerFaceColor',cScale(i,:),'MarkerEdgeColor','none','MarkerFaceAlpha',tempAlpha);
% end
% xlim([0.9 1.1]);
% ylim([tempRange(1) tempRange(numel(tempRange))]);
% set(gca,'TickLength',[0 0 ],'ytick',tempRange(1):0.4:tempRange(numel(tempRange)));
% set(gcf,'PaperUnits','inches','PaperPosition',[0 0 1 3]);
% print(strcat(defaultDir,expChoice,'_Color bar.png'),'-dpng');
% print(strcat(defaultDir,expChoice,'_Color bar.pdf'),'-dpdf');
% close();
