(*Mathematica:: Version 8.0      *)
(*:Context:: "mfit`"                     *)
(*:Title::  MFit            		*)
(*:Summary::   mfitgene                 	*)
(*:References::                         *)
(*Date  2013 Stuart G. Baker   *)

(*Input
FitCoreSet           mfitcore.m*)


BeginPackage["mfitkey`", "mfitplot`","mfitcore`"]  


MFitKey::usage="MFitKey[dataset,{splitQ,newfitQ,newbestQ},showdatainput,
           {newnameCORE,storednameCORE},{zbest,znext},par]"           

Clear[MFitKey]



Begin["Private`"]	  



  MFitKey[dataset_,{splitQ_,exportplotQ_,newfitQ_,newbestQ_},showdatainput_,modelsetQ_,kslope_,
     {newnameCORE_,storednameCORE_},{zbest_,znext_},par_]:=
  Module[{modelset,xname,yname,plotname,dataname,data,geneid,genename,times,
               numgenes,
               set,newname,storedname, posvec,
               zstatmat,zmodmat,zfuncmat,zmsemat,
               zbestmat,zbestset,zbest0,znext0,splitQ0,vec,modlist},
     (*SUPPRESS OVERFLOW UNDERFLOW MESSAGE*) 
         Off[General::ovfl];       
         Off[General::unfl];                  
     (*basic quantities*)
         {data,geneid,genename,times,xname,yname,plotname,dataname}=dataset;
         numgenes=Length[data];
   (*ADJUST NAME FOR MODELSET*)      
    If[modelsetQ=="full+poly", nameM="P"];
    If[modelsetQ=="full", 	       nameM="F"];
    If[modelsetQ=="reduced", nameM="R"];
    
         
 (*NEW FIT AND NEW BEST CALCULATION*)
    If[newfitQ==True, 
                  newname=StringJoin[newnameCORE,nameM,dataname,".m"];
                  Print["DATASET: ",dataname, "  export initial fit result to ",newname];
                   set=FitCoreSet[#,dataset,par,splitQ,modelsetQ,kslope]&  /@ Range[numgenes];
             Export[newname,set];
             (*compute best fitting model*)   
	        {posvec,zstatmat,zmodmat,zfuncmat,zparmat,zmsemat}=Transpose[set];
	         zbestmat=FindBestModel[zstatmat,zmodmat,zfuncmat,zparmat,zmsemat,#,
	                      zbest,znext,modelsetQ,geneid]& /@ Range[numgenes];
	      newnamebestQ=StringJoin[newnameCORE,nameM,dataname,"best",".m"];
	     Print["                export best fit result to ",newnamebestQ];
             zbestset={zbestmat,zbest,znext,splitQ};  
             Export[newnamebestQ,zbestset]];
  
  (*STORED FIT AND NEW BEST CALCULATION-- not generally used*)   
       If[newfitQ==False && newbestQ==True,
            storedname=StringJoin[storednameCORE,nameM,dataname,".m"];
            set=Quiet@Drop[ReadList[storedname],1][[1]];
           newnamebestQ=StringJoin[storednameCORE,dataname,"best",".m"];
           Print["DATASET: ",dataname, " import initial fit from ",storedname, " export best fit to ",newnamebestQ];
          (*compute best model*)   
           {posvec,zstatmat,zmodmat,zfuncmat,zparmat,zmsemat}=Transpose[set];
          zbestmat=FindBestModel[zstatmat,zmodmat,zfuncmat,zparmat,zmsemat,#,zbest,
           znext,modelsetQ,geneid]& 
                      /@ Range[numgenes];
           zbestset={zbestmat,zbest,znext,splitQ};  
          Export[newnamebestQ,zbestset]]; 
 
 (*STORED FIT AND STORED BEST FIT CALCULATION*)   
       If[newfitQ==False && newbestQ==False,
            storednamebest=StringJoin[storednameCORE,nameM,dataname,"best",".m"];
       	     zbestset=Quiet@Drop[ReadList[storednamebest],1][[1]];
       	    If[showdatainput,
       	    Print["DATASET: ",dataname,"    import best fit result from ",storednamebest]]];
   
   Return[zbestset]]





(*---------------------------------Find Best Model---------------------------------------*)
(*NOTE:zmod =  model that is more specific 
              zstat= goodness of fit statistic = sqrt (mse)/ (max-min)
              zfunc = functional form: fla, lin, sig, dbs, gen, poly3, poly5, poly7 
              zpar = parameter estimates for zfunc
              zmse = MSE of model *)

FindBestModel[zstatmat_,zmodmat_,zfuncmat_,zparmat_,zmsemat_,pos_,zbest_,znext_,modelsetQ_,geneid_]:=
Module[{res},
If[modelsetQ=="full+poly",
  res= FindBestModelF[zstatmat,zmodmat,zfuncmat,zparmat,zmsemat,pos,zbest,znext,geneid]];
If[modelsetQ=="reduced" || modelsetQ=="full", 
   res= FindBestModelR[zstatmat,zmodmat,zfuncmat,zparmat,zmsemat,pos,zbest,znext,geneid]];
Return[res]]


FindBestModelF[zstatmat_,zmodmat_,zfuncmat_,zparmat_,zmsemat_,pos_,zbest_,znext_,geneid_]:=
Module[{ztype,zmod,zstat,zfunc,zpar,zmse,
              zmodA,zstatA,zfuncA,zparA,zmseA,
              zmodB,zstatB,zfuncB,zparB,zmseB,
              zstatvecP,zparvecP,zstatC,posC,zmodC,zfuncC,zparC,zmseC,
               zstatvecR,zparvecR,zmodvecR,
               zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,res,
              posA,posAx,posB},
(*select position*)
        zmodvec=zmodmat[[pos]];
        zstatvec=zstatmat[[pos]];
        zfuncvec=zfuncmat[[pos]];
        zparvec=zparmat[[pos]]; 
        zmsevec=zmsemat[[pos]]; 
 (*polynomial*)
        zstatvecP=Take[zstatvec,-3];
         zparvecP=Take[zparvec,-3];
         zmodvecP=Take[zmodvec,-3];
         zstatC=Min[zstatvecP]; 
         posC= Position[zstatvec,zstatC][[1,1]];
          zmodC=zmodvec[[posC]];
        zfuncC=zfuncvec[[posC]];
          zparC=zparvec[[posC]];
          zmseC=zmsevec[[posC]];
    (*biologically relevant; drop last three corresponding to polynomials*)
         zstatvecR=Drop[zstatvec,-3];
         zparvecR=Drop[zparvec,-3];
         zmodvecR=Drop[zmodvec,-3];
    (*minimum of all models excluding polynomial*)
            zstatA=Min@zstatvecR;
            posA= Position[zstatvec,zstatA][[1,1]];
            zmodA=zmodvec[[posA]];
            zfuncA=zfuncvec[[posA]];
            zparA=zparvec[[posA]];
            zmseA=zmsevec[[posA]];
      (*minimum of all models with fewer parameters*) 
             (*fewest paramters less than zstatA + znstatnex*) 
            zstatB=ComputeZStatB[zstatvec,zstatA,znext,posA];           
             posB= Position[zstatvec,zstatB] [[1,1]];
            zmodB=zmodvec[[posB]];
             zfuncB=zfuncvec[[posB]];
             zparB=zparvec[[posB]];
             zmseB=zmsevec[[posB]];
   (*poor fit  for all models*)
      If[zstatA> zbest && zstatC >zbest, 
                    ztype="D";  
                    zstat=100;
                    zmod="poorfit";
                    zfunc="poorfit";
                    zpar={0};
                    zmse=100];
      (*poor fit biologically relevant; good fit polynomial*)
      If[zstatA> zbest && zstatC <=zbest, 
                    ztype="C";  
                    zmod=zmodC;
                    zstat=zstatC;
                    zfunc=zfuncC;
                    zpar=zparC;
                    zmse=zmseC];                
    If[zstatA <= zbest  && zstatB <= zstatA+znext,    
             ztype="B"; 
               zmod=zmodB;
               zstat=zstatB;
                zfunc=zfuncB;
               zpar=zparB;
              zmse=zmseB];
       If[zstatA <= zbest  && zstatB > zstatA+znext,
            ztype="A"; 
             zmod=zmodA;
	     zstat=zstatA;
	      zfunc=zfuncA;
	      zpar=zparA;
              zmse=zmseA];
  (*output*)
     res={ztype,zmod,zstat,zfunc,zpar,zmse,
              zmodA,zstatA,zfuncA,zparA,zmseA,
              zmodB,zstatB,zfuncB,zparB,zmseB,
               zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,pos};
      Return[res]];

FindBestModelR[zstatmat_,zmodmat_,zfuncmat_,zparmat_,zmsemat_,pos_,zbest_,znext_,geneid_]:=
Module[{ztype,zmod,zstat,zfunc,zpar,zmse,
              zmodA,zstatA,zfuncA,zparA,zmseA,
              zmodB,zstatB,zfuncB,zparB,zmseB,
              zstatvecP,zparvecP,zstatC,posC,zmodC,zfuncC,zparC,zmseC,
               zstatvecR,zparvecR,zmodvecR,
               zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,res,
              posA,posAx,posB},
(*select position*)
        zmodvec=zmodmat[[pos]];
        zstatvec=zstatmat[[pos]];
        zfuncvec=zfuncmat[[pos]];
        zparvec=zparmat[[pos]]; 
        zmsevec=zmsemat[[pos]]; 
    (*minimum of all models excluding polynomial*)
           zstatA=Min@zstatvec;
           posA= Position[zstatvec,zstatA][[1,1]];
            zmodA=zmodvec[[posA]];
            zfuncA=zfuncvec[[posA]];
            zparA=zparvec[[posA]];
            zmseA=zmsevec[[posA]];
      (*fewest paramters less than zstatA + znstatnex*) 
            zstatB=ComputeZStatB[zstatvec,zstatA,znext,posA];
             posB= Position[zstatvec,zstatB] [[1,1]];
             zmodB=zmodvec[[posB]];
             zfuncB=zfuncvec[[posB]];
             zparB=zparvec[[posB]];
             zmseB=zmsevec[[posB]];
   (*poor fit  for all models*)
      If[zstatA> zbest, 
                    ztype="C";  
                    zstat=100;
                    zmod="poorfit";
                    zfunc="poorfit";
                    zpar={0};
                    zmse=100];
    If[zstatA <= zbest  && zstatB <= zstatA+znext,    
             ztype="B"; 
               zmod=zmodB;
               zstat=zstatB;
                zfunc=zfuncB;
               zpar=zparB;
              zmse=zmseB];
       If[zstatA <= zbest  && zstatB > zstatA+znext,
            ztype="A"; 
             zmod=zmodA;
	     zstat=zstatA;
	      zfunc=zfuncA;
	      zpar=zparA;
              zmse=zmseA]; 
  (*output*)
     res={ztype,zmod,zstat,zfunc,zpar,zmse,
              zmodA,zstatA,zfuncA,zparA,zmseA,
              zmodB,zstatB,zfuncB,zparB,zmseB,
               zmodvec,zstatvec,zfuncvec,zparvec,zmsevec,pos};
      Return[res]];


ComputeZStatB[zstatvec_,zstatA_,znext_,posA_]:=
Module[{CurveSetB,CurveSubsetB,len,zstatB,posAX},
           (*if ties take first *)
              posAX=Max[posA-1,1];
              CurveSetB=Take[zstatvec,posAX];
              CurveSubsetB=Select[CurveSetB,(#<=zstatA+znext)&];
	       len=Length[CurveSubsetB];
	       If[len==0,  zstatB=CurveSetB[[-1]]];
	       If[len>0,    zstatB=CurveSubsetB[[1]]];
         Return[zstatB]]





End[] 
EndPackage[]


