(*Mathematica:: Version 8.0      *)
(*:Context:: "mfit`"                     *)
(*:Title::  mfitpair            		*)
(*:Summary:: fitting algorithm for gene sequences *)
(*:References::  Finding biologically relevant response patterns in high-throughput experiments*)
(*Date 2014   Stuart G. Baker   *)
(*Input
ComputeAIC			mfitsupport.m*)


BeginPackage["mfitdiagnostic`","mfitsupport`"]  

PlotRPEAIC::usage="PlotRPEAIC[modelset,zbestmat1,zbestmat2,n]"

Clear[PlotRPEAIC]



Begin["Private`"]	  


PlotRPEAIC[modelset_,zbestmat1_,zbestmat2_,n_]:=
  Module[{resx,ztype1,ztype2,zstatDEL1,zstatDEL2,zaicDEL1,zaicDEL2,pos,
   times,len,range,mat,matx,matA,vecA,matB,vecB, xname,yname,plotcore, plot,plotnamex,
  xlist,ylist,xmin,xmax,ymin,ymax,x0,x1,hline,vline,yA,yB,textA,textB,pointA,pointB,tri,triA,triB},
 Print[" "];
 Print["Plot change in RPE versus change in AIC"];
 (*number of points for AIC computation*)
     (*number of genes*) 
        len=Length@zbestmat1;
        range=Range[len];
     (*results for all data*)
          mat=GenMatRPEAIC[modelset,#,zbestmat1,zbestmat2,n]& /@ range;
     (*results for good fits*)      
          matx=Select[mat,(#[[1]]==2)&];
     (*re-arrange by each gene*) 
       {resx,ztype1,ztype2,zstatDEL1,zstatDEL2,zaicDEL1,zaicDEL2,pos}=Transpose[matx];
         posx=Join[pos,pos];     
          ztype=Join[ztype1,ztype2];
          zstatDEL=Join[zstatDEL1,zstatDEL2] 100;
          zaicDEL=Join[zaicDEL1,zaicDEL2];
       (*check by printing*)
          match=Transpose[{posx,ztype,zstatDEL,zaicDEL}];
	    colname={"pos","ztype","RPEDIF","AICDIF"};
           (* Print@TableForm[match,TableHeadings->{None,colname}]; *)   
      (*separate by type A and type B fits*)  
         tri=Transpose[{zaicDEL,zstatDEL,ztype}];
	    triA=Select[tri,(#[[-1]]=="A")&];
	    triB=Select[tri,(#[[-1]]=="B")&];
	   vecA=Take[#,2]& /@ triA; 
	   vecB=Take[#,2]& /@ triB; 
   lenA=Length[vecA];
   lenB=Length[vecB];
   If[lenA >0 && lenB >0,
  (*plot core*) 
        plotcore=ListPlot[{vecA,vecB},PlotStyle->{Red,Darker[Green]}];
    (*point labels*)
      {xlist,ylist}=Transpose@Join[vecA,vecB];
       xmax=Max[xlist];
     xmin=Min[xlist];
      ymax=Max[ylist];
      ymin=Min[ylist];
      x0=xmin+.04 (xmax-xmin);
      x1=xmin+.02 (xmax-xmin);
      yA=ymin+.95 (ymax-ymin);
      yB=ymin+.85 (ymax-ymin);
      textA=Graphics@{Text[Style["Curve A selected",FontSize->10],{x0,yA},{-1,0}]};
      textB=Graphics@{Text[Style["Curve B selected",FontSize->10],{x0,yB},{-1,0}]};    
      pointA=Graphics@{PointSize[.01],Red,Point[{x1,yA}]};
      pointB=Graphics@{PointSize[.01],Darker[Green],Point[{x1,yB}]};    
  (*frame labels*)   
     yname=DisplayForm[Subscript["RPE","A"]-Subscript["RPE","B"]];   
     xname=DisplayForm[Subscript["AIC","A"]-Subscript["AIC","B"]];  
   (*lines*)
     hline=Graphics@Line[{{xmin,5},{xmax,5}}];
     vline=Graphics@Line[{{0,ymin},{0,ymax}}];
   (*plot*)
   ticksy={{0,"0%"},{5,"5%"},{10,"10%"},{15,"15%"},{20,"20%"},{25,"25%"},{30,"30%"}};
    plot=Show[plotcore,textA,textB,pointA, pointB, hline,vline,Frame->True,
           FrameTicks->{Automatic,ticksy,None,None},
            FrameLabel->{xname,yname, None,None},AspectRatio->1,LabelStyle->Medium];
    Print[plot];
    plotnamex="mfitFIG-RPE-AIC.eps";
    Export[plotnamex,plot];
   Print["exporting plot named ",plotnamex],
   (*else*)
   Print["too few points for AIC RPE plot"]];
    Return[Null]]


GenPointsRPEAIC[mat_]:=
 Module[{resx,ztype1,ztype2,zstatDEL1,zstatDEL2,zaicDEL1,zaicDEL2,zstatDEL,zaicDEL,pair,pos0},
    {resx,ztype1,ztype2,zstatDEL1,zstatDEL2,zaicDEL1,zaicDEL2,pos0}=Transpose[mat];
    zstatDEL=Join[zstatDEL1,zstatDEL2];
    zaicDEL=Join[zaicDEL1,zaicDEL2];
    pair=Transpose[{zaicDEL,zstatDEL}]; 
 Return[pair]]  


GenMatRPEAIC[modelset_,pos_,zbestmat1_,zbestmat2_,n_]:=
Module[{zbestvec1,zbestvec2, 
ztype1,zmod1,zstatA1,zstatB1,zparA1,zparB1,zmseA1,zmseB1,
ztype2,zmod2,zstatA2,zstatB2,zparA2,zparB2,zmseA2,zmseB2,
zstatDEL1,zstatDEL2,zaicDEL1,zaicDEL2, res,resx},
     (*vector for position fo gene*)
           zbestvec1=zbestmat1[[pos]];
           zbestvec2=zbestmat2[[pos]];
       {ztype1,zmod1,zstatA1,zstatB1,zparA1,zparB1,zmseA1,zmseB1}=SelectParameters[zbestvec1];
       {ztype2,zmod2,zstatA2,zstatB2,zparA2,zparB2,zmseA2,zmseB2}=SelectParameters[zbestvec2];
        testmodQ=MemberQ[modelset,zmod1] && MemberQ[modelset,zmod2];
   (*AIC*)    
      zaicA1=ComputeAIC[zparA1,zmseA1,n];
      zaicA2=ComputeAIC[zparA2,zmseA2,n];
      zaicB1=ComputeAIC[zparB1,zmseB1,n];
      zaicB2=ComputeAIC[zparB2,zmseB2,n];
  (*AIC DIF*)    
      zaicDEL1=zaicB1-zaicA1;
      zaicDEL2=zaicB2-zaicA2;
  (*RPE DIF*)    
      zstatDEL1=zstatB1-zstatA1;
      zstatDEL2=zstatB2-zstatA2;
  (*If either best or next is selected in both*)
      If[(ztype1=="A" || ztype1=="B") && (ztype2=="B" || ztype2=="A") && testmodQ, resx=2]; 
       If[(ztype1=="A" || ztype1=="B") && (ztype2=="B" || ztype2=="A") && (testmodQ==False), 
            Print["MODEL EXCLUDED"];    resx=2]; 
      If[(ztype1=="A" || ztype1=="B") && (ztype2=="C" || ztype2=="D"),  resx=1]; 
      If[(ztype1=="C" || ztype1=="D") && (ztype2=="B" || ztype2=="A"), resx=1];   
      If[(ztype1=="C" || ztype1=="D") && (ztype2=="C" || ztype2=="D"), resx=0];   
  (*output*)   
     res={resx,ztype1,ztype2,zstatDEL1,zstatDEL2,zaicDEL1,zaicDEL2,pos};
Return[res]]


  
  


SelectParameters[zbestvec_]:=
Module[{ztype,zmod,zstat,zfunc,zpar,zmse,
              zmodA,zstatA,zfuncA,zparA,zmseA,
              zmodB,zstatB,zfuncB,zparB,zmseB,
               zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,pos},
{ztype,zmod,zstat,zfunc,zpar,zmse,
              zmodA,zstatA,zfuncA,zparA,zmseA,
              zmodB,zstatB,zfuncB,zparB,zmseB,
               zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,pos}=zbestvec;
Return[{ztype,zmod,zstatA,zstatB,zparA,zparB,zmseA,zmseB}]]




End[] 
EndPackage[]