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

(*Input
SplitData			mfitcore.m
GenFunc			mfitcore.m*)


BeginPackage["mfitform`","mfitcore`"]  

FormPair::usage="FormPair[posvec,zbestmat1,zbestmat2,dataset1,dataset2,par,splitQ,kform]"
FormPairVec::usage="FormPairVec[pos,zbestmat1,zbestmat2,dataset1,dataset2,par,splitQ,kform]"


Clear[FormPair,FormPairVec]

Begin["Private`"]	  


FormPair[posvec_,zbestmat1_,zbestmat2_,dataset1_,dataset2_,par_,splitQ_,kform_]:=
Module[{res}, 
 res=FormPairVec[#,zbestmat1,zbestmat2,dataset1,dataset2,par,splitQ,kform]&  /@ posvec;
Return[res]]


FormPairVec[pos_,zbestmat1_,zbestmat2_,dataset1_,dataset2_,par_,splitQ_,kform_]:=
Module[{zbestvec1,zbestvec2, 
           xrange1,func1,difx1,dify1,mod1,
           xrange2,func2,difx2,dify2,mod2,
           s,zmod,zmodfull,difMx,difMy,vec1,vec2,
           yA,yB,yM, xA,xB,xMID,hC,hM,htype,res},
  (*basics*)
       zbestvec1=zbestmat1[[pos]];
       zbestvec2=zbestmat2[[pos]];
       {func1,xrange1,difx1,dify1,zmod1}=GenFuncVec[zbestvec1,dataset1,par,splitQ];
       {func2,xrange2,difx2,dify2,zmod2}=GenFuncVec[zbestvec2,dataset2,par,splitQ];
       s=par[[-1]];
       zmod=zmod1; (*same model in both data sets*)
      difMx=Max[difx1,difx2];
       difMy=Max[dify1,dify2];
  (*FLAT*)
     If[zmod=="flat", 
         hC=0;
         hM=Abs[func2-func1]];
  (*LINEAR*)
   If[zmod=="lineU" || zmod=="lineD", 
       (*ymin and ymax is vec*)
       vec1=FormLin[func1,xrange1,s];
      vec2=FormLin[func2,xrange2,s];
       {yA, yB}=Abs[vec1-vec2];
       hC=0;
      hM=Mean[{yA,yB}]/difMy];      
    If[MemberQ[{"sigmoidD","sigmoidU","hockeyD","hockeyU",
           "stepD","stepU","step+D","step+U","transitionD","transitionU"},zmod],
      (*xmidpoint, ylevelA, ylevelB*)
       vec1=FormSig[func1,xrange1,s];
       vec2=FormSig[func2,xrange2,s];
      {xMID, yA, yB}=Abs[vec1-vec2];
       hC=xMID/difMx;
      hM=Mean[{yA,yB}]/difMy];    

  If[zmod=="impulseU" || zmod=="impulse+U", 
     vec1=FormImpulseU[func1,xrange1,s];
     vec2=FormImpulseU[func2,xrange2,s];
     {xM, yA, yB,yM}=Abs[vec1-vec2];
     hC=xM/difMx;
     hM=Mean[{yA,yB,yM}]/difMy];

  If[zmod=="impulseD"|| zmod=="impulse+D", 
  (*xM at impulse, yA at flat, yB at other flat, yM at impulse*)
     vec1=FormImpulseD[func1,xrange1,s];
     vec2=FormImpulseD[func2,xrange2,s];
    {xM, yA, yB,yM}=Abs[vec1-vec2];
     hC=xM/difMx;
     hM=Mean[{yA,yB,yM}]/difMy];
     
     
    (*POLY*)
   If[zmod=="polynomial" ,
      vec1=FormPoly[func1,xrange1,s];
      vec2=FormPoly[func2,xrange2,s];
       {yA, yB}=Abs[vec1-vec2];
       hC=0;
      hM=Mean[{yA,yB}]/difMy];        
             
     
 (*output*)
   If[hC >= kform && hM <  kform, htype=1];  (*HC*)
  If[hC <   kform && hM >= kform, htype=2];  (*HM*)
  If[hC >= kform && hM >= kform, htype=3]; (*both*)
  If[hC <   kform && hM <   kform, htype=0];  (*same*)
   res={zmod,htype, hC,hM,pos};

 Return[res]]
 
 
  
  GenFuncVec[zbestvec_,dataset_,par_,splitQ_]:=
    Module[{ztype,zmod,zstat,zfunc,zpar,zmse,
	        zmodA,zstatA,zfuncA,zparA,zmseA,
	        zmodB,zstatB,zfuncB,zparB,zmseB,
	         zmodvecX,zstatvecX,zfuncvecX,zparvecX,zmsevecX,pos,
               data,geneid,genename,times,xname,yname,plotname,dataname,
               y,x,s,func,
                xmin,xmax,ymin,ymax,difx,dify,xrange},
  (*INPUTS*)  
             {ztype,zmod,zstat,zfunc,zpar,zmse,
	        zmodA,zstatA,zfuncA,zparA,zmseA,
	        zmodB,zstatB,zfuncB,zparB,zmseB,
	         zmodvecX,zstatvecX,zfuncvecX,zparvecX,zmsevecX,pos}= zbestvec;
           {data,geneid,genename,times,xname,yname,plotname,dataname}=dataset;
       (*DATA*)
  	        y= data[[pos]];
  	        x= times;  	      
                s=par[[-1]];
                func=GenFunc[{zfunc,zpar},s];   
          (*Min and Max*)
                 xmin=Min[times];
                  xmax=Max[times];
 	          ymin=Min[y];
 	         ymax=Max[y];
 	     	   difx= xmax-xmin;
                 dify=ymax-ymin;
                 xrange={xmin,xmax};
        Return[{func,xrange,difx,dify,zmod}]]

 
(*----------------------------------------------FORMS------------------------------------------*)

FormPoly[func_,xrange_,s_]:=
Module[{ylist, yslope, yA,yB,res},
         {ylist, yslope}=GenYList[func,xrange,s];
         yA=Min@ylist;
          yB=Max@ylist;
         res={yA,yB};
   Return[res]]


FormLin[func_,xrange_,s_]:=
Module[{ylist, yslope, yA,yB,res},
         {ylist, yslope}=GenYList[func,xrange,s];
         yA=ylist[[1]];
          yB=ylist[[-1]];
         res={yA,yB};
   Return[res]]
     
FormSig[func_,xrange_,s_]:=
Module[{ylist,dify,yTAR,ylistA,ylistB,yA,yB,xA,xB,len,xMID,res},
         (*predicted values*)
         {ylist, yslope}=GenYList[func,xrange,s];
          yslopemax=Max[yslope];
          xMID=Position[yslope, yslopemax][[1,1]];
          yA=ylist[[1]];
          yB=ylist[[-1]];
          res={xMID,yA,yB};
   Return[res]]
     
 FormImpulseD[func_,xrange_,s_]:=
 Module[{ylist,yA,yB,xA,xB,yM,xM,res,xMID},
      {ylist, yslope}=GenYList[func,xrange,s];
      yM=Min[ylist];
        xM=Position[ylist,yM][[1,1]];
    {xMID, yA,yB}=FormSig[func,xrange,s];
         res={xM,yA,yB,yM};
      Return[res]]
 
 FormImpulseU[func_,xrange_,s_]:=
  Module[{ylist,ysplit,yA,yB,xA,xB,yM,xM,res,xMID},
         {ylist, yslope}=GenYList[func,xrange,s];
         yM=Max[ylist];
           xM=Position[ylist,yM][[1,1]];
       {xMID, yA,yB}=FormSig[func,xrange,s];
        res={xM,yA,yB,yM};
    Return[res]]
         
  GenYList[func_,{xmin_,xmax_},s_]:=
  Module[{dif,mintime,maxtime,ylist,yslope},
         dif=xmax-xmin;
          ylist=Flatten@Table[func, {s, xmin, xmax, dif/40}];
          yslope=Abs[Drop[ylist,-1]-Drop[ylist,1]];
 Return[{ylist,yslope}]];


 
 

 
End[] 
EndPackage[]



