fit_dms.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcbfm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_dmsurge(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *x2;
  double       *ymeas;
  double       *ysim;
  double       *w;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Fitting of the sum of surge functions and delay time",
  "to time-activity curves (TACs).",
  " ",
  "Function:",
  " ",
  "  f(x) = Sum of Ai*(x-dt)*exp(-Ki*(x-dt)) where i=1..N",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -n=<value>",
  "     Maximum number of surge functions (1<=n<=10); by default 5.",
  " -kMax=<value>",
  "     Maximum K value (kmax>kmin). By default 10.0.",
  " -kMin=<value>",
  "     Minimum K value (>0). By default, kMax*1.0E-10.",
  " -bfnr=<value>",
  "     Number of basis functions, used in determination of initial parameter",
  "     estimates and number of exponentials for nonlinear fitting;",
  "     default is 500.",
  " -dt=<min,max,step> | -dt=delay",
  "     The range and step size of precision of time delay estimation, or",
  "     fixed delay time.",
  " -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.",
  " -nllsq | -llsq",
  "     Non-linear LSQ fit (default), or only linear LSQ fit; In case of NLLSQ,",
  "     LLSQ is still performed first to obtain initial parameter values.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted TACs are written in specified TAC file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Function parameters are written the format determined by file name extension.",
  "PET time frames are used in fitting, if available in TAC file.",
  " ",
  "See also: fit_feng, fit_suri, fit_gvar, fit2dat, tacframe, fitdt",
  " ",
  "Keywords: TAC, curve fitting",
  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;
  char    tacfile[FILENAME_MAX], parfile[FILENAME_MAX], svgfile[FILENAME_MAX], fitfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  int     bfNr=500;
  int     maxSurgeNr=10;
  double  kMin=nan("");  
  double  kMax=10.0;  
  double  dtRange[3];
  int     nllsq=1; // whether to do NLLSQ after LLSQ
//  unsigned int model=MF_DMSURGE; // fid="dmsurge"
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=svgfile[0]=fitfile[0]=(char)0;
  for(int i=0; i<3; i++) dtRange[i]=nan("");
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    char *cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(!*cptr) continue;
    if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      if(atoiCheck(cptr+2, &maxSurgeNr)==0 && maxSurgeNr>0 && maxSurgeNr<=10) continue;
      fprintf(stderr, "Error: invalid option for number of surge functions '%s'.\n", argv[ai]);
      return(1);
    } else if(strncasecmp(cptr, "BFNR=", 5)==0) {
      if(atoiCheck(cptr+5, &bfNr)==0 && bfNr>5) continue;
    } else if(strncasecmp(cptr, "KMIN=", 5)==0) {
      kMin=atofVerified(cptr+5); if(kMin>0) continue;
    } else if(strncasecmp(cptr, "KMAX=", 5)==0) {
      kMax=atofVerified(cptr+5); if(kMax>0.0) continue;
    } else if(strncasecmp(cptr, "DT=", 3)==0) {
      int n=atofList(cptr+3, ",", dtRange, 3);
      if(n==1) {
        dtRange[1]=dtRange[0]; dtRange[2]=0.0; continue;
      } else if(n==3) {
        if(dtRange[0]<dtRange[1] && dtRange[2]<0.2*(dtRange[1]-dtRange[0])) continue;
      }
      fprintf(stderr, "Error: invalid delay time option '%s'.\n", argv[ai]);
      return(1);
    } else if(strcasecmp(cptr, "NLLSQ")==0) {
      nllsq=1; continue;
    } else if(strcasecmp(cptr, "LLSQ")==0) {
      nllsq=0; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-3;
  
  /* 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);
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  /* Did we get all the information that we need? */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* Check k range, and set kMin to default, unless user provided */
  if(!isfinite(kMin)) kMin=kMax*1.0E-10;
  else if(!(kMax>kMin)) {
    fprintf(stderr, "Error: invalid k limits.\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("tacfile := %s\n", tacfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(fitfile[0]) printf("fitfile := %s\n", fitfile);
    printf("nllsq := %d\n", nllsq);
    printf("maxSurgeNr := %d\n", maxSurgeNr);
    printf("bfNr := %d\n", bfNr);
    printf("kmin := %g\n", kMin);
    printf("kmax := %g\n", kMax);
    if(isfinite(dtRange[0]))
      printf("dtMin := %g\ndtMax := %g\ndtStep :=%g\n", dtRange[0], dtRange[1], dtRange[2]);
    printf("weights := %d\n", weights);
    fflush(stdout);
  }
 
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("isframe := %d\n", tac.isframe);
  }
  if(tac.tacNr<1) {
    fprintf(stderr, "Error: file contains no data.\n");
    tacFree(&tac); return(2);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Check and set weights */
  if(weights==0) {
    if(!tacIsWeighted(&tac)) {
      tac.weighting=WEIGHTING_OFF;
      for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
    } 
  } else if(weights==1) {
    tac.weighting=WEIGHTING_OFF;
    for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
  } else if(weights==2) {
    int ret=tacWByFreq(&tac, ISOTOPE_UNKNOWN, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(2);
    }
  }
 
 
  /*
   *  Get and check data range
   */
  double ymin, ymax, xmin, xmax;
  tacYRange(&tac, 0, &ymin, &ymax, NULL, NULL, NULL, NULL);
  tacXRange(&tac, &xmin, &xmax);
  if(verbose>2) {
    printf("  fitrange_minv := %g\n", ymin);
    printf("  fitrange_maxv := %g\n", ymax);
    printf("  fitrange_minx := %g\n", xmin);
    printf("  fitrange_maxx := %g\n", xmax);
  }
  if(!(xmax>xmin)) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac); return(2);
  }
  /* Get the first smallest sample distance (delta x) that is > 0 */
  double sdist=xmax;
  for(int i=1; i<tac.sampleNr; i++) {
    double dx=tac.x[i]-tac.x[i-1];
    if(dx>0.0 && dx<sdist) sdist=dx;
  }
  if(verbose>1) printf("sdist := %g\n", sdist);
  /* Set delay time range if necessary */
  if(!isfinite(dtRange[0])) {
    dtRange[1]=0.1*(xmax-xmin);
    dtRange[0]=-dtRange[1];
    dtRange[2]=0.01*(dtRange[1]-dtRange[0]);
//    dtRange[2]=0.001*(dtRange[1]-dtRange[0]);
    if(verbose>1)
      printf("dtMin := %g\ndtMax := %g\ndtStep :=%g\n", dtRange[0], dtRange[1], dtRange[2]);
  }
 
 
  /*
   *  Prepare the room for fitted TTACs, if requested
   */
  TAC ftac; tacInit(&ftac);
  if(fitfile[0] || svgfile[0]) {
    if(verbose>1) printf("allocating space for fitted TTACs\n");
    if(tacDuplicate(&tac, &ftac)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot allocate space for fitted TACs.\n");
      tacFree(&tac);
      return(3);
    }
  }
 
 
  /*
   *  Set PAR structure for the results
   */
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, 1+2*maxSurgeNr, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); tacFree(&ftac); return(4);
  }
  par.tacNr=tac.tacNr;
  par.parNr=1+2*maxSurgeNr;
  par.format=PAR_FORMAT_TSV_UK;
  /* Set parameter names and units */
  sprintf(par.n[0].name, "dT");
  par.n[0].unit=tac.tunit;
  for(int i=1; i<par.parNr; i+=2) {
    sprintf(par.n[i].name, "a%d", 1+i/2);
    sprintf(par.n[i+1].name, "k%d", 1+i/2);
    par.n[i].unit=tac.cunit;
    par.n[i+1].unit=unitInverse(tac.tunit);
  }
  int maxParNr=0;
 
 
  /*
   *  Spectral analysis to find good initial parameters for non-linear fitting for each TAC
   */
  for(int ci=0; ci<tac.tacNr; ci++) {
    if(verbose>1 && tac.tacNr>1) printf("linear fitting of TAC %s\n", tac.c[ci].name);
 
    /* First Spectral analysis */
    double pa[bfNr], pk[bfNr], yfit[tac.sampleNr], dt=nan("");
    int ret=0;
    if(tac.isframe)
      ret=spectralDMSurge(tac.x1, tac.x2, tac.c[ci].y, tac.w, tac.sampleNr, kMin, kMax, bfNr, 
                         dtRange[0], dtRange[1], dtRange[2], 
                         pk, pa, &dt, yfit, &status);
    else
      ret=spectralDMSurge(tac.x, NULL, tac.c[ci].y, tac.w, tac.sampleNr, kMin, kMax, bfNr, 
                         dtRange[0], dtRange[1], dtRange[2], 
                         pk, pa, &dt, yfit, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
    }
    if(verbose>2) {
      printf("delta_t := %g\n", dt);
      printf("solutions for each k:\n");
      for(int bi=0; bi<bfNr; bi++) printf("\t%g\t%g\n", pk[bi], pa[bi]);
    }
    /* Spectral analysis refined */
    if(verbose>2) printf("refining parameters of Spectral analysis\n");
    double kMin2, kMax2, dtRange2[3];
    if(spectralKRange(pk, pa, bfNr, &kMin2, &kMax2, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
    }
    if(kMax2<6.0*kMin2) {kMax2*=3.0; kMin2*=0.75; if(kMin2<kMin) kMin2=kMin;}
    dtRange2[0]=dt-2.0*dtRange[2]; if(dtRange2[0]<dtRange[0]) dtRange2[0]=dtRange[0];
    dtRange2[1]=dt+2.0*dtRange[2]; if(dtRange2[1]>dtRange[1]) dtRange2[1]=dtRange[1];
    dtRange2[2]=0.02*(dtRange2[1]-dtRange2[0]);
    if(verbose>2) {
      printf(" refined kMin := %g\n refined kMax := %g\n", kMin2, kMax2);
      printf(" refined delay range and step size : %g %g %g\n", dtRange2[0], dtRange2[1], dtRange2[2]);
    }
    if(tac.isframe)
      ret=spectralDMSurge(tac.x1, tac.x2, tac.c[ci].y, tac.w, tac.sampleNr, kMin2, kMax2, bfNr, 
                         dtRange2[0], dtRange2[1], dtRange2[2], 
                         pk, pa, &dt, yfit, &status);
    else
      ret=spectralDMSurge(tac.x, NULL, tac.c[ci].y, tac.w, tac.sampleNr, kMin2, kMax2, bfNr, 
                         dtRange2[0], dtRange2[1], dtRange2[2], 
                         pk, pa, &dt, yfit, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
    }
    if(verbose>2) {
      printf("refined delta_t := %g\n", dt);
      printf("refined solutions for each k:\n");
      for(int bi=0; bi<bfNr; bi++) printf("\t%g\t%g\n", pk[bi], pa[bi]);
    }
    if(ftac.tacNr>0)
      for(int i=0; i<tac.sampleNr; i++) ftac.c[ci].y[i]=yfit[i];
 
    /* Get max A parameter to be used as basis of upper limit */
    double aMax=pa[doubleMaxIndex(pa, bfNr)];
    if(verbose>3) printf("aMax := %g\n", aMax);
 
    /* Calculate how many zero-separated basis functions there really were */
    int fbfNr=spectralBFNr(pk, pa, bfNr);
    if(verbose>2) printf("fbfNr := %d\n", fbfNr);
    if(fbfNr<1) {
      fprintf(stderr, "Error: cannot estimate initial parameters.\n");
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
    }
 
    /* Set the nr of Surge functions, at most maxSurgeNr */
    int surgeNr=fbfNr;
    if(surgeNr>maxSurgeNr) surgeNr=maxSurgeNr;
    if(verbose>1) printf("surgeNr := %d\n", surgeNr);
    int localParNr=1+2*surgeNr;
    if(localParNr>maxParNr) maxParNr=localParNr;
 
    /* Extract parameters for specified number of Surge functions at most */
    double init_k[surgeNr], init_a[surgeNr];
    ret=spectralBFExtract(pk, pa, bfNr, init_k, init_a, surgeNr);
    if(ret<1) {
      fprintf(stderr, "Error: cannot estimate initial parameters.\n");
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
    }
    if(verbose>2) {
      printf("\n\tClustered a and k values\n");
      for(int i=0; i<surgeNr; i++) printf("\t%g\t%g\n", init_a[i], init_k[i]);
      printf("\n"); fflush(stdout);
    }
 
    /* Set parameters in PAR structure */
    par.r[ci].model=modelCodeIndex("dmsurge");
    //printf("model=%d\nold_id=%d\n", par.r[ci].model, modelOldId(par.r[ci].model));
    par.r[ci].dataNr=tacWSampleNr(&tac);
    par.r[ci].start=xmin; 
    par.r[ci].end=xmax;
    par.r[ci].wss=0.0;
    par.r[ci].fitNr=2*surgeNr; if(dtRange[2]>0.0) par.r[ci].fitNr++;
    par.r[ci].p[0]=dt;
    for(int i=0, pi=1; i<surgeNr && pi<par.parNr; i++, pi+=2) par.r[ci].p[pi]=init_a[i];
    for(int i=0, pi=2; i<surgeNr && pi<par.parNr; i++, pi+=2) par.r[ci].p[pi]=init_k[i];
    for(int i=localParNr; i<par.parNr; i++) par.r[ci].p[i]=0.0;
 
    /* If not refined with NLLSQ, calculate fitted TAC and WSS with final LLSQ parameters */
    if(nllsq==0) {
      if(ftac.isframe)
        ret=mfEvalFrameY("dmsurge", localParNr, par.r[ci].p, ftac.sampleNr, ftac.x1, ftac.x2, yfit, 0);
      else
        ret=mfEvalY("dmsurge", localParNr, par.r[ci].p, ftac.sampleNr, ftac.x, yfit, 0);
      if(ret) {
        fprintf(stderr, "Error: cannot calculate fitted TAC.\n");
        tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
      }
      for(int i=0; i<tac.sampleNr; i++) {
        double v=tac.c[ci].y[i]-yfit[i];
        par.r[ci].wss+=tac.w[i]*v*v;
      }
 
      /* Copy fitted TAC if needed later */
      if(ftac.tacNr>0) for(int i=0; i<tac.sampleNr; i++) ftac.c[ci].y[i]=yfit[i];
    }
 
 
 
    /* If needed, calculate fitted TAC with final LLSQ parameters */
    if(nllsq==0 && ftac.tacNr>0) {
      if(ftac.isframe)
        ret=mfEvalFrameY("dmsurge", localParNr, par.r[ci].p, ftac.sampleNr, ftac.x1, ftac.x2, ftac.c[ci].y, 0);
      else
        ret=mfEvalY("dmsurge", localParNr, par.r[ci].p, ftac.sampleNr, ftac.x, ftac.c[ci].y, 0);
      if(ret) {
        fprintf(stderr, "Error: cannot calculate fitted TAC.\n");
        tacFree(&tac); tacFree(&ftac); parFree(&par); return(6);
      }
    }
  }
  /* Set parameter number to the actual max among TACs */
  par.parNr=maxParNr;
 
 
  /*
   *  Non-linear fitting of each TAC using initial estimates from above, if requested.
   *  Delay time is constrained to the estimate from LLSQ.
   */
  if(nllsq!=0) {
    for(int ci=0; ci<tac.tacNr; ci++) {
      if(verbose>1 && tac.tacNr>1) printf("non-linear fitting of TAC %s\n", tac.c[ci].name);
      double final_wss=0.0;
 
      /* Prepare for LLSQ optimization */
      NLOPT ipar; nloptInit(&ipar); // dt, a[0], k[0], a[1], k[1], a[2], ...
      if(nloptAllocate(&ipar, par.parNr)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); tacFree(&ftac); parFree(&par); return(8);
      }
      ipar.totalNr=par.parNr;
      for(unsigned int i=0; i<ipar.totalNr; i++) ipar.xfull[i]=par.r[ci].p[i];
      /* deltaT */
      ipar.xlower[0]=ipar.xupper[0]=ipar.xfull[0];
      ipar.xdelta[0]=0.0;
      ipar.xtol[0]=0.0;
      for(unsigned int i=1; i<ipar.totalNr; i++) {
        if(i%2) { // A
          ipar.xlower[i]=0.0;
          ipar.xupper[i]=100.*ipar.xfull[i];
          ipar.xdelta[i]=0.01*ipar.xfull[i];
          ipar.xtol[i]=0.0001*ipar.xfull[i];
        } else { // k
          if(ipar.xfull[i]<kMin) ipar.xfull[i]=kMin;
          ipar.xlower[i]=0.01*ipar.xfull[i]; //if(ipar.xlower[i]<kMin) ipar.xlower[i]=kMin;
          ipar.xupper[i]=10.*ipar.xfull[i];
          ipar.xdelta[i]=0.01*ipar.xfull[i];
          ipar.xtol[i]=0.0001*ipar.xfull[i];
        } 
      }
      if(verbose>2) {
        printf("Initial parameters for optimization:\n");
        nloptWrite(&ipar, stdout); fflush(stdout);
      }
      /* Get nr of fixed parameters before deltas are modified between fittings */
      unsigned int fixedNr=nloptFixedNr(&ipar);
      /* Set function */
      ipar._fun=func_dmsurge;
      /* Set data pointers for the fit */
      FITDATA fitdata;
      fitdata.n=tac.sampleNr;
      if(!tac.isframe) {fitdata.x=tac.x; fitdata.x2=NULL;} else {fitdata.x=tac.x1; fitdata.x2=tac.x2;}
      fitdata.ymeas=tac.c[ci].y;
      fitdata.ysim=ftac.c[ci].y;
      fitdata.w=tac.w;
      ipar.fundata=&fitdata;
      double initial_wss=func_dmsurge(ipar.totalNr, ipar.xfull, ipar.fundata);
      if(verbose>3)
        printf(" initial_wss := %g\n", initial_wss);
      if(!isfinite(initial_wss)) {
        fprintf(stderr, "Error: invalid initial parameter guess.\n");
        tacFree(&tac); tacFree(&ftac); parFree(&par); nloptFree(&ipar); return(8);
      }
 
      /* Nonlinear optimization: Simplex */
      if(verbose>2) {printf("starting nonlinear optimization\n"); fflush(stdout);}
      if(verbose>3) printf("1st non-linear optimization\n");
      if(nloptSimplex(&ipar, 0, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); tacFree(&ftac); parFree(&par); nloptFree(&ipar); return(8);
      }
      final_wss=func_dmsurge(ipar.totalNr, ipar.xfull, ipar.fundata);
      if(verbose>3) nloptWrite(&ipar, stdout);
      if(verbose>6) {
        printf("measured and fitted TAC:\n");
        for(int i=0; i<tac.sampleNr; i++)
          printf("\t%g\t%g\n", tac.c[ci].y[i], ftac.c[ci].y[i]);
      }
      if(verbose>3) printf(" final_wss := %g\n", final_wss);
      /* Simplex again, with new deltas and smaller tolerance */
      for(unsigned int i=1; i<ipar.totalNr; i++) {
        ipar.xdelta[i]*=0.093;
        ipar.xtol[i]*=0.01;
      }
      if(verbose>3) printf("2nd non-linear optimization\n");
      if(nloptSimplex(&ipar, 0, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        tacFree(&tac); tacFree(&ftac); parFree(&par); nloptFree(&ipar); return(8);
      }
      final_wss=func_dmsurge(ipar.totalNr, ipar.xfull, ipar.fundata);
      if(verbose>3) nloptWrite(&ipar, stdout);
      if(verbose>6) {
        printf("measured and fitted TAC:\n");
        for(int i=0; i<tac.sampleNr; i++)
          printf("\t%g\t%g\n", tac.c[ci].y[i], ftac.c[ci].y[i]);
      }
      if(verbose>3) printf(" final_wss := %g\n", final_wss);
 
      /* Copy results into PAR structure */
      for(unsigned int i=0; i<ipar.totalNr; i++)
        par.r[ci].p[i]=ipar.xfull[i];
      par.r[ci].wss=final_wss;
      par.r[ci].fitNr=ipar.totalNr-fixedNr;
 
      nloptFree(&ipar);
    }
  }
 
 
  /* Print and save function parameters */
  if(verbose>0) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 0, NULL);
  if(parfile[0]) {
    /* Set file format */
    par.format=parFormatFromExtension(parfile);
    if(verbose>2) printf("parameter file format := %s\n", parFormattxt(par.format));
    if(par.format==PAR_FORMAT_UNKNOWN) par.format=PAR_FORMAT_FIT;
    /* set header contents */
    {
      iftPut(&par.h, "datafile", tacfile, 0, NULL);
      {
      char buf[256];
      time_t t=time(NULL);
      iftPut(&par.h, "analysis_time", ctime_r_int(&t, buf), 0, NULL);
      tpcProgramName(argv[0], 1, 1, buf, 256);
      iftPut(&par.h, "program", buf, 0, NULL);
      }
    }
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    FILE *fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&ftac); parFree(&par); return(12);
    }
    if(verbose>0) printf("function parameters saved in %s\n", parfile);
  }
 
 
  /*
   *  SVG plot of fitted and original data
   */
  if(svgfile[0]) {
 
    if(verbose>1) printf("saving SVG plot\n");
    if(tacPlotFitSVG(&tac, &ftac, "", 0.0, nan(""), 0.0, nan(""), svgfile, NULL)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot plot fitted data.\n");
      tacFree(&tac); tacFree(&ftac); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("Plots written in %s.\n", svgfile);
  }
 
  /*
   *  Save fitted TTACs
   */
  if(fitfile[0]) {
    if(verbose>1) printf("writing %s\n", fitfile);
    FILE *fp; fp=fopen(fitfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", fitfile);
      tacFree(&tac); tacFree(&ftac); parFree(&par);
      return(31);
    }
    int ret=tacWrite(&ftac, fp, TAC_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&ftac); parFree(&par);
      return(32);
    }
    if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
  }
 
  tacFree(&tac); tacFree(&ftac); parFree(&par);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_dmsurge(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  /* Calculate the curve values and weighted SS */
  if(d->x2!=NULL) { // PET frame start and end times available
    if(mfEvalFrameY("dmsurge", parNr, p, d->n, d->x, d->x2, d->ysim, 0)) return(nan(""));
  } else { 
    if(mfEvalY("dmsurge", parNr, p, d->n, d->x, d->ysim, 0)) return(nan(""));
  }
  double wss=0.0;
  for(unsigned i=0; i<d->n; i++) {
    double v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/