function [out] = SingClassROC(inStructFusion2,Label,varargin)
% % function [out] = SingClassROC(inStructFusion2,Label,varargin)

% FUNCTION 3
% This is the third function (3 out of 3) for one-class classification that
% classifies the test sample by selecting a SRD threshold and a window size
% from a Reciever Operatero Characteristic (ROC) curve
% - The window size is determined by the maximum area under the ROC curve
% - The SRD threshold is determined by finding the closes point to [0,1]
% on the ROC curve


%{

Inputs:
    inStructGenFus2 : struct file generated from the output of SingClassFusion.m
    Label: A cell array to put a label on you class of interest.
           Ex: If your class of interest is chicken, then you would enter
           Label as: Label{1} = 'Chicken';
Outputs:
      Figures:
              - Area under the curve
              - ROC curves for all windows
              - ROC curves for max and optimal window only
              - Single class classifiction result for max and optimal
                windows (at their respective SRD thesholds)
              - Distribution of SRD ranks for all test samples at max and
                optimal window
     
     Table containing the following information at max and optimal threshold:
              - Window
              - SRD Rank Threshold
              - Area Under the Curve
              - Sensitivity
              - Specificity
              - Accuracy


/// Created by Brett Brownfield & Tony Lemos August 1, 2018 (Last edited on 7/11/19)
/// Department of Chemistry, Idaho State University, Pocatello, Idaho 83209 USA

%}


if nargin < 2
    Label{1} = 'Class of Interest';
end
% Obtaining information from inputStruct
inputStruct=inStructFusion2.input;
inputStruct=rmfield(inputStruct, 'sig');

outfun.input=inputStruct;
nC = size(inputStruct.Xin,2);

% Obtaining results from different threshold based on the SRD ranking

% Threshold in intervals of 1
Thresh = 0:1:100;
lengThres = length(Thresh);
MaxWind = inputStruct.tuneMax(1); % All Threshold will be the same

% Preallocation
valClass = cell(1,lengThres);
ytrue =  cell(1,lengThres);
FOM = cell(1,lengThres);
seOA = nan(MaxWind,lengThres);
spOA = nan(MaxWind,lengThres);
oneMinspOA = nan(MaxWind,lengThres);
AreaUnCurv = nan(1,MaxWind);
DistOpt = nan(1,lengThres);
DistMax = nan(1,lengThres);

% Obtaining the classifiction perforance values at every threshold
for i = 1:lengThres
    % Performing classification at the given threshold 
    [valClass{i}, ytrue{i},valRanking]=classSingClassify(inStructFusion2.fusOutSing, Thresh(i), inputStruct);
    
    outfun.valClass=valClass{i};
    outfun.ytrue=ytrue{i};
    
    [FOM{i}] = SingClassFOM(outfun);
    
    % Obtaining Specificity and Sensitivity (and transforming it to be between 0 and 1)
    seOA(:,i) = FOM{i}.stdFOM.seOA./100;
    spOA(:,i) = FOM{i}.stdFOM.spOA./100;
    oneMinspOA(:,i) = 1-spOA(:,i);
end





% ROC curves
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

% Re-arranging seOA and spOA so that the area under the curve is not negative
seOA2 = fliplr(seOA);
oneMinspOA2 = fliplr(oneMinspOA);

% Plotting the ROC at every window
figure,  hold, xlabel('1 - Specificity'), ylabel('Sensitivity'), title('ROC: All Windows'), grid on
for i = 1:MaxWind
    plot(oneMinspOA2(i,:),seOA2(i,:),'-','linewidth',2)
    AreaUnCurv(i) = trapz(oneMinspOA2(i,:),seOA2(i,:)); % Finding the area under the curve
end

OptWind = find(max(AreaUnCurv)==AreaUnCurv);

% Removing other windows that give you the same results
if length(OptWind)>1, OptWind(1:end-1) = []; end

% Plotting the area under the curve across all windows
h = figure;
hStairs = stairs(AreaUnCurv,'Linewidth',2);
xlabel('Window Size'), title('AUC'), ylabel('Area'), xlim([1 MaxWind]), grid on;


% Putting points on the optimal and maximum window
cursorMode = datacursormode(h);
hDatatip = cursorMode.createDatatip(hStairs);
pos1 = [OptWind AreaUnCurv(OptWind) 0];
set(hDatatip,'Position',pos1);

% Setting the second data tip
cursorMode = datacursormode(h);
hDatatip = cursorMode.createDatatip(hStairs);
pos2 = [MaxWind AreaUnCurv(MaxWind) 0];
set(hDatatip,'Position',pos2);

% Plotting the optimal and max window ROC curve
figure, hold, grid on, xlabel('1 - Specificity'), ylabel('Sensitivity'), title('ROC')
plot(oneMinspOA(OptWind,:),seOA(OptWind,:),'-*','linewidth',2)
plot(oneMinspOA(MaxWind,:),seOA(MaxWind,:),'-*','linewidth',2)
legend({'Optimal Window','Max Window'},'location','best')


%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~




% Obtaining optimal threshold by calculating the distance between he ideal ROC
% point (0,1) and observed ROC points (Using Euclidean Distance)
xPoint = [0; 1];
for i = 1:lengThres
    x = [oneMinspOA(OptWind,i); seOA(OptWind,i)];
    x2 = [oneMinspOA(MaxWind,i); seOA(MaxWind,i)];
    DistOpt(i) = sqrt((xPoint-x)'*(xPoint-x));
    DistMax(i) = sqrt((xPoint-x2)'*(xPoint-x2));
end

OptThres = find(min(DistOpt) == DistOpt); % Optimal Threshold for optimal window
MaxThres = find(min(DistMax) == DistMax); % Max Threshold for maximum window

% In case there is more than one threshold that yield the same distance
%  - If so, choose the largest threshold
if length(MaxThres)>1, MaxThres(1:end-1) = []; end
if length(OptThres)>1,  OptThres(1:end-1) = []; end

% Plotting the the chosen thresholds on the ROC Curve
plot(oneMinspOA(OptWind,OptThres),seOA(OptWind,OptThres),'pb','linewidth',3,'markersize',.5)
plot(oneMinspOA(MaxWind,MaxThres),seOA(MaxWind,MaxThres),'pr','linewidth',3,'markersize',.5)
legend('Optimal Window','Max Window','Optimal Threshold','Max Threshold')




% Plotting Classification peformance and SRD rank distributions
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Optimal threshold
SingClassPlots(FOM{1,OptThres}.stdFOM.seOA,FOM{1,OptThres}.stdFOM.spOA,FOM{1,OptThres}.stdFOM.accOA,...
    Thresh(OptThres),OptWind,MaxWind,nC,valRanking,Label)

% Maximum threshold
SingClassPlots(FOM{1,MaxThres}.stdFOM.seOA,FOM{1,MaxThres}.stdFOM.spOA,FOM{1,MaxThres}.stdFOM.accOA,...
    Thresh(MaxThres),MaxWind,MaxWind,nC,valRanking,Label)
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~




% Obtaining value to make the final table
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Optimal_Threshold(1,1) = OptWind;
Max_Threshold(1,1) = MaxWind;

Optimal_Threshold(2,1) = Thresh(OptThres);
Max_Threshold(2,1) = Thresh(MaxThres);

Optimal_Threshold(3,1) = AreaUnCurv(OptWind);
Max_Threshold(3,1) = AreaUnCurv(MaxWind);

Optimal_Threshold(4,1) = FOM{1,OptThres}.stdFOM.seOA(OptWind);
Max_Threshold(4,1) = FOM{1,MaxThres}.stdFOM.seOA(MaxWind);

Optimal_Threshold(5,1) = FOM{1,OptThres}.stdFOM.spOA(OptWind);
Max_Threshold(5,1) = FOM{1,MaxThres}.stdFOM.spOA(MaxWind);

Optimal_Threshold(6,1) = FOM{1,OptThres}.stdFOM.accOA(OptWind);
Max_Threshold(6,1) = FOM{1,MaxThres}.stdFOM.accOA(MaxWind);

T = table(Optimal_Threshold,Max_Threshold);

T.Properties.RowNames = {'Window','SRD Rank Threshold','Area Under the Curve','Sensitivity','Specificity','Accuracy'};

out.ClassificationTable = T;




end

function  SingClassPlots(se,sp,acc,Threshold,UsedWind,MaxWind,nC,valRanking,Label)
% This function plot the single class classificiation pefromance (accuracy,
% sensitivity and specificity) as well as the SRD rank distribution for the
% test samples

h = figure;
hold,
plot(se,'b','LineWidth',4)
plot(sp,'g','LineWidth',4)
hPlot = plot(acc,'o-r','linewidth',4,'MarkerFaceColor','r','MarkerSize',6);
legend('Sensitivity', 'Specificity', 'Accuracy', 'Location','best')
title([Label{1},' as Target class: SRD Rank Threshold ', num2str(Threshold)])
grid on
a=gca;
a.YLim = [0 100];
a.XLim = [1 MaxWind];
xlabel('Window size')
dim = [.35 0 .17 .19]; 
str = ['Optimal Window Size: ',num2str(UsedWind)];
annotation('Textbox',dim,'String',str,'FitBoxToText','on')
a.FontSize = 14;

% Putting point on the optimal window
cursorMode = datacursormode(h);
hDatatip = cursorMode.createDatatip(hPlot);

pos = [UsedWind acc(UsedWind) 0];
set(hDatatip,'Position',pos);


%   Obtaining the SRD ranking distribution 

figure, hold
for i = 1:nC
    SIG = cell2mat(valRanking(:,i));
    h = histogram(SIG(:,UsedWind));   % Only getting sigma values for the used Window
    h.BinWidth =  (100)/32;
end
ax = gca;
ax.FontSize = 18;
a = ax.YLim;


plot([Threshold Threshold], a, '-r', 'linewidth',3)
legend({'Target','Non-target',['SRD Rank ',sprintf('%4.2f ',Threshold)]},'Location','best')
title(['Distribution of SRD rankings [',num2str(UsedWind),' window] ']); grid on
xlabel('SRD Rankings')

end

function [valClass, yTrue,valRanking]=classSingClassify(fusOutSing, Thresh, inputStruct)
% This funcion classifies each sample using the sigma threshold

struct2var(inputStruct);
nC=size(Xin,2);

nsamp = [nOutPure, nInNonPure];

if size(Xin,2) > length(nsamp)
    extra  =  (size(Xin,2) - length(nsamp));
    nsamp = [nsamp repmat(nsamp(end),1,extra)];
end

valClass=cell(nCV,nC);
valRanking=cell(nCV,nC); % SRD ranking of val
yTrue=logical([1 0]);

for i=1:nCV
    
    for ii=1:nC
        
        for j=1:nsamp(ii)
            
            for k=1:min(tuneMax)
                
                thres = Thresh; % SRD Threshold
               
                %Classifying Sample
                if ii==1
                    rankVal=fusOutSing{i,ii}(j,k).norm(end);% only works when srdcv=false  
                    valRanking{i,ii}(j,k) = (rankVal);
                else
                    rankVal=fusOutSing{i,ii}(j,k).norm(end); % only works when srdcv=false
                    valRanking{i,ii}(j,k) = rankVal;
                end
                
                if rankVal<thres
                    valClass{i,ii}(j,k)=false;
                else
                    valClass{i,ii}(j,k)=true;
                end %if
            end
        end % nsamp
    end % nC
end

end

function outfun= SingClassFOM(inStructFusion)
% Thisfunction calulates the figure of merits 
% of the one-class classification 
 


%{
Input:
    inStructFusion: the struct generated by SingClassClassify.m

Output:
    input: The inputs used to run one-class classification (see
           SingClassfusIn.m for more)
    stdFOM:(standard figure of merits) Contains the standard figure of
           merits for each cross validation split at each tuning parameter 
           window
        acc - Accuarcy for each cross validation split at each window
        se - Sensitivity for each cross validation split at each window
        sp - Specificity for each cross validation split at each window
        accOA - Overall accuracy (mean of accuracy across each cross 
                validation split)
        seOA - Overall sensitivity (mean of sensitivity across each cross
               valdiation split
        spOA - Overall specificity (mean of the specificity across each
               cross validation split
    
    outT: Shows the true postives, true negatives, etc. for each cross
          validation split at each tuning parameter window.
            TP: The number of true postives 
            TN: The number of true negatives
            FP: The number of false positives
            FN: The number of false negatives
%}
struct2var(inStructFusion);
inputStruct=inStructFusion.input;
struct2var(inputStruct);
[~, nC]=size(Xin);


% Checking if samples are classified correctly and calculating figures of merit for classification
[stdFOM, outT]=classcheckSingMrt(valClass, ytrue, nC, min(tuneMax), nCV);

% Saving outputs
outfun.input=inputStruct;
outfun.stdFOM=stdFOM;
outfun.outT=outT;

end

function [stdFOM, outT]=classcheckSingMrt(valClass, ytrue, nC, tMM, nCV)
% This function calculates the accuracy, senstivity, and selectivity for
% classification

if size(valClass,2) > 2
    reshp = size(valClass,2);   
    for h = 1:nCV
        for hh = 3:reshp
            valClass{h,2} = [valClass{h,2}; valClass{h,hh}];
        end
    end
    nC = 2;
end
%preallocation
TP=zeros(tMM, nCV);
TN=zeros(tMM, nCV);
FP=zeros(tMM, nCV);
FN=zeros(tMM, nCV);

for i=1:nCV
    for ii=1:nC
        [nsamp]=size(valClass{i,ii},1);
        for j=1:nsamp
            for k=1:tMM
                if valClass{i, ii}(j, k)==ytrue(ii) && ytrue(ii)==1
                    TP(k, i)=TP(k, i)+1; % 1 TP for correct classifcation
                elseif valClass{i,ii}(j, k)==ytrue(ii) && ytrue(ii)==0
                    TN(k, i)= TN(k, i)+1; % 1 True negative for correctly idenfity sample not in class
                elseif valClass{i,ii}(j, k)~=ytrue(ii)
                    if valClass{i,ii}(j,k)==0 && ytrue(ii)==1 % lives in Base but classified as not
                        FN(k, i)= FN(k, i)+1; % 1 FN for the class where it is supposed to be classified as
                    elseif valClass{i,ii}(j,k)==1 && ytrue(ii)==0 % Doesn't live in Base but classified into Base
                        FP(k, i)= FP(k, i)+1; % 1 FP in the class where it was incorrectly classified
                    end
                else
                    error('You have a serious problem')
                end
            end
        end
    end
end

outT.TP=TP;
outT.TN=TN;
outT.FP=FP;
outT.FN=FN;

acc=(TP+TN)./(TP+TN+FP+FN)*100; %accuracy
se=(TP)./(TP+FN)*100; %sensetivity
sp=(TN)./(TN+FP)*100; %specificity
stdFOM.acc=acc;
stdFOM.se=se;
stdFOM.sp=sp;
%Overall (OA) figure of merits across ALL classes
TPsum=sum(TP,2);
TNsum=sum(TN,2);
FPsum=sum(FP,2);
FNsum=sum(FN,2);
stdFOM.accOA=(TPsum+TNsum)./(TPsum+TNsum+FPsum+FNsum)*100;
stdFOM.seOA=(TPsum)./(TPsum+FNsum)*100;
stdFOM.spOA=(TNsum)./(TNsum+FPsum)*100;



end % function
