function gangliaCorrelation(focusfolder)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Trains the logistic regression and the SVM for the focus frame 
% identification
% Inputs:
%     - focusfolder: directory with the focus correlation data
% Output: (none)
%     - 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clearvars -except focusfolder
close all
load(focusfolder);
rng(3)
medlength = 21;
part = 1:15; 

for nn = 1:length(medlength)
    for fn = 1:15
        scorecorr = [];
        scoresnr = [];
        scoressim = [];
        GTm = [];
        scoretestcorr = [];
        scoretestsnr = [];
        scoretestssim = [];
        GTmtest = [];
        for k=1:15
            vcorr = focus_layers(k).correlation;
            ssim = focus_layers(k).ssim;
            snr = focus_layers(k).snr;
            quadloc = focus_layers(k).quadrant;
            focusTable = focus_layers(k).focusTable;
            
            GT = zeros(size(vcorr,2), 4);
            for period = 1:size(focusTable, 1)
                for q = 1:4
                    if focusTable(period, q+2) == 1
                        GT(focusTable(period, 1):focusTable(period, 2), q) = 1;
                    end
                end
            end
    
            nbins = 40;
            bins = linspace(-2,2,nbins+1);
            hgram = zeros(4,2,nbins);
            
            for mask = 1:size(vcorr, 1)
                if ~isnan(quadloc(mask)) && quadloc(mask) ~= 0
                    if all(k ~= part(:, fn))
                        scorecorr = cat(2, scorecorr, movmedian(normalize(vcorr(mask,:),'center','median'),medlength(nn),'omitnan'));
                        scoresnr = cat(2, scoresnr, movmedian(normalize(snr(mask,:),'center','median', 'scale', 'std'),medlength(nn),'omitnan'));
                        scoressim = cat(2, scoressim, movmedian(normalize(ssim(mask,:),'center','median'),medlength(nn),'omitnan'));
                        GTm = cat(1, GTm, GT(:, quadloc(mask)));
                    else
                        scoretestcorr = cat(2, scoretestcorr, movmedian(normalize(vcorr(mask,:),'center','median'),medlength(nn),'omitnan'));
                        scoretestsnr = cat(2, scoretestsnr, movmedian(normalize(snr(mask,:),'center','median', 'scale', 'std'),medlength(nn),'omitnan'));
                        scoretestssim = cat(2, scoretestssim, movmedian(normalize(ssim(mask,:),'center','median'),medlength(nn),'omitnan'));
                        GTmtest = cat(1, GTmtest, GT(:, quadloc(mask)));
                    end
                end
            end
            
        end
        
        X = cat(2, scoressim(:), scorecorr(:));
        Y = GTm(:);
    
        Pind = find(Y == 1);
        Nind = find(Y == 0);
        rNind = randsample(Nind, length(Pind));
        tInd = cat(1, Pind, rNind);
        mdl = fitglm(X(tInd, :),Y(tInd),'Distribution','binomial');
        LRscores = mdl.Fitted.Probability;
        
        Xtest = cat(2, scoretestssim(:), scoretestcorr(:));
        Ytest = GTmtest(:);
    
        ypredtest = predict(mdl,Xtest);
        LogRmodel{fn, nn} = mdl; 
    
        subplot(2, 1, 1)
        [fpr,tpr,T,AUC,opt] = perfcurve(Y(tInd),LRscores, 1, 'Cost', [0 1; 1 0]);
        plot(fpr,tpr,'k-',opt(1),opt(2),'ro')
        hold on
        [fprt,tprt,Tt,AUCt,optt] = perfcurve(Ytest(:),ypredtest,1,'Cost',[0 1; 1 0]);
        plot(fprt,tprt,'r-',optt(1),optt(2),'ko')
        axis equal
        title('Logistic Regression')
        xlabel('False positive rate')
        ylabel('True positive rate')
        axis([0, 1, 0, 1])
        legend({sprintf("AUC: %.2f\nThreshold: %.2f",AUC,T((fpr==opt(1))&(tpr==opt(2)))),sprintf("AUC: %.2f\nThreshold: %.2f",AUCt,Tt((fprt==optt(1))&(tprt==optt(2))))},'Location','SouthEast')
        subplot(2, 1, 2)
        plot(Ytest, 'k')
        hold on
        plot(ypredtest, 'r')
        yline(T((fpr==opt(1))&(tpr==opt(2))))
        saveY{fn} = Ytest;
        savepred{fn, nn} =ypredtest;
        threshold(fn, nn) = T((fpr==opt(1))&(tpr==opt(2)));
        saveX{fn} = Xtest;
        TP(fn, nn) = sum(ismember(find(Ytest), find(ypredtest > threshold(fn))));
        TN(fn, nn) = sum(ismember(find(Ytest == 0), find(ypredtest < threshold(fn))));
        FP(fn, nn) = sum(~ismember(find(ypredtest > threshold(fn)), find(Ytest)));
        FN(fn, nn) = sum(~ismember(find(ypredtest < threshold(fn)), find(Ytest == 0)));
        P(fn, nn) = length(find(Ytest));
        N(fn, nn) = length(find(Ytest == 0));
        Acc(fn, nn) = (TP(fn, nn)+TN(fn, nn))/(P(fn, nn)+N(fn, nn));
        phi(fn, nn) = (TP(fn, nn)*TN(fn, nn)-FP(fn, nn)*FN(fn, nn))/sqrt((TP(fn, nn)+FP(fn, nn))*(TP(fn, nn)+FN(fn, nn))*(TN(fn, nn)+FP(fn, nn))*(TN(fn, nn)+FN(fn, nn)));
    
        clf
    end
end

bar([1 2 3 4], [mean(phi) mean(Acc) mean(TN./N) mean(TP./P)], 'FaceColor', [0.8 0.8 0.8])
hold on
plot(ones(size(phi)), phi, 'ko')
errorbar(1, mean(phi), std(phi), 'k', 'LineWidth', 2)
plot(ones(size(Acc))*2, Acc, 'ko')
errorbar(2, mean(Acc), std(Acc), 'k', 'LineWidth', 2)
plot(ones(size(Acc))*3, TN./N, 'ko')
errorbar(3, mean(TN./N), std(TN./N), 'k', 'LineWidth', 2)
plot(ones(size(Acc))*4, TP./P, 'ko')
errorbar(4, mean(TP./P), std(TP./P), 'k', 'LineWidth', 2)
xticklabels({'Matthews correlation coefficient', 'Accuracy', 'True negative rate', 'True positive rate'})
ylabel('Value')


save('20220829_segmentationcorrelation_quadrant_additionalmetrics_trainedmodel.mat', 'LogRmodel', ...
    'savepred', 'saveY', 'threshold', 'saveX', 'TP', 'TN', 'FP', 'FN', 'P', 'N', 'Acc', 'phi');

pause

clearvars -except focusfolder
close all

load(focusfolder);
rng(3)
medlength = 21;

scorecorr = [];
scoresnr = [];
scoressim = [];
GTm = [];

scorediff = [];
scoreIOU = [];
nn =1;

for k=1:15
    vcorr = focus_layers(k).correlation;
    ssim = focus_layers(k).ssim;
    snr = focus_layers(k).snr;
    quadloc = focus_layers(k).quadrant;
    focusTable = focus_layers(k).focusTable;
    sdiff = focus_layers(k).diff;
    sIOU = focus_layers(k).IOU;
    order = [2 1 3 4];
    GT = zeros(size(vcorr,2), 4);
    for period = 1:size(focusTable, 1)
        for q = 1:4
            if focusTable(period, q+2) == 1
                GT(focusTable(period, 1):focusTable(period, 2), q) = 1;
            end
        end
    end
    for mask = 1:size(vcorr, 1)
        if ~isnan(quadloc(mask)) && quadloc(mask) ~= 0
                scorecorr = cat(2, scorecorr, movmedian(normalize(vcorr(mask,:),'center','median'),medlength(nn),'omitnan'));
                scoresnr = cat(2, scoresnr, movmedian(normalize(snr(mask,:),'zscore'),medlength(nn),'omitnan'));
                scoressim = cat(2, scoressim, movmedian(normalize(ssim(mask,:),'center','median'),medlength(nn),'omitnan'));
                scorediff = cat(2, scorediff, movmedian(normalize(sdiff(mask,:)./sIOU(mask,:),'zscore'),medlength(nn),'omitnan'));
                scoreIOU = cat(2, scoreIOU, movmedian(normalize(sIOU(mask,:),'center','median'),medlength(nn),'omitnan'));
                GTm = cat(1, GTm, GT(:, quadloc(mask)));
        end
    end  
end
    
    
X = cat(2,scoressim(:), scorecorr(:));
Y = GTm(:);

Pind = find(Y == 1);
Nind = find(Y == 0);
rNind = randsample(Nind, length(Pind));
tInd = cat(1, Pind, rNind);
mdl = fitglm(X(tInd, :),Y(tInd),'Distribution','binomial');
LRscores = mdl.Fitted.Probability;


[fpr,tpr,T,AUC,opt] = perfcurve(Y(tInd),LRscores, 1, 'Cost', [0 1; 1 0]);
plot(fpr,tpr,'k-',opt(1),opt(2),'ro')
hold on
axis equal
title('Logistic Regression')
xlabel('False positive rate')
ylabel('True positive rate')
axis([0, 1, 0, 1])
legend({sprintf("AUC: %.2f\nThreshold: %.2f",AUC,T((fpr==opt(1))&(tpr==opt(2))))},'Location','SouthEast')
clf