
(*:Mathematica:: Version 7  *)
(*:Context: "sr`"           *)
(*:Swirl-and-Ripple         *)
(*: Version 1               *)
(*Input	 swirlplot.m    PlotTwoGenesForSim
                  swirlroc.m   ComputeROCTestSample*)
(*:Author: 2010 Stuart G. Baker *)



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


PlotSim2D::usage="PlotSim[{x,y,k},n,type]"
GenSimData::usage="GenSimData[name,scenario]"


Clear[PlotSim2D,GenSimData]


Begin["Private`"]

  
  GenSimData[name_,scenario_,n_]:=
  Module[{data0a,data0b,data0c,data0d,data0,
			 data1a,data1b,data1c,data1d,data1, cname,namex,
			 data,genename,gename0,genename1,datalist},
	p=2000; (*number of non-informative genes*)
	Print["simulated data"];
	{data0a,data1a}=GenSimDataNonInf[p,n];
	{data0b,data1b}=GenSimDataTwoNormals[n];
 	If[scenario==1,
           genename1={"A","B"};
	   data0=Join[data0a,data0b];
	   data1=Join[data1a,data1b]];
	If[scenario==2,
           genename1={"A","B","C","D"}; 
           {data0c,data1c}=GenSimDataOneBin[n,{0,5},{3,1},{.3,.7}];
	   data0=Join[data0a,data0b,data0c];
	   data1=Join[data1a,data1b,data1c]];
	(*data*)   
        data={data0,data1};
	genename0=string[#]& /@ Range[p];
	genename=Join[genename0,genename1];
        namex=StringJoin[name," n=",ToString[n]];
	cname={namex,"class 0","class 1"};
	datalist={data,genename,cname};
	Return[datalist]]


string[x_]:=StringJoin["NonInfGene",ToString[x]]

(*------------------------------specific distributions----------------------*)

 GenSimDataNonInf[p_,n_]:= GenSimDataNonInf[p,n,{0,5}]

 GenSimDataNonInf[p_,n_,{mean_,se_}]:=
  Module[{data0,data1,i},
        Print["  ",p,"  noninformative genes with mean se for each class = ",{mean,se}];
	data0 = Table[RanNormList[mean,se,n], {i, 1, p}];
	data1 = Table[RanNormList[mean,se,n], {i, 1, p}];
 Return[{data0,data1}]]


GenSimDataTwoNormals[n_]:=GenSimDataTwoNormals[n,{0,5},{2,1}]

 GenSimDataTwoNormals[n_,{mean0_,se0_},{mean1_,se1_}]:=
  Module[{x0A,x0B,x1A,x1B,data0,data1,data,i},
       Print["  1 gene A  mean and se for class 0 ",{mean0, se0},"     class 1 ",{mean1, se1}];
       Print["  1 gene B  mean and se for class 0 ",{mean0, se0},"     class 1 ",{mean1, se1}];
        x0A = Table[RanNorm[mean0,se0], {i, 1, n}];
        x0B = Table[RanNorm[mean0,se0], {i, 1, n}];
         x1A = Table[RanNorm[mean1,se1], {i, 1, n}];
        x1B = Table[RanNorm[mean1,se1], {i, 1, n}];
        data0={x0A,x0B};
        data1={x1A,x1B};
	data={data0,data1};		   
 Return[data]]
    
    
   GenSimDataOneBin[n_]:= GenSimDataOneBin[n,{0,5},{3,1},{.3,.7}];

GenSimDataOneBin[n_,{mean0_,se0_},{mean1_,se1_},{p0_,p1_}]:=
  Module[{x0A,x0B,x1A,x1B,data0,data1,data,i},
 	Print["  1 gene A  mean and se for class 0 ",{mean0, se0},"   class 1 ",{mean1, se1}];
       Print["  1 gene B  binary with probability of 1 for class 0 ",p0,"  class 1 "p1];
       x0A = Table[RanNorm[mean0,se0], {i, 1, n}];
        x0B = RanBinX[p0,n];
	x1A = Table[RanNorm[mean1,se1], {i, 1, n}];
        x1B = RanBinX[p1,n];
        data0={x0A,x0B};
        data1={x1A,x1B};
	data={data0,data1};		   
 Return[data]]



RanNormList[mean_,se_,n_]:=
Module[{i},
  res=Table[RanNorm[mean,se], {i, 1, n}];
 Return[res]]




(*-------------------Plot  Simulated Data----------------------------------*)

 
PlotSim2D[{x_,y_,k_},numspecimens_,type_]:=
Module[{data,res,genenames,fset,
     lineR1,lineR2,lineS1,lineS2,
     textR1,textR2,textS1,textS2,
     rocR1,rocS1,rocR2,rocS2,
     plotR1,plotR2,plotS1,plotS2 ,
     plot2DR1,plot2DR2,plot2DS1,plot2DS2 ,
     plotROC1,plotROC2,plotdata,plot},
   {data,res,genenames}=GenSimData2[numspecimens,type];
   fset=N[Range[9]/10];
   rocR1=ComputeROCTestSample[data,res,{"Ripple",1},fset][[-1]];
   rocR2=ComputeROCTestSample[data,res,{"Ripple",2},fset][[-1]];
    rocS1=ComputeROCTestSample[data,res,{"Swirl",1},fset][[-1]];
   rocS2=ComputeROCTestSample[data,res,{"Swirl",2},fset][[-1]];
   {plotR1,lineR1,textR1}=PlotROCX[rocR1,fset,"Ripple",1];
   {plotR2,lineR2,textR2}=PlotROCX[rocR2,fset,"Ripple",2]; 
   {plotS1,lineS1,textS1}=PlotROCX[rocS1,fset,"Swirl",1];
    {plotS2,lineS2,textS2}=PlotROCX[rocS2,fset,"Swirl",2];
    plotROC1=Show[plotR1,plotS1,lineR1,textR1,lineS1,textS1];
    plotROC2=Show[plotR2,plotS2,lineR2,textR2,lineS2,textS2];
    {plot2DS1,plot2DS2,plot2DR1,plot2DR2}=PlotTwoGenesForSim[{x,y,k},{data,res,genenames}];
    plotset={{plot2DR1,plot2DS1,plotROC1},
                 {plot2DR2,plot2DS2,plotROC2}};
    plot=GraphicsGrid[plotset];
  Return[plot]]



GenSimData2[n_,type_:1]:=
Module[{data,genenames,datalist,data0,data1,xm0,xm1,vm0,vm1,posIN,posOUT,tlist,res,out},
        SeedRandom[18];
        If[type==1,  {data0,data1}=GenSimDataTwoNormals[n]];
        If[type==2,  {data0,data1}=GenSimDataOneBin[n]];  
         genenames={"Gene A","Gene B"};
         xm0=Mean[#]& /@ data0;
	 xm1=Mean[#]& /@ data1;
	 vm0=Var[#]& /@ data0;
	 vm1=Var[#]& /@ data1;
	 posIN={1,2};
	 posOUT={0,0};
	 tlist={1,1};
	 res={posIN,posOUT,tlist,xm0,xm1,vm0,vm1};
	 data={data0,data1};
	 out={data,res,genenames};
  Return[out]]




(*---------Support functions------------------------------*) 

 
 PlotROCX[tset_,fset_,score_,distance_]:=
 Module[{ticknamex,tickname,plotroc,roc,h,line,text,name},
   If[score=="Ripple",color=Hue[.1],h=.4];
   If[score=="Swirl", color=RGBColor[0,0,1],h=.2];	 
   If[distance==1, name="D = 1"];
   If[distance==2, name="D = 2"]; 
    roc={{0,0}} ~Join~ Transpose[{fset,tset}];
    ticknamex={.2,.4,.6,.8,1};
    tickname={ticknamex,ticknamex};
    plotroc=ListPlot[roc,Joined->True,
   namex=StringJoin["ROC curve   ",name];
    PlotStyle->{color},PlotLabel->namex,AxesLabel->{"FPR","TPR"},AxesOrigin->{0,0},	Ticks->tickname];
     line=ListPlot[{{.8,h},{.9,h}},Joined->True,PlotStyle->color];
     text =Graphics@Text[score,{.8,h},{1,0}];    
  Return[{plotroc,line,text}]]




 
RanBinX[p_,n_]:= 
Module[{ranuni,bin,i},
    ranuni=Table[Random[],{i,1,n}];
    bin=Sign[1+Sign[ranuni-p]];
 Return[bin]]


RanNorm[mu_,sigma_]:=
   N[mu + sigma Sqrt[-2 Log[Random[]]] Cos[2Pi Random[]]]


 RanUnitNorm[n_]:=Table[RanNorm[0,1],{n}]
 

 
 End[] 
EndPackage[]


