fit_gvar.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 func_gvar(int parNr, double *p, void*);
/*****************************************************************************/
double *x, *ymeas, *yfit, *w;
int dataNr=0, parNr=6;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of Gamma variate bolus plus recirculation function",
  "with time delay term to plasma and tissue time-activity curves (TACs).",
  " ",
  "Function:",
  "         | 0                                         , t < p1   ",
  "  f(t) = |                                                      ",
  "         | p2*((t-p1)^p3)*exp(-(t-p1)/p4)            , t >= p1  ",
  "         | + p5*(1-exp(-(t-p1)/p4))*exp(-(t-p1)/p6)             ",
  " ",
  "Usage: @P [Options] tacfile [fitfile]",
  " ",
  "Options:",
  " -recirc=<Y|n>",
  "     Fit recirculation terms p5 and p6 (y, default), or fit only",
  "     gamma variate function setting p5=p6=0 (n).",
  " -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.",
  " -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.",
  " -stoptime=<Fit end time>",
  "     All data with sample time > stoptime is excluded from the fit.",
  " -thr=<Threshold%>",
  "     If sample value after the peak is less than Threshold% of the peak,",
  "     then that sample and all samples after that are omitted from the fit.",
  "     Hint: try to set threshold close to 100, e.g. 99%, to get a good",
  "     estimate of tracer appearance time.",
  " -auto",
  "     Automatically sets the fit range from 0 to peak + at least one",
  "     sample more, until WSS/(sum of weights) gets worse.",
  " -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 and",
  "concentrations.",
  "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 (p1, p2, ...) of the fitted function to the fit file.",
  " ",
  "See also: fit2dat, fit_sinf, fit_exp, fit_ratf, fit_hiad",
  " ",
  "Keywords: curve fitting, TAC, 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;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  double   def_pmin[MAX_PARAMETERS], def_pmax[MAX_PARAMETERS];
  char     datfile[FILENAME_MAX], fitfile[FILENAME_MAX], parfile[FILENAME_MAX],
           svgfile[FILENAME_MAX], limfile[FILENAME_MAX];
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  int      autoend=0; // 0=no, 1=yes
  int      recirc=1;  // recirculation terms: 0=no, 1=yes
  double   threshold=0.0, stoptime=-1.0;
  char    *cptr;
  char     progname[256];
  int      ret;
 
 
  /* 
   *  Set defaults for parameter constraints; these will probably
   *  be refined before fitting.
   */
  /* delay */ def_pmin[0]=0.0;        def_pmax[0]=100.0;
  /* p2    */ def_pmin[1]=1.0E-06;    def_pmax[1]=1.0E+06;
  /* p3    */ def_pmin[2]=1.0E-03;    def_pmax[2]=2.0;
  /* p4    */ def_pmin[3]=1.0E-01;    def_pmax[3]=1.0E+03;
  /* p5    */ def_pmin[4]=1.0E-06;    def_pmax[4]=1.0E+06;
  /* p6    */ def_pmin[5]=1.0;        def_pmax[5]=1.0E+06;
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=fitfile[0]=parfile[0]=svgfile[0]=limfile[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(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "LIM=", 4)==0) {
      strlcpy(limfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "STOPTIME=", 9)==0) {
      if(atof_with_check(cptr+9, &stoptime)==0 && stoptime>0.0) continue;
    } else if(strncasecmp(cptr, "AUTO", 4)==0) {
      autoend=1; continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      if(atof_with_check(cptr+4, &threshold)==0) {
        if(threshold<0.0) threshold=0.0;
        else threshold/=100.0; // percentage to fraction
        if(threshold<1.0) continue;
      }
    } else if(strncasecmp(cptr, "RECIRC=", 7)==0) {
      if(strncasecmp(cptr+7, "YES", 1)==0) {recirc=1; continue;}
      else if(strncasecmp(cptr+7, "NO", 1)==0) {recirc=0; continue;}
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break;
  
  /* Change parameter number if NO recirculation terms */
  if(recirc==0) {
    parNr-=2;
    def_pmax[4]=def_pmin[4]=def_pmax[5]=def_pmin[5]=0.0;
  }
 
  /* 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(datfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {strlcpy(parfile, argv[ai++], FILENAME_MAX);}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'\n", argv[ai]);
    return(1);
  }
  if(!parfile[0]) strcpy(parfile, "stdout");
 
  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 only filename for parameter constraints was given, then write one
     with default contents, and exit */
  if(limfile[0] && !datfile[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);
    char buf[32];
    for(int i=0; i<parNr; i++) {
      sprintf(buf, "p%d_lower", 1+i);
      iftPutDouble(&ift, buf, def_pmin[i], NULL, 0);
      sprintf(buf, "p%d_upper", 1+i);
      iftPutDouble(&ift, buf, def_pmax[i], 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(!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) {
    if(limfile[0]) printf("limfile := %s\n", limfile);
    printf("datfile := %s\n", datfile);
    printf("parfile := %s\n", parfile);
    if(fitfile[0]) printf("fitfile := %s\n", fitfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    //printf("recirc := %d\n", recirc); // only later
    printf("weights := %d\n", weights);
    printf("threshold := %g\n", threshold);
    printf("autoend := %d\n", autoend);
    printf("stoptime := %g\n", stoptime);
    fflush(stdout);
  }
 
 
  /*
   *  Read model parameter upper and lower limits
   *  if file for that was given
   */
  if(limfile[0]) {
    if(verbose>1) printf("reading %s\n", limfile);
    IFT ift; iftInit(&ift);
    ret=iftRead(&ift, limfile, 1, 0);
    if(ret) {
      fprintf(stderr, "Error in reading '%s': %s\n", limfile, ift.status);
      iftEmpty(&ift); return(9);
    }
    if(verbose>10) iftWrite(&ift, "stdout", 0);
    double v;
    char buf[32];
    for(int i=0; i<6; i++) {
      sprintf(buf, "p%d_lower", 1+i);
      if(iftGetDoubleValue(&ift, 0, buf, &v, 0)>=0) def_pmin[i]=v;
      sprintf(buf, "p%d_upper", 1+i);
      if(iftGetDoubleValue(&ift, 0, buf, &v, 0)>=0) def_pmax[i]=v;
    }
    iftEmpty(&ift);
    /* Check that these are ok */
    int i, n=0;
    for(i=0, ret=0; i<parNr; i++) {
      if(def_pmin[i]<0.0) ret++;
      if(def_pmax[i]<def_pmin[i]) ret++;
      if(def_pmax[i]>def_pmin[i]) n++;
    }
    if(ret) {
      fprintf(stderr, "Error: invalid parameter constraints.\n");
      return(9);
    }
    if(n==0) {
      fprintf(stderr, "Error: no model parameters left free for fitting.\n");
      return(9);
    }
    /* Check if recirculation is turned off with limits */
    if(parNr==6) {
      if(def_pmax[4]<1.0E-10 && def_pmin[4]<1.0E-10) {
        def_pmax[4]=def_pmin[4]=def_pmax[5]=def_pmin[5]=0.0;
        recirc=0; parNr-=2;
      }
    }
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("recirc := %d\n", recirc);
    printf("parNr := %d\n", parNr);
    fflush(stdout);
  }
  if(verbose>2) {
    for(int i=0; i<parNr; i++) {
      printf("def_pmin[%d] := %g\n", i, def_pmin[i]);
      printf("def_pmax[%d] := %g\n", i, def_pmax[i]);
    }
    fflush(stdout);
  }
 
 
  /*
   *  Read data
   */
  if(verbose>1) {printf("reading %s\n", datfile); fflush(stdout);}
  DFT dft; dftInit(&dft);
  if(dftRead(datfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", datfile, dfterrmsg);
    fflush(stderr);
    return(2);
  }
  //int origDataNr=dft.frameNr;
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
  if(verbose>100) dftPrint(&dft);
  /* Apply specified fit stoptime if necessary */
  if(stoptime>0.0) {
    int fi;
    for(fi=1; fi<dft.frameNr; fi++) if(dft.x[fi]>stoptime) break;
    dft.frameNr=fi;
  }
  if((dft.frameNr<4 && recirc==0) || (dft.frameNr<6 && recirc==1)) {
    fprintf(stderr, "Error: not enough samples for decent fitting.\n");
    dftEmpty(&dft); fflush(stderr); return(3);
  }
  /* Check for NA's (data after stoptime is now 'removed') */
  if(dft_nr_of_NA(&dft) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", datfile);
    dftEmpty(&dft); fflush(stderr); return(3);
  }
 
  /* Get min and max X and Y */
  double tstart, tstop, miny, maxy;
  ret=dftMinMax(&dft, &tstart, &tstop, &miny, &maxy);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); return(3);
  }
  if(tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    dftEmpty(&dft); fflush(stderr); return(3);
  }
  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(int i=0; i<dft.frameNr; i++) dft.w[i]=1.0;}
  } else if(weights==1) {
    dft.isweight=1; for(int i=0; i<dft.frameNr; i++) dft.w[i]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&dft);
  }
 
 
 
  /*
   *  Allocate memory for fits
   */
  if(verbose>1) printf("allocating memory for fits.\n");
  FIT fit; fitInit(&fit);
  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;
  strlcpy(fit.datafile, datfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 256);
  strcpy(fit.unit, dft.unit);
  fit.time=time(NULL);
  for(int i=0; i<dft.voiNr; i++) {
    fit.voi[i].type=MF_GAMMAVR; // 1403
    fit.voi[i].parNr=6;
    fit.voi[i].start=tstart; fit.voi[i].end=tstop;
    fit.voi[i].dataNr=dft.frameNr;
  }
 
 
 
  /*
   *  Fit one TAC at a time
   */
  if(verbose>0) {printf("fitting\n"); fflush(stdout);}
  for(int 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;
    dataNr=dft.frameNr;
    /* Set the initial values for parameters */
    fit.voi[ri].wss=3.402823e+38;
 
    /* Get the peak value, time, and sample index */
    double peakx, peaky;
    int peaki;
    ret=dftMinMaxTAC(&dft, ri, NULL, &peakx, NULL, &peaky, NULL, NULL, NULL, &peaki);
    if(ret) {
      fprintf(stderr, "Error: invalid TAC data for fitting.\n");
      dftEmpty(&dft); fitEmpty(&fit); return(5);
    }
    if(verbose>4) printf(" peakx := %g\n peaky := %g\n", peakx, peaky);
 
    /* If required, apply threshold to remove points after the peak */
    if(threshold>0.0) {
      /* 'Remove' samples starting from a post-peak value that is smaller than threshold */
      double a=threshold*peaky;
      int i=peaki+1;
      for(; i<dft.frameNr; i++) if((ymeas[i])<a) {dataNr=i; break;}
      if(dft.timetype==DFT_TIME_STARTEND) fit.voi[ri].end=dft.x2[dataNr-1];
      else fit.voi[ri].end=dft.x[dataNr-1];
      fit.voi[ri].dataNr=dataNr;
    }
 
 
    /* Set parameter constraints */
    for(int pi=0; pi<6; pi++) {pmin[pi]=def_pmin[pi]; pmax[pi]=def_pmax[pi];}
    /* If user did not specifically set limits, then refine the limits based on TAC */
    if(!limfile[0]) {
      if(verbose>3) printf("  refining limits\n");
      /* delay */
      pmin[0]=-0.05*fabs(peakx);
      pmax[0]=0.95*fabs(peakx);
      /* gaussian variate scale factor */
      pmax[1]=10.0*peaky; if(pmin[1]>0.1*pmax[1]) pmin[1]=pmax[1]*1.0E-06;
      /* recirculation scale factor */
      if(parNr>4) {
        pmax[4]=10.0*peaky; if(pmin[4]>0.1*pmax[4]) pmin[4]=pmax[4]*1.0E-06;
      }
    }
    if(verbose>2) {
      printf("  constraints\n");
      for(int pi=0; pi<parNr; pi++) {
        printf("  pmin[%d] := %g\n", pi, pmin[pi]);
        printf("  pmax[%d] := %g\n", pi, pmax[pi]);
      }
    }
 
 
 
    /* fit */
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=1;
    int neighNr=15;
    int samNr=1000;
    int tgoNr=3;
    if(autoend==0) {
      ret=tgo(pmin, pmax, func_gvar, NULL, parNr, neighNr, &fit.voi[ri].wss, 
              fit.voi[ri].p, samNr, tgoNr, verbose-8);
      if(ret) {
        fprintf(stderr, "Error: error %d in TGO.\n", ret);
        dftEmpty(&dft); fitEmpty(&fit);
        return(6);
      }
    } else {
      double wsum=0.0;
      int saveDataNr=dataNr;
      /* Try to fit using different dataNr */
      dataNr=peaki+2; // peak sample plus one
      double lastRelWSS=3.402823e+38;
      while(dataNr<=saveDataNr) {
        ret=tgo(pmin, pmax, func_gvar, NULL, parNr, neighNr, &fit.voi[ri].wss, 
                fit.voi[ri].p, samNr, tgoNr, verbose-10);
        if(ret) {
          fprintf(stderr, "Error: error %d in TGO.\n", ret);
          dftEmpty(&dft); fitEmpty(&fit); return(6);
        }
        /* Compare to previous fit */
        wsum=doubleSum(w, (unsigned int)dataNr);
        if((fit.voi[ri].wss/wsum)>lastRelWSS) { // worse, stop
          /* Previous fit was better, but since it was not saved, fit it again */
          dataNr--;
          ret=tgo(pmin, pmax, func_gvar, NULL, parNr, neighNr, &fit.voi[ri].wss, 
                  fit.voi[ri].p, samNr, tgoNr, verbose-10);
          if(ret) {
            fprintf(stderr, "Error: error %d in TGO.\n", ret);
            dftEmpty(&dft); fitEmpty(&fit); return(6);
          } 
          break;
        }
        lastRelWSS=fit.voi[ri].wss/wsum;
        /* Try with one more sample */
        dataNr++;
      }
 
      if(dft.timetype==DFT_TIME_STARTEND) fit.voi[ri].end=dft.x2[dataNr-1];
      else fit.voi[ri].end=dft.x[dataNr-1];
 
    }
 
 
    /* Print measured and fitted curve */
    if(verbose>5) {
      printf("\n  Frame   Time        Value       Fitted     Weight \n");
      for(int i=0; i<dft.frameNr; i++)
        printf("  %02d   %8.3f  %12.4f  %12.4f  %8.4f\n",
          i+1, x[i], ymeas[i], yfit[i], w[i]);
    }
 
    /* 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
 
    /* Set final sample nr into results, not counting samples with zero weight */
    int n=0;
    for(int i=0; i<dataNr; i++) if(w[i]>1.0E-10) n++;
    fit.voi[ri].dataNr=n;
 
  } /* next TAC */
  if(verbose>0) fprintf(stdout, "\n");
 
 
  /*
   *  Print fit results on screen
   */
  if(verbose>0 && fit.voiNr<=100) fitWrite(&fit, "stdout");
 
 
  /*
   *  Save fit results, if file name was given by user
   */
  if(strcasecmp(parfile, "stdout")!=0) {
    if(verbose>1) printf("saving results in %s\n", parfile);
    if(fitWrite(&fit, parfile)) {
      fprintf(stderr, "Error in writing '%s': %s\n", parfile, fiterrmsg);
      dftEmpty(&dft); fitEmpty(&fit); return(11);
    }
  }
 
 
  /*
   *  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: cannot allocate memory for fitted curves.\n");
      dftEmpty(&dft); fitEmpty(&fit); dftEmpty(&adft);
      return(13);
    }
    //for(fi=0; fi<adft.frameNr; fi++) adft.w[fi]=1.0;
    for(int ri=0, ret=0; ri<adft.voiNr; ri++) {
      for(int 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, 32);
    if(strlen(adft.studynr)>1 && strlen(adft.studynr)<28) {
      strlcat(tmp, " ", 64);
      strlcat(tmp, adft.studynr, 64);
    }
    ret=plot_fitrange_svg(&dft, &adft, tmp, 0.0, nan(""), 0.0, nan(""),
                          svgfile, verbose-10);
    if(ret) {
      fprintf(stderr, "Error: cannot write '%s'.\n", 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(int ri=0, ret=0; ri<dft.voiNr; ri++) {
      for(int 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], 0, 0, progname, 32);
    snprintf(dft.comments, 128, "# program := %s\n", progname);
    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 func_gvar(int parNr, double *p, void *fdata)
{
  int i;
  double xt, et, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
 
  /* Calculate the curve and WSS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-pa[0];
    yfit[i]=0.0;
    if(xt>0.0) {
      et=exp(-xt/pa[3]);
      if(pa[1]>0.0) yfit[i]+=pa[1]*pow(xt, pa[2])*et;
      if(parNr>5 && pa[4]>0.0) yfit[i]+=pa[4]*(1.0-et)*exp(-xt/pa[5]);
    }
    /* Calculate weighted SS */
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  if(!isfinite(wss)) wss=nan(""); 
  wss*=penalty;
  if(0) printf("  wss=%g penalty=%g\n", wss, penalty);
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/