function evokedActivityAnalyses(flyData)
%% Settings

    % Protocol
    evokedActivityWindows = [3601:4800 8401:9600 13201:14400];
    fliesPerOdor = 3;
    stimulationTrains = 3;
    repetitionsPerTrain = 8;
    odorTime = 51:65;           % 1.5 s after odor response onset
    odorNames = {
                 '1-pentanol';
                 '2-heptanone';
                 '3-octanol';
                 '4-methylcyclohexanol';
                 'amyl acetate';
                 'benzaldehyde';
                 'ethyl acetate';
                 'ethyl valerate';
                                        };

    % Glomeruli Per Fly
    glomeruliPerFly = cumsum(flyData.glomeruliPerFly);
    glomeruliPerOdor = 0;
    for i=1:length(odorNames)
        glomeruliPerOdor = cat(1,glomeruliPerOdor,glomeruliPerFly(fliesPerOdor*i));
    end
    glomeruliPerFly = cat(1,0,glomeruliPerFly);
    
    % Analyses
    correlationAnalysis = 1;
    EuclideanDistanceAnalysis = 1;
    glomerularResponsesAnalysis = 1;
    relativeResponseChangeAnalysis = 1;
    sparsenessAnalysis = 1;


%% Data

    odorTimeSeries = flyData.glomeruliDFFTimeSeries(:,evokedActivityWindows);
    odorTimeSeries = reshape(odorTimeSeries,size(odorTimeSeries,1),150,repetitionsPerTrain*stimulationTrains);
    
    
%% Odor correlations

if correlationAnalysis == 1

    % Correlation matrices per fly
    R = [];
    for i=1:length(glomeruliPerFly)-1
        x = odorTimeSeries(glomeruliPerFly(i)+1:glomeruliPerFly(i+1),:,:);
        x = squeeze(mean(x(:,odorTime,:),2));
        R = cat(3,R,corrcoef(x));
        clear('x')
    end
    figure, set(gcf,'position',[80 80 1150 720])
        for i=1:length(glomeruliPerFly)-1
            subplot(4,6,i)
                imagesc(R(:,:,i))
                colormap(jet), caxis([0.7 1])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                title(['Fly # ',num2str(i)])
        end
        
    % Correlation matrices - average per odor (3 flies per odor)
    figure, set(gcf,'position',[40 125 1375 635])
        idxMat = cat(3,[1:8; 1:8; 1:8; 9:16; 9:16; 9:16; 17:24; 17:24; 17:24],[1:8; 9:16; 17:24; 1:8; 9:16; 17:24; 1:8; 9:16; 17:24]);
        subPlotID1 = [1 2 3 17 18 19 33 34 35];
        subPlotID2 = [65 66 67 81 82 83 97 98 99];
        colorAxisLimits = [0.8 0.8 0.8 0.6 0.8 0.6 0.7 0.8];
        for i=1:length(odorNames)
            x = mean(R(:,:,(i*3)-2:i*3),3);
            for j=1:9
                if i<5
                    subplot(7,16,subPlotID1(j)+4*(i-1))
                else
                    subplot(7,16,subPlotID2(j)+4*(i-5))
                end
                imagesc(x(idxMat(j,:,1),idxMat(j,:,2))), colormap(jet), caxis([colorAxisLimits(i) 1])
                set(gca,'Box','off','TickDir','out','XTickLabel',[],'YTickLabel',[])
                if j==2
                    title(odorNames(i),'Fontsize',12,'Fontweight','b')
                end
            end
        end
    % Stats on odor correlations - grouped by odor (3 flies per odor)
    figure, set(gcf,'position',[110 75 1355 720])
    disp('     ')
    disp('     Linear Mixed-Effects Model on Odor Correlations (grouped by odor):')
    for i=1:8
        subplot(2,4,i)
            X = R(:,:,(1:3)+3*(i-1));
            Y = cat(3,ones(size(R,1),size(R,2)),2*ones(size(R,1),size(R,2)),3*ones(size(R,1),size(R,2)));
            matIdx = [1:8;9:16;17:24];
            x = [];
            y = [];
            for j=1:3
                a = X(matIdx(j,:),matIdx(j,:),:);
                b = Y(matIdx(j,:),matIdx(j,:),:);
                idx = tril(ones(size(a,1),size(a,2)),-1);
                a = reshape(a(logical(idx)),[],1);
                b = reshape(b(logical(idx)),[],1);
                x = cat(2,x,a);
                y = cat(2,y,b);
                clear('a','b')
            end
            matIdx = cat(3,[9:16; 17:24; 17:24],[1:8; 1:8; 9:16]);  % trains 1-2, trains 1-3, trains 2-3
            z = [];
            w = [];
            for j=1:3
                a = X(matIdx(j,:,1),matIdx(j,:,2),:);
                b = Y(matIdx(j,:,1),matIdx(j,:,2),:);
                a = reshape(a(logical(idx)),[],1);
                b = reshape(b(logical(idx)),[],1);
                z = cat(2,z,a);
                w = cat(2,w,b);
                clear('a','b')
            end
            tbl1 = table([x(:,1);x(:,2)],[zeros(length(x),1);ones(length(x),1)],[y(:,1);y(:,2)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl2 = table([x(:,1);x(:,3)],[zeros(length(x),1);ones(length(x),1)],[y(:,1);y(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl3 = table([x(:,2);x(:,3)],[zeros(length(x),1);ones(length(x),1)],[y(:,2);y(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl4 = table([z(:,1);z(:,2)],[zeros(length(z),1);ones(length(z),1)],[w(:,1);w(:,2)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl5 = table([z(:,1);z(:,3)],[zeros(length(z),1);ones(length(z),1)],[w(:,1);w(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl6 = table([z(:,2);z(:,3)],[zeros(length(z),1);ones(length(z),1)],[w(:,2);w(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            lme1 = fitlme(tbl1,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme2 = fitlme(tbl2,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme3 = fitlme(tbl3,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme4 = fitlme(tbl4,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme5 = fitlme(tbl5,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme6 = fitlme(tbl6,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            p1 = lme1.Coefficients.pValue(2);
            p2 = lme2.Coefficients.pValue(2);
            p3 = lme3.Coefficients.pValue(2);
            p4 = lme4.Coefficients.pValue(2);
            p5 = lme5.Coefficients.pValue(2);
            p6 = lme6.Coefficients.pValue(2);
            hold on
                bar(1:3,mean(x,1),'FaceColor','w')
                errorbar(1:3,mean(x,1),std(x,[],1)/sqrt(size(x,1)),'Color','k','LineStyle','none','CapSize',15)
                bar(4:6,mean(z,1),'FaceColor','w')
                errorbar(4:6,mean(z,1),std(z,[],1)/sqrt(size(z,1)),'Color','k','LineStyle','none','CapSize',15)
            hold off
            xlim([0 7])
            ylim([0.7 1])
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1) 
            xlabel('train','Fontsize',12,'Fontweight','b')
            ylabel('correlation','Fontsize',12,'Fontweight','b')
            title(odorNames(i),'Fontsize',12,'Fontweight','b')
            disp(['         ',odorNames{i},':'])
            disp('             within-train correlations:')
            disp(['                 train 1 vs train 2: p-value = ',num2str(p1),'.'])
            disp(['                 train 1 vs train 3: p-value = ',num2str(p2),'.'])
            disp(['                 train 2 vs train 3: p-value = ',num2str(p3),'.'])
            disp('             across-trains correlations:')
            disp(['                 trains 1-2 vs trains 1-3: p-value = ',num2str(p4),'.'])
            disp(['                 trains 1-2 vs trains 2-3: p-value = ',num2str(p5),'.'])
            disp(['                 trains 1-3 vs trains 2-3: p-value = ',num2str(p6),'.'])
    end
    disp('     ')
        
    % Average correlations across all flies
    figure, set(gcf,'position',[180 150 600 600])
        x = mean(R,3);
        idxMat = cat(3,[1:8; 1:8; 1:8; 9:16; 9:16; 9:16; 17:24; 17:24; 17:24],[1:8; 9:16; 17:24; 1:8; 9:16; 17:24; 1:8; 9:16; 17:24]);
        for i=1:9
            subplot(3,3,i)
                imagesc(x(idxMat(i,:,1),idxMat(i,:,2))), colormap(jet), caxis([0.7 1])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        end
        clear('x','idxMat')
    % Quantifications of correlations
    idxMat = cat(3,[1:8; 9:16; 17:24],[1:8; 9:16; 17:24]);
    Y1 = [];
    flyLabels1 = [];
    for i=1:3
        x = R(idxMat(i,:,1),idxMat(i,:,2),:);
        y1 = [];
        for j=1:size(x,3)
            a = x(:,:,j);
            idx = tril(ones(size(a)),-1);
            a = reshape(a(logical(idx)),[],1);
            y1 = cat(1,y1,a);
            if i==1
                flyLabels1 = cat(1,flyLabels1,j*ones(length(a),1));
            end
            clear('a','idx')
        end
        Y1 = cat(2,Y1,y1);
        clear('x','y1')
    end
    idxMat = cat(3,[9:16; 17:24; 17:24],[1:8; 1:8; 9:16]);
    Y2 = [];
    flyLabels2 =[];
    for i=1:3
        x = R(idxMat(i,:,1),idxMat(i,:,2),:);
        y2 = [];
        for j=1:size(x,3)
            a = x(:,:,j);
            y2 = cat(1,y2,a(:));
            if i==1
                flyLabels2 = cat(1,flyLabels2,j*ones(length(a(:)),1));
            end
            clear('a')
        end
        Y2 = cat(2,Y2,y2);
        clear('x','y2')
    end
    odorLabels1 = [];
    odorLabels2 = [];
    for i=1:length(odorNames)
        x = length(find(flyLabels1>3*(i-1) & flyLabels1<=3*i));
        odorLabels1 = cat(1,odorLabels1,i*ones(x,1));
        x = length(find(flyLabels2>3*(i-1) & flyLabels2<=3*i));
        odorLabels2 = cat(1,odorLabels2,i*ones(x,1));
    end
    % Bar plots
    figure, set(gcf, 'position',[390 100 600 600])
        hold on
            bar(1:3,mean(Y1,1),'FaceColor','none')
            errorbar(1:3,mean(Y1,1),std(Y1,[],1)/sqrt(size(Y1,1)),'Color','k','LineStyle','none','CapSize',20)
            bar(4:6,mean(Y2,1),'FaceColor','none')
            errorbar(4:6,mean(Y2,1),std(Y2,[],1)/sqrt(size(Y2,1)),'Color','k','LineStyle','none','CapSize',20)
        hold off
        ylim([0.7 1.0])
        set(gca,'Box','off','TickDir','out','XLim',[0 7],'Fontsize',10,'Fontweight','b','LineWidth',1)
        ylabel('correlation','Fontsize',12,'Fontweight','b')
        title('responses correlations','Fontsize',12,'Fontweight','b')
        clear('p1','p2','p3','p4','p5','p6')
    % Stats
        pValues = NaN(6,1);
        % Tables
        tbl1 = table([Y1(:,1);Y1(:,2)],[zeros(size(Y1,1),1);ones(size(Y1,1),1)],[odorLabels1;odorLabels1],[flyLabels1;flyLabels1],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl2 = table([Y1(:,1);Y1(:,3)],[zeros(size(Y1,1),1);ones(size(Y1,1),1)],[odorLabels1;odorLabels1],[flyLabels1;flyLabels1],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl3 = table([Y1(:,2);Y1(:,3)],[zeros(size(Y1,1),1);ones(size(Y1,1),1)],[odorLabels1;odorLabels1],[flyLabels1;flyLabels1],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl4 = table([Y2(:,1);Y2(:,2)],[zeros(size(Y2,1),1);ones(size(Y2,1),1)],[odorLabels2;odorLabels2],[flyLabels2;flyLabels2],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl5 = table([Y2(:,2);Y2(:,3)],[zeros(size(Y2,1),1);ones(size(Y2,1),1)],[odorLabels2;odorLabels2],[flyLabels2;flyLabels2],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl6 = table([Y2(:,1);Y2(:,3)],[zeros(size(Y2,1),1);ones(size(Y2,1),1)],[odorLabels2;odorLabels2],[flyLabels2;flyLabels2],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        % Linear Mixed-Effects Model
        lme1 = fitlme(tbl1,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme2 = fitlme(tbl2,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme3 = fitlme(tbl3,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme4 = fitlme(tbl4,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme5 = fitlme(tbl5,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme6 = fitlme(tbl6,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        % p Values
        pValues(1) = lme1.Coefficients.pValue(2);
        pValues(2) = lme2.Coefficients.pValue(2);
        pValues(3) = lme3.Coefficients.pValue(2);
        pValues(4) = lme4.Coefficients.pValue(2);
        pValues(5) = lme5.Coefficients.pValue(2);
        pValues(6) = lme6.Coefficients.pValue(2);
        clear('tbl1','tbl2','tbl3','tbl4','tbl5','tbl6','lme1','lme2','lme3','lme4','lme5','lme6')
        % Display p Values
        disp('     ')
        disp('     Linear Mixed-Effects Model on Odor Correlations (across all flies): ')
        disp('         Intra-Train Correlations: ')
        disp(['             Train 1 vs Train 2, p-value = ',num2str(pValues(1)),'.'])
        disp(['             Train 1 vs Train 3, p-value = ',num2str(pValues(2)),'.'])
        disp(['             Train 2 vs Train 3, p-value = ',num2str(pValues(3)),'.'])
        disp('         Inter-Train Correlations: ')
        disp(['             Trains 1-2 vs Trains 1-3, p-value = ',num2str(pValues(4)),'.'])
        disp(['             Trains 1-3 vs Trains 2-3, p-value = ',num2str(pValues(5)),'.'])
        disp(['             Trains 1-2 vs Trains 2-3, p-value = ',num2str(pValues(6)),'.'])
        disp('     ')
    clear('Y1','Y2')
    
end

    
%% Euclidean distances

if EuclideanDistanceAnalysis == 1

    % EuD per Fly
    EuclideanDistances = [];
    for i=1:length(glomeruliPerFly)-1
        x = odorTimeSeries(glomeruliPerFly(i)+1:glomeruliPerFly(i+1),:,:);
        x = squeeze(mean(x(:,odorTime,:),2));
        x = pdist(x'); x = squareform(x);
        EuclideanDistances = cat(3,EuclideanDistances,x);
        clear('x')
    end
    figure, set(gcf,'position',[730 80 1150 720])
        for i=1:size(EuclideanDistances,3)
            subplot(4,6,i)
                imagesc(EuclideanDistances(:,:,i)), colormap(jet), caxis([0 50])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                title(['Fly # ',num2str(i)],'Fontsize',12,'Fontweight','b')
        end
        
    % EuD - average per odor (3 flies per odor)
    figure, set(gcf,'position',[690 125 1375 635])
        idxMat = cat(3,[1:8; 1:8; 1:8; 9:16; 9:16; 9:16; 17:24; 17:24; 17:24],[1:8; 9:16; 17:24; 1:8; 9:16; 17:24; 1:8; 9:16; 17:24]);
        subPlotID1 = [1 2 3 17 18 19 33 34 35];
        subPlotID2 = [65 66 67 81 82 83 97 98 99];
        colorAxisLimits = [40 60 60 50 60 70 60 60];
        for i=1:length(odorNames)
            x = mean(EuclideanDistances(:,:,(i*3)-2:i*3),3);
            for j=1:9
                if i<5
                    subplot(7,16,subPlotID1(j)+4*(i-1))
                else
                    subplot(7,16,subPlotID2(j)+4*(i-5))
                end
                imagesc(x(idxMat(j,:,1),idxMat(j,:,2))), colormap(jet), caxis([0 colorAxisLimits(i)])
                set(gca,'Box','off','TickDir','out','XTickLabel',[],'YTickLabel',[])
                if j==2
                    title(odorNames(i),'Fontsize',12,'Fontweight','b')
                end
            end
        end
    % Stats on EuD - grouped by odor (3 flies per odor)
    figure, set(gcf,'position',[510 75 1355 720])
    disp('     ')
    disp('     Linear Mixed-Effects Model on Odor Euclidean Distances (grouped by odor):')
    for i=1:8
        subplot(2,4,i)
            X = EuclideanDistances(:,:,(1:3)+3*(i-1));
            Y = cat(3,ones(size(R,1),size(R,2)),2*ones(size(R,1),size(R,2)),3*ones(size(R,1),size(R,2)));
            matIdx = [1:8;9:16;17:24];
            x = [];
            y = [];
            for j=1:3
                a = X(matIdx(j,:),matIdx(j,:),:);
                b = Y(matIdx(j,:),matIdx(j,:),:);
                idx = tril(ones(size(a,1),size(a,2)),-1);
                a = reshape(a(logical(idx)),[],1);
                b = reshape(b(logical(idx)),[],1);
                x = cat(2,x,a);
                y = cat(2,y,b);
                clear('a','b')
            end
            matIdx = cat(3,[9:16; 17:24; 17:24],[1:8; 1:8; 9:16]);  % trains 1-2, trains 1-3, trains 2-3
            z = [];
            w = [];
            for j=1:3
                a = X(matIdx(j,:,1),matIdx(j,:,2),:);
                b = Y(matIdx(j,:,1),matIdx(j,:,2),:);
                a = reshape(a(logical(idx)),[],1);
                b = reshape(b(logical(idx)),[],1);
                z = cat(2,z,a);
                w = cat(2,w,b);
                clear('a','b')
            end
            tbl1 = table([x(:,1);x(:,2)],[zeros(length(x),1);ones(length(x),1)],[y(:,1);y(:,2)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl2 = table([x(:,1);x(:,3)],[zeros(length(x),1);ones(length(x),1)],[y(:,1);y(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl3 = table([x(:,2);x(:,3)],[zeros(length(x),1);ones(length(x),1)],[y(:,2);y(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl4 = table([z(:,1);z(:,2)],[zeros(length(z),1);ones(length(z),1)],[w(:,1);w(:,2)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl5 = table([z(:,1);z(:,3)],[zeros(length(z),1);ones(length(z),1)],[w(:,1);w(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl6 = table([z(:,2);z(:,3)],[zeros(length(z),1);ones(length(z),1)],[w(:,2);w(:,3)],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            lme1 = fitlme(tbl1,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme2 = fitlme(tbl2,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme3 = fitlme(tbl3,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme4 = fitlme(tbl4,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme5 = fitlme(tbl5,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme6 = fitlme(tbl6,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            p1 = lme1.Coefficients.pValue(2);
            p2 = lme2.Coefficients.pValue(2);
            p3 = lme3.Coefficients.pValue(2);
            p4 = lme4.Coefficients.pValue(2);
            p5 = lme5.Coefficients.pValue(2);
            p6 = lme6.Coefficients.pValue(2);
            hold on
                bar(1:3,mean(x,1),'FaceColor','w')
                errorbar(1:3,mean(x,1),std(x,[],1)/sqrt(size(x,1)),'Color','k','LineStyle','none','CapSize',15)
                bar(4:6,mean(z,1),'FaceColor','w')
                errorbar(4:6,mean(z,1),std(z,[],1)/sqrt(size(z,1)),'Color','k','LineStyle','none','CapSize',15)
            hold off
            xlim([0 7])
            ylim([0 50])
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1) 
            xlabel('train','Fontsize',12,'Fontweight','b')
            ylabel('Euclidean distance','Fontsize',12,'Fontweight','b')
            title(odorNames(i),'Fontsize',12,'Fontweight','b')
            disp(['         ',odorNames{i},':'])
            disp('             within-train Euclidean distances:')
            disp(['                 train 1 vs train 2: p-value = ',num2str(p1),'.'])
            disp(['                 train 1 vs train 3: p-value = ',num2str(p2),'.'])
            disp(['                 train 2 vs train 3: p-value = ',num2str(p3),'.'])
            disp('             across-trains Euclidean distances:')
            disp(['                 trains 1-2 vs trains 1-3: p-value = ',num2str(p4),'.'])
            disp(['                 trains 1-2 vs trains 2-3: p-value = ',num2str(p5),'.'])
            disp(['                 trains 1-3 vs trains 2-3: p-value = ',num2str(p6),'.'])
    end
    disp('     ')
        
    % Average Euclidean distances across all flies
    figure, set(gcf,'position',[830 150 600 600])
        x = mean(EuclideanDistances,3);
        idxMat = cat(3,[1:8; 1:8; 1:8; 9:16; 9:16; 9:16; 17:24; 17:24; 17:24],[1:8; 9:16; 17:24; 1:8; 9:16; 17:24; 1:8; 9:16; 17:24]);
        for i=1:9
            subplot(3,3,i)
                imagesc(x(idxMat(i,:,1),idxMat(i,:,2))), colormap(jet), caxis([0 50])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        end
        clear('x','idxMat')
    % Quantifications of Euclidean distances
    idxMat = cat(3,[1:8; 9:16; 17:24],[1:8; 9:16; 17:24]);
    Y1 = [];
    flyLabels1 = [];
    for i=1:3
        x = EuclideanDistances(idxMat(i,:,1),idxMat(i,:,2),:);
        y1 = [];
        for j=1:size(x,3)
            a = x(:,:,j);
            idx = tril(ones(size(a)),-1);
            a = reshape(a(logical(idx)),[],1);
            y1 = cat(1,y1,a);
            if i==1
                flyLabels1 = cat(1,flyLabels1,j*ones(length(a),1));
            end
            clear('a','idx')
        end
        Y1 = cat(2,Y1,y1);
        clear('x','y1')
    end
    idxMat = cat(3,[9:16; 17:24; 17:24],[1:8; 1:8; 9:16]);
    Y2 = [];
    flyLabels2 = [];
    for i=1:3
        x = EuclideanDistances(idxMat(i,:,1),idxMat(i,:,2),:);
        y2 = [];
        for j=1:size(x,3)
            a = x(:,:,j);
            y2 = cat(1,y2,a(:));
            if i==1
                flyLabels2 = cat(1,flyLabels2,j*ones(length(a(:)),1));
            end
            clear('a')
        end
        Y2 = cat(2,Y2,y2);
        clear('x','y2')
    end
    odorLabels1 = [];
    odorLabels2 = [];
    for i=1:length(odorNames)
        x = length(find(flyLabels1>3*(i-1) & flyLabels1<=3*i));
        odorLabels1 = cat(1,odorLabels1,i*ones(x,1));
        x = length(find(flyLabels2>3*(i-1) & flyLabels2<=3*i));
        odorLabels2 = cat(1,odorLabels2,i*ones(x,1));
    end
    figure, set(gcf, 'position',[1040 100 600 600])
        hold on
            bar(1:3,mean(Y1,1),'FaceColor','w')
            errorbar(1:3,mean(Y1,1),std(Y1,[],1)/sqrt(size(Y1,1)),'Color','k','LineStyle','none','CapSize',20)
            bar(4:6,mean(Y2,1),'FaceColor','w')
            errorbar(4:6,mean(Y2,1),std(Y2,[],1)/sqrt(size(Y2,1)),'Color','k','LineStyle','none','CapSize',20)
        hold off
        ylim([0 50])
        set(gca,'Box','off','TickDir','out','XLim',[0 7],'Fontsize',10,'Fontweight','b','LineWidth',1)
        ylabel('Euclidean distance','Fontsize',12,'Fontweight','b')
        title('odor response, Euclidean distance','Fontsize',12,'Fontweight','b')
        clear('p1','p2','p3','p4','p5','p6')
    % Stats
        pValues = NaN(6,1);
        % Tables
        tbl1 = table([Y1(:,1);Y1(:,2)],[zeros(size(Y1,1),1);ones(size(Y1,1),1)],[odorLabels1;odorLabels1],[flyLabels1;flyLabels1],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl2 = table([Y1(:,1);Y1(:,3)],[zeros(size(Y1,1),1);ones(size(Y1,1),1)],[odorLabels1;odorLabels1],[flyLabels1;flyLabels1],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl3 = table([Y1(:,2);Y1(:,3)],[zeros(size(Y1,1),1);ones(size(Y1,1),1)],[odorLabels1;odorLabels1],[flyLabels1;flyLabels1],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl4 = table([Y2(:,1);Y2(:,2)],[zeros(size(Y2,1),1);ones(size(Y2,1),1)],[odorLabels2;odorLabels2],[flyLabels2;flyLabels2],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl5 = table([Y2(:,2);Y2(:,3)],[zeros(size(Y2,1),1);ones(size(Y2,1),1)],[odorLabels2;odorLabels2],[flyLabels2;flyLabels2],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl6 = table([Y2(:,1);Y2(:,3)],[zeros(size(Y2,1),1);ones(size(Y2,1),1)],[odorLabels2;odorLabels2],[flyLabels2;flyLabels2],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        % Linear Mixed-Effect Model
        lme1 = fitlme(tbl1,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme2 = fitlme(tbl2,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme3 = fitlme(tbl3,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme4 = fitlme(tbl4,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme5 = fitlme(tbl5,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme6 = fitlme(tbl6,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        % p Values
        pValues(1) = lme1.Coefficients.pValue(2);
        pValues(2) = lme2.Coefficients.pValue(2);
        pValues(3) = lme3.Coefficients.pValue(2);
        pValues(4) = lme4.Coefficients.pValue(2);
        pValues(5) = lme5.Coefficients.pValue(2);
        pValues(6) = lme6.Coefficients.pValue(2);
        clear('tbl1','tbl2','tbl3','tbl4','tbl5','tbl6','lme1','lme2','lme3','lme4','lme5','lme6')
        % Display p Values
        disp('     ')
        disp('     Linear Mixed-Effects Model on Odor Euclidean Distances (across all flies): ')
        disp('         Intra-Train Euclidean Distances: ')
        disp(['             Train 1 vs Train 2, p-value = ',num2str(pValues(1)),'.'])
        disp(['             Train 1 vs Train 3, p-value = ',num2str(pValues(2)),'.'])
        disp(['             Train 2 vs Train 3, p-value = ',num2str(pValues(3)),'.'])
        disp('         Inter-Train Euclidean Distances: ')
        disp(['             Trains 1-2 vs Trains 1-3, p-value = ',num2str(pValues(4)),'.'])
        disp(['             Trains 1-3 vs Trains 2-3, p-value = ',num2str(pValues(5)),'.'])
        disp(['             Trains 1-2 vs Trains 2-3, p-value = ',num2str(pValues(6)),'.'])
        disp('     ')
    clear('Y1','Y2')
    
end

    
%% Glomerular responses per odor (3 flies per odor)

if glomerularResponsesAnalysis == 1

    % Glomerular odor responses - grouped by odor (3 flies per odor)
    figure, set(gcf,'position',[40 125 1375 635])
        for i=1:length(odorNames)
            x = odorTimeSeries(glomeruliPerOdor(i)+1:glomeruliPerOdor(i+1),:,:);
            x = squeeze(mean(x(:,odorTime,:),2));
            [~,order] = sort(mean(x,2),'descend');
            x = x(order,:);
                for j=1:3
                    subplot(2,16,j+4*(i-1))
                        imagesc(x(:,(1:8)+8*(j-1))), colormap(jet), caxis([0 25])
                        set(gca,'Box','off','TickDir','out','XTick',[1 8],'Fontsize',10,'Fontweight','b','LineWidth',1)    % 'YTick',[1 size(x,1)],
                        if j==2
                            title(odorNames(i),'Fontsize',12,'Fontweight','b')
                        end
                        if j>1
                            set(gca,'YTickLabel',[])
                        end
                end
                clear('x')
        end
    % Heatmap of glomerular odor responses
    figure, set(gcf,'position',[85 180 800 500])
        averageOdorResponses = [];
        for i=1:length(glomeruliPerFly)-1
            x = odorTimeSeries(glomeruliPerFly(i)+1:glomeruliPerFly(i+1),:,:);
            x = squeeze(mean(x(:,odorTime,:),2));
            averageOdorResponses = cat(1,averageOdorResponses,x);
            clear('x')
        end
        x = mean(averageOdorResponses,2);
        [~,I] = sort(x,'descend');
        subplot(1,3,1)
            imagesc(averageOdorResponses(I,1:8)), caxis([0 20]), colormap(jet)
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        subplot(1,3,2)
            imagesc(averageOdorResponses(I,9:16)), caxis([0 20]), colormap(jet)
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
            title('Glomerular responses','Fontsize',12,'Fontweight','b')
        subplot(1,3,3)
            imagesc(averageOdorResponses(I,17:24)), caxis([0 20]), colormap(jet)
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        clear('x','I')
    % Histograms of glomerular odor responses
    figure, set(gcf,'position',[700 150 600 550])
        resolution = -6:4:22;
        probAxis = [0 0.4];
        subplot(1,3,1)
            x = averageOdorResponses(:,1:8); x = x(:);
            histogram(x,resolution,'Normalization','probability','FaceColor',[0.8 0.8 0.8]);
            clear('x')
            ylim(probAxis)
            set(gca,'Box','off','TickDir','out','XLim',[-8 26],'Fontsize',10,'Fontweight','b','LineWidth',1)
            xlabel('odor response {\Delta}F/F','Fontsize',12,'Fontweight','b')
            ylabel('probability','Fontsize',12,'Fontweight','b')
        subplot(1,3,2)
            x = averageOdorResponses(:,9:16); x = x(:);
            histogram(x,resolution,'Normalization','probability','FaceColor',[0.4 0.4 0.4]);
            clear('x')
            ylim(probAxis)
            set(gca,'Box','off','TickDir','out','XLim',[-8 26],'Fontsize',10,'Fontweight','b','LineWidth',1)
            title('odor response','Fontsize',12,'Fontweight','b')
        subplot(1,3,3)
            x = averageOdorResponses(:,17:24); x = x(:);
            histogram(x,resolution,'Normalization','probability','FaceColor',[0 0 0]);
            clear('x')
            ylim(probAxis)
            set(gca,'Box','off','TickDir','out','XLim',[-8 26],'Fontsize',10,'Fontweight','b','LineWidth',1)
        clear('resolution','probAxis')

    % Glomerular response change from 1st stimulus
    responseChangeFromFirstStim = odorTimeSeries - odorTimeSeries(:,:,1);
    responseChangeFromFirstStim = squeeze(mean(responseChangeFromFirstStim(:,odorTime,:),2));
    % Heatmaps of glomerular odor response change - grouped by odor (3 flies per odor)
    R = [zeros(1,66) linspace(0,1,17) ones(1,17)]'; G = [linspace(1,0,33) zeros(1,50) linspace(0,1,17)]'; B = [ones(1,33) linspace(1,0,33) zeros(1,34)]';
    cmap = [R G B];
    figure, set(gcf,'position',[500 125 1375 635])
    for i=1:length(odorNames)
        x = responseChangeFromFirstStim(glomeruliPerOdor(i)+1:glomeruliPerOdor(i+1),:);
        [~,order] = sort(mean(x(:,2:end),2),'descend');
        x = x(order,:);
            for j=1:3
                subplot(2,16,j+4*(i-1))
                    imagesc(x(:,(1:8)+8*(j-1))), colormap(cmap), caxis([-10 5])
                    set(gca,'Box','off','TickDir','out','XTick',[1 8],'Fontsize',10,'Fontweight','b','LineWidth',1)    % 'YTick',[1 size(x,1)],
                    if j==2
                        title(odorNames(i),'Fontsize',12,'Fontweight','b')
                    end
                    if j>1
                        set(gca,'YTickLabel',[])
                    end
            end
            clear('x')
    end
    % Heatmap of glomerular response change
    R = [zeros(1,66) linspace(0,1,17) ones(1,17)]'; G = [linspace(1,0,33) zeros(1,50) linspace(0,1,17)]'; B = [ones(1,33) linspace(1,0,33) zeros(1,34)]';
    cmap = [R G B];
    colorAxisLimits = [-10 5];
    figure, set(gcf,'position',[530 185 800 500])
        x = mean(responseChangeFromFirstStim(:,2:end),2);
        [~,I] = sort(x,'descend');
        subplot(1,3,1)
            imagesc(responseChangeFromFirstStim(I,1:8)), caxis(colorAxisLimits), colormap(cmap)
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        subplot(1,3,2)
            imagesc(responseChangeFromFirstStim(I,9:16)), caxis(colorAxisLimits), colormap(cmap)
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
            title('Glomerular responses change','Fontsize',12,'Fontweight','b')
        subplot(1,3,3)
            imagesc(responseChangeFromFirstStim(I,17:24)), caxis(colorAxisLimits), colormap(cmap)
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        clear('x','I')
    % Histograms of glomerular response change
    figure, set(gcf,'position',[1250 160 600 550])
        resolution = -10:2:6;
        probAxis = [0 0.3];
        subplot(1,3,1)
            x = responseChangeFromFirstStim(:,2:8); x = x(:);
            histogram(x,resolution,'Normalization','probability','FaceColor',[0.8 0.8 0.8]);
            clear('x')
            ylim(probAxis)
            set(gca,'Box','off','TickDir','out','XLim',[-10 6],'Fontsize',10,'Fontweight','b','LineWidth',1)
            xlabel('response change {\Delta}F/F','Fontsize',12,'Fontweight','b')
            ylabel('probability','Fontsize',12,'Fontweight','b')
        subplot(1,3,2)
            x = responseChangeFromFirstStim(:,9:16); x = x(:);
            histogram(x,resolution,'Normalization','probability','FaceColor',[0.4 0.4 0.4]);
            clear('x')
            ylim(probAxis)
            set(gca,'Box','off','TickDir','out','XLim',[-10 6],'Fontsize',10,'Fontweight','b','LineWidth',1)
            title('odor response change','Fontsize',12,'Fontweight','b')
        subplot(1,3,3)
            x = responseChangeFromFirstStim(:,17:24); x = x(:);
            histogram(x,resolution,'Normalization','probability','FaceColor',[0 0 0]);
            clear('x')
            ylim(probAxis)
            set(gca,'Box','off','TickDir','out','XLim',[-10 6],'Fontsize',10,'Fontweight','b','LineWidth',1)
        clear('resolution','probAxis')  
    % CDF of glomerular reponse change
    figure, set(gcf,'position',[650 190 600 500])
        colors = colormap(jet(24));
        hold on
            for i=2:24
                cdfplot(responseChangeFromFirstStim(:,i)); hline = findobj(gcf,'type','line'); set(hline(1),'Color',colors(i,:))
            end
        hold off
        xlim([-20 10])
        set(gca,'TickDir','out','Fontsize',10,'Fontweight','b','Box','off'), grid off
        title('')
        xlabel ('DFF change','Color','k','Fontsize',12,'Fontweight','b')
        ylabel ('CDF(DFF change)','Color','k','Fontsize',12,'Fontweight','b')   
    % Linear Mixed-Effects Model
    odorLabels = NaN(1268,1);
    for i=1:8
        odorLabels(glomeruliPerOdor(i)+1:glomeruliPerOdor(i+1)) = i;
    end
    flyLabels = NaN(1268,1); 
    for i=1:24
        flyLabels(glomeruliPerFly(i)+1:glomeruliPerFly(i+1)) = i;
    end
    DFFChangeStats = ones(size(responseChangeFromFirstStim,2),size(responseChangeFromFirstStim,2));
    for i=2:size(responseChangeFromFirstStim,2)
        for j=i:size(responseChangeFromFirstStim,2)
            tbl = table([responseChangeFromFirstStim(:,i);responseChangeFromFirstStim(:,j)],...
                  [zeros(size(responseChangeFromFirstStim,1),1);ones(size(responseChangeFromFirstStim,1),1)],...
                  [odorLabels;odorLabels],[flyLabels;flyLabels],...
                  'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
            lme = fitlme(tbl,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
            p = lme.Coefficients.pValue(2);
            DFFChangeStats(i,j) = p;
            DFFChangeStats(j,i) = p;
            clear('p')
        end
    end
    cmap = jet; cmap(end,:) = [1 1 1];
    figure, set(gcf,'position',[1265 135 600 600])
        idxMat = cat(3,[1:8; 1:8; 1:8; 9:16; 9:16; 9:16; 17:24; 17:24; 17:24],[1:8; 9:16; 17:24; 1:8; 9:16; 17:24; 1:8; 9:16; 17:24]);
        for i=1:9
            subplot(3,3,i)
                imagesc(DFFChangeStats(idxMat(i,:,1),idxMat(i,:,2))), colormap(cmap), caxis([0 0.05])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                if i==2
                    title('DFF change stats (mixed-effect model)','Fontsize',12,'Fontweight','b')
                end
        end
        clear('idxMat')
        set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        
end

        
%% Relative Change in Odor Responses

if relativeResponseChangeAnalysis == 1

    method = 'standardDeviation'; % standardDeviation, or percentile

    % Odor responses per fly
    splitResponsesPerFly = cell(length(glomeruliPerFly)-1,4);
    flyID = cell(length(glomeruliPerFly)-1,4);
    for i=1:length(glomeruliPerFly)-1
        x = odorTimeSeries(glomeruliPerFly(i)+1:glomeruliPerFly(i+1),:,:);
        x = squeeze(mean(x(:,odorTime,:),2));
        y = mean(x,2);
        switch method
            case 'standardDeviation'
                xMean = mean(y); xStd = std(y);
                idx1 = y <= xMean-1*xStd;
                idx2 = y >  xMean-1*xStd & y <= xMean+1*xStd;
                idx3 = y >  xMean+1*xStd;
            case 'percentile'
                idx1 = y <= prctile(y,25);
                idx2 = y > prctile(y,25) & y <= prctile(y,75);
                idx3 = y > prctile(y,75);
        end
        % Response Change Fraction
        x = (x - mean(x(:,1:8),2)) ./ mean(x(:,1:8),2);
        splitResponsesPerFly{i,1} = x(idx1,:);
        splitResponsesPerFly{i,2} = x(idx2,:);
        splitResponsesPerFly{i,3} = x(idx3,:);
        flyID{i,1} = length(find(idx1));
        flyID{i,2} = length(find(idx2));
        flyID{i,3} = length(find(idx3));
        clear('x','y','xMean','xStd','idx1','idx2','idx3')
    end

    % Plot Relative Change in Odor Responses
    colors = [0 0 1;0 0.85 0;1 0 0];
    figure, set(gcf,'position',[50 100 600 750])
        hold on
            for i=1:3
                x = [];
                for j=1:length(glomeruliPerFly)-1
                    x = cat(1,x,splitResponsesPerFly{j,i});
                end
                xMean = mean(x,1);
                xStd = std(x,[],1);
                plot(1:8,xMean(1:8),'Color',colors(i,:))
                plot(1:8,xMean(1:8)-xStd(1:8)/sqrt(size(x,1)),'Color',colors(i,:))
                plot(1:8,xMean(1:8)+xStd(1:8)/sqrt(size(x,1)),'Color',colors(i,:))
                plot(10:17,xMean(9:16),'Color',colors(i,:))
                plot(10:17,xMean(9:16)-xStd(9:16)/sqrt(size(x,1)),'Color',colors(i,:))
                plot(10:17,xMean(9:16)+xStd(9:16)/sqrt(size(x,1)),'Color',colors(i,:))
                plot(19:26,xMean(17:24),'Color',colors(i,:))
                plot(19:26,xMean(17:24)-xStd(17:24)/sqrt(size(x,1)),'Color',colors(i,:))
                plot(19:26,xMean(17:24)+xStd(17:24)/sqrt(size(x,1)),'Color',colors(i,:))
                text(1,-0.6-0.05*(i-1),[num2str(size(x,1)),' glomeruli'],'Color',colors(i,:))
                clear('x','xMean','xStd')
            end
        hold off
        xlim([0 27])
        ylim([-0.8 0.4])
        set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        xlabel('stim. #')
        ylabel('odor response change ({\Delta}FF)')
        
    % Stats
    pValuesMatrix = ones(3,3,4);
    for i=1:3
        odorRespChange = [];
        for j=1:length(glomeruliPerFly)-1
            odorRespChange = cat(1,odorRespChange,splitResponsesPerFly{j,i});
        end
        flyLabels = [];
        for j=1:length(glomeruliPerFly)-1
            flyLabels = cat(1,flyLabels,j*ones(flyID{j,i},1));
        end
        odorLabels = [];
        for k=1:length(odorNames)
            x = length(find(flyLabels>3*(k-1) & flyLabels<=3*k));
            odorLabels = cat(1,odorLabels,k*ones(x,1));
            clear('x')
        end
        flyLabels = repmat(flyLabels,[1 8]);
        odorLabels = repmat(odorLabels,[1 8]);
        % Data
        X1 = odorRespChange(:,1:8);
        X2 = odorRespChange(:,9:16);
        X3 = odorRespChange(:,17:24);
        X1 = reshape(X1,[],1);
        X2 = reshape(X2,[],1);
        X3 = reshape(X3,[],1);
        flyLabels = reshape(flyLabels,[],1);
        odorLabels = reshape(odorLabels,[],1);
        % Tables
        tbl1 = table([X1;X2],[zeros(length(X1),1);ones(length(X2),1)],[odorLabels;odorLabels],[flyLabels;flyLabels],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl2 = table([X1;X3],[zeros(length(X1),1);ones(length(X2),1)],[odorLabels;odorLabels],[flyLabels;flyLabels],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        tbl3 = table([X2;X3],[zeros(length(X1),1);ones(length(X2),1)],[odorLabels;odorLabels],[flyLabels;flyLabels],...
               'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
        % Linear Mixed-Effects Model
        pValues = NaN(3,1);
        lme1 = fitlme(tbl1,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme2 = fitlme(tbl2,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        lme3 = fitlme(tbl3,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
        % p Values
        pValues(1) = lme1.Coefficients.pValue(2);
        pValues(2) = lme2.Coefficients.pValue(2);
        pValues(3) = lme3.Coefficients.pValue(2);
        pValuesMatrix(1,2,i) = pValues(1); pValuesMatrix(2,1,i) = pValues(1);
        pValuesMatrix(1,3,i) = pValues(2); pValuesMatrix(3,1,i) = pValues(2);
        pValuesMatrix(2,3,i) = pValues(3); pValuesMatrix(3,2,i) = pValues(3);
        clear('X1','X2','X3','odorRespChange','odorLabels','flyLabels','tbl1','tbl2','tbl3','lme1','lme2','lme3')
        % Display p Values
        if i==1
            disp('     ')
            disp('     Linear Mixed-Effects Model on Relative Change in Odor Responses: ')
        end
        if i==1
            disp('         Weak responses (< mean - std): ')
        elseif i==2
            disp('         Medium responses (> mean - std & < mean + std): ')
        elseif i==3
            disp('         Strong responses (> mean + std): ')
        end
        disp(['             Train 1 vs Train 2, p-value = ',num2str(pValues(1)),'.'])
        disp(['             Train 1 vs Train 3, p-value = ',num2str(pValues(2)),'.'])
        disp(['             Train 2 vs Train 3, p-value = ',num2str(pValues(3)),'.'])
        if i==3
            disp('     ')
        end
    end
    % p-value Matrices
    figure, set(gcf,'position',[860 90 880 240])
        cmap = jet; cmap(end,:) = [1 1 1];
        for i=1:3
            subplot(1,3,i)
                imagesc(pValuesMatrix(:,:,i)), colormap(cmap), caxis([0 0.05])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                if i==1
                    title('weak')
                elseif i==2
                    title('medium')
                elseif i==3
                    title('strong')
                end
        end
    % Bar plots
    figure, set(gcf,'position',[685 430 1250 420])
        barColors = [0.8 0.8 0.8;0.4 0.4 0.4;0 0 0];
        matIdx = [1:8;9:16;17:24];
        for i=1:3
            odorRespChange = [];
            for j=1:length(glomeruliPerFly)-1
                odorRespChange = cat(1,odorRespChange,splitResponsesPerFly{j,i});
            end
            subplot(1,3,i)
                hold on
                    for j=1:3
                        bar(j,mean(odorRespChange(:,matIdx(j,:)),'all'),'FaceColor',barColors(j,:));
                        errorbar(j,mean(odorRespChange(:,matIdx(j,:)),'all'),...
                                 std(odorRespChange(:,matIdx(j,:)),[],'all')/sqrt(numel(odorRespChange(:,matIdx(j,:)))),...
                                 'k','CapSize',20)
                    end
                hold off
                xlim([0 4])
                ylim([-0.6 0.1])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                xlabel('train')
                ylabel('response change fraction')
                if i==1
                    title('weak')
                elseif i==2
                    title('medium')
                elseif i==3
                    title('strong')
                end
        end
        
end


%% Population sparseness

if sparsenessAnalysis == 1

    populationSparseness = [];
    for k=1:repetitionsPerTrain*stimulationTrains
        % Population Sparseness per Odor
        odorSparseness = [];
        if k==1
            flyLabels = [];
        end
        for j=1:length(glomeruliPerFly)-1
            x = squeeze(mean(odorTimeSeries(glomeruliPerFly(j)+1:glomeruliPerFly(j+1),odorTime,k),2));
            x = x';
                repetitionSparseness = zeros(size(x,1),1);
            for i=1:size(x,1)
                vector = x(i,:);
                vector = vector - min(vector);
                vector = vector / max(vector);
                a = sum(vector/length(vector))^2;
                b = sum((vector.^2)/length(vector));
                if b==0
                    repetitionSparseness(i,1) = 0;
                else
                    repetitionSparseness(i,1) = (1-a/b) / (1-1/length(vector));
                end
                clear('vector','a','b')
            end
                odorSparseness = cat(1,odorSparseness,repetitionSparseness);
                if k==1
                    flyLabels = cat(1,flyLabels,j*ones(length(repetitionSparseness),1));
                end
                clear('x','repetitionSparseness')
        end
        populationSparseness = cat(2,populationSparseness,odorSparseness);
        clear('odorSparseness')
    end    
    odorLabels = [];
    for i=1:length(odorNames)
        x = length(find(flyLabels>3*(i-1) & flyLabels<=3*i));
        odorLabels = cat(1,odorLabels,i*ones(x,1));
    end
    
    % Plot Population Sparseness
    figure, set(gcf,'position',[160 100 600 750])
        hold on
            errorbar(1:8,mean(populationSparseness(:,1:8),1),std(populationSparseness(:,1:8),[],1)/sqrt(size(populationSparseness(:,1:8),1)),'Color','k','CapSize',12.5)
            errorbar(10:17,mean(populationSparseness(:,9:16),1),std(populationSparseness(:,9:16),[],1)/sqrt(size(populationSparseness(:,9:16),1)),'Color','k','CapSize',12.5)
            errorbar(19:26,mean(populationSparseness(:,17:24),1),std(populationSparseness(:,17:24),[],1)/sqrt(size(populationSparseness(:,17:24),1)),'Color','k','CapSize',12.5)
        hold off
        ylim([0.2 0.4])
        xlim([0 27])
        set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
        xlabel('stim. #','Fontsize',12,'Fontweight','b')
        ylabel('sparseness','Fontsize',12,'Fontweight','b')
        title('population sparseness','Fontsize',12,'Fontweight','b')
        
    % Linear Mixed-Effects Model on Population Sparseness 
    sparsenessStats = NaN(size(populationSparseness,2),size(populationSparseness,2));
    for i=1:size(populationSparseness,2)
        for j=i:size(populationSparseness,2)
            tbl = table([populationSparseness(:,i);populationSparseness(:,j)],...
                  [zeros(size(populationSparseness,1),1);ones(size(populationSparseness,1),1)],...
                  [odorLabels;odorLabels],[flyLabels;flyLabels],...
                  'VariableNames',{'responses','responseFlag','odorLabels','flyLabels'});
            lme = fitlme(tbl,'responses~responseFlag+(1|odorLabels)+(1|flyLabels)','FitMethod','REML');
            p = lme.Coefficients.pValue(2);
            sparsenessStats(i,j) = p;
            sparsenessStats(j,i) = p;
            clear('p')
        end
    end
    figure, set(gcf,'position',[415 115 600 600])
        cmap = jet; cmap(end,:) = [1 1 1];
        idxMat = cat(3,[1:8; 1:8; 1:8; 9:16; 9:16; 9:16; 17:24; 17:24; 17:24],[1:8; 9:16; 17:24; 1:8; 9:16; 17:24; 1:8; 9:16; 17:24]);
        for i=1:9
            subplot(3,3,i)
                imagesc(sparsenessStats(idxMat(i,:,1),idxMat(i,:,2))), caxis([0 0.05])
                colormap(cmap)
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                if i==2
                    title('sparseness stats (mixed-effect model)','Fontsize',12,'Fontweight','b')
                end
        end
        clear('idxMat')
        set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)    
    
    % Population Sparseness per Odor
    figure, set(gcf,'position',[500 75 1355 720])
        yLimMat = [0.15 0.45; 0.20 0.45; 0.20 0.45; 0.10 0.40; 0.15 0.40; 0.15 0.45; 0.15 0.50; 0.25 0.55];
        for i=1:8
            subplot(2,4,i)
                idx = (1:3)+3*(i-1);
                hold on
                    errorbar(1:8,mean(populationSparseness(idx,1:8),1),std(populationSparseness(idx,1:8),[],1)/sqrt(size(populationSparseness(idx,1:8),1)),'Color','k','CapSize',7.5)
                    errorbar(10:17,mean(populationSparseness(idx,9:16),1),std(populationSparseness(idx,9:16),[],1)/sqrt(size(populationSparseness(idx,9:16),1)),'Color','k','CapSize',7.5)
                    errorbar(19:26,mean(populationSparseness(idx,17:24),1),std(populationSparseness(idx,17:24),[],1)/sqrt(size(populationSparseness(idx,17:24),1)),'Color','k','CapSize',7.5)
                hold off
                ylim(yLimMat(i,:))
                xlim([0 27])
                set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1)
                if i>4
                    xlabel('stim. #','Fontsize',12,'Fontweight','b')
                end
                if i==1 || i==5
                    ylabel('sparseness','Fontsize',12,'Fontweight','b')
                end
                title(odorNames(i),'Fontsize',12,'Fontweight','b')
        end
    % Stats on Population Sparseness per Odor
    figure, set(gcf,'position',[540 75 1355 720])
    disp('     ')
    disp('     Linear Mixed-Effects Model on Population Sparseness across Trains of Stimulation:')
    for i=1:8
        subplot(2,4,i)
            idx = (1:3)+3*(i-1);
            x1 = populationSparseness(idx,1:8);    x1 = x1(:);
            x2 = populationSparseness(idx,9:16);   x2 = x2(:);
            x3 = populationSparseness(idx,17:24);  x3 = x3(:);
            y = repmat(flyLabels(idx),[1 8]);       y = y(:);
            tbl1 = table([x1;x2],[zeros(length(x1),1);ones(length(x2),1)],[y;y],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl2 = table([x1;x3],[zeros(length(x1),1);ones(length(x3),1)],[y;y],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            tbl3 = table([x2;x3],[zeros(length(x2),1);ones(length(x3),1)],[y;y],...
                  'VariableNames',{'responses','responseFlag','flyLabels'});
            lme1 = fitlme(tbl1,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme2 = fitlme(tbl2,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            lme3 = fitlme(tbl3,'responses~responseFlag+(1|flyLabels)','FitMethod','REML');
            p1 = lme1.Coefficients.pValue(2);
            p2 = lme2.Coefficients.pValue(2);
            p3 = lme3.Coefficients.pValue(2);
            hold on
                bar(1:3,mean([x1 x2 x3],1),'FaceColor','w')
                errorbar(1:3,mean([x1 x2 x3],1),std([x1 x2 x3],[],1)/sqrt(size([x1 x2 x3],1)),'Color','k','LineStyle','none','CapSize',20)
            hold off
            xlim([0 4])
            ylim([0.1 0.5])
            set(gca,'Box','off','TickDir','out','Fontsize',10,'Fontweight','b','LineWidth',1) 
            xlabel('train','Fontsize',12,'Fontweight','b')
            ylabel('sparseness','Fontsize',12,'Fontweight','b')
            title(odorNames(i),'Fontsize',12,'Fontweight','b')
            disp(['         ',odorNames{i},':'])
            disp(['             train 1 vs train 2: p-value = ',num2str(p1),'.'])
            disp(['             train 1 vs train 3: p-value = ',num2str(p2),'.'])
            disp(['             train 2 vs train 3: p-value = ',num2str(p3),'.'])

    end
    disp('     ')
        
end


