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


BeginPackage["mfitcore`"]  

FitCoreSet::usage="FitCoreSet[pos,dataset,par,splitQ,kslope]"
GenFunc::usage="GenFunc[zfunc,zpar,s]"
SplitData::usage="SplitData[pos,dataset,splitQ]"


Clear[FitCoreSet,GenFunc,SplitData]

Begin["Private`"]	  

(*NOTE: zstat = fitting statistic
              zmod = particular function upward or downward
              zfunc = function 
              zpar =   parameter estimates
              zmse =  mse *)


FitCoreSet[pos_,dataset_,par_,splitQ_,modelsetQ_,kslope_]:=
Module[{plotname,xname,yname, dataname,data,geneid,genename,times,
        len,rangep,printposvec,fracdone,temp,
      zstatvec,ztarvec,zsevec,zmsevec,set,modelset,
      res,pair,genenamepos,geneidpos},
  (*Data*)
     {data,geneid,genename,times,xname,yname,plotname,dataname}=dataset;
  (*REPORT FRACTION COMPLETE IF MORE THAN 100 GENES*)
     len=Length[data];
    If[len>100,
       (*fractions of numbers of genes*)
           rangep={0} ~Join~ N[Range[9]/10];   
           printposvec=Quantile[Range[len], #]& /@ rangep;
        (*print at fractions of genes completed*)
           If[MemberQ[printposvec,pos]  && Length[data]>20,
              fracdone=Round[10 pos/len]/10//N;
              NotebookDelete[temp];
              temp=PrintTemporary["  fraction complete = ", fracdone]];
      If[pos==len,   NotebookDelete[temp]]];
(*FIT MODEL*)
      set=FitCore[pos,dataset, par,splitQ,modelsetQ,kslope];
     {zstatvec,zmodvec,zfuncvec,zparvec,zmsevec}= Transpose[set];
(*OUTPUT*)  
  res={pos,zstatvec,zmodvec,zfuncvec,zparvec,zmsevec};
   Return[res]]
  
  
FitCore[pos_,dataset_,{a_,b_,c_,d_,e_,f_,g_,h_,s_},splitQ_,modelsetQ_,kslope_]:=
Module[{xname,yname, plotname,dataname,data,geneid,genename,times,
             x,y, dataXY,      dataXYtest,dataXYfit,xtest,ytest,xfit,yfit,
          resfla,reslin,ressigD,ressigU,resdbsD,resdbsU,resgenD,resgenU,
          zmodlin,zstatlin,zparsig,zpardbs,zstatsig,  res},          
  (*DATA *)
   {data,geneid,genename,times,xname,yname,plotname,dataname}=dataset;
      y= data[[pos]];
      x=times;
     {dataXYfit, dataXYtest}=SplitData[pos,dataset,splitQ];
       {xfit,yfit}=Transpose[dataXYfit];
      {xtest,ytest}=Transpose[dataXYtest];

(*FULL +POLYNOMIAL SET*)
 If[modelsetQ=="full+poly",
    resfla=FitFlat[dataXYfit,xtest,ytest,y,s];
    reslin=FitLin[dataXYfit,xtest,ytest,y,s];
    zmodlin=reslin[[2]]; 
If[zmodlin=="lineD",
     ressig=FitSigD[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s}];
     resdbs=FitDbsD[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig];
     resgen=FitGenD[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,s},reslin,ressig,resdbs]];
If[zmodlin=="lineU",
     ressig=FitSigU[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s}];
     resdbs=FitDbsU[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig];
     resgen=FitGenU[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,s},reslin,ressig,resdbs]];
 respoly3=FitPoly3[dataXYfit,xtest,ytest,y,times,{a,b,c,d,s}];    
 respoly5=FitPoly5[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s}];    
 respoly7=FitPoly7[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,h,s}];    
 res={resfla,reslin,ressig,resdbs,resgen,respoly3,respoly5,respoly7}];
 
 (*FULL*)
  If[modelsetQ=="full",
     resfla=FitFlat[dataXYfit,xtest,ytest,y,s];
     reslin=FitLin[dataXYfit,xtest,ytest,y,s];
     zmodlin=reslin[[2]]; 
 If[zmodlin=="lineD",
      ressig=FitSigD[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s}];
      resdbs=FitDbsD[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig];
      resgen=FitGenD[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,s},reslin,ressig,resdbs]];
 If[zmodlin=="lineU",
      ressig=FitSigU[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s}];
      resdbs=FitDbsU[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig];
      resgen=FitGenU[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,s},reslin,ressig,resdbs]];
 res={resfla,reslin,ressig,resdbs,resgen}];
 
(*REDUCED SET*)
 If[modelsetQ=="reduced",
    resfla=FitFlat[dataXYfit,xtest,ytest,y,s];
    reslin=FitLin[dataXYfit,xtest,ytest,y,s];
    zmodlin=reslin[[2]]; 
If[zmodlin=="lineD",
     ressig=FitSigD[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s}];
     resdbs=FitDbsD[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig]];
 If[zmodlin=="lineU",
     ressig=FitSigU[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s}];
     resdbs=FitDbsU[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig]];
  res={resfla,reslin,ressig,resdbs}];
 Return[res]]  


(*-----------------------------------FIT FUNCTIONS-----------------------*)

FitFlat[dataXYfit_,xtest_,ytest_,y_,s_]:=
 Module[{fitfla,ypredfla,zstatfla,zparfla,zshapfla,zmsefla,resfla},
    fitfla=Quiet@LinearModelFit[dataXYfit, 1, s];
       ypredfla=fitfla[#]& /@ xtest;
       {zstatfla,zmsefla}=GoodnessOfFit[ytest,ypredfla,y];
       zmodfla="flat";
       zparfla=fitfla["BestFitParameters"];
        zfuncfla="fla";
        resfla={zstatfla,zmodfla,zfuncfla,zparfla,zmsefla};
Return[resfla]]  

FitLin[dataXYfit_,xtest_,ytest_,y_,s_]:=
 Module[{fitlin,ypredlin,zstatlin,zparlin,zshaplin,zmselin,reslin},
       fitlin=Quiet@LinearModelFit[dataXYfit, s, s];
       ypredlin=fitlin[#]& /@ xtest;
       {zstatlin,zmselin}=GoodnessOfFit[ytest,ypredlin,y];
       zmodlin=ModLIN[ypredlin];
       zparlin=fitlin["BestFitParameters"];
       zfunclin="lin";
       reslin={zstatlin,zmodlin,zfunclin,zparlin,zmselin};
   Return[reslin]]     

(*--------------------------------------------sigmoid------------------------------*)

FitSigD[dataXYfit_,xtest_,ytest_,y_,times_,kslope_,{a_,b_,c_,d_,s_}]:=
       FitSig[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s},{"sigD",sigDf}]


FitSigU[dataXYfit_,xtest_,ytest_,y_,times_,kslope_,{a_,b_,c_,d_,s_}]:=
       FitSig[dataXYfit,xtest,ytest,y,times,kslope,{a,b,c,d,s},{"sigU",sigUf}]

FitSig[dataXYfit_,xtest_,ytest_,y_,times_,kslope_,{a_,b_,c_,d_,s_},{shapename_,sigf_}]:=
Module[{fitsig,ypredsig,ytimessig,zparsig,zstatsig,zfuncsig,zmsesig,ressig},
              fitsig=Quiet@NonlinearModelFit[dataXYfit, sigf[{a,b,c,d,s}], {a, b,c,d}, s];
              ypredsig=fitsig[#]& /@ xtest;
              ytimessig=fitsig[#]& /@ times;
              zparsig={a,b,c,d}/.fitsig["BestFitParameters"];
              {zstatsig,zmsesig}=GoodnessOfFit[ytest,ypredsig,y];
              zfuncsig=shapename;
              zmodsig=ModSIG[ytimessig,times,kslope];    
              ressig={zstatsig,zmodsig,zfuncsig,zparsig,zmsesig};
Return[ressig]]


(*--------------------------double sigmoid--------------------------------*)

FitDbsD[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,s_},reslin_,ressig_]:=
       FitDbs[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig,{"dbsD",dbsDf}]

FitDbsU[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,s_},reslin_,ressig_]:=
       FitDbs[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,s},reslin,ressig,{"dbsU",dbsUf}]


FitDbs[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,s_},reslin_,ressig_,{shapename_,dbsf_}]:=
Module[{zstatlin,zstatsig,zparsig,a0,b0,c0,d0,fitdbs,ypreddbs,ytimesdbs,zstatdbs,zfuncdbs,zmsedbs,zfitdbs},
     zstatlin=reslin[[1]];
     zstatsig=ressig[[1]];
     zparsig=ressig[[4]]; 
     If[zstatsig<zstatlin,
                {a0,b0,c0,d0}=zparsig;
                 fitdbs=Quiet@NonlinearModelFit[dataXYfit, dbsf[{a,b,c,d,e,f,s}], 
       	                           {{a,a0},{b,(a0+c0)/2},{c,c0},{d,d0},{e,0},{f,0}}, s];    
                ypreddbs=fitdbs[#]& /@ xtest;
               ytimesdbs=fitdbs[#]& /@  times;
               zpardbs={a,b,c,d,e,f}/.fitdbs["BestFitParameters"];
	       {zstatdbs,zmsedbs}=GoodnessOfFit[ytest,ypreddbs,y];                
               zfuncdbs=shapename;
               zmoddbs=ModDBS[ytimesdbs,zfuncdbs];            
               zfitdbs=fitdbs;
               resdbs={zstatdbs,zmoddbs,zfuncdbs,zpardbs,zmsedbs}];
      If[zstatsig>=zstatlin,resdbs=reslin];                  
Return[resdbs]]


(*----------------generalized double sigmoid--------------------------------------*)

FitGenD[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,g_,s_},reslin_,ressig_,resdbs_]:=
       FitGen[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,s},reslin,ressig,resdbs,{"genD",genDf}]

FitGenU[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,g_,s_},reslin_,ressig_,resdbs_]:=
       FitGen[dataXYfit,xtest,ytest,y,times,{a,b,c,d,e,f,g,s},reslin,ressig,resdbs,{"genU",genUf}]


FitGen[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,g_,s_},reslin_,ressig_,resdbs_,
                 {shapename_,genf_}]:=
Module[{zstatlin,zstatsig,zstatdbs,zparsig,zpardbs,a0,b0,c0,d0,e0,f0, 
              fitgen,ypredgen,ytimesgen,zstatgen,zfuncgen,zmsegen,zfitgen},  
              zstatlin=reslin[[1]];
              zstatsig=ressig[[1]];
              zstatdbs=resdbs[[1]];
              zpardbs=resdbs[[4]]; 
              (*needed to avoid double sigmoid that reduces to linear*)
               If[zstatdbs<zstatlin,
                    {a0,b0,c0,d0,e0,f0}=zpardbs;
                     fitgen=Quiet@NonlinearModelFit[dataXYfit, genf[{a,b,c,d,e,f,g,s}], 
                                   {{a,a0},{b,b0},{c,c0},{d,d0},{e,e0},{f,f0},{g,1,.5,2}}, s];    
                    ypredgen=fitgen[#]& /@ xtest;
                     ytimesgen=fitgen[#]& /@ times;
                     zpargen={a,b,c,d,e,f,g}/.fitgen["BestFitParameters"];
      	             {zstatgen,zmsegen}=GoodnessOfFit[ytest,ypredgen,y];              
                     zfuncgen=shapename;
                     zmodgen=ModGEN[ytimesgen,zfuncgen];            
                     zfitgen=fitgen;                         
                    resgen={zstatgen,zmodgen,zfuncgen,zpargen,zmsegen}];
               If[zstatdbs>=zstatlin, resgen=resdbs];
Return[resgen]]



FitPoly3[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,s_}]:=
Module[{fitpoly,ypredpoly,ytimespoly,zparpoly,zstatpoly,zfuncpoly,zmsepoly,respoly},
              fitpoly=Quiet@NonlinearModelFit[dataXYfit, poly3f[{a,b,c,d,s}], {a, b,c,d}, s];
              ypredpoly=fitpoly[#]& /@ xtest;
              ytimespoly=fitpoly[#]& /@ times;
              zparpoly={a,b,c,d}/.fitpoly["BestFitParameters"];
              {zstatpoly,zmsepoly}=GoodnessOfFit[ytest,ypredpoly,y];
              zfuncpoly="poly3";
              zmodpoly="polynomial3";    
              respoly={zstatpoly,zmodpoly,zfuncpoly,zparpoly,zmsepoly};
Return[respoly]]

FitPoly5[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,s_}]:=
Module[{fitpoly,ypredpoly,ytimespoly,zparpoly,zstatpoly,zfuncpoly,zmsepoly,respoly},
              fitpoly=Quiet@NonlinearModelFit[dataXYfit, poly5f[{a,b,c,d,e,f,s}], {a, b,c,d,e,f}, s];
              ypredpoly=fitpoly[#]& /@ xtest;
              ytimespoly=fitpoly[#]& /@ times;
              zparpoly={a,b,c,d,e,f}/.fitpoly["BestFitParameters"];
              {zstatpoly,zmsepoly}=GoodnessOfFit[ytest,ypredpoly,y];
              zfuncpoly="poly5";
              zmodpoly="polynomial5";    
              respoly={zstatpoly,zmodpoly,zfuncpoly,zparpoly,zmsepoly};
Return[respoly]]

FitPoly7[dataXYfit_,xtest_,ytest_,y_,times_,{a_,b_,c_,d_,e_,f_,g_,h_,s_}]:=
Module[{fitpoly,ypredpoly,ytimespoly,zparpoly,zstatpoly,zfuncpoly,zmsepoly,respoly},
              fitpoly=Quiet@NonlinearModelFit[dataXYfit, poly7f[{a,b,c,d,e,f,g,h,s}], {a, b,c,d,e,f,g,h}, s];
              ypredpoly=fitpoly[#]& /@ xtest;
              ytimespoly=fitpoly[#]& /@ times;
              zparpoly={a,b,c,d,e,f,g,h}/.fitpoly["BestFitParameters"];
              {zstatpoly,zmsepoly}=GoodnessOfFit[ytest,ypredpoly,y];
              zfuncpoly="poly7";
              zmodpoly="polynomial7";    
              respoly={zstatpoly,zmodpoly,zfuncpoly,zparpoly,zmsepoly};
Return[respoly]]



 (*-------------------------------FUNCTIONS---------------------------------*)
 
 linf[{a_,b_,s_}]:=a + b s
  
  
 sigUf[{a_,b_,c_,d_,s_}]:=  a + (c-a) Expit[ b (s-d) ]
 sigDf[{a_,b_,c_,d_,s_}]:=  a + (c-a) Expit[ b (d-s) ]
 
 dbsUf[{a_,b_, c_,d_, e_, f_, s_}] := 
             (1/b)  (a + (b-a)  Expit[ e (s-d)]) (c+ (b-c) Expit[ -e (s-f) ] )   
 dbsDf[{a_,b_, c_,d_, e_, f_, s_}] := 
            (1/b)  (a  + (b-a)  Expit[ e (d-s)])  (c+ (b-c) Expit[ -e (f-s) ] )


genUf[{a_,b_, c_,d_, e_, f_, g_,s_}] := 
             (1/b)  (a + (b-a)  Expit[ e (s-d)]) (c+ (b-c) Expit[ -e g (s-f) ] )   
 genDf[{a_,b_, c_,d_, e_, f_, g_,s_}] := 
            (1/b)  (a  + (b-a)  Expit[ e (d-s)])  (c+ (b-c) Expit[ -e g (f-s) ] )

poly3f[{a_,b_, c_,d_,s_}] := a+ b s + c s^ +d s^3 

poly5f[{a_,b_, c_,d_, e_, f_,s_}] := a+ b s + c s^ +d s^3 + e s^4 + f s^5  

poly7f[{a_,b_, c_,d_, e_, f_, g_,h_,s_}] := a+ b s + c s^ +d s^3 + e s^4 + f s^5 + g s^6 + h s^7 

Expit[s_]:= Exp[s]/(1+Exp[s])

(*-----------------------------------------------MODEL FILTERS----------------------------------*)


ModLIN[ypred_]:=
Module[{zmod},
          If[ypred[[1]] > ypred[[-1]], zmod="lineD",zmod="lineU"]; 
Return[zmod]]



ModSIG[ytimessig_,times_,kslope_]:=
 Module[{ ymax,ymin,xmax,xmin,kslopex,first,last,yslope,yderiv,yderivmin,
     yslopeavg,yderivmax,zmod,yslopemin,type},
              (*overall max and min*)
                    ymax=Max@ytimessig;
                    ymin=Min@ ytimessig;
               (*first and last response*)
                     first=ytimessig[[1]];
                     last=ytimessig[[-1]];        
             (*compute second derivative to distinguish sigmoid from hockey stick*)        
               xslope=Drop[times,1]-Drop[times,-1];
               yslope=(Drop[ytimessig,1]-Drop[ytimessig,-1])/xslope;
               yslopefirst=Abs@yslope[[1]];
               yslopelast=Abs@yslope[[-1]];
               xmax=Max@times;
               xmin=Min@times;
               (*adjust for scale*)
               yslopeavg=((ymax-ymin)/(xmax-xmin));
               kslopex= kslope;  
             If[yslopefirst >kslopex && yslopelast <=kslopex, type="hoc"];
             If[yslopefirst <=kslopex && yslopelast > kslopex, type="hoc"];
             If[yslopefirst > kslopex && yslopelast > kslopex, type="mid"];
             If[yslopefirst <= kslopex && yslopelast <= kslopex, type="sig"];
             (* Print["  kslope,yslopefirst, yslopelast ",  {type,kslopex,yslopefirst,yslopelast}];	*)
             (*final assigment*)
             zmod=Which[
                    type=="sig" && last > first,     "sigmoidU",
                    type=="sig" && last  < first,     "sigmoidD",
                    type=="mid" && last > first,    "transitionU",
	           type=="mid" && last  < first,    "transitionD",
                    type=="hoc" && last > first,    "hockeyU",
                    type=="hoc" && last  < first,    "hockeyD",
                   True, 					    "sigmoidX"];
           Return[zmod]]




 ModDBS[ytimesdbs_,zfuncdbs_]:=
 Module[{maxQ,minQ,monoQ,monoQD,monoQU,zmod},
                {minQ,maxQ,monoQU,monoQD}=ModComp[ytimesdbs];
          If[zfuncdbs=="lin",  zmod=												   "impulseL"];
         If[zfuncdbs=!="lin",  zmod=Which[
              maxQ==True,                                    							                    "impulseU",
              minQ==True, 			                  								            "impulseD",
	      maxQ==False && minQ ==False && monoQD==True && monoQU==False,             "stepD",
	      maxQ==False && minQ ==False && monoQU==True && monoQD==False,             "stepU",
		       True, 														     "impulseX"]] ;
           Return[zmod]]

  
ModGEN[ytimesgen_,zfuncgen_]:=
 Module[{maxQ,minQ,monoQ,monoQD,monoQU,zmop},
                {minQ,maxQ,monoQU,monoQD}=ModComp[ytimesgen];
                If[zfuncgen=="lin",      zmod=										     "impulse+L"];
                 If[zfuncgen=!="lin",    zmod=Which[
                       maxQ==True,                                     							               "impulse+U",
                       minQ==True, 			                   							               "impulse+D",
		       maxQ==False && minQ ==False &&  monoQD==True && monoQU==False,      "step+D",
		       maxQ==False && minQ ==False && monoQU==True && monoQD==False,       "step+U",
		       True, 														        "impulse+X"]] ;
           Return[zmod]]

  
ModComp[ytimesdbs_]:=
Module[{  max,min,first,last,dmax,dmin,
                  maxQ,minQ,monoQ,monoQD,monoQU},
      (*overall max and min*)
                    max=Max@ytimesdbs;
                    min=Min@ ytimesdbs;
               (*first and last*)
                     first=ytimesdbs[[1]];
                     last=ytimesdbs[[-1]];        
             (*monotonically changing*)  
                 monoQD=(max==first &&  min==last);
                 monoQU=(min==first && max==last);
                 monoQ=monoQD || monoQU;
            (*difference between*)  
                dmax=Min[Abs[first-max],Abs[last-max]];
                dmin=Min[Abs[first-min],Abs[last-min]];
            (*max or min not at first or last indicating impulse*)  
                maxQ=dmax>dmin && max > first && max > last;
                minQ=dmax <=dmin && min < first && min < last;
          Return[{minQ,maxQ,monoQU,monoQD}]]

 
 

(*----------------------------------------Goodness of Fit-----------------------*)

GoodnessOfFit[ytest_,ypred_,y_]:=
Module[{error,mean,zmse,ztar,ztarx,zstat,zse,zseM},
       error=ytest-ypred;
       zmse=Mean[error^2];
       zseM=Sqrt@zmse;
       ztar= Max[y]-Min[y];
        If[ztar==0, Print["warning flat = 0 range"]; ztarx=.01, ztarx=ztar];
       zstat= zseM/ztarx;
   Return[{zstat,zseM}]]

(*----------------------------------------Split Data-----------------------*)

SplitData[pos_,dataset_,splitQ_]:=
Module[{data,geneid,genename,times,xname,yname,plotname,dataname,
           splittype,splitfittest,
             y,x,len,k, rfit,rtest,dataXYtest,dataXYfit},
 (*Data*)     
  {data,geneid,genename,times,xname,yname,plotname,dataname}=dataset;
   y= data[[pos]];
   x=times;
    len=Length[times];
    {splittype,splitfittest}=splitQ;
    If[splittype=="alternating",
     If[OddQ[len], k=(len-1)/2, k=len/2];
      rtest=2 Range[k];
      rfit= Complement[Range[len],rtest]];
    If[splittype=="custom",  
       rfit= splitfittest;
      rtest= Complement[Range[len],splitfittest]];
    (*DATA SPLIT*)
    dataXY=Transpose[{x,y}];
      dataXYtest=dataXY[[rtest]];
      dataXYfit=dataXY[[rfit]];
Return[{dataXYfit,dataXYtest}]]



(*----------------------------Generate Functions------------------*)

 GenFunc[{zfunc_,zpar_},s_]:=
  Module[{zpars,func},
   zpars=Join[zpar,{s}];
  If[zfunc=="fla",        func=zpar];
   If[zfunc=="lin",       func=linf[zpars]];
   If[zfunc=="sigU",   func=sigUf[zpars]];
   If[zfunc=="sigD",   func=sigDf[zpars]];
   If[zfunc=="dbsU",  func=dbsUf[zpars]];
    If[zfunc=="dbsD",  func=dbsDf[zpars]];
   If[zfunc=="genU",   func=genUf[zpars]];
    If[zfunc=="genD",  func=genDf[zpars]];
    If[zfunc=="poly3",    func=poly3f[zpars]];
    If[zfunc=="poly5",    func=poly5f[zpars]];
    If[zfunc=="poly7",    func=poly7f[zpars]];
   Return[func]]



End[] 
EndPackage[]


