fit_feng.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the exponential based function suggested by",
  "Feng et al (1993) to the PET plasma or blood time-activity curve (TAC).",
  " ",
  "Function:",
  "  when t<=dt : f(t) = 0 ",
  "  when t>dt  : f(t) = (p1*(t-dt)-p3-p5)*exp(p2*(t-dt)) ",
  "                     + p3*exp(p4*(t-dt)) + p5*exp(p6*(t-dt))",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -lim[=<filename>]",
  "     Specify the constraints for function parameters;",
  "     This file with default values can be created by giving this",
  "     option as the only command-line argument to this program.",
  "     Without file name the default values are printed on screen.",
  " -integ[ral]",
  "     Given datafile contains the AUC of TAC as function of time.",
  " -wf",
  "     Weight by sampling interval.",
//  " -ext[ended] or -simpl[ified]", // deprecated
//  "     Apply extended model (one extra exponential term with two more",
//  "     parameters) or simplified model (p5 and p6 fixed to zero).",
  " -n=<2|3|4|A>",
  "     The nr of exponential functions; A=determined automatically.",
  "     Default is 3; changing the default may lead to a bad fit.",
  " -ss",
  "     Function approaches steady level; last exponent term is a constant,",
  "     i.e. either p6 or p8 is 0.",
  " -delay=<time>",
  "     Delay time (dt) is constrained to specified value (>=0).",
  " -MRL",
  "     Error is returned if MRL check is not passed.",
  " -fast or -safe",
  "     Speed up the fitting but increase the chance of failure, or",
  "     increase the reliability at the cost of computing time",
  " -res=<Filename>",
  "     Fitted parameters are also written in result file format.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -svg1=<Filename>",
  "     Initial part of fitted and measured TACs are plotted in SVG file",
  " -svg2=<Filename>",
  "     Lower part of fitted and measured TACs are plotted in SVG file",
  " -fit=<Filename>",
  "     Fitted TACs are written in TAC format.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC file must contain at least two columns, sample times in the first,",
  "and concentrations of TACs in the following columns.",
  "To obtain a good fit, TAC should be corrected for physical decay and any",
  "circulating metabolites.",
  "Function parameters (p1, ..., dt) will be written in the parfile.",
  " ",
  "References:",
  "1. Feng D, Huang S-C, Wang X. Models for computer simulation studies of input",
  "   functions for tracer kinetic modeling with positron emission tomography.",
  "   Int J Biomed Comput. 1993;32:95-110.",
  " ",
  "See also: fit2dat, fit_sinf, fit_exp, fit_dexp, fit_ratf, extrapol",
  " ",
  "Keywords: curve fitting, input, modelling, simulation",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
enum parameter {
  F_A1, F_A2, F_A3, F_A4, F_L1, F_L2, F_L3, F_L4, F_DT
};
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
int dataNr=0;
int parNr=7;
double *x, *ymeas, *yfit, *w, *ytmp;
double func_normal(int parNr, double *p, void*);
double func_integ(int parNr, double *p, void*);
double func_normal_direct(int parNr, double *p, void*);
double func_integ_direct(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  int      ret;
  char     tacfile[FILENAME_MAX], fitfile[FILENAME_MAX], limfile[FILENAME_MAX],
           parfile[FILENAME_MAX], resfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  char     svgfile1[FILENAME_MAX], svgfile2[FILENAME_MAX];
  int      dat_type=0; /* 0=normal; 1=integral */
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  int      iter_nr=2; // TGO iteration number 
  int      steadystate=0;
  int      MRL_check=0; // 0=off, 1=on
  int      sumn=3; // Nr of exponentials: 0=automatic, otherwise 2-4.
  double   fixed_delay=nan("");
  int      pTransform=1; // 0=no transform, 1=transforms.
 
#ifdef MINGW
  // Use Unix/Linux default of two-digit exponents in MinGW on Windows
  _set_output_format(_TWO_DIGIT_EXPONENT);
#endif
 
  /* Set parameter initial values and constraints */
  double def_pmin[MAX_PARAMETERS], def_pmax[MAX_PARAMETERS];
  def_pmin[F_A1]=0.0;       def_pmax[F_A1]=1.0E+09;
  def_pmin[F_A2]=0.0;       def_pmax[F_A2]=1.0E+09;
  def_pmin[F_A3]=0.0;       def_pmax[F_A3]=1.0E+09;
  def_pmin[F_A4]=0.0;       def_pmax[F_A4]=1.0E+09;
  def_pmin[F_L1]=-5.0;      def_pmax[F_L1]=-0.001;
  def_pmin[F_L2]=-5.0;      def_pmax[F_L2]=-0.001;
  def_pmin[F_L3]=-0.5;      def_pmax[F_L3]=0.0;
  def_pmin[F_L4]=-0.01;     def_pmax[F_L4]=0.0;
  def_pmin[F_DT]=0.0;       def_pmax[F_DT]=5.0;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=fitfile[0]=limfile[0]=parfile[0]=resfile[0]=(char)0;
  svgfile[0]=svgfile1[0]=svgfile2[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "INTEGRAL", 5)==0 || strcasecmp(cptr, "AUC")==0) {
      dat_type=1; continue;
    } else if(strcasecmp(cptr, "MRL")==0) {
      MRL_check=1; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      if(strlcpy(svgfile, cptr+4, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "SVG1=", 5)==0) {
      if(strlcpy(svgfile1, cptr+5, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "SVG2=", 5)==0) {
      if(strlcpy(svgfile2, cptr+5, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "RES=", 4)==0 && strlen(cptr)>4) {
      strlcpy(resfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "EXTENDED", 1)==0) {
      sumn=4; continue;
    } else if(strncasecmp(cptr, "SIMPLIFIED", 5)==0) {
      sumn=2; continue;
    } else if(strcasecmp(cptr, "SS")==0) {
      steadystate=1; continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      cptr+=2; if(strcasecmp(cptr, "A")==0) {sumn=0; continue;}
      sumn=atoi(cptr); if(sumn>=2 && sumn<=4) continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      if(!atof_with_check(cptr+6, &fixed_delay) && fixed_delay>=0.0) continue;
    } else if(strncasecmp(cptr, "LIM=", 4)==0 && strlen(cptr)>4) {
      strlcpy(limfile, cptr+4, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "LIM")==0) {
      strcpy(limfile, "stdout"); continue;
    } else if(strcasecmp(cptr, "SAFE")==0) {
      iter_nr=5; continue;
    } else if(strcasecmp(cptr, "NORMAL")==0) {
      iter_nr=2; continue;
    } else if(strcasecmp(cptr, "FAST")==0) {
      iter_nr=1; continue;
    }
    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(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(parfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  if(!parfile[0]) strcpy(parfile, "stdout");
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    if(limfile[0]) printf("limfile := %s\n", limfile);
    printf("tacfile := %s\n", tacfile);
    printf("parfile := %s\n", parfile);
    printf("fitfile := %s\n", fitfile);
    printf("resfile := %s\n", resfile);
    printf("svgfile := %s\n", svgfile);
    printf("svgfile1 := %s\n", svgfile1);
    printf("svgfile2 := %s\n", svgfile2);
    printf("MRL_check := %d\n", MRL_check);
    printf("weights := %d\n", weights);
    if(sumn>0) printf("sumn := %d\n", sumn);
    printf("dat_type := %d\n", dat_type);
    if(!isnan(fixed_delay)) printf("fixed_delay := %g\n", fixed_delay);
    printf("tgo_iter_nr := %d\n", iter_nr);
  }
 
 
  /* If only file name for parameter constraints was given, then write one
     with default contents, and exit */
  if(limfile[0] && !tacfile[0]) {
    /* Check that initial value file does not exist */
    if(strcasecmp(limfile, "stdout")!=0 && access(limfile, 0) != -1) {
      fprintf(stderr, "Error: parameter constraint file %s exists.\n", limfile);
      return(9);
    }
    if(verbose>1 && strcasecmp(limfile, "stdout")!=0) 
      printf("writing parameter constraints file\n");
    /* Create parameter file */
    IFT ift; iftInit(&ift);
    iftPutDouble(&ift, "A1_lower", def_pmin[F_A1], NULL, 0);
    iftPutDouble(&ift, "A1_upper", def_pmax[F_A1], NULL, 0);
    iftPutDouble(&ift, "L1_lower", def_pmin[F_L1], NULL, 0);
    iftPutDouble(&ift, "L1_upper", def_pmax[F_L1], NULL, 0);
    iftPutDouble(&ift, "A2_lower", def_pmin[F_A2], NULL, 0);
    iftPutDouble(&ift, "A2_upper", def_pmax[F_A2], NULL, 0);
    iftPutDouble(&ift, "L2_lower", def_pmin[F_L2], NULL, 0);
    iftPutDouble(&ift, "L2_upper", def_pmax[F_L2], NULL, 0);
    iftPutDouble(&ift, "A3_lower", def_pmin[F_A3], NULL, 0);
    iftPutDouble(&ift, "A3_upper", def_pmax[F_A3], NULL, 0);
    iftPutDouble(&ift, "L3_lower", def_pmin[F_L3], NULL, 0);
    iftPutDouble(&ift, "L3_upper", def_pmax[F_L3], NULL, 0);
    iftPutDouble(&ift, "A4_lower", def_pmin[F_A4], NULL, 0);
    iftPutDouble(&ift, "A4_upper", def_pmax[F_A4], NULL, 0);
    iftPutDouble(&ift, "L4_lower", def_pmin[F_L4], NULL, 0);
    iftPutDouble(&ift, "L4_upper", def_pmax[F_L4], NULL, 0);
    iftPutDouble(&ift, "DT_lower", def_pmin[F_DT], NULL, 0);
    iftPutDouble(&ift, "DT_upper", def_pmax[F_DT], NULL, 0);
    ret=iftWrite(&ift, limfile, 0);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", limfile, ift.status);
      iftEmpty(&ift); return(9);
    }
    if(strcasecmp(limfile, "stdout")!=0)
      fprintf(stdout, "Parameter file %s with initial values written.\n", limfile);
    iftEmpty(&ift); return(0);
  }
 
 
  /* Is something missing? */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /*
   *  Read model parameter upper and lower limits if file for that was given
   */
  if(limfile[0]) {
    pTransform=0;
    IFT ift; iftInit(&ift);
    if(verbose>1) printf("reading %s\n", limfile);
    ret=iftRead(&ift, limfile, 1, 0);
    if(ret) {
      fprintf(stderr, "Error in reading '%s': %s\n", limfile, ift.status);
      return(9);
    }
    if(verbose>10) iftWrite(&ift, "stdout", 0);
    int n=0;
    if(iftGetDoubleValue(&ift, 0, "A1_lower", &def_pmin[F_A1], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A1_upper", &def_pmax[F_A1], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A2_lower", &def_pmin[F_A2], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A2_upper", &def_pmax[F_A2], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A3_lower", &def_pmin[F_A3], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A3_upper", &def_pmax[F_A3], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A4_lower", &def_pmin[F_A4], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "A4_upper", &def_pmax[F_A4], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L1_lower", &def_pmin[F_L1], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L1_upper", &def_pmax[F_L1], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L2_lower", &def_pmin[F_L2], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L2_upper", &def_pmax[F_L2], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L3_lower", &def_pmin[F_L3], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L3_upper", &def_pmax[F_L3], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L4_lower", &def_pmin[F_L4], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "L4_upper", &def_pmax[F_L4], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "DT_lower", &def_pmin[F_DT], 0)>=0) n++;
    if(iftGetDoubleValue(&ift, 0, "DT_upper", &def_pmax[F_DT], 0)>=0) n++;
    iftEmpty(&ift);
    if(n==0) {fprintf(stderr, "Error: invalid parameter file.\n"); return(9);}
    /* Change the order if necessary */
    for(int pi=0; pi<parNr; pi++) if(def_pmax[pi]<def_pmin[pi]) {
      double v=def_pmax[pi]; def_pmax[pi]=def_pmin[pi]; def_pmin[pi]=v;
    }
    /* Check that not all are constrained */
    n=parFreeNr(9, def_pmin, def_pmax);
    if(n==0) {
      fprintf(stderr, "Error: no model parameters left free for fitting.\n");
      return(9);
    }
    /* What about delay time? */
    if(isnan(fixed_delay)) {
      /* not fixed with option, thus check if fixed with lim file */
      if(fabs(def_pmax[F_DT]-def_pmin[F_DT])<1.0E-10) fixed_delay=def_pmin[F_DT];
    } else {
      /* fixed with option, which overwrites lim file */
      def_pmax[F_DT]=def_pmin[F_DT]=fixed_delay;
    }
    if(verbose>1) {
      printf("constraints :=");
      for(int pi=0; pi<9; pi++) printf(" [%g,%g]", def_pmin[pi], def_pmax[pi]);
      printf("\n");
    }
  }
 
 
  /*
   *  Read data file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  DFT tac; dftInit(&tac);
  if(dftRead(tacfile, &tac)) {
    fprintf(stderr, "Error in reading '%s': %s\n", tacfile, dfterrmsg);
    return(2);
  }
  if(verbose>1) printf("checking %s\n", tacfile);
  if(verbose>1 && tac.voiNr>1) printf("tacNr := %d\n", tac.voiNr);
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&tac);
  if(verbose>50) dftPrint(&tac);
 
  /* Get min and max X and Y */
  double tstart, tstop, miny, maxy;
  ret=dftMinMax(&tac, &tstart, &tstop, &miny, &maxy);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  if(tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  if(tstart<0.0 && verbose>0) fprintf(stderr, "Warning: negative x value(s).\n");
  if(miny<0.0 && verbose>0) fprintf(stderr, "Warning: negative y value(s).\n");
  if(verbose>1) {
    printf("x_range := %g - %g\n", tstart, tstop);
    printf("sample_number := %d\n", tac.frameNr);
  }
 
  /* Check for NA's */
  if(dft_nr_of_NA(&tac) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  /* Remove negative or zero values from the end of data, if only one TAC; */
  /* this is often needed to process a TAC measured using blood pump */
  if(tac.voiNr==1) {
    int i; for(i=tac.frameNr-1; i>=0; i--) if(tac.voi[0].y[i]>0.0) break;
    tac.frameNr=i+1;
  }
  /* Check frameNr */
  if(tac.frameNr<8) {
    fprintf(stderr, "Error: %s does not contain enough valid samples.\n", tacfile);
    dftEmpty(&tac); return(2);
  }
 
  /* Set study number, if not yet available */
  if(!tac.studynr[0]) studynr_from_fname(tacfile, tac.studynr);
 
  /* Check and set weights */
  if(weights==0) {
    if(tac.isweight==0) {
      tac.isweight=1; for(int i=0; i<tac.frameNr; i++) tac.w[i]=1.0;}
  } else if(weights==1) {
    tac.isweight=1; for(int i=0; i<tac.frameNr; i++) tac.w[i]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&tac);
  }
  if(verbose>20) {
    printf("data_weights := %g", tac.w[0]);
    for(int i=1; i<tac.frameNr; i++) printf(", %g", tac.w[i]);
    printf("\n");
  }
 
  /* Get the number of samples with weight > 0 */
  int fitDataNr=dftWSampleNr(&tac);
  if(verbose>2) printf("samples_in_fit := %d\n", fitDataNr);
 
 
  /*
   *  Allocate memory for the fit parameters
   */
  if(verbose>1) printf("allocating memory for fit parameters.\n");
  FIT fit; fitInit(&fit);
  ret=fit_allocate_with_dft(&fit, &tac);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for fit parameters.\n");
    dftEmpty(&tac); return(3);
  }
  tpcProgramName(argv[0], 1, 1, fit.program, 1024);
  strlcpy(fit.datafile, tacfile, FILENAME_MAX);
  fit.time=time(NULL);
  for(int ri=0; ri<tac.voiNr; ri++) {
    if(sumn==2) fit.voi[ri].type=1312; 
    else if(sumn==4) fit.voi[ri].type=1314;
    else fit.voi[ri].type=1313;
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=fitDataNr;
  }
 
 
  /*
   *  Fit each TAC
   */
  if(verbose>0) {printf("fitting\n"); fflush(stdout);}
  double (*func)(int parNr, double *p, void*);
  if(pTransform) {
    if(dat_type==0) func=func_normal; else func=func_integ;
  } else {
    if(dat_type==0) func=func_normal_direct; else func=func_integ_direct;
  }
  TGO_LOCAL_INSIDE=0;
  if(pTransform) TGO_SQUARED_TRANSF=0; else TGO_SQUARED_TRANSF=1;
 
  /* Fit one TAC at a time */
  for(int ri=0; ri<tac.voiNr; ri++) {
    if(tac.voiNr>1 && verbose>0) {printf("."); fflush(stdout);}
 
    /* Set data pointers */
    x=tac.x; ymeas=tac.voi[ri].y; yfit=tac.voi[ri].y2; w=tac.w;
    dataNr=tac.frameNr;
 
    double aic=1.0E+100;
 
    double xmax=nan(""), ymax=nan("");
    int imax=0;
    double dt=fixed_delay;
    double lambda1=nan(""), a1=nan("");
 
    if(!limfile[0]) { // Limits not given by user
      /* Determine reasonable guess for parameters of the surge function and delay:
         L1=1/(xmax-dt) , A1=ymax*L1*exp(1)  */
      if(dat_type==0) { // TAC
        ret=dftMinMaxTAC(&tac, ri, NULL, &xmax, NULL, &ymax, NULL, NULL, NULL, &imax);
      } else { // TAC AUC: max value is the highest slope
        int n=dataNr/5; if(n<2) n=2;
        double xstart=0.0; if(!isnan(fixed_delay)) xstart=fixed_delay;
        ret=highest_slope_after(x, ymeas, dataNr, n, xstart, &ymax, NULL, NULL, &xmax);
        /* Search the sample index for max */
        double xdmin=1.0E+10; imax=0;
        for(int i=0; i<dataNr; i++) {
          double xd=fabs(x[i]-xmax); if(xd<xdmin) {imax=i; xdmin=xd;}
        }
      }
      if(ret!=0) {
        fprintf(stderr, "Error: invalid TAC data in %s\n", tacfile);
        dftEmpty(&tac); fitEmpty(&fit); return(2);
      }
      if(verbose>2) printf("  maxy=%g max_x=%g max_i=%d\n", ymax, xmax, imax);
      if(isnan(fixed_delay)) {
        /* How long time (ht) before the peak do we have to go to get y<ymax/50 ? */
        int i; for(i=imax-1; i>0; i--) if(ymeas[i]<0.5*ymax) break;
        /* assume that activity start time dt=xmax-2*ht */
        double ht=xmax-x[i];
        dt=xmax-2.0*ht; if(!(dt>0.0)) dt=0.0;
        if(verbose>2) printf("  guessed dt=%g (ht=%g)\n", dt, ht);
      }
      /* Guess for the surge function */
      lambda1=1.0/(xmax-dt);
      a1=ymax*lambda1*M_E;
      if(verbose>2) printf("  a1=%g lambda1=%g\n", a1, lambda1);
    }
 
    /* Fit with surge function and one exponential function, if requested */
    if(sumn==0 || sumn==2) {
      parNr=5; //if(isnan(fixed_delay)) parNr++;
      /* Limits for the surge and exp functions */
      if(limfile[0]) {
        pmin[0]=def_pmin[F_A1]; pmax[0]=def_pmax[F_A1];
        pmin[1]=def_pmin[F_L1]; pmax[1]=def_pmax[F_L1];
        pmin[2]=def_pmin[F_A2]; pmax[2]=def_pmax[F_A2];
        pmin[3]=def_pmin[F_L2]; pmax[3]=def_pmax[F_L2];
        pmin[4]=def_pmin[F_DT]; pmax[4]=def_pmax[F_DT];
      } else {
        pmin[0]=0.25*a1; pmax[0]=20.0*a1;
        pmin[1]=0.25*lambda1; pmax[1]=10.0*lambda1;
        pmin[2]=0.0000001*a1; pmax[2]=50.0*a1;
        pmin[3]=0.0000001*lambda1; pmax[3]=20.0*lambda1;
        /* Limits for delay time */
        if(!isnan(fixed_delay)) {
          pmin[parNr-1]=pmax[parNr-1]=fixed_delay;
        } else {
          pmin[parNr-1]=0.25*dt; pmax[parNr-1]=0.5*(dt+xmax);
        }
      }
      /* Fit */
      int tgo_nr=150*parNr;
      int neigh_nr=3*parNr;
      double wss, p[MAX_PARAMETERS];
      for(int i=0; i<MAX_PARAMETERS; i++) p[i]=0.0;
      ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &wss, p, tgo_nr, iter_nr, verbose-5);
      if(ret) {
        fprintf(stderr, "Error %d in TGO.\n", ret);
        dftEmpty(&tac); fitEmpty(&fit); return(4);
      }
      if(verbose>2) {
        printf("fitted parameters:\n");
        for(int i=0; i<parNr; i++)
          printf("\tp%d\t%g\t(%g - %g)\n", 1+i, p[i], pmin[i], pmax[i]);
        printf("\tWSS := %g\n", wss);
      }
      /* Check the MRL */
      int mrl_passed=1;
      if(MRL_check) {
        if(verbose>1) printf("\tChecking the MRL.\n");
        int m=mrl_between_tacs(yfit, ymeas, dataNr);
        if(sumn==2) {
          if(m>3 && m>dataNr/3) {
            fprintf(stderr, "Error: bad fit.\n");
            fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
            dftEmpty(&tac); fitEmpty(&fit);
            return(7);
          } else if(m>2 && m>dataNr/4) {
            fprintf(stderr, "\tWarning: bad fit.\n");
            fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
          } else if(verbose>0) {
            printf("\tMRL test passed.\n");
            if(verbose>1) fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
          }
        }
        if(m>3 && m>dataNr/3) mrl_passed=0;
      }
      if(mrl_passed) {
        /* Force parameter values inside limits, as they are in the simulation */ 
        double a;
        modelCheckParameters(parNr, pmin, pmax, p, p, &a);
        /* Calculate AIC */
        aic=aicSS(wss/a, fitDataNr, parFreeNr(parNr, pmin, pmax));
        if(verbose>2) printf("\tAIC := %g\n", aic);
        /* Copy parameter values into results data structure */
        for(int i=0; i<parNr; i++) fit.voi[ri].p[i]=p[i];
        fit.voi[ri].wss=wss/a; // remove penalty from reported WSS
        fit.voi[ri].parNr=parNr;
      }
    }
 
    /* Fit with surge function and two exponential functions, if requested */
    if(sumn==0 || sumn==3) {
      parNr=7; //if(isnan(fixed_delay)) parNr++;
      /* Limits for the surge and exp functions */
      if(limfile[0]) {
        pmin[0]=def_pmin[F_A1]; pmax[0]=def_pmax[F_A1];
        pmin[1]=def_pmin[F_L1]; pmax[1]=def_pmax[F_L1];
        pmin[2]=def_pmin[F_A2]; pmax[2]=def_pmax[F_A2];
        pmin[3]=def_pmin[F_L2]; pmax[3]=def_pmax[F_L2];
        pmin[4]=def_pmin[F_A3]; pmax[4]=def_pmax[F_A3];
        pmin[5]=def_pmin[F_L3]; pmax[5]=def_pmax[F_L3];
        pmin[6]=def_pmin[F_DT]; pmax[6]=def_pmax[F_DT];
      } else {
        pmin[0]=0.5*a1; pmax[0]=10.0*a1;
        pmin[1]=0.5*lambda1; pmax[1]=5.0*lambda1;
        pmin[2]=0.00001*a1; pmax[2]=20.0*a1;
        pmin[3]=0.01*lambda1; pmax[3]=10.0*lambda1;
        pmin[4]=0.00001*ymax; pmax[4]=0.9*ymax;
        pmin[5]=0.0; pmax[5]=0.75; if(steadystate) pmax[5]=0.0;
        /* Limits for delay time */
        if(!isnan(fixed_delay)) {
          pmin[parNr-1]=pmax[parNr-1]=fixed_delay;
        } else {
          pmin[parNr-1]=0.25*dt; pmax[parNr-1]=0.5*(dt+xmax);
        }
      }
      /* Fit */
      int tgo_nr=250*parNr;
      int neigh_nr=1+2*parNr;
      double wss, p[MAX_PARAMETERS];
      for(int i=0; i<MAX_PARAMETERS; i++) p[i]=0.0;
      ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &wss, p, tgo_nr, iter_nr, verbose-5);
      if(ret) {
        fprintf(stderr, "Error %d in TGO.\n", ret);
        dftEmpty(&tac); fitEmpty(&fit); return(4);
      }
      if(verbose>2) {
        printf("fitted parameters:\n");
        for(int i=0; i<parNr; i++)
          printf("\tp%d\t%g\t(%g - %g)\n", 1+i, p[i], pmin[i], pmax[i]);
        printf("\tWSS := %g\n", wss);
      }
      /* Check the MRL */
      int mrl_passed=1;
      if(MRL_check) {
        if(verbose>1) printf("\tChecking the MRL.\n");
        int m=mrl_between_tacs(yfit, ymeas, dataNr);
        if(sumn==3) {
          if(m>3 && m>dataNr/3) {
            fprintf(stderr, "Error: bad fit.\n");
            fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
            dftEmpty(&tac); fitEmpty(&fit);
            return(7);
          } else if(m>2 && m>dataNr/4) {
            fprintf(stderr, "\tWarning: bad fit.\n");
            fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
          } else if(verbose>0) {
            printf("\tMRL test passed.\n");
            if(verbose>1) fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
          }
        }
        if(m>3 && m>dataNr/3) mrl_passed=0;
      }
      if(mrl_passed) {
        /* Force parameter values inside limits, as they are in the simulation */ 
        double a;
        modelCheckParameters(parNr, pmin, pmax, p, p, &a);
        /* Calculate AIC */
        double laic=aicSS(wss/a, fitDataNr, parFreeNr(parNr, pmin, pmax));
        if(verbose>2) printf("\tAIC := %g\n", laic);
        if(sumn==3 || laic<aic) {
          aic=laic;
          /* Copy parameter values into results data structure */
          for(int i=0; i<parNr; i++) fit.voi[ri].p[i]=p[i];
          fit.voi[ri].wss=wss/a; // remove penalty from reported WSS
          fit.voi[ri].parNr=parNr;
        }
      }
    }
 
    /* Fit with surge function and three exponential functions, if requested */
    if(sumn==0 || sumn==4) {
      parNr=9; //if(isnan(fixed_delay)) parNr++;
      /* Limits for the surge and exp functions */
      if(limfile[0]) {
        pmin[0]=def_pmin[F_A1]; pmax[0]=def_pmax[F_A1];
        pmin[1]=def_pmin[F_L1]; pmax[1]=def_pmax[F_L1];
        pmin[2]=def_pmin[F_A2]; pmax[2]=def_pmax[F_A2];
        pmin[3]=def_pmin[F_L2]; pmax[3]=def_pmax[F_L2];
        pmin[4]=def_pmin[F_A3]; pmax[4]=def_pmax[F_A3];
        pmin[5]=def_pmin[F_L3]; pmax[5]=def_pmax[F_L3];
        pmin[6]=def_pmin[F_A4]; pmax[6]=def_pmax[F_A4];
        pmin[7]=def_pmin[F_L4]; pmax[7]=def_pmax[F_L4];
        pmin[8]=def_pmin[F_DT]; pmax[8]=def_pmax[F_DT];
      } else {
        if(sumn==4 || aic>1.0E+10) {
          pmin[0]=0.5*a1; pmax[0]=10.0*a1;
          pmin[1]=0.5*lambda1; pmax[1]=10.0*lambda1;
          pmin[2]=0.00001*a1; pmax[2]=20.0*a1;
          pmin[3]=0.01*lambda1; pmax[3]=10.0*lambda1;
          pmin[4]=0.00001*ymax; pmax[4]=0.9*ymax;
          pmin[5]=0.0001; pmax[5]=1.25;
          pmin[6]=0.0; pmax[6]=0.9*ymax;
          pmin[7]=0.0; pmax[7]=0.50; if(steadystate) pmax[7]=0.0;
          /* Limits for delay time */
          if(!isnan(fixed_delay)) {
            pmin[parNr-1]=pmax[parNr-1]=fixed_delay;
          } else {
            pmin[parNr-1]=0.25*dt; pmax[parNr-1]=0.5*(dt+xmax);
          }
        } else {
          pmin[0]=0.5*fit.voi[ri].p[0]; pmax[0]=1.2*fit.voi[ri].p[0];
          pmin[1]=0.9*fit.voi[ri].p[1]; pmax[1]=2.5*fit.voi[ri].p[1];
          pmin[2]=0.5*fit.voi[ri].p[2]; pmax[2]=1.2*fit.voi[ri].p[2];
          pmin[3]=0.9*fit.voi[ri].p[3]; pmax[3]=1.5*fit.voi[ri].p[3];
          pmin[4]=0.9*fit.voi[ri].p[4]; pmax[4]=1.2*fit.voi[ri].p[4];
          pmin[5]=0.01; pmax[5]=1.10;
          pmin[6]=0.0; pmax[6]=0.5*ymax;
          pmin[7]=0.0; pmax[7]=0.25; if(steadystate) pmax[7]=0.0;
          /* Limits for delay time */
          if(!isnan(fixed_delay)) {
            pmin[parNr-1]=pmax[parNr-1]=fixed_delay;
          } else {
            pmin[parNr-1]=0.9*fit.voi[ri].p[6]; pmax[parNr-1]=1.1*fit.voi[ri].p[6];
          }
        }
      }
      /* Fit */
      int tgo_nr=300*parNr;
      int neigh_nr=2+2*parNr;
      double wss, p[MAX_PARAMETERS];
      for(int i=0; i<MAX_PARAMETERS; i++) p[i]=0.0;
      ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &wss, p, tgo_nr, iter_nr, verbose-5);
      if(ret) {
        fprintf(stderr, "Error %d in TGO.\n", ret);
        dftEmpty(&tac); fitEmpty(&fit); return(4);
      }
      if(verbose>2) {
        printf("fitted parameters:\n");
        for(int i=0; i<parNr; i++)
          printf("\tp%d\t%g\t(%g - %g)\n", 1+i, p[i], pmin[i], pmax[i]);
        printf("\tWSS := %g\n", wss);
      }
      /* Check the MRL */
      int mrl_passed=1;
      if(MRL_check) {
        if(verbose>1) printf("\tChecking the MRL.\n");
        int m=mrl_between_tacs(yfit, ymeas, dataNr);
        if(sumn==4) {
          if(m>3 && m>dataNr/3) {
            fprintf(stderr, "Error: bad fit.\n");
            fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
            dftEmpty(&tac); fitEmpty(&fit);
            return(7);
          } else if(m>2 && m>dataNr/4) {
            fprintf(stderr, "\tWarning: bad fit.\n");
            fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
          } else if(verbose>0) {
            printf("\tMRL test passed.\n");
            if(verbose>1) fprintf(stderr, "\tMRL := %d / %d\n", m, dataNr);
          }
        }
        if(m>3 && m>dataNr/3) mrl_passed=0;
      }
      if(mrl_passed) {
        /* Force parameter values inside limits, as they are in the simulation */ 
        double a;
        modelCheckParameters(parNr, pmin, pmax, p, p, &a);
        /* Calculate AIC */
        double laic=aicSS(wss/a, fitDataNr, parFreeNr(parNr, pmin, pmax));
        if(verbose>2) printf("\tAIC := %g\n", laic);
        if(sumn==4 || laic<aic) {
          aic=laic;
          /* Copy parameter values into results data structure */
          for(int i=0; i<parNr; i++) fit.voi[ri].p[i]=p[i];
          fit.voi[ri].wss=wss/a; // remove penalty from reported WSS
          fit.voi[ri].parNr=parNr;
        }
      }
    }
 
    if(aic>1.0E+10) {
      fprintf(stderr, "Error: fitting not successful.\n");
      dftEmpty(&tac); fitEmpty(&fit); return(5);
    }
 
    if(pTransform) {
      /* Separate exponential terms to independent values */
      for(int pi=5; pi<fit.voi[ri].parNr; pi+=2) fit.voi[ri].p[pi]*=fit.voi[ri].p[pi-2];
      /* Switch the sign of exponential terms */
      for(int pi=1; pi<fit.voi[ri].parNr; pi+=2) fit.voi[ri].p[pi]*=-1.0;
    }
 
    /* Compute fitted TAC again with final sumn, if sumn was not predetermined */
    if(sumn==0) {
      if(verbose>3) printf("\tparNr := %d\n", fit.voi[ri].parNr);
      int n=fit.voi[ri].parNr-1; //if(!isnan(fixed_delay)) n--;
      n/=2; if(verbose>2) printf("\tsumn := %d\n", n);
      if(n==2) fit.voi[ri].type=1312;
      else if(n==3) fit.voi[ri].type=1313;
      else fit.voi[ri].type=1314;
      fitEvaltac(&fit.voi[ri], x, yfit, tac.frameNr);
    }
    if(verbose>3) {
      printf("Sample\tX\tY\tYfit\tWeight\n");
      for(int i=0; i<dataNr; i++)
        printf("%d\t%9.2e\t%9.2e\t%9.2e\t%7.1e\n", 1+i, x[i], ymeas[i], yfit[i], w[i]);
    }
 
  } // next TAC
  if(tac.voiNr>1 && verbose>0) {printf("\n"); fflush(stdout);}
 
 
  /*
   *  Save fit results
   */
  if(verbose>1) printf("saving results in %s\n", parfile);
  ret=fitWrite(&fit, parfile); if(ret) {
    fprintf(stderr, "Error in writing '%s': %s\n", parfile, fiterrmsg);
    dftEmpty(&tac); fitEmpty(&fit); return(11);
  }
  if(verbose>0) {printf("Function parameters written in %s\n", parfile); fflush(stdout);}
 
  /* Also allocate memory for result file, if needed */
  if(resfile[0]) {
    if(verbose>1) {printf("allocating memory for results.\n"); fflush(stdout);}
    RES res; resInit(&res);
    char buf[256];
    ret=fitToResult(&fit, &res, buf);
    if(ret!=0) {
      fprintf(stderr, "Error in making results: %s\n", buf); fflush(stderr);
      dftEmpty(&tac); fitEmpty(&fit); resEmpty(&res); return(12);
    }
    int sumn=0;
    for(int i=0; i<fit.voiNr; i++) if(fit.voi[i].parNr>sumn) sumn=fit.voi[i].parNr;
    sumn--; sumn/=2; if(verbose>3) {printf("  sumn=%d\n", sumn); fflush(stdout);}
    int n=0; 
    strcpy(res.parname[n++], "Function");
    for(int i=1; i<=sumn; i++) {
      sprintf(res.parname[n++], "A%d", i);
      sprintf(res.parname[n++], "L%d", i);
    }
    strcpy(res.parname[n++], "dT");
    strcpy(res.parname[n++], "WSS");
    strcpy(res.program, fit.program);
    sprintf(res.datarange, "%g - %g", tstart, tstop);
    strcpy(res.studynr, tac.studynr);
    /* Move dT parameter into correct place in case parameter number is variable */
    for(int ri=0; ri<fit.voiNr; ri++) if(fit.voi[ri].parNr!=res.parNr-2) {
      int imv=res.parNr-2-fit.voi[ri].parNr;
      int i2=res.parNr-2;
      int i1=i2-imv;
      res.voi[ri].parameter[i2]=res.voi[ri].parameter[i1];
      res.voi[ri].parameter[i1]=0.0;
    }
    if(verbose>1) resPrint(&res);
    if(verbose>1) {printf("saving results in %s\n", resfile); fflush(stdout);}
    if(resWrite(&res, resfile, verbose-3)) {
      fprintf(stderr, "Error in writing '%s': %s\n", resfile, reserrmsg);
      dftEmpty(&tac); fitEmpty(&fit); resEmpty(&res); return(11);
    }
    resEmpty(&res);
    if(verbose>1) {printf("Function parameters written in %s\n", resfile); fflush(stdout);}
  }
 
 
 
  /*
   *  Plotting of fitted TACs
   */
  if(svgfile[0] || svgfile1[0] || svgfile2[0]) {
 
    if(verbose>1) printf("calculating fitted curve at automatically generated sample times\n");
    DFT adft; dftInit(&adft);
    ret=dftAutointerpolate(&tac, &adft, 1.02*tac.x[tac.frameNr-1], verbose-10);
    if(ret) {
      fprintf(stderr, "Error %d in memory allocation for fitted curves.\n", ret);
      dftEmpty(&tac); fitEmpty(&fit); dftEmpty(&adft);
      return(13);
    }
    if(dat_type==1) for(int ri=0; ri<adft.voiNr; ri++) 
      ret=fitIntegralEvaltac(&fit.voi[ri], adft.x, adft.voi[ri].y, adft.frameNr);
    else for(int ri=0; ri<adft.voiNr; ri++)
      ret=fitEvaltac(&fit.voi[ri], adft.x, adft.voi[ri].y, adft.frameNr);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted curve for '%s'.\n", svgfile);
      dftEmpty(&tac); dftEmpty(&adft); fitEmpty(&fit);
      return(14);
    }
    /* Save SVG plot of fitted and original data */
    char tmp[256];
    tpcProgramName(argv[0], 0, 0, tmp, 128); strcat(tmp, " ");
    if(strlen(adft.studynr)>1) strcat(tmp, adft.studynr);
    if(svgfile[0]) {
      if(verbose>1) printf("writing %s\n", svgfile);
      ret=plot_fitrange_svg(&tac, &adft, tmp, 0.0, nan(""), 0.0, nan(""), svgfile, verbose-10);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
        dftEmpty(&adft); dftEmpty(&tac); fitEmpty(&fit);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile);
    }
 
    double tmax=tac.x[tac.frameNr/2];
    if(dat_type==0 && (svgfile1[0] || svgfile2[0])) 
      dftMinMaxTAC(&tac, 0, NULL, &tmax, NULL, NULL, NULL, NULL, NULL, NULL);
 
    /* Save SVG plot of first part of fitted and original data */
    if(svgfile1[0]) {
      if(verbose>1) printf("writing %s\n", svgfile1);
      ret=plot_fitrange_svg(&tac, &adft, tmp, 0.0, 2.*tmax, 0.0, nan(""), svgfile1, verbose-10);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile1);
        dftEmpty(&adft); dftEmpty(&tac); fitEmpty(&fit);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile1);
    }
    /* Save SVG plot of lower part of fitted and original data */
    if(svgfile2[0]) {
      if(verbose>1) printf("writing %s\n", svgfile2);
      ret=plot_fitrange_svg(&tac, &adft, tmp, 0.8*tmax+0.2*tac.x[tac.frameNr-1],
                            nan(""), 0.0, nan(""), svgfile2, verbose-10);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile2);
        dftEmpty(&adft); dftEmpty(&tac); fitEmpty(&fit);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile2);
    }
    dftEmpty(&adft);
  }
 
 
 
  /*
   *  Save fitted TAC with original sample times;
   *  this destroys the original sample values
   */
  if(fitfile[0]) {
    if(verbose>1) printf("writing fitted TAC(s)\n");
 
    if(dat_type==1) for(int ri=0; ri<tac.voiNr; ri++) 
      ret=fitIntegralEvaltac(&fit.voi[ri], tac.x, tac.voi[ri].y, tac.frameNr);
    else for(int ri=0; ri<tac.voiNr; ri++)
      ret=fitEvaltac(&fit.voi[ri], tac.x, tac.voi[ri].y, tac.frameNr);
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted curve for '%s'.\n", fitfile);
      dftEmpty(&tac); fitEmpty(&fit);
      return(17);
    }
    tac.isweight=0;
 
    ret=dftWrite(&tac, fitfile);
    if(ret) {
      fprintf(stderr, "Error: cannot write '%s'.\n", fitfile);
      dftEmpty(&tac); fitEmpty(&fit);
      return(18);
    }
  }
 
 
  dftEmpty(&tac); fitEmpty(&fit);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
 
double func_normal(int parNr, double *p, void *fdata)
{
  int i;
  double A1, A2, A3, A4, L1, L2, L3, L4;
  double dt, xt, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  double e1, e2, e3, e4;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  if(fdata) {/*prevent compiler warning*/}
  if(parNr==5 || parNr==7 || parNr==9) dt=pa[parNr-1]; else dt=0.0;
  A1=pa[0]; A2=pa[2]; 
  if(parNr>=6) A3=pa[4]; else A3=0.0; 
  if(parNr>=8) A4=pa[6]; else A4=0.0;
  L1=pa[1]; L2=pa[3]; 
  if(parNr>=6) L3=L2*pa[5]; else L3=0.0;
  if(parNr>=8) L4=L3*pa[7]; else L4=0.0;
  L1*=-1.0; L2*=-1.0; L3*=-1.0; L4*=-1.0;
 
  /* Check that D(0)>0 */
  double d= A1 + A2*L2 + A3*L3 + A4*L4 - L1*(A2+A3+A4);
  if(d<=0.0) penalty*=100.0; // penalty*=(1.0-v)*1.0E+03;
 
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-dt;
    if(xt<=0.0)
      yfit[i]=0.0;
    else {
      e1=exp(L1*xt);
      e2=exp(L2*xt);
      e3=exp(L3*xt);
      e4=exp(L4*xt);
      // evaluate the function value at xt
      yfit[i]=(A1*xt -A2 - A3 - A4)*e1 + A2*e2 + A3*e3 + A4*e4;
    }
    double v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  wss*=penalty;
  return(wss);
}
 
double func_integ(int parNr, double *p, void *fdata)
{
  int i;
  double A1, A2, A3, A4, L1, L2, L3, L4;
  double dt, xt, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  double e1, e2, e3, e4;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  if(fdata) {/*prevent compiler warning*/}
  if(parNr==5 || parNr==7 || parNr==9) dt=pa[parNr-1]; else dt=0.0;
  A1=pa[0]; A2=pa[2]; 
  if(parNr>=6) A3=pa[4]; else A3=0.0; 
  if(parNr>=8) A4=pa[6]; else A4=0.0;
  L1=pa[1]; L2=pa[3]; 
  if(parNr>=6) L3=L2*pa[5]; else L3=0.0;
  if(parNr>=8) L4=L3*pa[7]; else L4=0.0;
  L1*=-1.0; L2*=-1.0; L3*=-1.0; L4*=-1.0;
 
  /* Check that D(0)>0 */
  double d= A1 + A2*L2 + A3*L3 + A4*L4 - L1*(A2+A3+A4);
  if(d<=0.0) penalty*=100.0; //penalty*=(1.0-v)*1.0E+03;
 
  /* Calculate the integral curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-dt;
    yfit[i]=0.0;
    if(xt>0.0) {
      // evaluate the function value at xt
      e1=exp(L1*xt);
      if(L1!=0.0) {
        yfit[i]+=(A2/L1)*(1.0-e1) + (A3/L1)*(1.0-e1) + (A4/L1)*(1.0-e1);
        yfit[i]+=(A1/(L1*L1))*(1.0 + e1*(L1*xt-1.0)); 
      }
      e2=exp(L2*xt);
      if(L2!=0.0)
        yfit[i]+=(A2/L2)*(e2-1.0); 
      e3=exp(L3*xt);
      if(L3!=0.0)
        yfit[i]+=(A3/L3)*(e3-1.0); 
      e4=exp(L4*xt);
      if(L4!=0.0)
        yfit[i]+=(A4/L4)*(e4-1.0);
    }
    double v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  wss*=penalty;
  return(wss);
}
 
/* Functions with parameters that are used directly without any transformations */
 
double func_normal_direct(int parNr, double *p, void *fdata)
{
  int i;
  double A1, A2, A3, A4, L1, L2, L3, L4;
  double dt, xt, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  double e1, e2, e3, e4;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  if(fdata) {/*prevent compiler warning*/}
  if(parNr==5 || parNr==7 || parNr==9) dt=pa[parNr-1]; else dt=0.0;
  A1=pa[0]; A2=pa[2]; 
  if(parNr>=6) A3=pa[4]; else A3=0.0; 
  if(parNr>=8) A4=pa[6]; else A4=0.0;
  L1=pa[1]; L2=pa[3]; 
  if(parNr>=6) L3=pa[5]; else L3=0.0;
  if(parNr>=8) L4=pa[7]; else L4=0.0;
 
  /* Check that D(0)>0 */
  double d= A1 + A2*L2 + A3*L3 + A4*L4 - L1*(A2+A3+A4);
  if(d<=0.0) penalty*=100.0; // penalty*=(1.0-v)*1.0E+03;
 
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-dt;
    if(xt<=0.0)
      yfit[i]=0.0;
    else {
      e1=exp(L1*xt);
      e2=exp(L2*xt);
      e3=exp(L3*xt);
      e4=exp(L4*xt);
      // evaluate the function value at xt
      yfit[i]=(A1*xt -A2 - A3 - A4)*e1 + A2*e2 + A3*e3 + A4*e4;
    }
    double v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  wss*=penalty;
  return(wss);
}
 
double func_integ_direct(int parNr, double *p, void *fdata)
{
  int i;
  double A1, A2, A3, A4, L1, L2, L3, L4;
  double dt, xt, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  double e1, e2, e3, e4;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  if(fdata) {/*prevent compiler warning*/}
  if(parNr==5 || parNr==7 || parNr==9) dt=pa[parNr-1]; else dt=0.0;
  A1=pa[0]; A2=pa[2]; 
  if(parNr>=6) A3=pa[4]; else A3=0.0; 
  if(parNr>=8) A4=pa[6]; else A4=0.0;
  L1=pa[1]; L2=pa[3]; 
  if(parNr>=6) L3=pa[5]; else L3=0.0;
  if(parNr>=8) L4=pa[7]; else L4=0.0;
 
  /* Check that D(0)>0 */
  double d= A1 + A2*L2 + A3*L3 + A4*L4 - L1*(A2+A3+A4);
  if(d<=0.0) penalty*=100.0; //penalty*=(1.0-v)*1.0E+03;
 
  /* Calculate the integral curve values and weighted SS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-dt;
    yfit[i]=0.0;
    if(xt>0.0) {
      // evaluate the function value at xt
      e1=exp(L1*xt);
      if(L1!=0.0) {
        yfit[i]+=(A2/L1)*(1.0-e1) + (A3/L1)*(1.0-e1) + (A4/L1)*(1.0-e1);
        yfit[i]+=(A1/(L1*L1))*(1.0 + e1*(L1*xt-1.0)); 
      }
      e2=exp(L2*xt);
      if(L2!=0.0)
        yfit[i]+=(A2/L2)*(e2-1.0); 
      e3=exp(L3*xt);
      if(L3!=0.0)
        yfit[i]+=(A3/L3)*(e3-1.0); 
      e4=exp(L4*xt);
      if(L4!=0.0)
        yfit[i]+=(A4/L4)*(e4-1.0);
    }
    double v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  wss*=penalty;
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/