(*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 2013   Stuart G. Baker   *)


(*Inputs
PlotRPEAIC		    mfitdiagnostic.m
GoodModelPairQ        mfitsupport.m
PosGoodPair:             mfitsupport.m
ShortModelName	     mfitsupport.m
FormPair			     mfitform.m
FormPairVec		     mfitform.m
ReportShapePairs	      mfitreport.m
PlotShapePairs           mfitplot.m*)
  
 

BeginPackage["mfitpair`","mfitdiagnostic`","mfitplot`","mfitreport`", "mfitcore`","mfitform`","mfitsupport`"]  

MFitPairReport::usage="MFitPairReport[{zbestmat1,zbestmat2,dataset1,dataset2,par,plitQ,exportplotQ}]"
           
Clear[MFitPairReport]


Begin["Private`"]	  


MFitPairReport[modelset_,{zbestmat1_,zbestmat2_,dataset1_,dataset2_,par_,splitQ_,
                         exportplotQ_,modelsetQ_,kform_},showsummary_,showreport_,showplot_,
                         showRPEAIC_,format_]:=
Module[{input,modpairvecS,plotname1,plotname2,times,n0,n},
      plotname1=dataset1[[-2]];
      plotname2=dataset2[[-2]];
       tallyvec=TallyGoodModelPairs[modelset,zbestmat1,zbestmat2,plotname1,plotname2];
      modpairvec=(Take[#,1][[1]])& /@ tallyvec;
      modpairvecS=ReOrder[modelset,modpairvec];
      input={zbestmat1,zbestmat2,dataset1,dataset2,par,splitQ,exportplotQ,kform};

    (*fitting diagnostics*)
    If[showRPEAIC,
        times=dataset1[[4]];
        n0=Length[times];
        If[OddQ[n0], n=(n0-1)/2, n=n0/2];
        PlotRPEAIC[modelset,zbestmat1,zbestmat2,n]];            
   
     If[showsummary,
       Print["OVERALL COUNTS"]
        If[modelsetQ=="full+poly",         ReportBRPolySummary[zbestmat1,zbestmat2,plotname1,plotname2]];
        If[modelsetQ=="reduced" || modelsetQ=="full",
            ReportBRSummary[zbestmat1,zbestmat2,plotname1,plotname2]];
        
        Print["TRANSITION TABLES"];  
            TransitionTable[modelset,tallyvec,plotname1,plotname2,modelsetQ];     
        Print["MODEL PAIR COUNTS"];
         ReportPairCounts[modelset,input,modpairvecS]];
   
   (*print out tables of pairs*)       
     If[showreport,ReportShapePairs[modelset,input,modpairvecS]];
   
   (*plot pairs*)
   If[showplot, PlotShapePairs[modelset,input,modpairvecS,format]];
 
Return[Null]]



  (*------------------------------REPORT SUMMARY GOOD FITS-----------------*)

ReportBRPolySummary[zbestmat1_,zbestmat2_,plotname1_,plotname2_]:=
Module[{ztypevec1,ztypevec2,zmodvec1,zmodvec2,
               ztypeBR1,zmodBR1,ztypeBR2,zmodBR2,
               zadd,zABset,zACset,zBCset,zCCset,numgenes,
               resBR,resFLA,resLIN,resSIG,resDBS,resGEN,resPOLY,
               zBRset,zFLAset,zLINset,zSIGset,zDBSset,zGENset,sPOLYset,
              mat, pmat, rowname, colname},
     ztypevec1=(Transpose[zbestmat1])[[1]];
     ztypevec2=(Transpose[zbestmat2])[[1]];
     zmodvec1=(Transpose[zbestmat1])[[2]];
    zmodvec2=(Transpose[zbestmat2])[[2]];
    (*sets*)   
        zABset={"A","B"};
        zACset={"A","C"};
        zBCset={"B","C"};
        zCCset={"C"};
        resAB=GenCounts[zABset,ztypevec1,ztypevec2];
        resAC=GenCounts[zACset,ztypevec1,ztypevec2];
        resBC=GenCounts[zBCset,ztypevec1,ztypevec2];
        resCC=GenCounts[zCcset,ztypevec1,ztypevec2];
   (*create table*)  
     mat={resAB, resAC,resBC,resCC};
    rowname={"AB","AC","BC","CC"};
    colname={plotname1,plotname2,"both"};
    Print["SUMMARY COUNTS: A B = biologically relevant curves C= polynomial" ];
    Print["    ",TableForm[mat,TableHeadings->{rowname,colname}]];
   numgenes=Length[zbestmat1];
   pmat=GenProb[#,numgenes]& /@ mat;
    Print["SUMMARY FRACTIONS OF ALL GENES A B = biologically relevant curves C= polynomial" ];
    Print["    ",TableForm[pmat,TableHeadings->{rowname,colname}]];
   Return[Null]]

GenProb[vec_,numgenes_]:=Round[(vec/ numgenes) 100]/100//N;

ReportBRSummary[zbestmat1_,zbestmat2_,plotname1_,plotname2_]:=
Module[{ztypevec1,ztypevec2,zmodvec1,zmodvec2,numgenes,
               ztypeBR1,zmodBR1,ztypeBR2,zmodBR2,
               zadd,
               resBR,resFLA,resLIN,resSIG,resDBS,resGEN,resPOLY,
               zBRset,zFLAset,zLINset,zSIGset,zDBSset,zGENset,sPOLYset,
              mat, rowname, colname},
     ztypevec1=(Transpose[zbestmat1])[[1]];
     ztypevec2=(Transpose[zbestmat2])[[1]];
     zmodvec1=(Transpose[zbestmat1])[[2]];
    zmodvec2=(Transpose[zbestmat2])[[2]];
    (*sets*)   
        zBRset={"A","B"};
        resBR=GenCounts[zBRset,ztypevec1,ztypevec2];
   (*create table*)  
     mat={resBR};
    rowname={"biologically relevant"};
    colname={plotname1,plotname2,"both"};
    Print["SUMMARY COUNTS: Biologically relevant"];
    Print["    ",TableForm[mat,TableHeadings->{rowname,colname}]];
       numgenes=Length[zbestmat1];
    Print["SUMMARY FRACTIONS OF ALL GENES: Biologically relevant"];
     pmat=GenProb[#,numgenes]& /@ mat;
    Print["    ",TableForm[pmat,TableHeadings->{rowname,colname}]];
  Return[Null]]



 GenCounts[ztypeset_,ztypevec1_,ztypevec2_]:=
 Module[{type1,type2,type12,len1,len2,len12,res},
     type1=Ind[ztypeset,#]& /@ ztypevec1;
    type2=Ind[ztypeset,#]& /@ ztypevec2;
     type12=IndBothQ[ztypeset,#]& /@ Transpose[{ztypevec1, ztypevec2}];
    len1=Plus@@type1;
    len2=Plus@@type2;
    len12=Plus@@type12;
   res={len1,len2,len12};
   Return[res]]

Ind[ztypeset_,ztype_]:=
Module[{res},
If[MemberQ[ztypeset,ztype], res=1,res=0];
Return[res]]

IndBothQ[ztypeset_,{ztype1_,ztype2_}]:=
Module[{res},
If[(MemberQ[{"A","B"}, ztype1] && MemberQ[ztypeset,ztype1] &&  MemberQ[{"A","B"}, ztype2]) || 
  (MemberQ[{"A","B"}, ztype2] && MemberQ[ztypeset,ztype2] &&  MemberQ[{"A","B"}, ztype1]), res=1,res=0];
Return[res]]



(*------------REPORT SUMMARY PAIR COUNTS-------------------*)	

TallyGoodModelPairs[modelset_,zbestmat1_,zbestmat2_,plotname1_,plotname2_]:=
Module[{range,trimat,goodmat,onevec,len,modvec1,modvec2,modpairset,tallyvec,modpairvec},
    (*identify good fits*)
     range=Range@Length[zbestmat1];
     trimat=GoodModelPairQ[modelset,#,zbestmat1,zbestmat2]& /@ range;
     goodmat=Select[trimat,(#[[1]]==1)&];
    len=Length[goodmat];
    If[len>0,
      {onevec,modvec1,modvec2}= Transpose[goodmat];
       modpairset=Transpose[{modvec1,modvec2}];
      tallyvec=Tally[modpairset]];
    If[len==0, 
       Print["no good fits"]; tallyvec={}];
 Return[tallyvec]]
 

  (*------------------------------REPORT  SUMMARY PAIR COUNTS--------------*)

ReportPairCounts[modelset_,{zbestmat1_,zbestmat2_,dataset1_,dataset2_,
  par_,splitQ_,exportplotQ_,kform_}, modpairvecS_]:=
  Module[{hmat,hmat1,htab,hmat2},
   Print["   heterochrony (HC) and heterometry (HM)"];
     hmat=SummaryPairCounts[modelset,zbestmat1,zbestmat2,dataset1,dataset2,
                 par,splitQ,kform,#]& /@ modpairvecS;
    hmat1=Select[hmat,(#[[3]]=!=0)&];
    If[Length[hmat1]>0,
    hmat2=Sort[hmat1, (#1[[3]] > #2[[3]])&];
    htab=TableForm[hmat2,TableHeadings->
    {None, {"model1","model2","total","similar","heterochrony","heterometry","both"}}];
    Print["          ",htab],
    Print["           None"]];
Return[Null]]


SummaryPairCounts[modelset_,zbestmat1_,zbestmat2_,dataset1_,dataset2_,par_,splitQ_,kform_,modpair_]:=
Module[{mod1,mod2,len,range,posvec,posvecFIT,lenx,set,modvec,htypevec,hCvec,hMvec,hvec},
   (*find good fits for the model pair*)
     {mod1,mod2}=modpair;
     len=Length@zbestmat1;
    range=Range[len];
    posvec=PosGoodPair[modelset,#,zbestmat1,zbestmat2,modpair]& /@ range;
    posvecFIT=Select[posvec,(#>0)&];
  If[mod1==mod2,
    set=FormPair[posvecFIT,zbestmat1,zbestmat2,dataset1,dataset2,par,splitQ,kform];
    {modvec,htypevec, hCvec,hMvec,posvec}= Transpose[set];
    hvec=CompHCHM[htypevec,modpair],
  (*else*)
    lenx=Length[posvecFIT];
   hvec={mod1,mod2, lenx,0,0,0,0}];	
  Return[hvec]]


CompHCHM[htypevec_,modpair_]:=
Module[{len0,len1,len2,len3,zmod,res,mod1,mod2},
len0=Length@Select[htypevec,(#==0)&];
len1=Length@Select[htypevec,(#==1)&];
len2=Length@Select[htypevec,(#==2)&];
len3=Length@Select[htypevec,(#==3)&];
(*pick one of same models in pair*)
lentot=len0+len1+len2+len3;
{mod1,mod2}=modpair;
res={mod1,mod2,lentot, len0,len1,len2,len3};
Return[res]];



(*----------------------ordering of model pairs------------------*)


ReOrder[modelset_,modpairvec_]:=
Module[{len,range,rule,res,modpairvec1,modpairvec2,modpairvec3},
   len=Length[modelset];
   range=Range[len];
   rule=MapThread[Rule,{modelset,range}];
   modpairvec1=ExtendOne[rule,#]& /@ modpairvec;
  modpairvec2=Sort[modpairvec1];
  modpairvec3 =Flatten[Drop[#,1]]& /@ modpairvec2;
Return[modpairvec3]]



ExtendOne[rule_,modpair_]:=
Module[{c,res,numpair},
  numpair={100, 1} (modpair/.rule);
   c=Apply[Plus,numpair];
  res={c,modpair};
Return[res]]


PrintSame[mat_]:=
Module[{mats,matx,colname,headings,posvecs},
  (*largest to smallest*)
           mats=Reverse@Sort[mat];  
        (*drop sorting variable and positiion variable*)
          matx=Transpose@Drop[Drop[Transpose[mats],1],-1];
       (*print table groups*)
          colname={"genename","genid","class","HC%","HM%"};
          headings={None,colname};
         PrintTableFormGroupCore[matx,10,headings];
        (*new sorting variable*)
          posvecs=Flatten@Take[Transpose[mats],-1];
     Return[posvecs]]
     
     
PrintDiffer[mat_]:=
Module[{mats,matx,colname,headings},
  (*largest to smallest*)
           mats=Reverse@Sort[mat];  
     (*drop sorting variable and positiion variable*)
          matx=Transpose@Drop[Drop[Transpose[mats],1],-1];
       (*print table groups*)
          colname={"genename","genid"};
          headings={None,colname};
          PrintTableFormGroupCore[matx,10,headings];
       (*new sorting variable*)
          posvecs=Flatten@Take[Transpose[mats],-1];
     Return[posvecs]]     
     
     
 PrintTableFormGroupCore[mat_,lenblock_,headings_]:=
 Module[{ncol,nrow,i,zerovec,numblocks,numtot,matpad,matset},
   {nrow,ncol}=Dimensions[mat];
   numblocks=Ceiling[nrow/lenblock];
   numtot=numblocks lenblock;
  zerovec={Table[0,{i,1,ncol}]};
   matpad=PadRight[mat,numtot,zerovec];
  matset=Partition[matpad,lenblock];
  PrintTableFormGroupOne[#,headings]& /@ matset;
  Return[Null]]
   

PrintTableFormGroupOne[mat_,headings_]:=
Module[{matx},
 matx=Select[mat,(#[[1]]=!=0)&];
 Print["            ",TableForm[matx,TableHeadings->headings]];
 Return[Null]]
 

 (*------------------------TransitionTables------------*)
 
TransitionTable[modelset_,tallyvec_,plotname1_,plotname2_,modelsetQ_]:=
 Module[{modelsetD,modelsetU,summaryname},
    If[modelsetQ=="full+poly" || modelsetQ=="full",
       modelsetD={"flat",   "lineD", "transitionD", "hockeyD", "sigmoidD", "impulseD", 
                              "stepD","impulse+D","step+D"};
       modelsetU={"flat","lineU", "transitionU", "hockeyU", "sigmoidU", "impulseU", 
                   "stepU", "impulse+U","step+U"}];
    If[modelsetQ=="reduced",
       modelsetD={"flat",   "lineD", "transitionD", "hockeyD", "sigmoidD", "impulseD", "stepD"};
       modelsetU={"flat","lineU", "transitionU", "hockeyU", "sigmoidU", "impulseU", "stepU"}];
       summaryname=StringJoin[{" rows: ", plotname1,";    columns: ", plotname2}];
      Print["Summary Table ",summaryname];
      TransitionTableCore["Down-Down",modelsetD,modelsetD,tallyvec];
      TransitionTableCore["Down-Up",modelsetD,modelsetU,tallyvec];
      TransitionTableCore["Up-Down",modelsetU,modelsetD,tallyvec];
     TransitionTableCore["Up-Up",modelsetU,modelsetU,tallyvec];
 Return[Null]]
 

TransitionTableCore[tablename_,modelset1_,modelset2_,tallyvec_]:=
Module[{len1,len2,mat,i,j,matlist,matsum,name,index,mat1,mat2,modname1,modname2,
rowtot,coltot,modnamevec1,modnamevec2,tallyvecS,
      tot, summaryname,summarytable},
     len1=Length[modelset1];
     len2=Length[modelset2];
     mat=Table[0,{i,1,len1},{j,1,len2}];
     tallyvecS=TallySelect[modelset1,modelset2,tallyvec];
     matlist=SummaryCore[modelset1,modelset2,mat,#]& /@ tallyvecS;
    mat=Apply[Plus,matlist];
    len=Length[mat];
    If[len==0,     Print["    ",tablename,"  total number in this table = ",0]];
    If[len>0,
     tot=Apply[Plus,Flatten[mat]];
    modnamevec1=ShortModelName[#]& /@ modelset1;          
    modnamevec2=ShortModelName[#]& /@ modelset2;          
    Print["    ",tablename,"  total number in this table = ",tot];
    summarytable=TableForm[mat,TableHeadings->{modnamevec1,modnamevec2},
    TableSpacing->{1,1}];
    Print["         ",summarytable]];
Return[Null]]

TallySelect[set1_,set2_,tallyvec_]:=
Module[{tallyvec1,tallyvec2},
  tallyvec1=TallySelectCore[set1,set2,#]& /@ tallyvec;
  tallyvec2=Select[tallyvec1,(#=!=0)&];
Return[tallyvec2]]
 
TallySelectCore[set1_,set2_,{{a_,b_},c_}]:=
Module[{test,res},
  test=MemberQ[set1,a] && MemberQ[set2,b];
  If[test==True, res={{a,b},c}, res=0];
Return[res]]

SummaryCore[modelset1_,modelset2_,mat_,{{mod1_,mod2_},c_}]:=
Module[{i,j,mat2},
  i=Position[modelset1,mod1][[1,1]];
  j=Position[modelset2,mod2][[1,1]];
  mat2=ReplacePart[mat,{i,j}->c];
 Return[mat2]]


 PrintTally[tallyvec_,plotname1_,plotname2_]:=
 Module[{tallyvec1,tallyvec2,tallyvec3,tallymat,tallymat1,colname},
      (*list counts first*)
       tallyvec1=Reverse[#]& /@ tallyvec;
       (*largest counts to smallest counts*)
         tallyvec2=Reverse@Sort[tallyvec1];
         tallymat= Flatten[#]& /@ tallyvec2;
         tallymat1=Reverse[#]& /@ tallymat;
    colname={plotname1,plotname2, "count"};
    Print@TableForm[tallymat1,TableHeadings->{None,colname}];
Return[Null]]



 
End[] 
EndPackage[]

(*
     {ztype,zmod,zstat,zfunc,zpar,zmse,
   	       zmodA,zstatA,zfuncA,zparA,zmseA,
   	        zmodB,zstatB,zfuncB,zparB,zmseB,
	     zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,pos}=zbestvec*)