fit_dexp.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_exp(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *ymeas;
  double       *ysim;
  double       *w;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the sum of 1-5 exponential functions to decaying",
  "time-activity curve (TAC).",
  " ",
  "Function:",
  " ",
  "  f(x) = a1*exp(k1*x) + a2*exp(k2*x) + a3*exp(k3*x) + ...",
  " ",
  ", where a's > 0 and k's < 0.",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -5 | -4 | -3 | -2 | -1 | -1BL",
  "     Fit to sum of specified number of exponential functions, or",
  "     monoexponential with baseline; by default determined automatically.",
  " -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.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -nr=<value>",
  "     Number of basis functions, used in determination of initial parameter",
  "     estimates and number of exponentials for nonlinear fitting;",
  "     default is 100.",
  " -peak | -peak1",
  "     Fit is started only at the TAC peak, or at the next sample after peak.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If TAC file contains more than one TAC, only the first is used.",
  "Function parameters are written in fit format.",
  " ",
  "See also: extrapol, fit_exp, fit_fexp, fit2dat, fit2auc, paucinf, tacln",
  " ",
  "Keywords: TAC, curve fitting, extrapolation",
  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];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  int     bfNr=100;
  int     expnr=0;
  int     peakMode=0; // 0=fit all, 1=fit from peak, 2=fit from peak+1
  int     ret;
  unsigned int model=0;
  int     baseline=0; // last exponential term is freely fitted (0) or set to zero (1)
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=svgfile[0]=(char)0;
  /* 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, "NR=", 3)==0) {
      if(atoiCheck(cptr+3, &bfNr)==0 && bfNr>5) continue;
    } else if(strcasecmp(cptr, "5")==0 || strcasecmp(cptr, "5e")==0) {
      expnr=5; model=modelCodeIndex("5exp"); continue;
    } else if(strcasecmp(cptr, "4")==0 || strcasecmp(cptr, "4e")==0) {
      expnr=4; model=modelCodeIndex("4exp"); continue;
    } else if(strcasecmp(cptr, "3")==0 || strcasecmp(cptr, "3e")==0) {
      expnr=3; model=modelCodeIndex("3exp"); continue;
    } else if(strcasecmp(cptr, "2")==0 || strcasecmp(cptr, "2e")==0) {
      expnr=2; model=modelCodeIndex("2exp"); continue;
    } else if(strcasecmp(cptr, "1")==0 || strcasecmp(cptr, "1e")==0) {
      expnr=1; model=modelCodeIndex("1exp"); continue;
    } else if(strcasecmp(cptr, "1BL")==0 || strcasecmp(cptr, "1eBL")==0) {
      expnr=2; model=modelCodeIndex("2exp"); baseline=1; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "PEAK")==0) {
      peakMode=1; continue;
    } else if(strcasecmp(cptr, "PEAK1")==0) {
      peakMode=2; 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);
  }
 
 
  /* 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(expnr>0) printf("expnr := %d\n", expnr);
    if(expnr==2) printf("baseline := %d\n", baseline);
    printf("bfNr := %d\n", bfNr);
    printf("weights := %d\n", weights);
    printf("peakMode := %d\n", peakMode);
    fflush(stdout);
  }
 
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  ret=tacRead(&tac, tacfile, &status);
  if(ret!=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.tacNr>1) {
    fprintf(stderr, "Warning: file contains %d TACs; using the first one.\n", tac.tacNr);
    tac.tacNr=1;
  }
  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);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples.\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) {
    ret=tacWByFreq(&tac, ISOTOPE_UNKNOWN, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(2);
    }
  }
 
 
  /* Get the index of TAC peak */
  {
    TAC ptac; tacInit(&ptac);
    /* Take average of any duplicate samples */
    if(tacMultipleSamples(&tac, 1, &ptac, verbose-2)!=0) {
      fprintf(stderr, "Error: cannot process possible duplicates.\n");
      tacFree(&tac); tacFree(&ptac); return(2);
    }
    if(ptac.sampleNr<3) {
      fprintf(stderr, "Error: too few sample times.\n");
      tacFree(&tac); tacFree(&ptac); return(2);
    }
    /* Find peak value and its sample index */
    double ymin, ymax; int imax;
    if(tacYRange(&ptac, 0, &ymin, &ymax, NULL, &imax, NULL, NULL)!=0) {
      fprintf(stderr, "Error: no maximum found for TAC.\n");
      tacFree(&tac); tacFree(&ptac); return(2);
    }
    if(verbose>2) {
      printf("  minv := %g\n", ymin);
      printf("  maxv := %g\n", ymax);
      printf("  maxx := %g\n", ptac.x[imax]);
    }
    if(imax>0 && peakMode==0 && verbose>0) {
      fprintf(stderr, "Warning: TAC does not start with peak.\n");
    }
    if(!((ymax-ymin)>0.0) || !(ymax>0.0)) {
      fprintf(stderr, "Error: invalid TAC value range.\n");
      tacFree(&tac); tacFree(&ptac); return(2);
    }
    if(peakMode>0) {
      if(imax>ptac.sampleNr-3) {
        fprintf(stderr, "Error: too few sample times after the peak.\n");
        tacFree(&tac); tacFree(&ptac); return(2);
      }
      /* Delete samples before the peak */
      double tstart=ptac.x[imax]; if(peakMode==2) tstart=ptac.x[imax+1];
      if(verbose>1) printf("starting fit from %g\n", tstart);
      int ret=tacExtractRange(&tac, &tac, tstart, 1.1*ptac.x[ptac.sampleNr-1]);
      if(ret!=TPCERROR_OK) {
        fprintf(stderr, "Error: cannot exclude pre-peak data.\n");
        tacFree(&tac); tacFree(&ptac); return(2);
      }
    }
    tacFree(&ptac);
  }
 
 
  /*
   *  Get the first sample distance (that is > 0)
   */
  double sdist=0.0;
  {
    for(int i=1; i<tac.sampleNr; i++) {
      sdist=tac.x[i]-tac.x[i-1];
      if(sdist>0.0) break;
    }
    if(!(sdist>0.0)) {
      fprintf(stderr, "Error: invalid sample times.\n");
      tacFree(&tac); return(2);
    }
  }
  if(verbose>1) {
    printf("sdist := %g\n", sdist);
  }
 
  if(verbose>1) printf("determining initial parameter estimates\n");
  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);
  }
 
  /*
   *  Spectral analysis to find better range of parameters
   */
  if(verbose>2) printf("estimating initial k value range\n");
  double kmin, kmax; // kmin must be >0
  kmin=1.0E-06/xmax; // e^(-k*xmax) still almost 1
  kmax=3.0/(xmin+sdist); // e^(-k*(xmin+sdist) is about 0.01
  if(kmax<100.0*kmin) kmax=100.0*kmin;
 
// kmax=-log(yrange*1.0E-06)/tac.x[0]; if(kmax>100.) kmax=100.;
 
 
  if(verbose>2) printf(" kmin := %g\n kmax := %g\n", kmin, kmax);
  double pa[bfNr], pk[bfNr], yfit[tac.sampleNr];
  ret=spectralDExp(tac.x, tac.c[0].y, tac.w, tac.sampleNr, kmin, kmax, bfNr,
                   pk, pa, yfit, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  if(verbose>2) {
    printf("solutions for each k:\n");
    for(int bi=0; bi<bfNr; bi++) printf("\t%g\t%g\n", pk[bi], pa[bi]);
  }
  if(verbose>5) {
    printf("measured and fitted TAC:\n");
    for(int i=0; i<tac.sampleNr; i++)
      printf("\t%g\t%g\n", tac.c[0].y[i], yfit[i]);
  }
 
  /* Set better k value range */
  if(!(pa[0]>0.0 && pa[bfNr-1]>0.0)) {
    if(verbose>2) printf("refine k value range\n");
    ret=spectralKRange(pk, pa, bfNr, &kmin, &kmax, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(3);
    }
    if(kmax<2.0*kmin) {kmax=2.0*kmin; kmin*=0.5;}
    if(verbose>2) printf(" refined kmin := %g\n kmax := %g\n", kmin, kmax);
    ret=spectralDExp(tac.x, tac.c[0].y, tac.w, tac.sampleNr, kmin, kmax, bfNr,
                     pk, pa, yfit, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); return(3);
    }
    if(verbose>4) {
      printf("solutions for each k:\n");
      for(int bi=0; bi<bfNr; bi++) printf("\t%g\t%g\n", pk[bi], pa[bi]);
    }
    if(verbose>3) {
      printf("measured and fitted TAC:\n");
      for(int i=0; i<tac.sampleNr; i++)
        printf("\t%g\t%g\n", tac.c[0].y[i], yfit[i]);
    }
  }
  if(verbose>2) {
    double v, init_wss=0.0;
    for(int i=0; i<tac.sampleNr; i++) {
      v=tac.c[0].y[i]-yfit[i];
      init_wss+=tac.w[i]*v*v;
    }
    printf(" initial_wss := %g\n", init_wss);
  }
 
  /* Get max a parameter to be used as basis of upper limit */
  double amax=pa[doubleMaxIndex(pa, bfNr)];
  if(verbose>3) printf("max_a := %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); return(3);
  }
 
 
  /* If user did not give the nr of exponentials, then set it here */
  if(expnr==0) {
    expnr=fbfNr;
    if(expnr>5) expnr=5; // there should be no need for more
    if(verbose>1) printf("expnr := %d\n", expnr);
    char s[8]; sprintf(s, "%dexp", expnr);
    model=modelCodeIndex(s);
  }
 
  /* Extract parameters for at most expnr exponentials */
  double init_k[expnr], init_a[expnr];
  ret=spectralBFExtract(pk, pa, bfNr, init_k, init_a, expnr);
  if(ret<1) {
    fprintf(stderr, "Error: cannot estimate initial parameters.\n");
    tacFree(&tac); return(3);
  }
  if(verbose>2) {
    printf("\n\tClustered a and k values\n");
    for(int i=0; i<expnr; i++) printf("\t%g\t%g\n", init_a[i], init_k[i]);
    printf("\n"); fflush(stdout);
  }
 
  /* Copy parameters as initial values for nonlinear optimization */
  NLOPT ipar; nloptInit(&ipar);
  ret=nloptAllocate(&ipar, 2*expnr); // a[0], k[0], a[1], k[1], a[2], ...
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  ipar.totalNr=2*expnr;
  for(int ei=0; ei<expnr; ei++) {
    /* A */
    int pi=2*ei;
    ipar.xfull[pi]=init_a[ei];
    ipar.xlower[pi]=0.0;
    ipar.xupper[pi]=10.0*init_a[ei];
    ipar.xdelta[pi]=0.05*ipar.xfull[pi];
    ipar.xtol[pi]=fabs(0.001*ipar.xdelta[pi]);
    /* k */
    pi++;
    ipar.xfull[pi]=init_k[ei];
    ipar.xlower[pi]=0.0;
    ipar.xupper[pi]=10.0*ipar.xfull[pi]; if(ipar.xupper[pi]<10.0) ipar.xupper[pi]=10.0;
    ipar.xdelta[pi]=0.05*ipar.xfull[pi];
    ipar.xtol[pi]=fabs(0.001*ipar.xdelta[pi]);
  }
  if(expnr==2 && baseline!=0) {
    int pi=3; if(fabs(ipar.xfull[1])<fabs(ipar.xfull[3])) pi=1;
    ipar.xfull[pi]=ipar.xlower[pi]=ipar.xupper[pi]=ipar.xdelta[pi]=ipar.xtol[pi]=0.0;
  }
  nloptRemoveEmpties(&ipar);
  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);
 
 
  /*
   *  Nonlinear optimization:
   *  Simplex
   */
  if(verbose>1) printf("starting nonlinear optimization\n");
  /* Set function */
  ipar._fun=func_exp;
  /* Set data pointers for the fit */
  FITDATA fitdata;
  fitdata.n=tac.sampleNr;
  fitdata.x=tac.x;
  fitdata.ymeas=tac.c[0].y;
  fitdata.ysim=yfit;
  fitdata.w=tac.w;
  ipar.fundata=&fitdata;
  if(verbose>2) printf("1st non-linear optimization\n");
  ret=nloptSimplex(&ipar, 0, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&ipar); return(4);
  }
  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[0].y[i], yfit[i]);
  }
  if(verbose>3) {
    double v, final_wss=0.0;
    for(int i=0; i<tac.sampleNr; i++) {
      v=tac.c[0].y[i]-yfit[i];
      final_wss+=tac.w[i]*v*v;
    }
    printf(" final_wss := %g\n", final_wss);
  }
  /* Simplex again, with new deltas */
  for(unsigned int i=0; i<ipar.totalNr; i++) {
    if(fabs(ipar.xdelta[i])<1.0E-100) continue;
    if(ipar.xfull[i]>1.0E-06) ipar.xdelta[i]=0.0103*ipar.xfull[i];
    else ipar.xdelta[i]*=0.0103;
  }
  if(verbose>2) printf("non-linear optimization, 2nd time\n");
  ret=nloptSimplex(&ipar, 0, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&ipar); return(4);
  }
  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[0].y[i], yfit[i]);
  }
  if(verbose>3) {
    double v, final_wss=0.0;
    for(int i=0; i<tac.sampleNr; i++) {
      v=tac.c[0].y[i]-yfit[i];
      final_wss+=tac.w[i]*v*v;
    }
    printf(" final_final_wss := %g\n", final_wss);
  }
  /* and Simplex again, with new deltas */
  for(unsigned int i=0; i<ipar.totalNr; i++) {
    ipar.xdelta[i]=0.093*ipar.xdelta[i];
  }
  if(verbose>2) printf("non-linear optimization, 3rd time\n");
  ret=nloptSimplex(&ipar, 0, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&ipar); return(4);
  }
  /* Call objective function again with the final parameters */
  (void)func_exp(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[0].y[i], yfit[i]);
  }
  double final_wss=0.0;
  {
    double v;
    for(int i=0; i<tac.sampleNr; i++) {
      v=tac.c[0].y[i]-yfit[i];
      final_wss+=tac.w[i]*v*v;
    }
    if(verbose>3) printf(" final_final_final_wss := %g\n", final_wss);
  }
 
 
  /*
   *  Copy function parameters into PAR struct for printing and saving
   */
  PAR par; parInit(&par);
  ret=parAllocateWithTAC(&par, &tac, ipar.totalNr, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&ipar); return(6);
  }
  par.tacNr=1; par.parNr=ipar.totalNr;
  par.format=PAR_FORMAT_TSV_UK;
  par.r[0].model=model;
  par.r[0].dataNr=tacWSampleNr(&tac);
  par.r[0].start=tac.x[0]; 
  par.r[0].end=tac.x[tac.sampleNr-1];
  par.r[0].wss=final_wss;
  par.r[0].fitNr=ipar.totalNr-fixedNr;
  /* Set parameter names and units */
  for(int i=0; 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);
  }
  /* Set parameters in opposite order, and write k's as negatives */
  for(unsigned int i=0; i<ipar.totalNr; i+=2) {
    par.r[0].p[par.parNr-2-i]=ipar.xfull[i];
    par.r[0].p[par.parNr-1-i]=-ipar.xfull[i+1];
  }
  /* Print and save */
  if(verbose>0) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, 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); nloptFree(&ipar); parFree(&par); return(11);
    }
    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); nloptFree(&ipar); parFree(&par); return(12);
    }
    if(verbose>0) printf("parameters saved in %s\n", parfile);
  }
 
 
  /*
   *  Plot measured and fitted data, if requested
   */
  if(svgfile[0]) {
    if(verbose>1) printf("plotting measured and fitted data\n");
    /* Prepare place for function TAC */
    int snr;
    snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end-par.r[0].start, 2, NULL);
    TAC ftac; tacInit(&ftac);
    ret=tacAllocate(&ftac, snr, 1);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); nloptFree(&ipar); parFree(&par); return(21);
    }
    ftac.tacNr=1; ftac.sampleNr=snr;
    ret=tacCopyHdr(&tac, &ftac);
    if(ret==TPCERROR_OK) ret=tacCopyTacchdr(tac.c, ftac.c);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); nloptFree(&ipar); parFree(&par); tacFree(&ftac);
      return(22);
    }
    snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end-par.r[0].start, 2, ftac.x);
    ftac.isframe=0;
    for(int i=0; i<ftac.sampleNr; i++) ftac.x[i]+=par.r[0].start;
    /* Compute function values at sample times */
    ret=mfEvalY(modelCode(par.r[0].model), par.parNr, par.r[0].p, 
                ftac.sampleNr, ftac.x, ftac.c[0].y, verbose-2);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate function values.\n");
      tacFree(&tac); nloptFree(&ipar); parFree(&par); tacFree(&ftac);
      return(23);
    }
    if(verbose>10) tacWrite(&ftac, stdout, TAC_FORMAT_UNKNOWN, 1, NULL);
    /* Plot */
    ret=tacPlotFitSVG(&tac, &ftac, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); nloptFree(&ipar); parFree(&par); tacFree(&ftac);
      return(24);
    }
    tacFree(&ftac);
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
  }
 
 
  nloptFree(&ipar);
  tacFree(&tac);
  parFree(&par);
 
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_exp(int parNr, double *p, void *fdata)
{
  unsigned int i, ei;
  double v, wss;
  FITDATA *d=(FITDATA*)fdata;
 
  /* Calculate the curve values and weighted SS */
  for(i=0, wss=0.0; i<d->n; i++) {
    d->ysim[i]=0.0;
    if(isnan(d->ymeas[i]) || isnan(d->x[i])) continue;
    for(ei=0; ei<(unsigned int)parNr/2; ei++) 
      d->ysim[i]+=p[2*ei]*exp(-p[2*ei+1]*d->x[i]);
    v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
  if(0) {
    for(i=0; i<(unsigned int)parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/