(*:Mathematica:: Version 7*)
(*:Context: "swirl`"         *)
(*:Swirl-and-Ripple       *)
(*Compute ROC and concave envelope*)
(*:Author: 2010 Stuart G. Baker *)



BeginPackage["swirlroc`","swirlsup`"]

PlotROCGoal1::usage="PlotROC[fset,tset,siftname]"
PlotROCGoal2::usage="PlotROC[fset,tset,siftname]"

ComputeROCTestSample::usage="ComputeROCTestSample[data,{posIN,posOUT,tlist,xm0,xm1,vm0,vm1},{score,distance},fset]"
ComputeROCTestBoot::usage="ComputeROCTestBoot[dataTEST,resTRAIN,{score,distance},fset,maxboot]"

ROCToRU::usage="ROCToRU[roc,pi]"


Clear[ComputeROCTestSample,PlotROC,ROCToRU,ComputeROCTestBoot] 


Begin["Private`"]


 (*-------------------PlotROC--------------*)

 PlotROCGoal1[fset_,tset_,siftname_,pi_,dname_,cutvec_]:=
	Module[{rmat,tlow,tupp,tmed,
	  rocmed0,roclow0,rocupp0,
	  rocmed,roclow,rocupp,plotroc,plotruc},
  (*compute TPR for FPR which are fixed*)
       tmat=Transpose[tset];
       tlow=LowerB[#]& /@tmat;
       tupp=UpperB[#]& /@tmat;
       tmed=Mean[#]& /@tmat;
   (*ROC curve includes (0,0) and (1,1)*)
      rocmed0=ROCAdd0011[fset,tmed];
      roclow0=ROCAdd0011[fset,tlow];
      rocupp0=ROCAdd0011[fset,tupp];
      plotroc=PlotROCCore[rocmed0,roclow0,rocupp0,siftname,dname];
   (*RU curve based on concave ROC*)
        rocmed=ConcaveROC@rocmed0;
        {fpr0,tpr0}=Transpose[rocmed0];
        {fpr,tpr}=Transpose[rocmed];
       pair0=Transpose[{fpr0,cutvec}];
      pair1=Select[pair0,MemberQ[fpr,#[[1]]]&];
      {fpr1,cutvec1}=Transpose[pair1];
      plotruc=PlotRUGoal1[rocmed,pi,siftname,dname,cutvec1];
    Return[{plotroc,plotruc}]]

PlotRUGoal1[rocmed_,pi_,siftname_,dname_,cutvec_]:=
	  Module[{lab,rumedpair,rulowpair,ruupppair,rurange,tickname,  plotrumedx,
   		plotrumed,plotrulow,plotruuup,plotru},
	       rumedpair=ROCToRU[rocmed,pi];	
	       {r,ru}=Transpose[rumedpair];
	       cutveclow=Drop[cutvec,-1];
	       cutvecupp=Drop[cutvec,1];
	       mat={ru,r,cutveclow,cutvecupp};
	       matr=Round[mat 100]/100//N;
	       Print["     ",TableForm[matr, TableHeadings->{{"RU","R","cut-low","cut-upp"},None}]];
	       rurange={Automatic,{0,1}};
    	       plotrumed=PlotRUCore[rumedpair,siftname,dname,rurange,{Black}];
	       plotru=Show[plotrumed];
	 Return[plotru]]




PlotROCGoal2[fset_,tset_,siftname_,pi_,dname_]:=
	Module[{rmat,tlow,tupp,tmed,
	  rocmed0,roclow0,rocupp0,
	  rocmed,roclow,rocupp,plotroc,plotruc},
  (*compute TPR for FPR which are fixed*)
       tmat=Transpose[tset];
       tlow=LowerB[#]& /@tmat;
       tupp=UpperB[#]& /@tmat;
       tmed=Mean[#]& /@tmat;
   (*ROC curve includes (0,0) and (1,1)*)
      rocmed0=ROCAdd0011[fset,tmed];
      roclow0=ROCAdd0011[fset,tlow];
      rocupp0=ROCAdd0011[fset,tupp];
      plotroc=PlotROCCore[rocmed0,roclow0,rocupp0,siftname,dname];
   (*RU curve based on concave ROC*)
        rocmed=ConcaveROC@rocmed0;
       plotruc=PlotRUGoal2[rocmed,pi,siftname,dname];
    Return[{plotroc,plotruc}]]



	PlotROCCore[rocmed_,roclow_,rocupp_,siftname_,dname_]:=
	Module[{lab,ticknamex,tickname,plotrocmed,plotroclow,plotrocupp,plotroc},
	  ticknamex={.2,.4,.6,.8,1};
	  lab=StringJoin["ROC: ",dname,": ",siftname];
	  tickname={ticknamex,ticknamex};
	  plotrocmed=ListPlot[rocmed,Joined->True,
	      PlotLabel->lab,
			PlotStyle->{Black},
			AxesLabel->{"FPR","TPR"},
			AxesOrigin->{0,0},
			Ticks->tickname];
            plotroclow=ListPlot[roclow,Joined->True,
     	      PlotStyle->{Black,Dashing[{.05,.05}]},
  	       PlotLabel->lab, 
			 AxesLabel->{"FPR","TPR"},
			 AxesOrigin->{0,0},
			 Ticks->tickname];
     plotrocupp=ListPlot[rocupp,Joined->True,
	       PlotStyle->{Black,Dashing[{.05,.05}]},
	       PlotLabel->lab,
			 AxesLabel->{"FPR","TPR"},
			 AxesOrigin->{0,0},
			 Ticks->tickname];
      plotroc=Show[plotrocmed,plotroclow,plotrocupp];
	  Return[plotroc]]

	
	  
	PlotRUGoal2[rocmed_,pi_,siftname_,dname_]:=
	  Module[{lab,rumedpair,rulowpair,ruupppair,rurange,tickname,  plotrumedx,
   		plotrumed,plotrulow,plotruuup,plotru},
	       rumedpair=ROCToRU[rocmed,pi];	
	       Print[rumedpair];
               rurange={Automatic,{0,1}};
    	       plotrumed=PlotRUCore[rumedpair,siftname,dname,rurange,{Black}];
	       plotru=Show[plotrumed];
	 Return[plotru]]


	 PlotRUCore[rupair_,siftname_,dname_,rurange_,style_]:=
	 Module[{len,plotru,lab},				 
           lab=StringJoin["RU: ",dname,": ",siftname];
	  len=Length@Transpose[rupair][[1]];
	  If[len>2,
	  plotru=ListPlot[rupair,Joined->True,
	      PlotLabel->lab,
			PlotRange->rurange,
			PlotStyle->{Black},
			AxesLabel->{"risk","RU"},
			AxesOrigin->{0,0}],
	 (*else*)
	   plotru0=ListPlot[{{0,0},{1.1,1.1}},
	 	      PlotLabel->lab,
	 			PlotRange->rurange,
	 			PlotStyle->{Black},
	 			AxesLabel->{"risk","RU"},
	 			AxesOrigin->{0,0}];
	plotru1=Graphics@Text["insufficent",{.3,.7},{-1,0}];
	plotru2=Graphics@Text["data",{.3,.4},{-1,0}];
	plotru=Show[plotru0,plotru1,plotru2]];
    Return[plotru]]



  LowerB[y_]:=Max[0,Mean[y] -1.96 Sqrt[Var[y]]]
 
  UpperB[y_]:=Min[1,Mean[y] + 1.96 Sqrt[Var[y]]]


  ROCAdd0011[fset_,tset_]:=
  Module[{roc1,roc2},
    roc1=Transpose[{fset,tset}];
    roc2=Join[{{0,0}},roc1,{{1,1}}];
  Return[roc2]]

  


 (*----------------ROCToRU-----------------*)


ROCToRU[roc_,pi_]:=
Module[{slope0,slope,z,r,f,t,t0,t1,w,rocz,x,n},
  (*slopes s1 s2 s3 correspond to risks r3, r2, r1*)
       slope0=ROCSlopeX[roc];
  (*orders risks as r1 r2 r3*)
       slope=Reverse[slope0];
 (*compute risks*)
	  z=slope pi / (1-pi);
          r= (z/(1+z));
  (*compute weights*)
     {f,t}=Transpose[roc];
      t0=Drop[t,-1];
      t1=Drop[t,1];
      w=Reverse[WeightComp[#,pi]& /@ Transpose[{t0,t1,Reverse[r]}]];
  (*----------COMPUTE ROC---------*)
      {fprx,tprx}=Transpose[roc];			   
   (*----------COMPUTE RU---------*)
  (*reverse order from ROC because risk is reversed *)
         fpr=Reverse[Drop[fprx,1]];
	  tpr=Reverse[Drop[tprx,1]];    
	   tri=Transpose[{r,fpr,tpr}];	
  (*different formulas depend on prevalence*)
     tri0=Select[tri,(#[[1]]<pi)&];
	  If[Length[tri0]==0,
	   Print["warning in computing RU curve"];
		{r0,fpr0,tpr0}={{0},{0},{0}},
	  {r0,fpr0,tpr0}=Transpose[tri0]];
  (*greater than prevalence*)  
	  tri1=Select[tri,(#[[1]]>=pi)&];
            {r1,fpr1,tpr1}=Transpose[tri1];
  (*slopes*)	 
  	  s0 =((1-pi)/pi)(r0 /(1-r0));
            s1 =((1-pi)/pi)(r1 /(1-r1));  
	  sc=Join[s0,s1];
  (*relative utilities*)
        ru0=(WeightRU0[#]& /@ Transpose[{fpr0,tpr0,s0}]); 
	 ru1=tpr1- s1 fpr1;		 
	 ruc=Join[ru0,ru1];
 (*risks*)
	  rc= Join[r0,r1];
 	  pairruc=Transpose[{rc,ruc}];
  Return[pairruc]]
  
	  
  WeightComp[{t0_,t1_,r_},pi_]:=
      Module[{w},
	   If[r==0,w=0,w=(t1-t0) pi /r];
	  Return[w]]
  
  WeightRU0[{fpr0_,tpr0_,s0_}]:=
	  Module[{ru0},
	  If[s0==0,   ru0=0,ru0=(1-fpr0)- (1- tpr0)/s0];
    Return[ru0]]

ROCSlopeX[roc_]:=
 Module[{pair2,res},
  pair2=Partition[roc,2,1];
  res=ROCSlopePairX[#]& /@ pair2;
 Return[res]]

ROCSlopePairX[x_]:=
 Module[{fpr0,tpr0,fpr1,tpr1,slope},
 {{fpr0,tpr0},{fpr1,tpr1}}=x;
 slope=(tpr1-tpr0)/(fpr1-fpr0);
Return[slope]]





(*-------------------------ROC Bootstrap for Goal 1------------------------------------*)
 
 ComputeROCTestBoot[dataTEST_,resTRAIN_,{score_,distance_},fset_,maxboot_]:=
  Module[{mat,i},
     mat=Table[ComputeROCTestBootOne[dataTEST,resTRAIN,{score,distance},fset],{i,1,maxboot}];
	  res=Transpose[mat];
	 Return[res]]

 ComputeROCTestBootOne[dataTEST_,resTRAIN_,{score_,distance_},fset_]:=
	Module[{x0,x1,p0,n0,p1,n1,ranvec0,ranvec1,x0R,x1R,dataTESTR,ROCsetTESTR},		 
          {x0,x1}=dataTEST;
			 {p0,n0}=Dimensions[x0];
			 {p1,n1}=Dimensions[x1];
		    ranvec0=RandomInteger[{1,n0},n0];
			 ranvec1=RandomInteger[{1,n1},n1];
			 x0R=SelectColumns[x0,ranvec0];
			 x1R=SelectColumns[x1,ranvec1];
			 dataTESTR={x0R,x1R};
			 ROCsetTESTR=ComputeROCTestSample[dataTESTR,resTRAIN,{score,distance},fset,False];
		Return[ROCsetTESTR]]


(*------------------------Compute ROC of Score in Test Sample------------------------------------------*)

ComputeROCTestSample[{x0_,x1_},{posIN_,posOUT_,tlistx_,xm0_,xm1_,vm0_,vm1_},{score_,distance_},fset_,show_:False]:=
Module[{xtestt0,xtest1,d00,d01,d10,d11,d00sq,d01sq,d10sq,d0sum,d1sum,r0,r1,plotroc,ROCset}, 
   (*selected genes from training sample*)
      xtest0=x0[[posIN]];
      xtest1=x1[[posIN]];
	(*distance measures of selected genes from training sample*)
		{d00sq,d01sq,d10sq,d11sq}=SumOfDistanceSquaredTest[{xtest0,xtest1},{xm0,xm1,vm0,vm1},distance];
	If[score=="Swirl",
	    d00=Sqrt[d00sq];
		 d01=Sqrt[d01sq];
		 d10=Sqrt[d10sq];
		 d11=Sqrt[d11sq];	
		 d0sum=ReplaceZero[d00+d01];
        	    d1sum=ReplaceZero[d10+d11];
        	    r0=d00/d0sum;	 (*probability classify x0 as 0 = true positive rate*)
                 r1=d10/d1sum]; (*probability classify x1 as 0 = false positive rate*)
	If[score=="Ripple",
        	 r0=d00sq-d01sq;
                 r1=d10sq-d11sq];
    	  ROCset=ComputeROC[r0,r1,fset];
Return[ROCset]]

SumOfDistanceSquaredTest[{x0_,x1_},{xc0_,xc1_,vm0_,vm1_},distance_:1]:=
Module[{p0,n0,p1,n1,v0,v1,vm,xc00,xc01,xc10,xc11,d00sq,d01sq,d10sq,d11sq},
 (*p is genes, n is subjects*)
   {p0,n0}=Dimensions[x0];
   {p1,n1}=Dimensions[x1];
(*variance*)
   If[distance==1,   vm=(vm0 (n0-1) + vm1 (n1-1) )/(n0+n1-2);   v0=vm;      v1=vm];	   
   If[distance==2,                                              v0=vm0;     v1=vm1];								  																	 
(*distance measures square root of sum over squares of genes: vector for specimens*)  
   d00sq= Apply[Plus,(x0-xc0)^2/v0];		    (*distance of x0 to centroid 0*)
   d01sq= Apply[Plus,(x0-xc1)^2/v1];		    (*distance of x0 to centroid 1*)
   d10sq= Apply[Plus,(x1-xc0)^2/v0];		    (*distance of x1 to centroid 0*)
   d11sq= Apply[Plus,(x1-xc1)^2/v1];    	    (*distance of x1 to centroid 1*)
Return[{d00sq,d01sq,d10sq,d11sq}]];


ComputeROC[pred0_,pred1_,fset_]:=
  Module[{kset,cutvec, f,plotRR,
         fprvec,tprvec,rocvec0,rocvec,rocvec1,rocvec2,
	 cutvecR,fprvecR,tprvecR,rocvec0R,rocvecR,rocvec1R,rocvec2R,			
	rocvec3,auc,aucR,ROCset},
    (*cutpoints*)
	       cutvec=Quantile[Join[pred0,pred1],fset];
      (*ROC points larger values classified as 1 *)
	    fprvec=rocx[pred0,cutvec];
            tprvec=rocx[pred1,cutvec];
	    rocvec0=Transpose[{fprvec,tprvec}];
 	    rocvec1=Join[{{1,1}},rocvec0,{{0,0}}];
 	    rocvec2=Reverse[rocvec1];
    (*computation of AUC, TPR's*)										   
            ROCset=GenROCSet[rocvec2,fset];  
 Return[ROCset]]				  



 rocx[pred_,cutvec_]:=(roci[pred,#]& /@ cutvec)/Length[pred]//N;

 roci[pred_,cut_]:=  Length@ Select[pred, (#>=cut)&]


 (*------------------------------AUC and TPR---------------------------------*)

  GenROCSet[rocvec1_,fset_]:=
  Module[{auc,t0,t1,fpr0,fpr1,pair},
    auc=AUC[rocvec1];
	 fpr0=fset[[1]];
	 fpr1=fset[[2]];
	 (*tpr for first fpr*)
	 t0=CompTPR[rocvec1,fpr0];
	 (*tpr for second fpr*)
	 t1=CompTPR[rocvec1,fpr1];
	 tset=CompTPR[rocvec1,#]& /@fset;
  Return[{auc,t0,t1,tset}]]



  CompTPR[rocvec_,f_]:=
  Module[{pairset,pair,f0,t0,f1,t1},
   pairset=Partition[rocvec,2,1];
   pair=Select[pairset,((#[[1,1]] < f) && (#[[2,1]] >= f))&];
   {{f0,t0},{f1,t1}}=pair[[1]];
	If[Abs[f1-f0]>0,
	 slope=(t1-t0)/(f1-f0);
    t= t0 + slope (f-f0),
	 t=t0; Print["ROC slope not defined ",pair]];
  Return[t]]

 
 
AUC[roc_]:=
Module[{fprvec,trpvec,tmvec,auc},
  (*fpr and tpr*)
     {fprvec,tprvec}=Transpose[roc];
  (*sum dif fpr x mean tpr*)
     wvec=dif[#]& /@ Partition[fprvec,2,1];
	  tmvec=Mean[#]& /@ Partition[tprvec,2,1];
     auc=Apply[Plus, wvec tmvec]; 
 Return[auc]]


  dif[{a_,b_}]:= b-a

 
(*-----------------------------------CONCAVE ROC FROM INITIAL ROC--------------------------------------------*)

ConcaveROC[list_]:= 
Module[{set,listnew,listnew1},
  (*select each point with highest slope starting from point (0,0) at position 1*)
  (*stops just before end of list*)
     {listnew,maxpox}=NestWhile[pick[list, #] &, {{{0, 0}}, 1}, (#[[2]] < Length[list]) &] ;
 (*avoids repeats e.g of (1,1)*)
     listnew1=UnionPairs[listnew];
Return[listnew1]]


 UnionPairs[list_]:=
  Module[{list1,list2,list3},
   (*sort by fpr*)  
      list1=Sort[list];
	(*Split groups by same lists {{.3,.7}},{{1,1},{1,1}}*)
      list2=Split[list1];
	(*Take first in each list*)
	   list3=#[[1]]& /@ list2;
  Return[list3]]


pick[list_,{listsofar_,possofar_}]:=
Module[{start,remlist,remslope,remfpr,remtpr,rempos,remset,remmax,
         range,maxslope,maxpos,maxfpr, maxtpr,listnew,res,startfpr,slopepairlist,
		starttpr,area}, 
 (*starting point on list*)
	    start=list[[possofar]];
	   {startfpr,starttpr}=start;
(*remainder of the list*)
         remlist=Drop[list,possofar];
         {remfpr,remtpr}=Transpose[remlist];
  (*compute pairs from starting value*)
     	 slopepairlist=Transpose[{(remtpr-starttpr),(remfpr-startfpr)}];  
  (*compute slope from starting value*)
  	     remslope=calcslope[slopepairlist];
  (*compute remainig points*)
         range=Range@Length[list];
        rempos=Drop[range,possofar];
         remset=Transpose[{remslope,rempos,remfpr,remtpr}];
        remmax=(Sort[remset])[[-1]];
      {maxslope,maxpos,maxfpr,maxtpr}=remmax;   
    listnew=Append[listsofar,{maxfpr,maxtpr}];
 Return[{listnew,maxpos}]]


calcslope[list_]:=calcslopei[#]& /@ list			
	
calcslopei[{num_,den_}]:= If[den==0,num+1000, num/den]; 




				  
				  
 End[] 
EndPackage[]


