(*:Mathematica:: Version 7 *)	   
(*:Context: "swirl`"       *)
(*:Swirl-and-Ripple        *)	
(*:Version  1              *)
(*Title: Illustrative plots of two genes and optimal distributionsy*)
(*Illustrative plots with two genes as well as formulas*)
(*llustrative plots of optimal distributions for Swirls and Ripples*)
(*Author: 2010 Stuart G. Baker *)

BeginPackage["swirlplotsym`","MultivariateStatistics`"]

PlotSym2D::usage="PlotSym2D[{x,y,k,a,b}] plots classification boundaries in various cases"
PlotSym2DOne::usage="PlotSym2DOne[{x,y,k,a,b},True] gives formulas for classification boundaries"
PlotOptDist::usage="PlotOptDist[q,{x,y},plot3DQ] gives examples of optimal distributions for Ripple and Swirl 
                    boundaries"

Clear[PlotSym2D,PlotSym2DOne,PlotOptDist] 

Begin["Private`"]

						
PlotSym2D[{x_,y_,k_,a_,b_},colortype_:"color"]:=
Module[{plot11,plot12, plot13, parset,plot21,plot22,plot23, plot},
 (*centroid for group 0 is 0,0*)
 (*centroid for group 1 is a,b*)
  parset={x,y,k,a,b};
  plot12=PlotSym2DOne[parset,{"Swirl",1},False,colortype];	 
  plot13=PlotSym2DOne[parset,{"Swirl",2},False,colortype];
  plot22=PlotSym2DOne[parset,{"Ripple",1},False,colortype]; 
  plot23=PlotSym2DOne[parset,{"Ripple",2},False,colortype]; 						
  plot=Show[GraphicsGrid[{{plot12,plot13},{plot22,plot23}}]];
Return[plot]]


PlotSym2DOne[{x_,y_,k_,a_,b_},{score_,distance_},show_,colortype_]:=
Module[{d0,d1,r,rx,yresk,core,xcrit,xcritx, range,
  w0x,w0y,w1x,w1y,
  v0x,v0y,v1x,v1y,
  dmax,h,ax,bx,
  plot0,plot0x,plot0y,plot1,plot1x,plot1y,
  xmin,xmax,
  yM,plotM, colorPoint},
(*variances specified w's are for plots*)
  If[distance==1,
    labvar="unit variance";
	 {w0x,w0y,w1x,w1y}={0,0,0,0}; 
    {v0x,v0y,v1x,v1y}={1,1,1,1}]; 
  If[distance==1,
    labvar="D=1";	 
	 {w0x,w0y,w1x,w1y}={1,4,1,4};
    {v0x,v0y,v1x,v1y}={1,4,1,4}];
 If[distance==2,
   labvar="D=2";	 
	{w0x,w0y,w1x,w1y}={1,1,1,4};
	{v0x,v0y,v1x,v1y}={1,1,1,4}];
(*scores formula*)
If[score=="Swirl",
    d0 = Sqrt[x^2 /v0x + y^2/v0y];
    d1 = Sqrt[(x - a)^2/v1x + (y - b)^2/v1y];
     r = d0 /(d0 + d1)];
If[score=="Ripple",
   d0 = x^2/v0x + y^2/v0y;
   d1 = (x - a)^2/v1x + (y - b)^2/v1y;
   r = d0 - d1];
If[score=="Swirl" && distance==1,  lab="Swirl,  D=1"];
If[score=="Swirl" && distance==2,  lab="Swirl,  D=2"];  
If[score=="Ripple" && distance==1,  lab="Ripple,  D=1"];
If[score=="Ripple" && distance==2,  lab="Ripple,  D=2"];  
(*----------Print Formulas--------------------------------------------------*)
	ShowFormula[{x,y,k},r,{score,distance},show]; 
 (*------------Set Up Values----------------------------------------------*)
 ax=4;
 bx=4;
 rule={a->ax,b->bx};
 rx=r/.rule;
range={Automatic,Automatic}; 
 yresk=y/.Solve[rx==k,y]//Simplify;  
If[score=="Swirl",
 (*solve for values of x such so y does not involve imaginary numbers*)
   core= yresk[[1,-1,2]];	 
(*critical valeus of x*)
  xcrit= x/.Solve[core^2 == 0, x] // Simplify; 
   xcritx=xcrit/.rule//Simplify];
(*------------plot confidence intervals around centroids---------------------------------*)
If[colortype=="color", colorPoint=RGBColor[0,0,1]];
If[colortype=="blackandwhite", colorPoint= RGBColor[0,0,0]];
 h=.4;
 plot0 = ListPlot[{{0, 0}, {ax, bx}}, Axes -> None, Frame -> True, 
    PlotRange -> {{-8, 8}, {-8, 8}}, PlotStyle -> {colorPoint,PointSize[.05]},
	 PlotRange->range]/.rule;					  
 plot0y = ListPlot[{{0, -1.96 w0y h }, {0,1.96 w0y h}}, Joined->True, Axes -> None, Frame -> True, 
    PlotRange -> {{-8, 8}, {-8, 8}}, PlotStyle -> {colorPoint},
	 PlotRange->range]/.rule;					  
plot0x = ListPlot[{{-1.96 w0x h, 0}, {1.96 w0x h,0}}, Joined->True, Axes -> None, Frame -> True, 
    PlotRange -> {{-8, 8}, {-8, 8}}, PlotStyle -> {colorPoint},
	 PlotRange->range]/.rule;					  
 plot1y = ListPlot[{{ax, bx-1.96 w1y h}, {ax,bx+1.96 w1y h}}, Joined->True, Axes -> None, Frame -> True, 
    PlotRange -> {{-8, 8}, {-8, 8}}, PlotStyle -> {colorPoint},
	 	 PlotRange->range]/.rule;				  
 plot1x = ListPlot[{{ax-1.96 w1x h, bx}, {ax+1.96 w1x h,bx}}, Joined->True, Axes -> None, Frame -> True, 
    PlotRange -> {{-8, 8}, {-8, 8}}, PlotStyle -> {colorPoint},
	 PlotRange->range]/.rule;
(*--------------Swirl Plot-----------------------------------------------------*)
If[score=="Swirl",
 	 yM=y/.Solve[rx==.5,y]//Simplify;
     plotM=Plot[yM,{x,-8,8},PlotStyle-> ColorQ[.5,.5,colortype],PlotRange->range];
 (*plots of centroid points*)
   cutvec={.2,.4,.6,.8};
	  cutvec={.1,.2,.3,.4,.6,.7,.8,.9};		   
   plotC0 = ListPlot[{{0, 0}}, Axes -> None, PlotRange->range,Frame -> True, 
       PlotStyle -> {colorPoint,PointSize[.03]}]; 
   plotsetx=Flatten@(plotswirlcore[{x,y,k},rx,xcritx,range,plotC0,#,colortype]& /@ cutvec);
   plotset=Join[{plot0,plot0x,plot0y,plot1x,plot1y},plotsetx,{plotM}]; 
   plot=Show[plotset,PlotRange->{{-8,8},{-8,8}},Frame->True,
	 FrameTicks->None, FrameLabel->{"Gene A","Gene B",lab,None},BaseStyle->Medium, AspectRatio->1]];	
(*--------------Ripple Plot-----------------------------------------------------*)
If[score=="Ripple", 
 dmax=60;
 cutvec={-.9,-.4,0,.3,.8} dmax;
 cutvec={-.9,-.6,-.4,0,.3,.6,.8} dmax;
   plotsetx=Flatten@(plotripplecore[{x,y,k},rx,range,#,colortype]& /@ cutvec);
   plotset=Join[{plot0,plot0x,plot0y,plot1x,plot1y},plotsetx]; 
   plot=Show[plotset,PlotRange->{{-8,8},{-8,8}},Frame->True,
	 FrameTicks->None, FrameLabel->{"Gene A","Gene B",lab,None},BaseStyle->Medium, AspectRatio->1]];
Return[plot]]



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



plotswirlcore[{x_,y_,k_},rx_,xcritx_,range_,plotC0_,k0_,colortype_]:=
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,.5,colortype],PlotRange->range];   
  y0b=y/.Solve[rx==k0,y][[2]]//Simplify;
  plot0b=Plot[y0b,{x,xcritk0[[1]],xcritk0[[2]]},PlotStyle->ColorQ[k0,.5,colortype],PlotRange->range],
  plot0a=plotC0; plot0b=plotC0];
 Return[{plot0a,plot0b}]]



ColorQ[k_,a_,colortype_]:=
Module[{res},
If[colortype=="blackandwhite",
 If[k<a,res={Black}];
 If[k==a,res={Black}];
 If[k>a,res={Black}]]; 
If[colortype=="color",
 If[k<a,res={RGBColor[1,0,0]}];
 If[k==a,res={RGBColor[0,0,0]}];
 If[k>a,res={RGBColor[0,1,0]}]];
Return[res]]



ShowFormula[{x_,y_,k_},r_,{score_,distance_},show_]:=
Module[{yreskQ,coreQ,xcritQ},
(*------------SYMBOLIC FORMULAS*)
If[score=="Swirl" && show==True,
 Print["Model ",{score,distance}];
 yreskQ=y/.Solve[r==k,y]//Simplify;
 coreQ= yreskQ[[1,-1,2]];	 
 xcritQ= x/.Solve[coreQ^2 == 0, x] // Simplify; 
 Print["Swirls Formula"];
 Print["one centroid is 0 0 and one is a, b"];
 Print["Solution to  r= d0/(d0+d1) ==k, y as a functinon of x"];
 Print[yreskQ];
 Print@InputForm[yreskQ];
 Print["yset=part of solution needed to be real"];
 Print[coreQ];
 Print["xset= range of x so that yset==0"];
 Print[xcritQ];
 Print["swirls plot y in region of xset"]]; 
If[score=="Ripple" && show==True,
  Print["Model ",{score,distance}];
  yreskQ=y/.Solve[r==k,y]//Simplify;
  Print["Line formula"];
  Print["one centroid is 0 0 and one is a, b"];
  Print["Solution to  r= d0^2 -d1^2 ==k, y as a functinon of x"];
  Print@InputForm[yresQ];	 
  Print[yreskQ]];
Return[Null]]



(*-------------------------Plots of Optimal Distributions----------------------------*)


PlotOptDist[{x_,y_},color_,plot3DQ_:False]:=PlotOptDist[{0,0},{3,9},{x,y},color]

PlotOptDist[{x0_,y0_},{x1_,y1_},{x_,y_},colortype_:"blackandwhite"]:=
Module[{k,plot1,plot2,plot3,plot4,plot1x,plot2x,plot3x,plot4x,contourlabels},
 plot3DQ=False;
 k=10;
 If[colortype=="blackandwhite",color="GrayTones"];
 If[colortype=="color",color="BlueGreenYellow"];
 contourlabels=Automatic;
 plot1=RippleNormal[{x0,y0},{x1,y1},{28,5},k,False,color,contourlabels];
 plot2=SwirlNormal[{x0,y0},{x1,y1},{28,5},k,False,color,contourlabels];
 plot3=SwirlNormal[{x0,y0},{x1,y1},{128,5},k,False,color,contourlabels];
 plot4=SwirlUniform[{x0,y0},{x1,y1},k,False,color,contourlabels];
 Print@GraphicsGrid[{plot1,plot2,plot3,plot4}];
 Return[Null]]

 RippleNormal[{x0_,y0_},{x1_,y1_},{var_,cov_},k_,plot3DQ_,color_,contourlabels_]:=
 Module[{g0,g1,plotg0,plotg1,plotf0,plotf1,plotres},
  g0=PDF[MultinormalDistribution[{x0,y0},{{var,cov},{cov,var}}],{x,y}];
  g1=PDF[MultinormalDistribution[{x1,y1},{{var,cov},{cov,var}}],{x,y}];
  plotg0=ContourPlot[g0,{x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color,
   ContourLabels -> contourlabels, PlotLabel->"Ripples class 0"];
  plotg1=ContourPlot[g1,{x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color,
     ContourLabels -> contourlabels, PlotLabel->"Ripples class 1"];
  plotres=GraphicsRow[{Null,plotg0,plotg1}];
  Print[plotres];
  If[plot3DQ==True,
   plotf0=Plot3D[g0,{x,x0-k,x1+k},{y,y0-k,y1+k}];
   plotf1=Plot3D[g1,{x,x0-k,x1+k},{y,y0-k,y1+k}];
   Print@GraphicsRow[{Null,plotf0,plotf1}]];	 
 Return[{Null,plotg0,plotg1}]]

SwirlNormal[{x0_,y0_},{x1_,y1_},{var_,cov_},k_,plot3DQ_,color_,contourlabels_]:=
Module[{d0,d1,p,xn,yn,f,z,w0,w1,range,plotN,plotN0,plotN1,plot3N,plot3N0,plot3N1,plotres},
    d0=Sqrt[(x-x0)^2 + (y-y0)^2];
    d1=Sqrt[(x-x1)^2 + (y-y1)^2];
    p=d0/(d0+d1);
  	 xn=(x0+x1)/2;
	 yn=(y0+y1)/2;
    f=PDF[MultinormalDistribution[{xn,yn},{{var,cov},{cov,var}}],{x,y}];
   (*NIntegrate sometimes reports a problem but gives an answer*)
   (*w0=NIntegrate[f (1-p), {x, -Infinity, Infinity},{y, -Infinity, Infinity}];
     w1=NIntegrate[f p, {x, -Infinity, Infinity},{y, -Infinity, Infinity}]; *)
   (*approximation to integral; result is similar to result from NIntegrate*)
   z=.3;
	w0=z^2 Plus@@Flatten@Table[f (1-p),{x,x0-2k,x1+2k,z},{y,y0-2k,y1+2k,z}];
   w1=z^2 Plus@@Flatten@Table[f p,{x,x0-2k,x1+2k,z},{y,y0-2k,y1+2k,z}];		
   range={{x0-k,x1+k},{y0-k,y1+k}};
   plotN=ContourPlot[f,               {x,x0-k,x1+k},{y,y0-k,y1+k},PlotRange->range,ColorFunction->color,
     ContourLabels -> contourlabels,PlotLabel->"Swirls background"];
   plotN0=ContourPlot[f (1- p) /w0,   {x,x0-k,x1+k},{y,y0-k,y1+k},PlotRange->range,ColorFunction->color, 
	  ContourLabels -> contourlabels,PlotLabel->"Swirls class 0"];
   plotN1=ContourPlot[f     p  /w1,   {x,x0-k,x1+k},{y,y0-k,y1+k},PlotRange->range,ColorFunction->color, 
	  ContourLabels -> contourlabels,PlotLabel->"Swirls class 1"];
   plotres=GraphicsRow[{plotN,plotN0,plotN1}];
	Print[plotres];
  If[plot3DQ==True,
    plot3N=Plot3D[f,             {x,x0-k,x1+k},{y,y0-k,y1+k},  ColorFunction->"GrayTones", PlotLabel->"background"];
    plot3N0=Plot3D[f (1-p)/ w0,  {x,x0-k,x1+k},{y,y0-k,y1+k},  ColorFunction->"GrayTones", PlotLabel->"class 0"];
    plot3N1=Plot3D[f     p/ w1,  {x,x0-k,x1+k},{y,y0-k,y1+k},  ColorFunction->"GrayTones", PlotLabel->"class 1"];
    Print@GraphicsRow[{plot3N,plot3N0,plot3N1}]];
  Return[{plotN,plotN0,plotN1}]]

SwirlUniform[{x0_,y0_},{x1_,y1_},k_,plot3DQ_,color_,contourlabels_]:=
Module[{d0,d1,p,g,plotU,s0,s1,plotU0,plotU1,plot3U,plot3U0,plot3U1,plotres},
    d0=Sqrt[(x-x0)^2 + (y-y0)^2];
    d1=Sqrt[(x-x1)^2 + (y-y1)^2];
    p=d0/(d0+d1);
 (*Uniform underlying distribution*)
   gsum=NIntegrate[1 ,{x,x0-k,x1+k},{y,y0-k,y1+k}];
	g=1/gsum//N;
   plotU=ContourPlot[g,{x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color,
   ContourLabels -> contourlabels, PlotLabel->"Swirls background"];
   s0=NIntegrate[g (1-p),{x,x0-k,x1+k},{y,y0-k,y1+k}];
   s1=NIntegrate[g  p,   {x,x0-k,x1+k},{y,y0-k,y1+k}];		  
   plotU0=ContourPlot[(1-p) g/ s0,  {x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color, 
	   ContourLabels -> contourlabels, PlotLabel->"Swirls class 0"];
   plotU1=ContourPlot[ p    g/ s1,  {x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color, 
	   ContourLabels->contourlabels,PlotLabel->"Swirls class 1"];
   plotres=GraphicsRow[{plotU,plotU0,plotU1}];
	Print[plotres];
   If[plot3DQ==True,
     plot3U=Plot3D[g,        {x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color, PlotLabel->"Swirls background"];
     plot3U0=Plot3D[g (1-p), {x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color, PlotLabel->"Swirls class 0"];
     plot3U1=Plot3D[g  p,    {x,x0-k,x1+k},{y,y0-k,y1+k},ColorFunction->color, PlotLabel->"Swirls class 1"];
    Print@GraphicsRow[{plot3U,plot3U0,plot3U1}]];	
  Return[{plotU,plotU0,plotU1}]]



End[] 

EndPackage[]


parset={x0,y0,k0,a0,b0};