%[A nnn nnp npn npp]=roc_ratios(nrm,w1,fa0,pic,shp)
%INPUT ARGUMENTS
%nrm=normalization_flag     nrm=0=raw Yout; nrm=1=z-scores
%w1=integration_window      w=2*w1+1; w1=0=no smoothing
%fa0=fixed_FA_for_each_AU   fa0=0=no fixed F.A.
%pic=graph_flag             0=no graphs; 1=plot ROC curves
%shp=subject_set            shp1=0 or 1 or 2 or 3 or shp=subject

%EG   [A nnn nnp npn npp]=roc_ratios(1,5,.1,0,1);
% zscore interval analysis on shipment 1, false alarm 10%, no graphs


%OUTPUT VARIABLES
%A = area under roc for each AU 
%npp=#hits=TP npn=#false_alarm=FP nnp=#misses=FN nnn=#true_neg=TN   ==>
%hit=npp./(npp+nnp); fa=npn./(nnn+npn); PPV=npp./(npp+npn); NPV=nnn./(nnn+nnp);

%calls window.m  make_X.m  rocA.m  rocA_pic.m

%MPLAB SPECIFIC LINES WHICH NEED TO BE CHANGED
%To make this run w/o mplab access, you need to modify the following
%
%remove lines "sub=".. through line "s=sub(c" .., and instead 
%variables s,sn and strtfrm should be input arguments to the function. 
%startframe will be 1 in future, so no need to remove frames before start
%sn is the 3-string version of s=list of subject numbers
%s and sn are used in: 1) "i=1:size(s,2)" the loop size
% 2) "eval(['load /home/gwen/4FACS_HS/rufacs_out/S',sn(i,:)]);" output_file
% 3) "X=[X;make_X(s(i),g,frm)]; % human labels"
%
%make_X.m not only receives input s, but it also reads au_fr.mat file
%which is mplab specific. This file could be copied to CERT, but there
%needs to be explanation of how to extend it for new subjects.

function [A nnn nnp npn npp]=roc_ratios(nrm,w1,fa0,pic,shp)

  %mplab-specific chunk:next 20 lines
sub=[10 12 14 15 16 17 18 21 23 24 28 30 31 32 45 46 47 49 51 56 60 70 71 ...
     75 77 79 80 86 89 90 94 103];%list of NSF subjects
sname(32,:)='103'; 
for i=1:31  sname(i,:)=['0' num2str(sub(i))];end
shp1=[2:6 11 12 15 18:20 22:25 28:32];%NSF shipment_1
shp2=setdiff([1:32],shp1);%NSF shipment_2
shp_old=[2:6 11 12 15 22:24 28];
startframe(shp2)=[4410 213 117 106 146 152 158 82 164 95 196 165];
startframe(shp1)=1;
if (size(shp,2)==1)
  if (shp==1) choose=shp1; 
  elseif (shp==2) choose=shp2;
  elseif (shp==3 choose=[1:32];
  elseif (shp==0) choose=shp_old; 
  else choose=find(sub==shp);
  end
else choose=find(ismember(sub,shp));
end
s=sub(choose);strtfrm=startframe(choose);sn=sname(choose,:);
  
  
g=[1 2 4 5 10 12 14 20]; %list of au's under study

Y=[];X=[];
for i=1:size(s,2)
  %load file where SVM output for this subject is written
    eval(['load /home/gwen/4FACS_HS/rufacs_out/S',sn(i,:)]);

    % remove frames before start frame
    frm=frm-strtfrm(i)+1;
    ffrm=find(frm>0);
    frm=frm(ffrm);
    Yout=Yout(ffrm,:);
   
  %concatenate variables for all subjects
   z=Yout;
  if (nrm==1)          % use z-scores
    z=zscore(z); 
  end
  if (w1>0)            % interval analysis
    z=window(z,w1,frm);
  end
  Y=[Y;z];             % system output
  X=[X;make_X(s(i),g,frm)]; % human labels
end

% area under ROC curves
 
   if(pic==0)                       % no graphs
     A=rocA(Y,X);
   else                             % with graphs
     A=rocA_pic(Y,X);
   end
 
 %fixed false alarm rate for each AU
 
if (fa0>0)             
  for em=1:size(g,2)
    x=X(:,em);
    nn=size(find(x==-1),1);
    fn=sort(Y(find(x==-1),em));
    th=fn(floor(nn*(1-fa0))); %threshold -> fa=fa0
    yb=sign(Y(:,em)-th);
  
    npp(em)=size(find(yb==1&x==1),1);
    npn(em)=size(find(yb==1&x==-1),1);
    nnp(em)=size(find(yb==-1&x==1),1);
    nnn(em)=size(find(yb==-1&x==-1),1); 
  end

else% fa0=0 no fixed F.A. maximize PPV.*NPV
  
  for em=1:size(g,2)
    x=X(:,em);
    np=size(find(x==1),1);
    if (np>0)
      th_hi=max(Y(find(x==-1),em));
      th_lo=min(Y(find(x==1),em));
      dth=(th_hi-th_lo)/100;  
      for i=1:101
	th=th_lo+(i-1)*dth;
	yb=sign(Y(:,em)-th);
	PPV(i)=size(find(yb==1&x==1),1)/size(find(yb==1),1);
	NPV(i)=size(find(yb==-1&x==-1),1)/size(find(yb==-1),1);
	HIT(i)=size(find(yb==1&x==1),1)/size(find(x==1),1);
	CFA(i)=size(find(yb==-1&x==-1),1)/size(find(x==-1),1);
      end
      L=PPV.*NPV.*HIT.*CFA;%4 reasonable threshold choice
 %     L=HIT.*PPV;
      i=find(L==max(L),1);
      th=th_lo+(i-1)*dth;
      yb=sign(Y(:,em)-th);
      
      npp(em)=size(find(yb==1&x==1),1);
      npn(em)=size(find(yb==1&x==-1),1);
      nnp(em)=size(find(yb==-1&x==1),1);
      nnn(em)=size(find(yb==-1&x==-1),1);      
    end
  end

end
