fit_sinf.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the exponential based function (Wagner, 1976)",
  "to the PET plasma or blood time-activity curve (TAC). This function can",
  "be used when PET tracer is introduced into venous blood as a short infusion.",
  "Function: f(t)=",
  "  if t<=Ta :",
  "     0",
  "  if t>Ta and t<Ta+Ti :",
  "     (1/Ti)*Sum[i=1..n, (Ai/Li)*(1-exp(-Li*(t-Ta)))]",
  "  if t>=Ta+Ti :",
  "     (1/Ti)*Sum[i=1..n, (Ai/Li)*(exp(-Li*(t-Ta-Ti)) - exp(-Li*(t-Ta)))]",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
/*" -lim=<filename>",
  "     Specify the constraints for model parameters;",
  "     This file with default values can be created by giving this",
  "     option as the only command-line argument to this program.",*/
  " -Ti=<infusion time>",
  "     Duration of tracer infusion; 0, if short bolus.",
  " -Ta=<appearance time>",
  "     Time when tracer concentration starts ascending.",
  " -tau1=<Dispersion time constant (s)>",
  " -tau2=<2nd dispersion time constant (s)>",
  "     Dispersion is added to function, also to saved curve and plot.",
  " -n=<1|2|3|4|5|A>",
  "     The nr of summed functions; A=determined automatically (default).",
  "     Changing the default may lead to a bad fit.",
  " -wf",
  "     Weight by sampling interval.",
  " -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.",
  " -svg1=<Filename>",
  "     Initial part of fitted and measured TACs are plotted in SVG file",
  " -svg2=<Filename>",
  "     Lower part of fitted and measured TACs are plotted in SVG file",
  " -fit=<Filename>",
  "     Fitted TACs are written in TAC format.",
  " -dcfit=<Filename>",
  "     Fitted dispersion corrected TAC is written in given file;",
  "     same output as with -fit when tau1 and tau2 are not set.", 
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC file must contain two columns, sample times (x) and concentrations (y);",
  "any additional columns are ignored. For a good fit, TAC should be corrected",
  "for physical decay and circulating metabolites.",
  "Function parameters (Ta, Ti, A1, L1, ...) will be written in the parfile.",
  " ",
  "References:",
  "1. Wagner JG. J Pharmacokin Biopharm. 1976;4:443-467.",
  "2. Kudomi N et al. Eur J Nucl Med Mol Imaging 2008;35:1899-1911.",
  " ",
  "See also: fit2dat, fit_feng, fit_exp, fit_ratf, extrapol, disp4dft",
  " ",
  "Keywords: curve fitting, input, modelling, simulation, dispersion",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
#define MAX_FUNC_NR 5
double tau1=0.0, tau2=0.0;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
int dataNr=0;
double *x, *ymeas, *yfit, *w, *ytmp;
double func331(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  int      ret;
  char    *cptr, progname[256];
  char     tacfile[FILENAME_MAX], fitfile[FILENAME_MAX], limfile[FILENAME_MAX],
           parfile[FILENAME_MAX], resfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  char     svgfile1[FILENAME_MAX], svgfile2[FILENAME_MAX], dcfile[FILENAME_MAX];
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  int      MRL_check=0; // 0=off, 1=on
  double   fixedTi, fixedTa;
  int      sumn=0; // Nr of summed functions: 0=automatic, otherwise 1-5.
  DFT      tac;
  double   (*func)(int parNr, double *p, void*);
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=fitfile[0]=dcfile[0]=limfile[0]=parfile[0]=resfile[0]=(char)0;
  svgfile[0]=svgfile1[0]=svgfile2[0]=(char)0;
  fixedTa=fixedTi=nan("");
  dftInit(&tac); //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, "MRL")==0) {
      MRL_check=1; continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0) {
      if(strlcpy(svgfile, cptr+4, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "SVG1=", 5)==0) {
      if(strlcpy(svgfile1, cptr+5, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "SVG2=", 5)==0) {
      if(strlcpy(svgfile2, cptr+5, FILENAME_MAX)>0) continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strlcpy(fitfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "DCFIT=", 6)==0 && strlen(cptr)>4) {
      strlcpy(dcfile, cptr+6, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "RES=", 4)==0 && strlen(cptr)>4) {
      strlcpy(resfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      cptr+=2; if(strcasecmp(cptr, "A")==0) {sumn=0; continue;}
      sumn=atoi(cptr); if(sumn>=1 && sumn<=MAX_FUNC_NR) continue;
    } else if(strncasecmp(cptr, "tau=", 4)==0) {
      if(!atof_with_check(cptr+4, &tau1) && tau1>=0.0) continue;
    } else if(strncasecmp(cptr, "tau1=", 5)==0) {
      if(!atof_with_check(cptr+5, &tau1) && tau1>=0.0) continue;
    } else if(strncasecmp(cptr, "tau2=", 5)==0) {
      if(!atof_with_check(cptr+5, &tau2) && tau2>=0.0) continue;
    } else if(strncasecmp(cptr, "Ta=", 3)==0) {
      if(!atof_with_check(cptr+3, &fixedTa) && fixedTa>=0.0) continue;
    } else if(strncasecmp(cptr, "Ti=", 3)==0) {
      if(!atof_with_check(cptr+3, &fixedTi) && fixedTi>=0.0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
 
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) {strlcpy(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(parfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  if(!parfile[0]) strcpy(parfile, "stdout");
 
  /* Is something missing? */
  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);
    printf("parfile := %s\n", parfile);
    printf("fitfile := %s\n", fitfile);
    printf("dcfile := %s\n", dcfile);
    printf("resfile := %s\n", resfile);
    printf("svgfile := %s\n", svgfile);
    printf("svgfile1 := %s\n", svgfile1);
    printf("svgfile2 := %s\n", svgfile2);
    printf("MRL_check := %d\n", MRL_check);
    printf("weights := %d\n", weights);
    if(sumn>0) printf("sumn := %d\n", sumn);
    if(!isnan(fixedTa)) printf("fixedTa := %g\n", fixedTa);
    if(!isnan(fixedTi)) printf("fixedTi := %g\n", fixedTi);
  }
 
 
 
  /*
   *  Read data file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  if(dftRead(tacfile, &tac)) {
    fprintf(stderr, "Error in reading '%s': %s\n", tacfile, dfterrmsg);
    return(2);
  }
  if(verbose>1) printf("checking %s\n", tacfile);
 
  /* Check voiNr */
  if(tac.voiNr>1) {
    fprintf(stderr, "Warning: %s contains more than one TAC; first is used.\n", tacfile);
    tac.voiNr=1;
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&tac);
  if(verbose>30) dftPrint(&tac);
 
  /* Get min and max X and Y */
  double tstart, tstop, miny, maxy;
  ret=dftMinMax(&tac, &tstart, &tstop, &miny, &maxy);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  if(tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  if(tstart<0.0) fprintf(stderr, "Warning: negative x value(s).\n");
  if(miny<0.0) fprintf(stderr, "Warning: negative y value(s).\n");
  if(verbose>1) {
    printf("x_range := %g - %g\n", tstart, tstop);
    printf("sample_number := %d\n", tac.frameNr);
  }
 
  /* Check for NA's */
  if(dft_nr_of_NA(&tac) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  /* Remove negative or zero values from the end of data, if only one TAC; */
  /* this is often needed to process a TAC measured using blood pump */
  if(tac.voiNr==1) {
    int i; for(i=tac.frameNr-1; i>=0; i--) if(tac.voi[0].y[i]>0.0) break;
    tac.frameNr=i+1;
  }
  /* Check frameNr */
  if(tac.frameNr<4) {
    fprintf(stderr, "Error: %s does not contain enough valid samples.\n", tacfile);
    dftEmpty(&tac); return(3);
  }
 
  /* Set study number, if not yet available */
  if(!tac.studynr[0]) studynr_from_fname(tacfile, tac.studynr);
  if(verbose>29) dftPrint(&tac);
 
  /* Check and set weights */
  if(weights==0) {
    if(tac.isweight==0) {
      tac.isweight=1; for(int i=0; i<tac.frameNr; i++) tac.w[i]=1.0;}
  } else if(weights==1) {
    tac.isweight=1; for(int i=0; i<tac.frameNr; i++) tac.w[i]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&tac);
  }
  if(verbose>2) {
    printf("data_weights := %g", tac.w[0]);
    for(int i=1; i<tac.frameNr; i++) printf(", %g", tac.w[i]);
    printf("\n");
  }
 
 
  /*
   *  Convert dispersion time constants to the same units as sample times
   */
  if(tau1>0.0 || tau2>0.0) {
    if(tac.timeunit==TUNIT_SEC) {
      if(verbose>3) printf("data and tau's are in the same units.\n");
    } else if(tac.timeunit==TUNIT_MIN) {
      tau1/=60.0; tau2/=60.0;
      if(verbose>2) printf("tau's are converted to minutes.\n");
    } else {
      fprintf(stderr, "Warning: tau's could not be converted to sample units.\n");
    }
  }
  if(verbose>1) {
    printf("tau1 := %g\n", tau1);
    printf("tau2 := %g\n", tau2);
  }
 
 
 
  /*
   *  Find the sample time and index of max concentration
   */
  double tmax; int imax;
  ret=dftMinMaxTAC(&tac, 0, NULL, &tmax, NULL, NULL, NULL, NULL, NULL, &imax);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
  if(verbose>1) {
    printf("peak_x := %g\n", tmax);
    printf("peak_y := %g\n", maxy);
    printf("peak_index := %d\n", imax);
  }
  if((tac.frameNr-imax)<2) {
    fprintf(stderr, "Error: invalid TAC shape in %s\n", tacfile);
    dftEmpty(&tac); return(2);
  }
 
 
 
  /*
   *  Determine the nr of sum functions and/or parameter constraints, unless
   *  already set by the user.
   */
  double p1[MAX_FUNC_NR], p2[MAX_FUNC_NR];
  int proposed_sumn;
  {
    int i, j, n;
    /* Make copy of the data starting from the peak, without any missing values */
    n=tac.frameNr-imax;
    double *x, *y;
    x=(double*)malloc(2*n*sizeof(double));
    if(x==NULL) {
      fprintf(stderr, "Error: out of memory.\n");
      dftEmpty(&tac); return(3);
    }
    y=x+n;
    for(i=imax, j=0; i<tac.frameNr && j<n; i++) {
      x[j]=tac.x[i]-tmax; y[j]=tac.voi[0].y[i];
      if(!isnan(x[j]) && !isnan(y[j])) j++;
    }
    n=j;
    if(n<2) {
      fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
      dftEmpty(&tac); free(x); return(3);
    }
    ret=fitExpDecayNNLS(x, y, n, -1.0, 1.0E-06, 5.0E+01, MAX_FUNC_NR, p1, p2, &n, verbose-2);
    if(ret || n<1) {
      fprintf(stderr, "Error: cannot determine constraints from the data.\n");
      dftEmpty(&tac); free(x); return(3);
    }
    free(x);
    if(verbose>1) {
      printf("Initial guesses for exp parameters:\n");
      for(i=0; i<n; i++) printf("\t%g\t%g\n", p1[i], p2[i]);
    }
    proposed_sumn=n;
    if(sumn==0) sumn=n; // user accepted automatic setting of sumn
  }
 
 
  /*
   *  Allocate memory for the fit parameters
   */
  if(verbose>1) printf("allocating memory for fit parameters.\n");
  FIT fit; fitInit(&fit);
  ret=fit_allocate_with_dft(&fit, &tac);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for fit parameters.\n");
    dftEmpty(&tac); return(4);
  }
  tpcProgramName(argv[0], 1, 1, fit.program, 1024);
  strlcpy(fit.datafile, tacfile, FILENAME_MAX);
  fit.time=time(NULL);
  fit.voi[0].type=331;
  fit.voi[0].start=tstart; fit.voi[0].end=tstop;
  fit.voi[0].dataNr=tac.frameNr;
 
 
  
  /*
   *  Set parameter constraints 
   */
  int parNr=2*sumn+2;
  fit.voi[0].parNr=parNr;
  /* Set Ta */
  if(isnan(fixedTa)) {
    pmin[0]=0.0; pmax[0]=tmax;
  } else {
    pmin[0]=pmax[0]=fixedTa;
  }
  /* Set Ti */
  if(isnan(fixedTi)) {
    pmin[1]=0.02*tmax; pmax[1]=1.5*tmax;
  } else {
    pmin[1]=pmax[1]=fixedTi;
  }
  /* Set the rest of the parameters */
  if(sumn==proposed_sumn) {
    /* User lets we set as many constraints as we feel right */
    int i, pi;
    for(i=0, pi=2; i<sumn; i++, pi+=2) {
      pmin[pi]=0.25*p1[i];      pmax[pi]=2.5*p1[i];
      //pmin[pi]=0.25*p1[i]*(-p2[i]);   pmax[pi]=2.5*p1[i]*(-p2[i])*tmax;
      pmin[pi+1]=0.15*(-p2[i]); pmax[pi+1]=2.0*(-p2[i]);
    }
    pmin[3]=1.0E-10;
  } else if(sumn>proposed_sumn) {
    /* User wants us to add more functions */
    int i, j, pi;
    for(i=0, pi=2; i<proposed_sumn; i++, pi+=2) {
      pmin[pi]=0.25*p1[i];      pmax[pi]=2.5*p1[i];
      pmin[pi+1]=0.15*(-p2[i]); pmax[pi+1]=2.0*(-p2[i]);
    }
    for(j=1; i<sumn; i++, j++) {
      pmin[pi]=0.0;       pmax[pi]=0.2*maxy*(double)j;
      pmin[pi+1]=1.0E-10; pmax[pi+1]=0.1*(double)j;
    }
  } else {
    /* User wants us to use fewer functions */
    int i, pi;
    for(i=0, pi=2; i<sumn; i++, pi+=2) {
      pmin[pi]=0.05*p1[i];      pmax[pi]=10.0*p1[i];
      pmin[pi+1]=0.05*(-p2[i]); pmax[pi+1]=5.0*(-p2[i]);
    }
    pmin[3]=1.0E-10;
  }
 
 
  
  /*
   *  Fit 
   */
  {
    if(verbose>2) printf("fitting\n");
    dataNr=tac.frameNr;
    x=tac.x; ymeas=tac.voi[0].y; yfit=tac.voi[0].y2; ytmp=tac.voi[0].y3; w=tac.w;
    int iter_nr=1; 
    int tgo_nr=400; 
    int neigh_nr=8;
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=0;
    func=func331;
    double wss, p[MAX_PARAMETERS];
    ret=tgo(pmin, pmax, func, NULL, parNr, neigh_nr, &wss, p, tgo_nr, iter_nr, verbose-5);
    if(ret) {
      fprintf(stderr, "Error %d in TGO.\n", ret);
      dftEmpty(&tac); fitEmpty(&fit);
      return(6);
    }
    if(verbose>2) {
      printf("fitted parameters:\n");
      for(int i=0; i<parNr; i++)
        printf("\tp%d\t%g\t(%g - %g)\n", 1+i, p[i], pmin[i], pmax[i]);
      printf("\tWSS\t:=\t%g\n", wss);
    }
    if(verbose>3) {
      printf("Sample\tX\tY\tYfit\tWeight\n");
      for(int i=0; i<dataNr; i++) printf("%d\t%9.2e\t%9.2e\t%9.2e\t%7.1e\n", 
                                         1+i, x[i], ymeas[i], yfit[i], w[i]);
    }
    /* Check the MRL */
    if(MRL_check) {
      if(verbose>1) printf("Checking the MRL.\n");
      int m=mrl_between_tacs(yfit, ymeas, dataNr);
      if(m>3 && m>dataNr/3) {
        fprintf(stderr, "Error: bad fit.\n");
        fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
        dftEmpty(&tac); fitEmpty(&fit);
        return(7);
      } else if(m>2 && m>dataNr/4) {
        fprintf(stderr, "Warning: bad fit.\n");
        fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
      } else if(verbose>0) {
        printf("MRL test passed.\n");
        if(verbose>1) fprintf(stderr, "MRL := %d / %d\n", m, dataNr);
      }
    }
    /* Force parameter values inside limits, as they are in the simulation */ 
    double a;
    modelCheckParameters(parNr, pmin, pmax, p, fit.voi[0].p, &a);
    fit.voi[0].wss=wss/a; // remove penalty from reported WSS
    for(int i=0; i<parNr; i++) fit.voi[0].p[i]=p[i];
  }
 
  /*
   *  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(&tac); 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");
    int n;
    RES res; resInit(&res);
    char buf[256];
    ret=fitToResult(&fit, &res, buf);
    if(ret!=0) {
      fprintf(stderr, "Error in making results: %s\n", buf);
      dftEmpty(&tac); fitEmpty(&fit); resEmpty(&res); return(12);
    }
    n=0; strcpy(res.parname[n++], "Function");
    if(1 /*type==331*/) {
      strcpy(res.parname[n++], "Ta");
      strcpy(res.parname[n++], "Ti");
      for(int i=1; i<=sumn; i++) {
        sprintf(res.parname[n++], "A%d", i);
        sprintf(res.parname[n++], "L%d", i);
      }
    }
    strcpy(res.parname[n++], "WSS");
    strcpy(res.program, fit.program);
    sprintf(res.datarange, "%g - %g", tstart, tstop);
    strcpy(res.studynr, tac.studynr);
    if(verbose>1) 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(&tac); fitEmpty(&fit); resEmpty(&res); return(11);
    }
    resEmpty(&res);
    if(verbose>1) printf("Function parameters written in %s\n", resfile);
  }
 
 
  
 
  /*
   *  Plotting of fitted TACs
   */
  if(svgfile[0] || svgfile1[0] || svgfile2[0]) {
 
    if(verbose>1)
      printf("calculating fitted curve at automatically generated sample times\n");
    DFT adft; dftInit(&adft);
    ret=dftAutointerpolate(&tac, &adft, 1.02*tac.x[tac.frameNr-1], verbose-10);
    if(ret) {
      fprintf(stderr, "Error %d in memory allocation for fitted curves.\n", ret);
      dftEmpty(&tac); fitEmpty(&fit); dftEmpty(&adft);
      return(13);
    }
    double a;
    for(int i=0; i<adft.frameNr; i++) {
      ret=fitEval(&fit.voi[0], adft.x[i], &a);
      if(ret!=0) break;
      adft.voi[0].y[i]=a;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted curve for '%s'.\n", svgfile);
      dftEmpty(&tac); dftEmpty(&adft); fitEmpty(&fit);
      return(14);
    }
    /* Add dispersion if required */
    if(tau1>0.0 || tau2>0.0) {
      ret=simDispersion(adft.x, adft.voi[0].y, adft.frameNr, tau1, tau2, adft.voi[0].y3);
      if(ret!=0) {
        fprintf(stderr, "Error: cannot add dispersion to fitted curve for '%s'.\n", svgfile);
        dftEmpty(&tac); dftEmpty(&adft); fitEmpty(&fit);
        return(14);
      }
    }
    /* Save SVG plot of fitted and original data */
    char tmp[256];
    tpcProgramName(argv[0], 0, 0, tmp, 128); strcat(tmp, " ");
    if(strlen(adft.studynr)>1) strcat(tmp, adft.studynr);
    if(svgfile[0]) {
      if(verbose>1) printf("writing %s\n", svgfile);
      ret=plot_fitrange_svg(&tac, &adft, tmp, 0.0, nan(""), 0.0, nan(""), svgfile, verbose-10);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
        dftEmpty(&adft); dftEmpty(&tac); fitEmpty(&fit);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile);
    }
    /* Save SVG plot of first part of fitted and original data */
    if(svgfile1[0]) {
      if(verbose>1) printf("writing %s\n", svgfile1);
      ret=plot_fitrange_svg(&tac, &adft, tmp, 0.0, 2.*tmax, 0.0, nan(""), svgfile1, verbose-10);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile1);
        dftEmpty(&adft); dftEmpty(&tac); fitEmpty(&fit);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile1);
    }
    /* Save SVG plot of lower part of fitted and original data */
    if(svgfile2[0]) {
      if(verbose>1) printf("writing %s\n", svgfile2);
      ret=plot_fitrange_svg(&tac, &adft, tmp, 0.8*tmax+0.2*tac.x[tac.frameNr-1],
                            nan(""), 0.0, nan(""), svgfile2, verbose-10);
      if(ret) {
        fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile2);
        dftEmpty(&adft); dftEmpty(&tac); fitEmpty(&fit);
        return(30+ret);
      }
      if(verbose>0) printf("Plots written in %s\n", svgfile2);
    }
    dftEmpty(&adft);
  }
 
 
 
  /*
   *  Save fitted TAC with original sample times;
   *  this destroys the original sample values
   */
  if(fitfile[0] || dcfile[0]) {
    if(verbose>1) printf("writing fitted TAC(s)\n");
 
    double a;
    for(int i=0; i<tac.frameNr; i++) {
      ret=fitEval(&fit.voi[0], tac.x[i], &a);
      if(ret!=0) break;
      tac.voi[0].y[i]=a;
    }
    if(ret!=0) {
      fprintf(stderr, "Error: cannot calculate fitted curve for '%s'.\n", fitfile);
      dftEmpty(&tac); fitEmpty(&fit);
      return(17);
    }
 
    if(dcfile[0]) {
      if(verbose>1) printf("writing %s\n", dcfile);
      ret=dftWrite(&tac, dcfile);
      if(ret) {
        fprintf(stderr, "Error: cannot write '%s'.\n", dcfile);
        dftEmpty(&tac); fitEmpty(&fit);
        return(18);
      }
    }
 
    if(fitfile[0]) {
      if(tau1>0.0 || tau2>0.0) {
        if(verbose>1) printf("  adding dispersion...\n");
        ret=simDispersion(tac.x, tac.voi[0].y, tac.frameNr, tau1, tau2, tac.voi[0].y3);
        if(ret!=0) {
          fprintf(stderr, "Error: cannot add dispersion to fitted curve for '%s'.\n", fitfile);
            dftEmpty(&tac); fitEmpty(&fit);
            return(17);
        }
      }
      if(verbose>1) printf("writing %s\n", fitfile);
      ret=dftWrite(&tac, fitfile);
      if(ret) {
        fprintf(stderr, "Error: cannot write '%s'.\n", fitfile);
        dftEmpty(&tac); fitEmpty(&fit);
        return(18);
      }
    }
  }
 
  dftEmpty(&tac); fitEmpty(&fit);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
 
double func331(int parNr, double *p, void *fdata)
{
  int fi, n, i;
  double d, wss=0.0, f, t, e1, e2;
  double Ta, Ti, one_per_Ti, A[MAX_FUNC_NR], L[MAX_FUNC_NR], AL[MAX_FUNC_NR];
  double pa[MAX_PARAMETERS], penalty=1.0;
  if(fdata) {/*prevent compiler warning*/}
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
 
  /* Interpretation of parameters */
  n=(parNr-2)/2;
  Ta=pa[0];
  Ti=pa[1]; if(Ti>0.0) one_per_Ti=1.0/Ti; else one_per_Ti=1.0;
  for(i=0; i<n; i++) {
    A[i]=pa[2*i+2];
    L[i]=pa[2*i+3];
  }
  for(i=0; i<n; i++) if(L[i]>1.0E-12) AL[i]=A[i]/L[i]; else AL[i]=A[i];
 
  /* Compute function value at each sample time */
  for(fi=0; fi<dataNr; fi++) {
    t=x[fi];
    if(t<=Ta) {
      yfit[fi]=0.0;
    } else if(t<Ta+Ti) {
      for(i=0, f=0.0; i<n; i++) {
        e2=-L[i]*(t-Ta);
        if(e2>-0.000001) e2=1.0; else if(e2<-50) e2=0.0; else e2=exp(e2);
        f+=AL[i]*(1.0-e2);
      }
      if(Ti>0.0) f*=one_per_Ti;
      yfit[fi]=f;
    } else {
      for(i=0, f=0.0; i<n; i++) {
        e1=-L[i]*(t-Ta-Ti);
        if(e1>-0.000001) e1=1.0; else if(e1<-50) e1=0.0; else e1=exp(e1);
        e2=-L[i]*(t-Ta);
        if(e2>-0.000001) e2=1.0; else if(e2<-50) e2=0.0; else e2=exp(e2);
        f+=AL[i]*(e1-e2);
      }
      if(Ti>0.0) f*=one_per_Ti;
      yfit[fi]=f;
    }
  }
 
  /* Add dispersion if required */
  if(tau1>0.0 || tau2>0.0)
    (void)simDispersion(x, yfit, dataNr, tau1, tau2, ytmp);
 
  /* Compute weighted SS */
  for(fi=0; fi<dataNr; fi++) {
    d=ymeas[fi]-yfit[fi]; if(isnan(ymeas[fi])) continue;
    //if(min_meas==MINMEAS_SS) d*=d; else d=fabs(d);
    wss+=w[fi]*d*d;
    //printf("yfit=%g ymeas=%g d=%g wss=%g\n", yfit[fi], ymeas[fi], d, wss);
  }
  wss*=penalty;
  //printf(" -> %g\n", wss);
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/