fit_hiad.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"
/*****************************************************************************/
 
/*****************************************************************************/
double func1821(int parNr, double *p, void*);
double func1801(int parNr, double *p, void*);
/*****************************************************************************/
double (*func)(int parNr, double *p, void*);
double *x, *ymeas, *yfit, *w; /* Local */
int dataNr=0, parNr=5;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double dlowtime;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the sum of Hill and Hill derivative function to",
  "the PET time-activity curves (TACs):",
  "f(x) = A * { [n*xt^(n-1)]/[h^n+xt^n] - [n*xt^(2*n-1)]/[h^n+xt^n]^2   ",
  "             + [k*xt^n]/[h^n+xt^n] }",
  ", where xt=x-dt and n>=1.",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -Hill",
  "     Hill function without the Hill derivative function is used;",
  "     f(x) = (A*xt^n)/(h^n+xt^n).",
  " -w1 All weights are set to 1.0 (no weighting); by default, weights in",
  "     data file are used, if available.",
  " -wf",
  "     Weight by sampling interval.",
  " -k=<value>",
  "     Parameter k is constrained to given value; setting k to zero causes",
  "     the function to approach zero.",
  " -n=<value>",
  "     Parameter n (HillSlope) is constrained to the given value;",
  "     n<1 is not allowed.",
  " -delay=<<value>|mean|median>",
  "     Delay time (dt) is constrained to specified value or to mean or median",
  "     of all TACs.",
/*" -dlow=<time>",
  "     Penalize Hill derivative higher than 1% of total tissue activity",
  "     after specified time",*/
  " -MRL",
  "     Error is returned if MRL check is not passed.",
  " -res=<Filename>",
  "     Fitted parameters are also written in result file format.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified 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 2 columns, sample times (x) and",
  "concentrations (y).",
  "Weights can be specified as usual if data is in DFT format.",
  " ",
  "Program writes the fit start and end times, nr of points, WSS,",
  "and parameters of the fitted function to the fit file.",
  " ",
  "See also: fit2dat, fit_sigm, fit_bpr, fit_ppf, fit_exp, fit_ratf",
  " ",
  "Keywords: curve fitting, Hill function",
  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, m, n, ret;
  int      type=1821; // 1821 or 1801
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  double   fixed_n=nan("");
  double   fixed_k=nan("");
  double   fixed_delay=nan("");
  int      MRL_check=0; // 0=off, 1=on
  int      fix_delay=0; // 0: not fixed; 1: fixed to given value;
                        // 2: fixed to mean; 3: fixed to median.
  char     datfile[FILENAME_MAX], fitfile[FILENAME_MAX], parfile[FILENAME_MAX],
           svgfile[FILENAME_MAX], resfile[FILENAME_MAX], *cptr, temp[512];
  char     progname[256];
  DFT      dft;
  FIT      fit;
  double   tstart, tstop, miny, maxy, maxt;
  int      tgo_nr, iter_nr, neigh_nr;
  double   def_pmin[MAX_PARAMETERS], def_pmax[MAX_PARAMETERS];
 
  /* 
   *  Set defaults for parameter constraints
   */
  dlowtime=nan("");
  /* delay */ def_pmin[0]=0.0;        def_pmax[0]=0.0; // set later
  /* A     */ def_pmin[1]=1.0E-06;    def_pmax[1]=1.0E+06; // set later
  /* h     */ def_pmin[2]=1.0E-06;    def_pmax[2]=200.0;  // set later
  /* n     */ def_pmin[3]=1.0;        def_pmax[3]=10.0;
  /* k     */ def_pmin[4]=0.0;        def_pmax[4]=20.0; // set later
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=fitfile[0]=parfile[0]=svgfile[0]=resfile[0]=(char)0;
  dftInit(&dft); fitInit(&fit);
  tpcProgramName(argv[0], 1, 1, progname, 256);
  /* 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(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strcasecmp(cptr, "HILL")==0) {
      type=1801; parNr=4; continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      if(!atof_with_check(cptr+2, &fixed_n) && fixed_n>0.0) {
        def_pmin[3]=def_pmax[3]=fixed_n; continue;
      }
    } else if(strncasecmp(cptr, "k=", 2)==0 && strlen(cptr)>2) {
      if(!atof_with_check(cptr+2, &fixed_k) && fixed_k>=0.0) {
        def_pmin[4]=def_pmax[4]=fixed_k; continue;
      }
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      cptr+=6;
      if(strcasecmp(cptr, "MEAN")==0) {fix_delay=2; continue;}
      else if(strncasecmp(cptr, "AVERAGE", 2)==0) {fix_delay=2; continue;}
      if(strcasecmp(cptr, "MEDIAN")==0) {fix_delay=3; continue;}
      if(!atof_with_check(cptr, &fixed_delay)) {
        fix_delay=1; def_pmin[0]=def_pmax[0]=fixed_delay;
        continue;
      }
    } else if(strncasecmp(cptr, "DLOW=", 5)==0) {
      if(!atof_with_check(cptr+5, &dlowtime) && dlowtime>0.0) continue;
    } else if(strcasecmp(cptr, "MRL")==0) {
      MRL_check=1; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strcpy(svgfile, cptr+4); continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strcpy(fitfile, cptr+4); continue;
    } else if(strncasecmp(cptr, "RES=", 4)==0 && strlen(cptr)>4) {
      strcpy(resfile, cptr+4); 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 */
  for(; ai<argc; ai++) {
    if(!datfile[0]) {strcpy(datfile, argv[ai]); continue;}
    else if(!parfile[0]) {strcpy(parfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'\n", argv[ai]);
    return(1);
  }
  if(!parfile[0]) strcpy(parfile, "stdout");
  if(type==1821) func=func1821; else func=func1801;
  
 
  /* Is something missing? */
  if(!datfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("parfile := %s\n", parfile);
    printf("fitfile := %s\n", fitfile);
    printf("resfile := %s\n", resfile);
    printf("svgfile := %s\n", svgfile);
    printf("MRL_check := %d\n", MRL_check);
    if(!isnan(fixed_n)) printf("fixed_n := %g\n", fixed_n);
    if(!isnan(fixed_k)) printf("fixed_k := %g\n", fixed_k);
    if(!isnan(dlowtime)) printf("dlowtime := %g\n", dlowtime);
    printf("fix_delay := %d\n", fix_delay);
    if(fix_delay==1) printf("fixed_delay := %g\n", fixed_delay);
    printf("type := %d\n", type);
    printf("weights := %d\n", weights);
  }
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", datfile);
  if(dftRead(datfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", datfile, dfterrmsg);
    return(2);
  }
  dataNr=dft.frameNr;
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
  if(verbose>30) dftPrint(&dft);
 
  /* Get min and max X and Y */
  ret=dftMinMax(&dft, &tstart, &tstop, &miny, &maxy);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
  if(tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
  if(tstart<0.0) fprintf(stderr, "Warning: negative x value(s).\n");
  if(miny<0.0) fprintf(stderr, "Warning: negative y value(s).\n");
 
  /* Check and set weights */
  if(weights==0) {
    if(dft.isweight==0) {
      dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;}
  } else if(weights==1) {
    dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&dft);
  }
  for(fi=0; fi<dft.frameNr; fi++) if(isnan(dft.voi[0].y[fi])) dft.w[fi]=0;
  if(verbose>2) {
    printf("data_weights := %g", dft.w[0]);
    for(fi=1; fi<dft.frameNr; fi++) printf(", %g", dft.w[fi]);
    printf("\n");
  }
 
  /*
   *  If only one TAC to fit, ignore delay fixed to mean or median
   */
  if(dft.voiNr==1 && (fix_delay==2 || fix_delay==3)) {
    fprintf(stderr, "Note: delay setting ignored because only one TAC to fit.\n");
    fix_delay=0;
  }
 
  /*
   *  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); return(4);
  }
  /* Set necessary information in fit data structure */
  fit.voiNr=dft.voiNr;
  strncpy(fit.datafile, datfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 1024);
  fit.time=time(NULL);
  for(ri=0; ri<dft.voiNr; ri++) {
    fit.voi[ri].type=type; fit.voi[ri].parNr=parNr;
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=dataNr;
  }
 
 
 
 
  /*
   *  Estimate mean or median delay time, to be fixed in second fits
   */
  if(fix_delay==2 || fix_delay==3) {
    double maxv;
    if(verbose>0) {printf("fitting common delay time\n"); fflush(stdout);}
    def_pmin[0]=def_pmax[0]=0.0;
    for(ri=0; ri<dft.voiNr; ri++) {
      if(verbose>1) {
        printf("%d: %s\n", ri+1, dft.voi[ri].name); fflush(stdout);}
 
      /* Set data pointers */
      x=dft.x; ymeas=dft.voi[ri].y; yfit=dft.voi[ri].y2; w=dft.w;
      /* Set the initial values for parameters */
      fit.voi[ri].wss=3.402823e+38;
 
      /* Set parameter constraints */
      for(pi=0; pi<parNr; pi++) {pmin[pi]=def_pmin[pi]; pmax[pi]=def_pmax[pi];}    
      /* some are based on TAC max value and time */
      ret=dftMinMaxTAC(&dft, ri, NULL, &maxt, NULL, &maxv, NULL, NULL, NULL, NULL);
      if(ret) {
        fprintf(stderr, "Error: invalid TAC data for fitting.\n");
        dftEmpty(&dft); fitEmpty(&fit); return(5);
      }
      if(verbose>4) printf(" maxt := %g\n maxv := %g\n", maxt, maxv);
      /* delay */
      pmin[0]=-0.05*fabs(tstop); def_pmin[0]+=pmin[0];
      pmax[0]=0.95*fabs(maxt); def_pmax[0]+=pmax[0];
      /* A     */
      pmin[1]=0.01*maxv;
      pmax[1]=50.0*maxv;
      /* h     */
      pmax[2]=0.2*fabs(maxt); //2.0*maxt;
      if(verbose>4) {
        printf("  constraints :=");
        for(pi=0; pi<parNr; pi++) printf(" [%g,%g]", pmin[pi], pmax[pi]);
        printf("\n");
      }
 
      /* fit */
      TGO_LOCAL_INSIDE=0;
      TGO_SQUARED_TRANSF=1;
      neigh_nr=6; iter_nr=1; tgo_nr=1000;
      ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &fit.voi[ri].wss, 
              fit.voi[ri].p, tgo_nr, iter_nr, verbose-5);
      if(ret) {
        fprintf(stderr, "Error %d in TGO.\n", ret);
        dftEmpty(&dft); fitEmpty(&fit);
        return(5);
      }
      /* Warn user about parameters that collide with limits */
      ret=modelCheckLimits(parNr, pmin, pmax, fit.voi[ri].p);
      if(ret!=0 && verbose>0) 
        fprintf(stderr, "Warning: fit collided with the limits.\n");
      if(verbose>2) {
        printf(" fitted parameters:");
        for(pi=0; pi<parNr; pi++) printf(" %g", fit.voi[ri].p[pi]);
        printf("\n wss := %g\n", fit.voi[ri].wss);
      }
      if(verbose==1 && dft.voiNr>1) {printf("."); fflush(stdout);}
    } // next TAC
    if(verbose==1) {printf("\n");}
 
    /* calculate median */
    double dfm[dft.voiNr], delay_mean, delay_median, delay_sd, v;
    for(ri=0; ri<dft.voiNr; ri++) dfm[ri]=fit.voi[ri].p[0];
    delay_mean=dmean(dfm, dft.voiNr, &delay_sd);
    delay_median=dmedian(dfm, dft.voiNr);
    if(verbose>1) {
      printf("mean_delay_time := %6.4f +- %6.4f\n", delay_mean, delay_sd);
      printf("median_delay_time := %.4f\n", delay_median);
    }
    if(fix_delay==2) v=delay_mean; else v=delay_median;
    def_pmin[0]/=(double)dft.voiNr; def_pmax[0]/=(double)dft.voiNr;
    if(v<def_pmin[0] || v>def_pmax[0]) {
      if(verbose>2) {
        printf("def_pmin[0] := %g\n", def_pmin[0]);
        printf("def_pmax[0] := %g\n", def_pmax[0]);
      }
      fprintf(stderr,
        "Error: delay time could not be fixed to regional mean or median\n");
      dftEmpty(&dft); fitEmpty(&fit); return(5);
    }
    def_pmin[0]=def_pmax[0]=v;
    if(verbose>0) printf("fixed_delay_time := %g\n", def_pmin[0]);
  }
 
 
 
  /*
   *  Fit one TAC at a time
   */
  if(verbose>0) {printf("fitting\n"); fflush(stdout);}
  for(ri=0; ri<dft.voiNr; ri++) {
    if(verbose>1) {printf("%d: %s\n", ri+1, dft.voi[ri].name); fflush(stdout);}
    /* Set data pointers */
    x=dft.x; ymeas=dft.voi[ri].y; yfit=dft.voi[ri].y2; w=dft.w;
    /* Set the initial values for parameters */
    fit.voi[ri].wss=3.402823e+38;
 
    /* Set parameter constraints */
    for(pi=0; pi<parNr; pi++) {pmin[pi]=def_pmin[pi]; pmax[pi]=def_pmax[pi];}    
    /* some are based on TAC max value and time */
    double maxv;
    ret=dftMinMaxTAC(&dft, ri, NULL, &maxt, NULL, &maxv, NULL, NULL, NULL, NULL);
    if(ret) {
      fprintf(stderr, "Error: invalid TAC data for fitting.\n");
      dftEmpty(&dft); fitEmpty(&fit); return(5);
    }
    if(verbose>4) printf(" maxt := %g\n maxv := %g\n", maxt, maxv);
    /* delay */
    if(fix_delay==0) { // if delay time was not fixed
      pmin[0]=-0.05*fabs(tstop);
      pmax[0]=0.95*fabs(maxt);
    }
    /* A     */
    pmin[1]=0.01*maxv;
    pmax[1]=50.0*maxv;
    /* h     */
    pmax[2]=0.3*dft.x[dft.frameNr-1];
    if(verbose>2) {
      printf("  constraints :=");
      for(pi=0; pi<parNr; pi++) printf(" [%g,%g]", pmin[pi], pmax[pi]);
      printf("\n");
    }
 
    /* fit */
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=1;
    neigh_nr=6; iter_nr=1; tgo_nr=1000;
    ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &fit.voi[ri].wss, 
            fit.voi[ri].p, tgo_nr, iter_nr, verbose-5);
    if(ret) {
      fprintf(stderr, "Error %d in TGO.\n", ret);
      dftEmpty(&dft); fitEmpty(&fit);
      return(5);
    }
    /* Warn user about parameters that collide with limits */
    ret=modelCheckLimits(parNr, pmin, pmax, fit.voi[ri].p);
    if(ret!=0 && verbose>0) 
      fprintf(stderr, "Warning: fit collided with the limits.\n");
    /* Force parameter values inside limits, as they are in the simulation */
    double a;
    modelCheckParameters(parNr, pmin, pmax,
                         fit.voi[ri].p, fit.voi[ri].p, &a);
    fit.voi[ri].wss/=a; // remove penalty from reported WSS
    if(verbose>2) {
      printf(" fitted parameters:");
      for(pi=0; pi<parNr; pi++) printf(" %g", fit.voi[ri].p[pi]);
      printf("\n wss := %g\n", fit.voi[ri].wss);
    }
    /* Calculate the maximum of Hill derivative part */
    if(verbose>1 && type!=1801) {
      double a, b, c, A, h, n, s, maxx1, maxx2;
      A=fit.voi[ri].p[1]; h=fit.voi[ri].p[2]; n=fit.voi[ri].p[3];
      a=2.0*n-n*n-1.0; b=-n*(n-1.0)*pow(h, n); c=n*(n-1.0)*pow(h, 2.0*n);
      s=sqrt(b*b - 4.0*a*c);
      maxx1=pow((-b-s)/a, 1.0/n)/pow(2.0, 1.0/n);
      maxx2=pow((-b+s)/a, 1.0/n)/pow(2.0, 1.0/n);
      double hn, xn, hnxn, hilld;
      hn=pow(h, n);
      xn=pow(maxx1, n); hnxn=hn+xn;
      hilld= A*n*pow(maxx1, n-1.0)/hnxn - n*pow(maxx1, 2.0*n-1.0)/(hnxn*hnxn);
      printf("maxx1=%g maxy1=%g\n", maxx1+fit.voi[ri].p[0], hilld);
      xn=pow(maxx2, n); hnxn=hn+xn;
      hilld= A*n*pow(maxx2, n-1.0)/hnxn - n*pow(maxx2, 2.0*n-1.0)/(hnxn*hnxn);
      printf("maxx2=%g maxy2=%g\n", maxx2+fit.voi[ri].p[0], hilld);
    }
 
    /* Check the MRL */
    if(verbose>1) printf("Checking the MRL.\n");
    m=mrl_between_tacs(yfit, ymeas, dataNr);
    if(m>3 && m>dataNr/3) {
      if(MRL_check!=0) {
        fprintf(stderr, "Error: bad fit.\n");
        fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
        dftEmpty(&dft); fitEmpty(&fit);
        return(7);
      }
    }
    if(m>2 && m>dataNr/4) {
      fprintf(stderr, "Warning: bad fit.\n");
      fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
    } else if(verbose>1) {
      printf("MRL test passed.\n");
      if(verbose>2) fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
    }
 
    if(verbose==1 && dft.voiNr>1) {printf("."); fflush(stdout);}
  } // next TAC
  if(verbose==1) {printf("\n");}
 
 
 
  /*
   *  Save fit results
   */
  if(verbose>1) printf("saving results in %s\n", parfile);
  ret=fitWrite(&fit, parfile); if(ret) {
    fprintf(stderr, "Error (%d) in writing '%s': %s\n", ret, parfile, fiterrmsg);
    dftEmpty(&dft); fitEmpty(&fit); return(11);
  }
  if(verbose>0) printf("Function parameters written in %s\n", parfile);
 
 
  /* Also allocate memory for result file, if needed */
  if(resfile[0]) {
    if(verbose>1) printf("allocating memory for results.\n");
    RES res;
    resInit(&res);
    ret=fitToResult(&fit, &res, temp);
    if(ret!=0) {
      fprintf(stderr, "Error in making results: %s\n", temp);
      dftEmpty(&dft); fitEmpty(&fit); resEmpty(&res); return(4);
    }
    n=0; strcpy(res.parname[n++], "Function");
    strcpy(res.parname[n++], "Delay");
    strcpy(res.parname[n++], "A");
    strcpy(res.parname[n++], "h");
    strcpy(res.parname[n++], "n");
    if(type==1821) {
      strcpy(res.parname[n++], "K");
    }
    strcpy(res.parname[n++], "WSS");
    strcpy(res.program, progname);
    sprintf(res.datarange, "%g - %g", tstart, tstop);
    strcpy(res.studynr, dft.studynr);
    if(verbose>0) resPrint(&res);
    if(verbose>1) printf("saving results in %s\n", resfile);
    if(resWrite(&res, resfile, verbose-3)) {
      fprintf(stderr, "Error in writing '%s': %s\n", resfile, reserrmsg);
      dftEmpty(&dft); fitEmpty(&fit); resEmpty(&res); return(11);
    }
    resEmpty(&res);
    if(verbose>1) printf("Function parameters written in %s\n", resfile);
  }
 
 
 
  /*
   *  Plot fitted curves
   */
  if(svgfile[0]) {
    if(verbose>1) printf("creating SVG plot\n");
    DFT adft;
    char tmp[FILENAME_MAX];
    double a;
 
    if(verbose>1)
      printf("calculating fitted curve at automatically generated sample times\n");
    dftInit(&adft);
    ret=dftAutointerpolate(&dft, &adft, 1.05*tstop, verbose-10);
    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;
    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 TAC 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, 128); strcat(tmp, " ");
    if(strlen(adft.studynr)>1) strcat(tmp, adft.studynr);
    ret=plot_fitrange_svg(&dft, &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(&dft); fitEmpty(&fit);
      return(30+ret);
    }
    if(verbose>0) printf("Plots written in %s\n", svgfile);
 
    dftEmpty(&adft);
  }
 
  /*
   *  Save fitted TACs (original data is lost here, so do this in the end)
   */
  if(fitfile[0]) {
    if(verbose>1) printf("saving fitted curves in %s\n", fitfile);
    double a;
    for(ri=0, ret=0; ri<dft.voiNr; ri++) {
      for(fi=0; fi<dft.frameNr; fi++) {
        ret=fitEval(&fit.voi[ri], dft.x[fi], &a); if(ret!=0) break;
        dft.voi[ri].y[fi]=a;
      }
      if(ret!=0) break;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted TACs for '%s'.\n", fitfile);
      dftEmpty(&dft); fitEmpty(&fit);
      return(21);
    }
    tpcProgramName(argv[0], 1, 0, temp, 128);
    sprintf(dft.comments, "# program := %s\n", temp);
    dft.isweight=0;
    if(dftWrite(&dft, fitfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", fitfile, dfterrmsg);
      dftEmpty(&dft); fitEmpty(&fit);
      return(22);
    }
    if(verbose>0) printf("Fitted TACs written in %s\n", fitfile);
  }
 
  /* Free memory */
  dftEmpty(&dft); fitEmpty(&fit);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func1821(int parNr, double *p, void *fdata)
{
  int i;
  double t, A, h, n, k,
         xt, xn, hn, hnxn, hill, hilld, v, wss;
  //double dmaxv, dmaxt;
  double pa[MAX_PARAMETERS], penalty=1.0;
  if(fdata) {}
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  t=pa[0]; A=pa[1]; h=pa[2]; n=pa[3]; k=pa[4];
 
  /* Calculate the curve values and weighted SS */
  hn=pow(h, n); //dmaxv=dmaxt=0.0;
  for(i=0, wss=0.0; i<dataNr; i++) if(w[i]>0.0) {
    xt=x[i]-t; if(xt<=0.0) {
      yfit[i]=0.0; hilld=0.0;
    } else {
      xn=pow(xt, n); hnxn=hn+xn;
      hilld= n*pow(xt, n-1.0)/hnxn - n*pow(xt, 2.0*n-1.0)/(hnxn*hnxn);
      hill=k*xn/hnxn;
      if(hilld>=0.0 && hilld<1.0E20 && hill>=0.0 && hill<1.0E20) {
        yfit[i]=A*(hilld+hill);
      } else 
        return(nan(""));
      //if(hilld>dmaxv) {dmaxv=hilld; dmaxt=x[i];}
    }
    v=yfit[i]-ymeas[i];
    /* Penalize Hill derivative > 0 after specified time by
       removing its contribution to fitted TAC */
    //if(x[i]>=dlowtime) v=fabs(v)+fabs(A*hilld)*(1.0+x[i]-dlowtime);
    //if(x[i]>=dlowtime) v-=A*hilld;
    v*=v;
    wss+=w[i]*v;
  }
 
 
  /* Penalize Hill derivative > 1% at specified time */
  if(dlowtime>0.0) {
    xt=dlowtime-t;
    xn=pow(xt, n); hnxn=hn+xn;
    hilld= n*pow(xt, n-1.0)/hnxn - n*pow(xt, 2.0*n-1.0)/(hnxn*hnxn);
    hill=k*xn/hnxn;
    v=hilld/(hill+hilld);
    v*=10000.0; v-=1.0; if(v<0.0) v=0.0;
    v+=1.0; v=v*v*v*v; penalty*=v;
    //v=hilld/(hill+hilld) - 0.01; if(v<0.0) v=0.0;
    //v+=1.0; v*=v; penalty*=v;
  }
 
  wss*=penalty;
 
  return(wss);
}

double func1801(int parNr, double *p, void *fdata)
{
  int i;
  double t, A, h, n,
         xt, xn, hn, hnxn, hill, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
  if(fdata) {}
 
  /* Check parameters against the constraints */
  modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Get parameters */
  t=pa[0]; A=pa[1]; h=pa[2]; n=pa[3];
 
  /* Calculate the curve values and weighted SS */
  hn=pow(h, n);
  for(i=0, wss=0.0; i<dataNr; i++) if(w[i]>0.0) {
    xt=x[i]-t; if(xt<=0.0) {
      yfit[i]=0.0;
    } else {
      xn=pow(xt, n); hnxn=hn+xn;
      hill=xn/hnxn;
      if(hill>=0.0 && hill<1.0E20) { 
        yfit[i]=A*hill;
      } else 
        return(nan(""));
    }
    v=yfit[i]-ymeas[i]; v*=v;
    wss+=w[i]*v;
  }
 
  wss*=penalty;
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/