fit_ppf.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
#define MAX_FUNC_NR 7
#define MAX_METAB_NR 3 
enum {MODEL_UNKNOWN, MODEL_PF, MODEL_HILL, MODEL_GCDF, MODEL_EP, MODEL_MPF, MODEL_MHILL};
int model=MODEL_UNKNOWN;
int func_par_nr[MAX_FUNC_NR]={0,5,5,5,6,5,5};
int func_type[MAX_FUNC_NR]={0,MF_EMAYER_PAR,MF_EHILL_PAR,MF_PPFIGAM,MF_PF_MU,MF_PF3M_PAR,MF_HILL3M_PAR};
double func_pf(int parNr, double *p, void*);
double func_hill(int parNr, double *p, void*);
double func_gcdf(int parNr, double *p, void*);
double func_ep(int parNr, double *p, void*);
double func_mpf(int parNr, double *p, void*);
double func_mhill(int parNr, double *p, void*);
double (*func_list[MAX_FUNC_NR])(int parNr, double *p, void*)=
  {NULL,func_pf,func_hill,func_gcdf,func_ep,func_mpf,func_mhill};
/*****************************************************************************/
#define MAX_MINMEAS_NR 4
enum {MINMEAS_UNKNOWN, MINMEAS_OLS, MINMEAS_LAD, MINMEAS_ODR};
int min_meas=MINMEAS_UNKNOWN;
/*****************************************************************************/
double *x, *ymeas[1+MAX_METAB_NR], *yfit[1+MAX_METAB_NR], *w; /* Local */
int origdataNr=0, fitdataNr=0, parNr=5, metabNr=0;
double frwgt[1+MAX_METAB_NR];
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double lastWSS=0.0;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Fits an empirical function to plasma parent (unchanged) tracer",
  "fraction curve (1), where the fractions are between 0 and 1.",
  " ",
  "Usage: @P [Options] fractionfile [fitfile]",
  " ",
  "Options:",
  " -model=<PF | HILL | GCDF | EP | MPF | MHILL>",
  "     Select the function (see descriptions below); default is HILL.",
  "     Use MPF or MHILL to fit functions to 2-3 metabolite fraction curves",
  "     instead of parent fraction curve.",
/* Currently in testing, therefore hidden.
  "     HILL2 can be used to save parameters for the 2nd Hill function",
  "     representation.",
*/
  " -min=<OLS|LAD|ODR>",
  "     Sum-of-squares (ordinal least squares, OLS) is minimized by default,",
  "     but optionally, sum of absolute deviations (LAD), or",
  "     sum-of-squares of Akaho distances (orthogonal distances, ODR) ",
  "     can be selected. ODR can be used only with PF and HILL.",
  " -delaymin=<value>",
  "     Set the lower limit for delay parameter.",
  " -mdelay=<separated | joint>",
  "     Delay parameter for each metabolite is fitted separately, or",
  "     all metabolites are assumed to share common delay time (default).",
  " -fix0",
  "     Parent fraction at time zero is fixed to zero sample, if available.",
  " -a=<value>, -b=<value>, -c=<value>, ...",
  "     Specified parameter is fixed to given (population mean) value.",
  "     Option -fix0 overrides other options when possible.",
  " -WC",
  "     The last data column contains sample weights.",
  " -W1",
  "     All weights are set to 1.0 (no weighting)",
  " -WF",
  "     Less frequent samples are given more weight.",
  " -WP=<weight>, -WM1=<weight>, -WM2=<weight>, -WM3=<weight>",
  "     Put additional or less weight to parent and/or metabolite fractions.", 
  " -ND",
  "     Some fractions are known to exceed 1, not divided by 100.",
  " -MRL",
  "     Error is returned if MRL check is not passed.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Fraction datafile must contain at least two columns: sample times (min) and",
  "fractions of parent tracer. Fractions of 1-3 metabolites can be given in",
  "additional columns. Weights column can be given as specified in",
  "DFT format (2) or with option -wc.",
  "Program writes the fit start and end times, nr of points, WSS,",
  "and parameters of the fitted function to the FIT file (3).",
  " ",
  "Available functions:",
  "PF - extended power function (1,4,5)",
  "  f(x) = {D^(-1/C) + (A*(x-E))^B }^-C , when x>E,",
  "  f(x) = D                            , when x<=E ,",
  "  where 0<A<=1, B>1.5, C>0, 0<D<=1, 0<=E.",
  "  With option -d=1 this is essentially the same function as proposed in (4)",
  "  or with options -b=2 -d=1 -e=0 the same as suggested in (5).", 
  "HILL - Extended Hill type function (1,6)",
  "  f(x) = D - {(D-A)(x-E)^B}/{C+(x-E)^B} , when x>E,",
  "  f(x) = D                              , when x<=E ,",
  "  where 0<=A<=D, 1<=B, 0<C, 0<D<=1, 0<=E.",
  "  With options -d=1 -e=0 this is the traditional Hill type function (6)",
  "  f(x) = 1 - {(1-A)x^B}/(C+x^B)",
  "  where C>0<=A<=D, 1<=B, 0<C, 0<D<=1, 0<=E.",
  "GCDF - Inverted Gamma CDF (7)",
  "  f(x) = A{1-B*gcdf(C(x-E),D)}, where 0<A<=1,0<=B<=1,C>0,D>0.",
  "EP - Exponential/Power function based on Mu et al (8)",
  "  f(x) = B*(C*exp(-D*(x-A)^3/(E+(x-A)^2)) + (1-C)*exp(-F*(x-A))), when x>A,",
  "  f(x) = B                                                      , when x<=A ,",
  "MPF - Extended power function is fitted to 1-3 metabolite fractions.",
  "MHILL - Extended Hill type function is fitted to 1-3 metabolite fractions.",
  " ",
  "References:",
  "1. Fitting the fractions of parent tracer in plasma.",
  "   https://www.turkupetcentre.net/petanalysis/input_parent_fitting.html",
  "2. https://www.turkupetcentre.net/petanalysis/format_tpc_dft.html",
  "3. https://www.turkupetcentre.net/petanalysis/format_tpc_fit.html",
  "4. Meyer PT, et al. J Cereb Blood Flow Metab. 2004;24(3):323-333.",
  "   https://doi.org/10.1097/01.wcb.0000110531.48786.9d",
  "5. Watabe H, et al. J Cereb Blood Flow Metab 2000;20:899-909.",
  "   https://doi.org/10.1097/00004647-200006000-00002",
  "6. Wu S, et al. J Nucl Med 2007;48:926-931.",
  "   https://doi.org/10.2967/jnumed.106.038075",
  "7. Naganawa M, et al. J Cereb Blood Flow Metab. 2014;34:1818-1825.",
  "   https://doi.org/10.1038/jcbfm.2014.150",
  "8. Mu L et al. EJNMMI Res. 2020;10:114.",
  "   https://doi.org/10.1186/s13550-020-00700-7",
  " ",
  "See also: fit2dat, metabcor, avgfract, fitedit, fit_fexp, fith2met, tac2svg",
  " ",
  "Keywords: input, plasma, modelling, simulation, metabolite correction",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int        ai, help=0, version=0, verbose=1;
  int        fi, pi, ri, type=0, ret, m, n;
  int        noDivide=0, last_col_is_weight=0;
  int        MRL_check=0; // 0=off, 1=on
  double     fixed_p[MAX_PARAMETERS]; // values <0 mean not fixed
  double     min_delay=-1.0;   // <0 not set
  int        fixed0=0; // fix fraction at time zero to measured value at time zero, if available
  int        delay_par_index=0; // Which parameter represents delay time
  int        iscale_par_index=0; // Which parameter represents scale
  int        ffile_type=1; // Parameters are saved in 1=usual or 2=Hill#2
                           // format; fitting is still the same. 
  int        mdelay_joint=1; // 0=separated, 1=joint
  char      *cptr, dfile[FILENAME_MAX], ffile[FILENAME_MAX],
             svgfile[FILENAME_MAX];
  DFT        dft;
  FIT        fit;
  double     a, b, tstart, tstop, maxfract, aic;
  // Weights: 0=left as it is in datafile, 1=set to 1, 2=sampling frequency
  int        weighting=0; 
 
#ifdef MINGW
  // Use Unix/Linux default of two-digit exponents in MinGW on Windows
  _set_output_format(_TWO_DIGIT_EXPONENT);
#endif
 
 
  /* Initially, assume that none of parameters is fixed */
  for(pi=0; pi<MAX_PARAMETERS; pi++) fixed_p[pi]=-1;
  /* By default, no additional weight for any of fractions */
  for(m=0; m<=MAX_METAB_NR; m++) frwgt[m]=1.0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&dft); fitInit(&fit);
  dfile[0]=ffile[0]=svgfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strncasecmp(cptr, "MODEL=", 6)==0) {
      cptr+=6;
      if(strncasecmp(cptr, "PF", 1)==0) {model=MODEL_PF; continue;}
      if(strcasecmp(cptr, "HILL2")==0) {model=MODEL_HILL; ffile_type=2; continue;}
      if(strncasecmp(cptr, "HILL", 1)==0) {model=MODEL_HILL; continue;}
      if(strncasecmp(cptr, "GCDF", 1)==0) {model=MODEL_GCDF; continue;}
      if(strcasecmp(cptr, "EP")==0) {model=MODEL_EP; continue;}
      if(strncasecmp(cptr, "MPF", 3)==0) {model=MODEL_MPF; continue;}
      if(strncasecmp(cptr, "MHILL", 3)==0) {model=MODEL_MHILL; continue;}
      fprintf(stderr, "Error: invalid model selection.\n"); return(1);
    } else if(strcasecmp(cptr, "ND")==0) {
      noDivide=1; continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weighting=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weighting=2; continue;
    } else if(strcasecmp(cptr, "WC")==0) {
      last_col_is_weight=1; continue;
    } else if(strncasecmp(cptr, "WP=", 3)==0) {
      frwgt[0]=atof_dpi(cptr+3); if(frwgt[0]>=0.0) continue;
    } else if(strncasecmp(cptr, "WM1=", 4)==0) {
      frwgt[1]=atof_dpi(cptr+4); if(frwgt[1]>=0.0) continue;
    } else if(strncasecmp(cptr, "WM2=", 4)==0) {
      frwgt[2]=atof_dpi(cptr+4); if(frwgt[2]>=0.0) continue;
    } else if(strncasecmp(cptr, "WM3=", 4)==0) {
      frwgt[3]=atof_dpi(cptr+4); if(frwgt[3]>=0.0) continue;
    } else if(strncasecmp(cptr, "A=", 2)==0) {
      fixed_p[0]=atof_dpi(cptr+2); if(fixed_p[0]>=0.0) continue;
    } else if(strncasecmp(cptr, "B=", 2)==0) {
      fixed_p[1]=atof_dpi(cptr+2); if(fixed_p[1]>=0.0) continue;
    } else if(strncasecmp(cptr, "C=", 2)==0) {
      fixed_p[2]=atof_dpi(cptr+2); if(fixed_p[2]>=0.0) continue;
    } else if(strncasecmp(cptr, "D=", 2)==0) {
      fixed_p[3]=atof_dpi(cptr+2); if(fixed_p[3]>=0.0) continue;
    } else if(strncasecmp(cptr, "E=", 2)==0) {
      fixed_p[4]=atof_dpi(cptr+2); if(fixed_p[4]>=0.0) continue;
    } else if(strncasecmp(cptr, "F=", 2)==0) {
      fixed_p[5]=atof_dpi(cptr+2); if(fixed_p[5]>=0.0) continue;
    } else if(strncasecmp(cptr, "DELAYMIN=", 9)==0) {
      min_delay=atof_dpi(cptr+9); if(min_delay>=0.0) continue;
    } else if(strcasecmp(cptr, "FIX0")==0 || strcasecmp(cptr, "FIXED0")==0) {
      fixed0=1; continue;
    } else if(strcasecmp(cptr, "MRL")==0) {
      MRL_check=1; continue;
    } else if(strncasecmp(cptr, "MDELAY=", 7)==0) {
      cptr+=7;
      if(strncasecmp(cptr, "JOINT", 1)==0) {mdelay_joint=1; continue;}
      if(strncasecmp(cptr, "SEPARATED", 1)==0) {mdelay_joint=0; continue;}
      fprintf(stderr, "Error: invalid mdelay selection.\n"); return(1);
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "MIN=", 4)==0) {
      cptr+=4;
      if(strcasecmp(cptr, "OLS")==0 || strcasecmp(cptr, "SS")==0 || strcasecmp(cptr, "WSS")==0) {
        min_meas=MINMEAS_OLS; continue;}
      if(strcasecmp(cptr, "LAD")==0 || strcasecmp(cptr, "ABS")==0 || strcasecmp(cptr, "WABS")==0) {
        min_meas=MINMEAS_LAD; continue;}
      if(strcasecmp(cptr, "ODR")==0) {
        min_meas=MINMEAS_ODR; continue;}
      fprintf(stderr, "Error: invalid minimization measure.\n"); return(1);
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(dfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ffile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!dfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Apply default model, if not selected by user */
  if(model==MODEL_UNKNOWN) model=MODEL_HILL;
  if(min_meas==MINMEAS_UNKNOWN) min_meas=MINMEAS_OLS;
  /* Check compatibility of options */
  if(min_meas==MINMEAS_ODR) {
    if(model!=MODEL_HILL && model!=MODEL_PF) {
      fprintf(stderr, "Warning: ODR is not available with selected model.\n");
      fprintf(stderr, "Note: using OLS instead of ODR.\n");
      min_meas=MINMEAS_OLS;
    }
  }
  /* Set function specific parameters */
  type=func_type[model];
  if(model==MODEL_HILL && ffile_type==2) type=MF_EHILL2_PAR;
  parNr=func_par_nr[model];
  delay_par_index=parNr-1;
  if(model!=MODEL_GCDF) iscale_par_index=parNr-2;
  if(model==MODEL_EP) {delay_par_index=0; iscale_par_index=1;}
  for(pi=0, n=0; pi<parNr; pi++) if(fixed_p[pi]>=0.0) n++;
  if((parNr-n)<1) {
    fprintf(stderr, "Error: too many parameters were fixed.\n");
    return(1);
  }
  if(fixed_p[delay_par_index]>=0.0 && mdelay_joint==0)
    fprintf(stderr, "Warning: option -mdelay is not effective with option -e\n");
  if(!ffile[0]) strcpy(ffile, "stdout");
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dfile := %s\n", dfile);
    printf("ffile := %s\n", ffile);
    printf("svgfile := %s\n", svgfile);
    printf("noDivide := %d\n", noDivide);
    printf("weighting := %d\n", weighting);
    printf("mdelay_joint := %d\n", mdelay_joint);
    printf("last_col_is_weight := %d\n", last_col_is_weight);
    printf("model := %d\n", model);
    printf("type := %d\n", type);
    printf("parNr := %d\n", parNr);
    printf("min_meas := %d\n", min_meas);
    for(pi=0; pi<parNr; pi++) if(fixed_p[pi]>=0.0) printf("fixed_p[%d] := %g\n", pi, fixed_p[pi]);
    if(min_delay>=0.0) printf("min_delay := %g\n", min_delay);
    printf("fixed0 := %d\n", fixed0);
    printf("MRL_check := %d\n", MRL_check);
    for(m=0; m<=MAX_METAB_NR; m++) printf("frwgt[%d] := %g\n", m, frwgt[m]);
    fflush(stdout);
  }
  if(verbose>2) CSV_TEST=verbose-2; else CSV_TEST=0;
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", dfile);
  if(dftRead(dfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dfile, dfterrmsg);
    return(2);
  }
  origdataNr=dft.frameNr;
 
  /* If required, set the last column as weights */
  if(last_col_is_weight!=0) {
    if(verbose>1) printf("reading weights from the last column.\n");
    if(dft.voiNr<2) {
      fprintf(stderr, "Error: no column for weights in %s\n", dfile);
      dftEmpty(&dft); return(2);
    }
    for(fi=0; fi<dft.frameNr; fi++) {
      dft.w[fi]=dft.voi[dft.voiNr-1].y[fi];
      if(isnan(dft.w[fi])) dft.w[fi]=0.0;
    }
    dft.isweight=1; dft.voiNr--;
  }
 
  if(verbose>1) {
    if(strlen(dft.studynr)>0 && strcmp(dft.studynr, ".")!=0)
      printf("study_number := %s\n", dft.studynr);
    printf("tacNr := %d\n", dft.voiNr);
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
 
  /* Make sure that times are in minutes */
  if(dft.timeunit==TUNIT_SEC) {
    if(verbose>1) printf("Sample times are converted to min\n");
    dftSec2min(&dft);
  }
  if(dft.timeunit==TUNIT_UNKNOWN) {
    dft.timeunit=TUNIT_MIN;
  }
 
  /* Ignore extra columns, if parent fraction is fitted */
  if(model==MODEL_PF || model==MODEL_HILL || model==MODEL_GCDF || model==MODEL_EP) {
    if(dft.voiNr>1) {
      fprintf(stderr, "Warning: extra columns in %s are ignored.\n", dfile);
      dft.voiNr=1;
    }
  } else { /* Metabolite fraction(s) will be fitted */
    if(dft.voiNr==1) {
      if(verbose>1) printf("calculating metabolite fractions\n");
      /* If no metabolite fractions were found, then add one */
      ret=dftAddmem(&dft, 1); if(ret!=0) {
        fprintf(stderr, "Error: cannot allocate memory for metabolite.\n");
        dftEmpty(&dft); return(2);
      }
      dft.voiNr++;
      for(fi=0; fi<dft.frameNr; fi++) dft.voi[1].y[fi]=1.0-dft.voi[0].y[fi];
    } else if(dft.voiNr>(1+MAX_METAB_NR)) {
      /* If too many metabolite fractions, then combine the last ones */
      for(fi=0; fi<dft.frameNr; fi++) for(ri=1+MAX_METAB_NR; ri<dft.voiNr; ri++)
      {
        if(isnan(dft.voi[ri].y[fi])) continue;
        if(isnan(dft.voi[MAX_METAB_NR].y[fi])) {
          dft.voi[MAX_METAB_NR].y[fi]=dft.voi[ri].y[fi]; continue;}
        dft.voi[MAX_METAB_NR].y[fi]+=dft.voi[ri].y[fi];
      }
      fprintf(stderr, "Warning: extra metabolite(s) are summed in Metab3.\n");
      dft.voiNr=1+MAX_METAB_NR;
    }
    /* Set the number of metabolite fraction curves to be fitted */
    metabNr=dft.voiNr-1; if(metabNr>MAX_METAB_NR) metabNr=MAX_METAB_NR;
    if(verbose>1) printf("metabNr := %d\n", metabNr);
    /* Set only later the fitted parameter nr based on the nr of metabolites */
  }
  /* Set the names for curves */
  if(strlen(dft.voi[0].name)==0) {
    strcpy(dft.voi[0].name, "Parent");
    strcpy(dft.voi[0].voiname, dft.voi[0].name);
  }
  for(ri=1; ri<dft.voiNr; ri++) if(strlen(dft.voi[ri].name)==0) {
    strcpy(dft.voi[ri].name, "Metab1");
    strcpy(dft.voi[ri].voiname, dft.voi[ri].name);
  }
  if(verbose>11) dftPrint(&dft);
  /* Set the number of fitted samples here; sort the data first if necessary */
  // not implemented, all data fitted
  fitdataNr=origdataNr;
  if(fitdataNr<2) { // More precise testing later
    fprintf(stderr, "Error: too few samples for fitting in %s\n", dfile);
    dftEmpty(&dft); return(2);
  }
  dft.frameNr=fitdataNr;
 
  /* Get start and end times */
  ret=dftMinMax(&dft, &tstart, &tstop, NULL, &maxfract);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", dfile);
    dftEmpty(&dft); return(2);
  }
  if(verbose>2) printf("max := %g\n", maxfract);
  /* Check that x values are >=0 and that there is range in x values */
  if(tstart<0.0 || !(tstop>tstart)) {
    fprintf(stderr, "Error: invalid sample times.\n");
    dftEmpty(&dft); return(2);
  }
  /* Check that y values are <=1 */
  if(maxfract>100.0 || maxfract<0.001) {
    fprintf(stderr, "Error: invalid fraction values.\n");
    dftEmpty(&dft); return(2);
  }
  if(noDivide==0 && maxfract>1.0) {
    fprintf(stderr, "Warning: converting percentages to fractions.\n");
    for(fi=0; fi<origdataNr; fi++) for(ri=0; ri<dft.voiNr; ri++)
      if(!isnan(dft.voi[ri].y[fi])) dft.voi[ri].y[fi]/=100.0;
    /* make sure that units are not percentages any more */
    dftUnitToDFT(&dft, CUNIT_UNITLESS);
  }
  if(verbose>9) dftPrint(&dft);
 
  /* Check and set weights */
  if(dft.isweight==0 || weighting==1) {
    dft.isweight=1;
    for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;
  }
  if(weighting==2) { // weighting by sample frequency
    dftWeightByFreq(&dft); // must be sorted by time before this
  }
  /* Set sample weights to zero, if parent fraction is NaN */
  for(fi=0; fi<dft.frameNr; fi++) if(isnan(dft.voi[0].y[fi])) dft.w[fi]=0;
  if(verbose>1) {
    printf("data_weights := %g", dft.w[0]);
    for(fi=1; fi<dft.frameNr; fi++) printf(", %g", dft.w[fi]);
    printf("\n");
  }
 
 
  /*
   *  Allocate memory for fits
   */
  if(verbose>1) printf("allocating memory for fits.\n");
  ret=fit_allocate_with_dft(&fit, &dft);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate space for fit results (%d).\n", ret);
    dftEmpty(&dft); fitEmpty(&fit); return(4);
  }
  /* Set necessary information in fit data structure */
  fit.voiNr=dft.voiNr;
  strncpy(fit.datafile, dfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 256);
  fit.time=time(NULL);
  for(fi=n=0; fi<fitdataNr; fi++) if(dft.x[fi]>=0.0 && !isnan(dft.voi[0].y[fi])) n++;
  for(ri=0; ri<fit.voiNr; ri++) {
    fit.voi[ri].type=type;
    fit.voi[ri].parNr=parNr;
    if(model==MODEL_MPF || model==MODEL_MHILL) fit.voi[ri].parNr*=metabNr; 
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=n;
  }
  if(model==MODEL_MHILL) {
    fit.voi[0].type=MF_HILL3M_PAR;
    fit.voi[1].type=MF_HILL3M_M1;
    if(metabNr>1) fit.voi[2].type=MF_HILL3M_M2;
    if(metabNr>2) fit.voi[3].type=MF_HILL3M_M3;
  } else if(model==MODEL_MPF) {
    fit.voi[0].type=MF_PF3M_PAR;
    fit.voi[1].type=MF_PF3M_M1;
    if(metabNr>1) fit.voi[2].type=MF_PF3M_M2;
    if(metabNr>2) fit.voi[3].type=MF_PF3M_M3;
  }
 
  /*
   *  Set parameter constraints
   */
  switch(model) {
    case MODEL_PF:
    case MODEL_MPF:
      pmin[0]=1.0e-20;  pmax[0]=1.0e+00;
      pmin[1]=1.5;      pmax[1]=5.0e+01;
      pmin[2]=1.0e-20;  pmax[2]=1.0e+00;
      pmin[3]=1.0e-20;  pmax[3]=1.0e+00;
      break;
    case MODEL_HILL:
    case MODEL_MHILL:
      pmin[0]=0.0;      pmax[0]=1.0e+00;
      pmin[1]=1.0;      pmax[1]=1.5e+01;
      pmin[2]=1.0e-06;  pmax[2]=1.0e+05;
      pmin[3]=1.0e-06;  pmax[3]=1.0e+00;
      break;
    case MODEL_GCDF:
      pmin[0]=0.0;      pmax[0]=1.0e+00;
      pmin[1]=0.1;      pmax[1]=1.0e+00;
      pmin[2]=1.0e-06;  pmax[2]=1.0e+01;
      pmin[3]=1.0e-06;  pmax[3]=1.0e+02;
      break;
    case MODEL_EP:
      pmin[0]=0.0;      pmax[0]=1.0; // delay
      pmin[1]=0.1;      pmax[1]=1.0; // max
      pmin[2]=0.001;    pmax[2]=1.0;
      pmin[3]=0.001;    pmax[3]=1.0;
      pmin[4]=1.0;      pmax[4]=1.0e+03;
      pmin[5]=1.0e-06;  pmax[5]=0.2;
      break;
  }
  /* Delay parameter limits; first the min */
  if(min_delay>=0) {
    if(min_delay>0.9*dft.x[fitdataNr-1]) min_delay=0.9*dft.x[fitdataNr-1];
    pmin[delay_par_index]=min_delay;
  } else {
    if(fixed_p[iscale_par_index]==1.0) { // fraction starts from 1.0
      a=0.0; // search last sample that is 1.0
      for(fi=0; fi<fitdataNr; fi++) if(dft.voi[0].y[fi]>=1.0 && dft.x[fi]>a) a=dft.x[fi];
      pmin[delay_par_index]=a;
    } else { // here initial fraction can be <1
      pmin[delay_par_index]=0.0;
    }
  }
  /* then the max delay */
  if(model==MODEL_PF || model==MODEL_HILL || model==MODEL_GCDF || model==MODEL_EP || mdelay_joint==1) {
    /* parent fraction is fitted, or if metabolite fraction(s) are fitted, then
       common delay time for all metabolites is assumed */ 
    if(fixed_p[iscale_par_index]==1.0) {
      // fractions start from 1.0/0.0
      a=1000.0; // search first sample time which is lower than 1.0
      for(fi=0; fi<fitdataNr; fi++) if(dft.voi[0].y[fi]<1.0 && dft.x[fi]<a) a=dft.x[fi];
      pmax[delay_par_index]=a;
    } else { // here initial fraction can be <1
      // max delay is set to the first sample time after max value
      for(fi=0, n=fitdataNr-1, a=0.0; fi<fitdataNr; fi++)
        if(dft.voi[0].y[fi]>=a) {a=dft.voi[0].y[fi]; n=fi;}
      if(n<fitdataNr-1) n++;
      pmax[delay_par_index]=dft.x[n];
    }
  } else {
    /* Metabolite fractions are fitted with separated delay times */
    /* Find highest and lowest parent fraction sample times,
       and set max delay time to the time where fraction has fallen to half
       of the difference */
    m=0; n=fitdataNr-1; a=dft.voi[0].y[m]; b=dft.voi[0].y[n];
    for(fi=0; fi<fitdataNr; fi++) {
      if(dft.voi[0].y[fi]>=a) {m=fi; a=dft.voi[0].y[m];}
      if(dft.voi[0].y[fi]<b) {n=fi; b=dft.voi[0].y[n];}
    }
    if(verbose>2) {
      printf("max parent fraction %g at %g\n", a, dft.x[m]);
      printf("min parent fraction %g at %g\n", b, dft.x[n]);
    }
    a=0.5*(a-b);
    for(fi=m; fi<=n; fi++) if(dft.voi[0].y[fi]<=a) break;
    if(verbose>2) printf("parent fraction dropped to %g at %g\n", a, dft.x[fi]);
    pmax[delay_par_index]=dft.x[fi];
  }
  if(min_delay>=0 && pmin[delay_par_index]==pmax[delay_par_index]) {
    pmax[delay_par_index]=1.5*pmin[delay_par_index];
    if(pmax[delay_par_index]>dft.x[fitdataNr-1]) pmax[delay_par_index]=dft.x[fitdataNr-1];
  }
  if(pmax[delay_par_index]<pmin[delay_par_index]) pmax[delay_par_index]=pmin[delay_par_index];
  /* Fix to user-defined values, if available */
  for(pi=0; pi<parNr; pi++) {if(fixed_p[pi]>=0.0) pmin[pi]=pmax[pi]=fixed_p[pi];}
  if(fixed0 && fabs(tstart)<1.0E-08) {
    /* There might be more than one zero samples, thus get the mean */
    double y0=0.0; int n=0;
    for(int i=0; i<fitdataNr; i++) if(fabs(dft.x[i])<1.0E-08) {y0+=dft.voi[0].y[i]; n++;}
    if(n>1) y0/=(double)n;
    if(!(y0>0.0)) {
      fprintf(stderr, "Warning: initial fraction is not fixed to sample zero.\n");
      fflush(stderr);
    } else {
      if(verbose>1) {printf("Note: initial fraction fixed to %g\n", y0); fflush(stdout);}
      pmin[iscale_par_index]=pmax[iscale_par_index]=y0;
      if(fixed_p[iscale_par_index]>0.0) {
        fprintf(stderr, "Warning: non-effective option -d for initial fraction.\n");
        fflush(stderr);
      }
    }
  }
  /* If metabolites are fitted, then fix the parameter limits and number */
  if(model==MODEL_MPF || model==MODEL_MHILL) {
    /* Copy the same limits to the rest of parameters */
    for(m=1; m<metabNr; m++) {
      for(pi=0; pi<parNr-1; pi++) { // other parameters except delay
        pmin[m*parNr+pi]=pmin[pi]; pmax[m*parNr+pi]=pmax[pi];}
      // delay difference compared to the first
      if(mdelay_joint==0) {
        a=pmax[pi]-pmin[pi];
        pmin[m*parNr+pi]=-0.5*a; pmax[m*parNr+pi]=+0.5*a;
      } else {
        pmin[m*parNr+pi]=0.0; pmax[m*parNr+pi]=0.0;
      }
    }
    /* Nr of fitted parameters can be set only now */
    parNr*=metabNr;
    if(verbose>2) printf("final_parNr := %d\n", parNr);
  } 
  if(verbose>1) {
    printf("Constraints for the fit:\n");
    for(pi=0; pi<parNr; pi++) printf(" limit[%d]:  %g - %g\n", pi+1, pmin[pi], pmax[pi]);
    fflush(stdout);
  }
 
 
  /*
   *  Check that sample number (not including NaNs) is at least one more
   *  than the number of actually fitted parameters.
   */
  for(fi=n=0; fi<fitdataNr; fi++) {
    if(dft.x[fi]<0.0) continue;
    if(!isnan(dft.voi[0].y[fi])) n++;
  }
  for(pi=m=0; pi<parNr; pi++) if(pmin[pi]<pmax[pi]) m++;
  if(verbose>1) printf("Comparison of nr of samples and params to fit: %d / %d\n", n, m); 
  if((n-1)<m) {
    fprintf(stderr, "Error: too few samples for fitting in %s\n", dfile);
    if(verbose>0) printf("n := %d\nm := %d\n", n, m);
    dftEmpty(&dft);  fitEmpty(&fit); return(2);
  }
 
 
  /*
   *  Fit (one fit only)
   */
  if(verbose>1) {printf("preparing to fit.\n"); fflush(stdout);}
  ri=0;
  /* Set data pointers */
  x=dft.x; w=dft.w; ymeas[0]=dft.voi[ri].y; yfit[0]=dft.voi[ri].y2;
  for(m=0; m<metabNr; m++) {
    ymeas[m+1]=dft.voi[m+1].y; yfit[m+1]=dft.voi[m+1].y2;
  }
  /* Set the initial value for WSS */
  fit.voi[ri].wss=3.402823e+38;
  /* fit */
  if(verbose>1) printf("fitting\n");
  TGO_LOCAL_INSIDE=1;
  TGO_LOCAL_OPT=0; // 0=Powell-Brent, 1=Bobyqa
  TGO_SQUARED_TRANSF=1; //if(model==MODEL_GCDF) TGO_SQUARED_TRANSF=0;
  int neighNr=6, iterNr=1, sampleNr=800;
  if(model==MODEL_MPF || model==MODEL_PF) {
    TGO_LOCAL_INSIDE=1;
    TGO_LOCAL_OPT=0;
    TGO_SQUARED_TRANSF=1;
    neighNr=5; iterNr=1; sampleNr=800;
  }
  ret=tgo(pmin, pmax, func_list[model], NULL, parNr, neighNr, &fit.voi[ri].wss,
          fit.voi[ri].p, sampleNr, iterNr, verbose-8);
  if(ret) {
    fprintf(stderr, "Error %d in TGO.\n", ret);
    dftEmpty(&dft); fitEmpty(&fit);
    return(4);
  }
  if(verbose>4) {
    printf("Results from tgo():\n");
    for(pi=0; pi<parNr; pi++) printf("  parameter[%d] := %g\n", pi, fit.voi[ri].p[pi]);
    printf("WSS := %g\n", fit.voi[ri].wss);
    fflush(stdout);
  }
  fit.voi[ri].wss=lastWSS; // remove penalty from reported WSS (do not use a)
  if(verbose>5) {printf("lastWSS := %g\n", lastWSS); fflush(stdout);}
  /* Correct fitted parameters to match constraints like inside the function */
  ri=0;
  (void)modelCheckParameters(parNr, pmin, pmax, fit.voi[ri].p, fit.voi[ri].p, &a);
  if(model==MODEL_HILL) { // Must be A<=D
    if(fit.voi[ri].p[0]>fit.voi[ri].p[3]) fit.voi[ri].p[0]=fit.voi[ri].p[3];
  } else if(model==MODEL_MHILL) {
    n=parNr/metabNr;
    for(m=0; m<metabNr; m++)
      if(fit.voi[ri].p[m*n]>fit.voi[ri].p[m*n+3])
        fit.voi[ri].p[m*n]=fit.voi[ri].p[m*n+3];
  }
  if(model==MODEL_MPF || model==MODEL_MHILL) { // delay can not be negative
    n=parNr/metabNr;
    for(m=1; m<metabNr; m++)
      if( (fit.voi[ri].p[n-1]+fit.voi[ri].p[(m+1)*n-1]) < 0.0)
        fit.voi[ri].p[(m+1)*n-1] = -fit.voi[ri].p[n-1];
  }
  /* Warn user about parameters that collide with limits */
  if(verbose>5) {
    printf("Corrected parameters:\n");
    for(pi=0; pi<parNr; pi++) printf("  parameter[%d] := %g\n", pi, fit.voi[ri].p[pi]);
  }
  if(verbose>1) {
    ret=modelCheckLimits(parNr, pmin, pmax, fit.voi[ri].p);
    if(ret==0) {if(verbose>2) fprintf(stdout, "ok\n");}
    else fprintf(stderr, "warning, fit collided with the limits.\n");
  }
  if(verbose>1) {fprintf(stdout, "WSS := %g\n", fit.voi[ri].wss); fflush(stdout);}
  /* Print measured and fitted curve */
  if(verbose>3) {
    printf("                  Measured  Fitted     Weight\n");
    for(fi=0; fi<fitdataNr; fi++)
      if(!isnan(ymeas[0][fi]))
        printf("  %2d: %9.4e   %9.4e %9.4e  %7.1e\n", fi+1, x[fi], ymeas[0][fi], yfit[0][fi], w[fi]);
      else
        printf("  %2d: %9.4e   %-9.9s %9.4e  %7.1e\n", fi+1, x[fi], ".", yfit[0][fi], w[fi]);
    printf(" fitted parameters:");
    for(pi=0; pi<parNr; pi++) printf(" %g", fit.voi[ri].p[pi]);
    printf("\n");
  }
  /* Check the MRL */
  if(verbose>1) printf("Checking the MRL.\n");
  fi=mrl_between_tacs(yfit[0], ymeas[0], fitdataNr);
  if(fi>3 && fi>fitdataNr/3) {
    if(MRL_check!=0) {
      fprintf(stderr, "Error: bad fit.\n");
      fprintf(stderr, "MRL := %d / %d\n", fi, fitdataNr);
      dftEmpty(&dft); fitEmpty(&fit);
      return(7);
    }
  }
  if(fi>2 && fi>fitdataNr/4) {
    fprintf(stderr, "Warning: bad fit.\n");
    fprintf(stderr, "MRL := %d / %d\n", fi, fitdataNr);
  } else if(verbose>1) {
    printf("MRL test passed.\n");
    if(verbose>2) fprintf(stderr, "MRL := %d / %d\n", fi, fitdataNr);
  }
  /* Calculate WSS for separate curves */
  if(verbose>1) {
    double wss[1+metabNr], f=0.0;
    for(ai=0; ai<=metabNr; ai++) {
      wss[ai]=0.0;
      for(fi=0; fi<fitdataNr; fi++) {
        a=dft.voi[ai].y[fi]-dft.voi[ai].y2[fi];
        wss[ai]+=w[fi]*a*a;
      }
      f+=wss[ai]; printf("WSS[%d] := %g\n", ai+1, wss[ai]);
    }
    printf("WSS[total] := %g\n", f);
  }
  /* Calculate AIC */
  if(verbose>1) {
    ri=0;
    for(ai=0, m=0; ai<=metabNr; ai++) if(frwgt[ai]>0.0)
      for(fi=0; fi<fitdataNr; fi++) if(w[fi]>0.0) m++;
    printf("nr_of_fitted_samples := %d\n", m);
    for(pi=0, n=0; pi<parNr; pi++) if(pmax[pi]>pmin[pi]) n++;
    printf("nr_of_fitted_parameters := %d\n", n);
    aic=aicSS(fit.voi[ri].wss, m, n);
    printf("AIC := %g\n", aic);
  }
  if(verbose>2) {
    printf("Fitted parameters:\n");
    ri=0;
    for(pi=0; pi<parNr; pi++) printf(" par[%d]:  %g\n", pi+1, fit.voi[ri].p[pi]);
    if(verbose>4) for(pi=0; pi<parNr; pi++)
      printf(" par[%d]:  %g  [%g, %g]\n", pi+1, fit.voi[ri].p[pi], pmin[pi], pmax[pi]);
  }
  /* In case of metabolite fitting, convert fitted delay parameter differences to actual delays */
  if(model==MODEL_MPF || model==MODEL_MHILL) {
    n=parNr/metabNr;
    ri=0;
    for(m=1; m<metabNr; m++) {
      fit.voi[ri].p[(m+1)*n-1]+=fit.voi[ri].p[n-1];
    }
  }
  /* Copy the parameters for metabolites, when necessary */
  for(ri=1; ri<fit.voiNr; ri++) {
    for(pi=0; pi<fit.voi[0].parNr; pi++) fit.voi[ri].p[pi]=fit.voi[0].p[pi];
    fit.voi[ri].wss=fit.voi[0].wss;
  }
  /* In case of Hill function when function format #2 is required divide a with d */
  if(model==MODEL_HILL && ffile_type==2) {
    fit.voi[0].p[0]/=fit.voi[0].p[3];
  }
 
  
  /*
   *  Save fit results
   */
  if((verbose>1) && strcmp(ffile, "stdout")!=0) printf("saving results in %s\n", ffile);
  ret=fitWrite(&fit, ffile); if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, ffile, fiterrmsg);
    dftEmpty(&dft); fitEmpty(&fit); return(11);
  }
  if(strcmp(ffile, "stdout")!=0 && verbose>0)
    printf("Function parameters written in %s\n", ffile);
 
 
  /*
   *  Plot fitted curves
   */
  if(svgfile[0]) {
    if(verbose>1) printf("creating SVG plot\n");
    DFT adft;
    char tmp[FILENAME_MAX];
 
    if(verbose>1) printf("calculating fitted curve at automatically generated sample times\n");
    dftInit(&adft);
    ret=dftAutointerpolate(&dft, &adft, 1.10*dft.x[dft.frameNr-1], verbose-7);
    if(ret) {
      fprintf(stderr, "Error %d in memory allocation for fitted curves.\n", ret);
      dftEmpty(&dft); fitEmpty(&fit); dftEmpty(&adft);
      return(13);
    }
    for(fi=0; fi<adft.frameNr; fi++) adft.w[fi]=1.0;
    x=adft.x; w=adft.w;
    for(ri=0; ri<adft.voiNr; ri++) {
      yfit[ri]=adft.voi[ri].y;
      ymeas[ri]=adft.voi[ri].y2;
    }
    // use library function
    for(ri=0, ret=0; ri<adft.voiNr; ri++) {
      for(fi=0; fi<adft.frameNr; fi++) {
        ret=fitEval(&fit.voi[ri], adft.x[fi], &a); if(ret!=0) break;
        adft.voi[ri].y[fi]=a;
      }
      if(ret!=0) break;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted curve for '%s'.\n", svgfile);
      dftEmpty(&dft); dftEmpty(&adft); fitEmpty(&fit);
      return(14);
    }
 
    /* Save SVG plot of fitted and original data */
    if(verbose>1) printf("writing %s\n", svgfile);
    tpcProgramName(argv[0], 0, 0, tmp, 64); strcat(tmp, " ");
    if(strlen(adft.studynr)>1) strcat(tmp, adft.studynr);
    ret=plot_fitrange_svg(&dft, &adft, tmp, 0.0, nan(""), 0.0, 1.0, svgfile, verbose-8);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
      dftEmpty(&adft); dftEmpty(&dft); fitEmpty(&fit);
      return(30+ret);
    }
    if(verbose>0) printf("Plots written in %s\n", svgfile);
 
    dftEmpty(&adft);
  }
 
  /* Free memory */
  dftEmpty(&dft); fitEmpty(&fit);
 
  return(0);
}
/*****************************************************************************/
 
/******************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_pf(int parNr, double *p, void *fdata)
{
  int i;
  double a, b, c, d, dt, v, f, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  a=pa[0]; b=pa[1]; c=pa[2]; d=pa[3]; dt=pa[4];
 
  /* Calculate the curve values and weighted SS */
  f=pow(d, -1.0/c);
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    if(x[i]<=dt) yfit[0][i]=d;
    else yfit[0][i]=pow(f+pow(a*(x[i]-dt), b), -c);
    v=yfit[0][i]-ymeas[0][i];
    if(min_meas==MINMEAS_LAD) v=fabs(v);
    else if(min_meas==MINMEAS_OLS) v*=v;
    else if(min_meas==MINMEAS_ODR) {
      /* Calculate function derivative at x */
      double df=0.0;
      if(x[i]>dt) df=-(b*c*pow(a*(x[i]-dt),b)/(x[i]-dt)) * pow(pow(d,-1.0/c)+pow(a*(x[i]-dt),b),-c-1.0);
      if(!isfinite(df)) df=0.0;
      /* Calculate Akaho distance */
      v*=v; v/=(1.0 + df*df);
    }
    wss+=w[i]*v;
  }
  lastWSS=wss; wss*=penalty;
  if(0) printf("%g %g %g %g %g -> %g\n", a, b, c, d, dt, wss);
  return(wss);
}
 
double func_hill(int parNr, double *p, void *fdata)
{
  int i;
  double a, b, c, d, dt, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  a=pa[0]; b=pa[1]; c=pa[2]; d=pa[3]; dt=pa[4];
  /* additional check: a<=d */
  if(a>d) {
    penalty+=100.0*(a/d);
    a=d; 
  }
 
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    if(x[i]<=dt)
      yfit[0][i]= d;
    else {
      double f=pow(x[i]-dt, b);
      yfit[0][i]= d + (a-d)*f/(c+f);
    }
    v=yfit[0][i]-ymeas[0][i];
    if(min_meas==MINMEAS_LAD) v=fabs(v);
    else if(min_meas==MINMEAS_OLS) v*=v;
    else if(min_meas==MINMEAS_ODR) {
      /* Calculate function derivative at x */
      double df=0.0; if(x[i]>dt) df=-b*c*(d-a)*pow(x[i]-dt,b-1.0)/pow(c+pow(x[i]-dt,b),2);
      if(!isfinite(df)) df=0.0;
      /* Calculate Akaho distance */
      v*=v; v/=(1.0 + df*df);
    }
    wss+=w[i]*v;
  }
  lastWSS=wss; wss*=penalty; // needed because additional penalty is possible
  //printf("   objfunc(): lastWSS=%g WSS=%g\n", lastWSS, wss);
  if(0) printf("%g %g %g %g %g -> %g\n", a, b, c, d, dt, wss);
  return(wss);
}
 
double func_gcdf(int parNr, double *p, void *fdata)
{
  int i;
  double a, b, c, d, dt, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  a=pa[0]; b=pa[1]; c=pa[2]; d=pa[3]; dt=pa[4];
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    if(x[i]<=dt)
      yfit[0][i]=a;
    else {
      yfit[0][i]= a*(1.0-b*igam(d, c*(x[i]-dt)));
    }
    v=yfit[0][i]-ymeas[0][i];
    if(min_meas==MINMEAS_LAD) v=fabs(v);
    else v*=v;
    wss+=w[i]*v;
    //printf("    x=%g ymeas=%g yfit=%g v=%g wss=%g\n", x[i], ymeas[0][i], yfit[0][i], v, wss);
  }
  lastWSS=wss; wss*=penalty;
  //printf("   objfunc(): lastWSS=%g WSS=%g\n", lastWSS, wss);
  if(0) printf("%g %g %g %g %g -> %g\n", a, b, c, d, dt, wss);
  return(wss);
}
 
double func_ep(int parNr, double *p, void *fdata) // Based on Mu et al 2020
{
  int i;
  double a, b, c, d, dt, s, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  dt=pa[0]; s=pa[1]; a=pa[2]; b=pa[3]; c=pa[4]; d=pa[5];
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    double xt=x[i]-dt;
    if(xt<=0.0)
      yfit[0][i]=s;
    else {
      yfit[0][i]= a*exp(-b*xt*xt*xt/(xt*xt+c)) + (1.0-a)*exp(-d*xt);
      yfit[0][i]*=s;
    }
    v=yfit[0][i]-ymeas[0][i];
    if(min_meas==MINMEAS_LAD) v=fabs(v);
    else v*=v;
    wss+=w[i]*v;
  }
  lastWSS=wss; wss*=penalty;
  return(wss);
}
 
double func_mpf(int parNr, double *p, void *fdata)
{
  int i, m, pn;
  double a[1+MAX_METAB_NR], b[1+MAX_METAB_NR], c[1+MAX_METAB_NR],
         d[1+MAX_METAB_NR], dt[1+MAX_METAB_NR];
  double pa[1+MAX_PARAMETERS], penalty=1.0;
  double mf[1+MAX_METAB_NR], g[1+MAX_METAB_NR], v, f, wss;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  if(fdata) {}
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  pn=parNr/metabNr;
  for(m=0; m<metabNr; m++) {
    a[m]=pa[m*pn+0]; b[m]=pa[m*pn+1]; c[m]=pa[m*pn+2];
    d[m]=pa[m*pn+3]; dt[m]=pa[m*pn+4]; if(m>0) dt[m]+=dt[0];
  }
  for(; m<MAX_METAB_NR; m++) a[m]=b[m]=c[m]=d[m]=dt[m]=0.0;
  /* additional check: delay time should be positive */
  for(m=1; m<metabNr; m++) if(dt[m]<0.0) {
    penalty+=100.0*(-dt[m]); dt[m]=0.0;}
 
  /* Calculate the curve values and weighted SS */
  for(m=0; m<metabNr; m++) g[m]=pow(d[m], -1.0/c[m]);
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    for(m=0; m<metabNr && m<=MAX_METAB_NR; m++) {
      if(x[i]<=dt[m]) mf[m]=1.0 - d[m];
      else mf[m]=1.0 - pow(g[m]+pow(a[m]*(x[i]-dt[m]), b[m]), -c[m]);
    }
    for(; m<MAX_METAB_NR; m++) mf[m]=0.0;
    yfit[0][i]=1.0;
    f=1.0-mf[0]*mf[1]-mf[0]*mf[2]-mf[1]*mf[2]+2.0*mf[0]*mf[1]*mf[2];
    yfit[1][i]=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/f; yfit[0][i]-=yfit[1][i];
    if(metabNr>1) {
      yfit[2][i]=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/f; yfit[0][i]-=yfit[2][i];
    }
    if(metabNr>2) {
      yfit[3][i]=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/f; yfit[0][i]-=yfit[3][i];
    }
    for(m=0; m<=metabNr && m<=MAX_METAB_NR; m++) {
      v=yfit[m][i]-ymeas[m][i];
      if(min_meas==MINMEAS_OLS) v*=v; else v=fabs(v);
      wss+=w[i]*frwgt[m]*v;
    }
  }
  lastWSS=wss; wss*=penalty; // needed because additional penalty is possible
  //if(!isfinite(wss)) wss=nan("");
  return(wss);
}
 
double func_mhill(int parNr, double *p, void *fdata)
{
  int i, m, pn;
  double a[1+MAX_METAB_NR], b[1+MAX_METAB_NR], c[1+MAX_METAB_NR],
         d[1+MAX_METAB_NR], dt[1+MAX_METAB_NR];
  double pa[1+MAX_PARAMETERS], penalty=1.0;
  double mf[1+MAX_METAB_NR], v, f, wss;
 
  /* Check parameters against the constraints */
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
  if(fdata) {}
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  pn=parNr/metabNr;
  for(m=0; m<metabNr; m++) {
    a[m]=pa[m*pn+0]; b[m]=pa[m*pn+1]; c[m]=pa[m*pn+2];
    d[m]=pa[m*pn+3]; dt[m]=pa[m*pn+4]; if(m>0) dt[m]+=dt[0];
    /* additional check: a<=d */
    if(a[m]>d[m]) {
      penalty+=100.0*(a[m]/d[m]);
      a[m]=d[m]; 
    }
  }
  for(; m<MAX_METAB_NR; m++) a[m]=b[m]=c[m]=d[m]=dt[m]=0.0;
  /* additional check: delay time should be positive */
  for(m=1; m<metabNr; m++) if(dt[m]<0.0) {
    penalty+=100.0*(-dt[m]); dt[m]=0.0;}
 
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<fitdataNr; i++) if(w[i]>0.0) {
    for(m=0; m<metabNr && m<=MAX_METAB_NR; m++) {
      if(x[i]<=dt[m])
        mf[m]=1.0 - d[m];
      else {
        f=pow(x[i]-dt[m], b[m]);
        mf[m]=1.0 - (d[m] + (a[m]-d[m])*f/(c[m]+f));
      }
    }
    for(; m<MAX_METAB_NR; m++) mf[m]=0.0;
    yfit[0][i]=1.0;
    f=1.0-mf[0]*mf[1]-mf[0]*mf[2]-mf[1]*mf[2]+2.0*mf[0]*mf[1]*mf[2];
    yfit[1][i]=mf[0]*(1.0-mf[1]-mf[2]+mf[1]*mf[2])/f; yfit[0][i]-=yfit[1][i];
    if(metabNr>1) {
      yfit[2][i]=mf[1]*(1.0-mf[0]-mf[2]+mf[0]*mf[2])/f; yfit[0][i]-=yfit[2][i];
    }
    if(metabNr>2) {
      yfit[3][i]=mf[2]*(1.0-mf[0]-mf[1]+mf[0]*mf[1])/f; yfit[0][i]-=yfit[3][i];
    }
    for(m=0; m<=metabNr && m<=MAX_METAB_NR; m++) {
      v=yfit[m][i]-ymeas[m][i];
      if(min_meas==MINMEAS_OLS) v*=v; else v=fabs(v);
      wss+=w[i]*frwgt[m]*v;
    }
  }
  lastWSS=wss; wss*=penalty; // needed because additional penalty is possible
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/