(*:Mathematica:: Version 7        *)
(*:Context: "swirl`"              *)
(*:Swirl-and-Ripple               *)
(*: Version 1                     *)
(*Export:  swirlcore.m   PlotTwoGenes
          swirlsim.m    PlotTwoGenesForSim*)	   
(*:Title:  Plotting with two genes*)
(*Author: 2010 Stuart G. Baker    *)


BeginPackage["swirlplot`"]

PlotTwoGenes::usage="PlotTwoGenes[dataTRAIN,dataTEST,parsym,resTRAIN]"
PlotTwoGenesForSim::usage="PlotTwoGenesForSim[{x,y,k},{data,res,genenames}]"
 										   	
Clear[PlotTwoGenes,PlotTwoGenesForSim]


Begin["Private`"]


PlotTwoGenes[dataTRAIN_,dataTEST_,parsym_,resTRAIN_,genenames_:0]:=
Module[{plots1TRAIN,plots2TRAIN, plotr1TRAIN,plotr2TRAIN,
	     minmaxvec,
	     plots1TEST,plots2TEST,plotr1TEST,plotr2TEST},
          (*same range for all data*)
	        minmaxvec=MinMax[dataTRAIN,dataTEST,resTRAIN];	 
	 (*training sample plots*)
	       {plots1TRAIN,plots2TRAIN,plotr1TRAIN,plotr2TRAIN}=PlotTwoGenesCore[parsym,dataTRAIN,resTRAIN,minmaxvec,genenames];
	 (*test sample plots*)
                {plots1TEST,plots2TEST,plotr1TEST,plotr2TEST}=PlotTwoGenesCore[parsym,dataTEST,resTRAIN, minmaxvec,genenames];
	  (*group by distance measure*)
	        plotset1=GraphicsGrid[{{plots1TRAIN,plots1TEST},{plotr1TRAIN,plotr1TEST}}]; 		 
		plotset2=GraphicsGrid[{{plots2TRAIN,plots2TEST},{plotr2TRAIN,plotr2TEST}}];
	 (*print*);
	     Print[plotset1];
	     Print[plotset2];
 Return[Null]]




PlotTwoGenesForSim[parsym_,{data_,res_,genenames_}]:=
	Module[{minmax0,resx}, 
           minmax0=MinMax[data,data,res];
	    resx=PlotTwoGenesCore[parsym,data,res,minmax0,genenames];
	 Return[resx]]


 (*---------------Core of plot-------------------*)

PlotTwoGenesCore[parsym_,data_,res_,minmaxvec_,genenames_]:=
Module[{plots1,plots2,
            colorS,colorR, plotr1,plotr2,
	     textS,lineS,textR,lineR, resx},
(*swirl plots*)
     colorS=RGBColor[1,0,1];
     plots1=PlotSwirlCore[parsym,data,res,1,colorS,minmaxvec,genenames,1];
     plots2=PlotSwirlCore[parsym,data,res,2,colorS,minmaxvec,genenames,2];
(*ripple plots*)
     colorR=RGBColor[0,1,1];
     plotr1=PlotRippleCore[parsym,data,res,1,colorR,minmaxvec,genenames,1];
     plotr2=PlotRippleCore[parsym,data,res,2,colorR,minmaxvec,genenames,2];
  (*all plots*)
     resx={plots1,plots2,plotr1,plotr2};
 Return[resx]]




(*-----------------------Swirl------------------------------------------------------------------------------*)						  

PlotSwirlCore[{x_,y_,k_},{xmat0_,xmat1_},{posvecIN_,posvecOUT_,tlistx_,xm0_,xm1_,vm0_,vm1_},
    distance_,color_,minmaxvec_,genenames_,dname_]:=
Module[{x0,y0,x1,y1, labx,laby,
        pair0,pair1,plotpair0,plotpair1,
        xmin,xmax,ymin,ymax,xdif,ydif,
		  xcen0,ycen0,xcen1,ycen1,
		  v0x,v0y, v1x,v1y,
		  d0,d1,rx,
		  num1, level1, core1,
		  num2, level2, core2,
		  xcritx,
		  plotC0,plotC1,name,
		  xvec,yvec,d0vec,d1vec,rxvec,kset,cutvec,
		  d00,d01,d10,d11,r0,r1},
	  (*names*)
	   If[Length[genenames]==0, 
		  {labx,laby}={"gene A","gene B"},
		  {labx,laby}=genenames];
    (*selected genes from training sample*)
         {x0,y0}=xmat0[[posvecIN]];
	  {x1,y1}=xmat1[[posvecIN]];
 (*plot points*)
	  pair0=Transpose[{x0,y0}];
	  pair1=Transpose[{x1,y1}];
	  plotpair0=ListPlot[pair0, PlotStyle -> {RGBColor[1,0,0],PointSize[.03]}]; 
	  plotpair1=ListPlot[pair1, PlotStyle -> {RGBColor[0,1,0],PointSize[.03]}];   
(*range + extra sides*)
	  {xmin,xmax,ymin,ymax}=minmaxvec;
	  range={{xmin,xmax},{ymin,ymax}};
(*key parameters*)	
    {xcen0,ycen0}=Flatten[xm0];
    {xcen1,ycen1}=Flatten[xm1];
 (*variance*)
	{v0x,v0y,v1x,v1y}=GenVarList[vm0,vm1,distance];
 (*SYMBOLIC DISTANCE*)	
    d0 = Sqrt[(x - xcen0)^2 /v0x + (y -ycen0)^2/v0y];
    d1 = Sqrt[(x - xcen1)^2 /v1x + (y -ycen1)^2/v1y];
     rx = d0 /(d0 + d1);
 (*FORMULA: Solve for y at equal distance in terms of ratio rx*)
     yresk=y/.Solve[rx==k,y]//Simplify;
(*CRITICAL REGION FOR X*)
    (*solve for values of x such so y does not involve imaginary numbers--find part in square root in numerator*)
	 (*there are two solutions for yresk*)
	  num1=Numerator@yresk[[1]];
	  level1=Level[num1,1];
	  core1=level1[[-1]];
	  num2=Numerator@yresk[[2]];
	  level2=Level[num2,1];
	  core2=level2[[-1]];
	  core={core1,core2};
 (*critical valeus of x, so part in square root is 0*)
     xcritx= x/.Solve[core^2 == 0, x] // Simplify; 
 (*plots of centroid points*)
       plotC0 = ListPlot[{{xcen0, ycen0}}, Axes -> None, PlotRange->range,Frame -> True, 
                     PlotStyle -> {RGBColor[0,0,1],PointSize[.04]}]; 
     plotC1 = ListPlot[{{xcen1, ycen1}}, Axes -> None, PlotRange->range,Frame -> True, 
           PlotStyle -> {RGBColor[0,0,1],PointSize[.04]}]; 
  (*generate cutpoints for data*)	
	  xvec=Join[x0,x1];
           yvec=Join[y0,y1];					   
	  d0vec = Sqrt[(xvec - xcen0)^2 /v0x + (yvec -ycen0)^2/v0y];
  	  d1vec = Sqrt[(xvec - xcen1)^2 /v1x + (yvec -ycen1)^2/v1y];
           rxvec = d0vec /(d0vec + d1vec);
	  kset=Range[9]/10//N;
         cutvec=Quantile[Sort[rxvec],kset];
(*plot of equal distance boundary*)
     yM=y/.Solve[rx==.5,y]//Simplify;
      plotM=Plot[yM,{x,xmin,xmax},PlotStyle->ColorQ[.5],PlotRange->range];
 (*plot of other distance boundaries based on cutpoints*)
     plotset=Flatten@plotswirlcore[{x,y,k},rx,xcritx,range,plotC0,#]& /@ cutvec;
 (*combine plots*)
     plotset1=Join[{plotpair0,plotpair1,plotC0,plotC1,plotM},plotset];
	  name=StringJoin["Swirls  D=",ToString[dname]];
	  plot= Show[plotset1,PlotRange->range,Frame->True, Axes->False,  AspectRatio->1,
           FrameLabel->{labx,laby,name,None},DisplayFunction->$DisplayFunction];
 Return[plot]]



plotswirlcore[{x_,y_,k_},rx_,xcritx_,range_,plotC0_,k0_]:=
Module[{xcrik0,y0a,plot0a,y0b,plot0b},
(*critical value at cutpoint*)
    xcritk0=xcritx/.{k->k0};
(*rule out imaginary solutions*)
   If[Im[xcritk0[[1]]]==0 && Im[xcritk0[[2]]]==0,
   y0a=y/.Solve[rx==k0,y][[1]]//Simplify;
   plot0a=Plot[y0a,{x,xcritk0[[1]],xcritk0[[2]]},PlotStyle->ColorQ[k0],PlotRange->range];   
   y0b=y/.Solve[rx==k0,y][[2]]//Simplify;
  plot0b=Plot[y0b,{x,xcritk0[[1]],xcritk0[[2]]},PlotStyle->ColorQ[k0],PlotRange->range],
  plot0a=plotC0; plot0b=plotC0];
Return[{plot0a,plot0b}]]




(*---------------------------------Ripple---------------------------------------------------------------*)


PlotRippleCore[{x_,y_,k_},{xmat0_,xmat1_},{posvecIN_,posvecOUT_,tlistx_,xm0_,xm1_,vm0_,vm1_},
     distance_,color_,minmaxvec_,genenames_,dname_]:=
Module[{x0,y0,x1,y1, xdif,ydif,labx,laby,
        pair0,pair1,plotpair0,plotpair1,
        xmin,xmax,ymin,ymax,
		  xcen0,ycen0,xcen1,ycen1,
		  v0x,v0y, v1x,v1y,
		  d0sq,d1sq,rx,
		  plotC0,plotC1,
		  xvec,yvec,d0sqvec,d1sqvec,rxvec,kset,cutvec,
		  d00,d01,d10,d11,r0,r1},
	  (*names*)
	   If[Length[genenames]==0, 
		  {labx,laby}={"gene A","gene B"},
		  {labx,laby}=genenames];
    (*selected genes from training sample*)
     {x0,y0}=xmat0[[posvecIN]];
	  {x1,y1}=xmat1[[posvecIN]];
	(*plot points*)
	  pair0=Transpose[{x0,y0}];
	  pair1=Transpose[{x1,y1}];
	  plotpair0=ListPlot[pair0, PlotStyle -> {RGBColor[1,0,0],PointSize[.03]}]; 
	  plotpair1=ListPlot[pair1, PlotStyle -> {RGBColor[0,1,0],PointSize[.03]}];   
 	(*range + extra sides*)
	  {xmin,xmax,ymin,ymax}=minmaxvec;
	  range={{xmin,xmax},{ymin,ymax}};
	(*key parameters*)	
    {xcen0,ycen0}=Flatten[xm0];
    {xcen1,ycen1}=Flatten[xm1];
 (*variance*)
	{v0x,v0y,v1x,v1y}=GenVarList[vm0,vm1,distance];
 (*SYMBOLIC DISTANCE*)	
    d0sq = (x - xcen0)^2 /v0x + (y -ycen0)^2/v0y;
    d1sq = (x - xcen1)^2 /v1x + (y -ycen1)^2/v1y;
    rx = d0sq-d1sq;
 (*FORMULA: Solve for y at equal distance in terms of ratio rx*)
     yresk=y/.Solve[rx==k,y]//Simplify;
 (*plots of centroid points*)
    plotC0 = ListPlot[{{xcen0, ycen0}}, Axes -> None, PlotRange->range,Frame -> True, 
        PlotStyle -> {RGBColor[0,0,1],PointSize[.04]}]; 
	 plotC1 = ListPlot[{{xcen1, ycen1}}, Axes -> None, PlotRange->range,Frame -> True, 
        PlotStyle -> {RGBColor[0,0,1],PointSize[.04]}]; 
 (*generate cutpoints for data*)	
	 xvec=Join[x0,x1];
    yvec=Join[y0,y1];					   
	 d0sqvec = (xvec - xcen0)^2 /v0x + (yvec -ycen0)^2/v0y;
  	 d1sqvec = (xvec - xcen1)^2 /v1x + (yvec -ycen1)^2/v1y;
    rxvec = d0sqvec - d1sqvec;
	 kset=Range[9]/10//N;
    cutvec=Quantile[Sort[rxvec],kset];
(*plot of equal distance boundary*)
    yM=y/.Solve[rx==.5,y]//Simplify;
    plotM=Plot[yM,{x,xmin,xmax},PlotStyle->ColorQ[.5],PlotRange->range];
 (*plot of other distance boundaries based on cutpoints*)
    plotset=plotripplecore[{x,y,k},rx,range,#]& /@ cutvec;
 (*combine plots*)
    plotset1=Join[{plotpair0,plotpair1,plotC0,plotC1,plotM},plotset];
	 name=StringJoin["Ripples  D=",ToString[dname]];
    plot= Show[plotset1,PlotRange->range,Frame->True,	  Axes->False, AspectRatio->1,
           FrameLabel->{labx,laby,name,None},DisplayFunction->$DisplayFunction];
  Return[plot]]
    


plotripplecore[{x_,y_,k_},rx_,range_,k0_]:=
Module[{y0,plot0,xmin,xmax,ymin,ymax},				 
{{xmin,xmax},{ymin,ymax}}=range;
y0=y/.Solve[rx==k0, y];
plot0=Plot[y0,{x,xmin,xmax},PlotStyle->ColorQ[k0],PlotRange->range];
Return[plot0]]


(*------------------functions for both------------------------*)

 GenVarList[vm0_,vm1_,distance_]:=
 Module[{n0,n1,
         w0x,w0y,w1x,w1y,
		v0x,v0y,v1x,v1y},
    n0=Length[vm0];
 n1=Length[vm1];
    If[distance==1,   vm=(vm0 (n0-1) + vm1 (n1-1) )/(n0+n1-2);   {w0x,w0y}=vm;     {w1x,w1y}=vm];	   
    If[distance==2,                                                                   {w0x,w0y}=vm0;    {w1x,w1y}=vm1];	
	(*normalize*) 
	 v0x=w0x/w0x;
	 v1x=w1x/w0x;
	 v0y=w0y/w0x;
	 v1y=w1y/w0x;	   
  	Return[{v0x,v0y,v1x,v1y}]];


ColorQ[k_]:=
Module[{res},
  If[k<.5,res=RGBColor[1,0,0]];
   If[k==.5,res=RGBColor[0,0,0]]; 
   If[k>.5,res=RGBColor[0,1,0]];
Return[res]]


(*----------For uniform plot ranges----------------------*)

MinMax[dataTRAIN_,dataTEST_,res_]:=
Module[{xminTRAIN,xmaxTRAIN,yminTRAIN,ymaxTRAIN,
       xminTEST,xmaxTEST,yminTEST,ymaxTEST,
		 xmin,xmax,ymin,ymax,xdif,ydif,minmaxvec}, 
 {xminTRAIN,xmaxTRAIN,yminTRAIN,ymaxTRAIN}=MinMaxOne[dataTRAIN,res];
 {xminTEST,xmaxTEST,yminTEST,ymaxTEST}=MinMaxOne[dataTEST,res];
  xmin=Min[xminTRAIN,xminTEST];
  xmax=Max[xmaxTRAIN,xmaxTEST];
  ymin=Min[yminTRAIN,yminTEST];
  ymax=Max[ymaxTRAIN,ymaxTEST];
  xdif=(xmax-xmin)/10;
  ydif=(ymax-ymin)/10;
  minmaxvec={xmin-xdif,xmax+xdif,ymin-ydif,ymax+ydif};
 Return[minmaxvec]]


MinMaxOne[data_,{posvecIN_,posvecOUT_,tlistx_,xm0_,xm1_,vm0_,vm1_}]:=
 Module[{x0,y0,x1,y1,xmin0,xmax0,ymin0,ymax0,xmin1,xmax1,ymin1,ymax2,xmin,xmax,ymin,ymax,xdif,ydif,range},
  {xmat0,xmat1}=data;
 { x0,y0}=xmat0[[posvecIN]];
  {x1,y1}=xmat1[[posvecIN]];
  xmin=Min@Join[x0,x1];
  xmax=Max@Join[x0,x1];
  ymin=Min@Join[y0,y1];
  ymax=Max@Join[y0,y1];
 Return[{xmin,xmax,ymin,ymax}]]







End[] 

EndPackage[]




