fit_blddr.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 "tpctacmod.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
#include "tpcfileutil.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_ipop(int parNr, double *p, void*);
double func_ifengm2(int parNr, double *p, void *);
double func_ifengm2s(int parNr, double *p, void *);
/*****************************************************************************/
 
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *ymeas;
  double       *ysim;
  double       *w;
  int          verbose;
} FITDATA;
const int maxParNr=7;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of AIF integral function to bladder time-activity",
  "curves (TACs). Bladder VOI must include all of the urine in the bladder when",
  "it is largest in the last time frame.",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -model=<M2S | M2 | LATEFDG>",
  "     Select the function; default is M2S.",
  "     M2S and M2 are the functions proposed by Feng et al (1993), including",
  "     two and three exponential terms, respectively.",
  "     LATEFDG is a simplified version, proposed by Phillips et al (1995),",
  "     which can only fit the late part of FDG input TAC, but is designed to",
  "     have two of its parameters fixed to population means.",
  " -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.",
  " -delay=<time>",
  "     Delay time (dT) is constrained to specified value (>=0); by default",
  "     delay time is fitted.",
  " -start=<time>",
  "     Samples before specified time are set to zero for the fit.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -m=<value>",
  "     Parameter m of LATEFDG function is constrained to given value.",
  " -n=<value>",
  "     Parameter n of LATEFDG function is constrained to given value.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: fit_sinf, fit_feng, fit2dat",
  " ",
  "Keywords: TAC, input, 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];
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  unsigned int model=0;
  double   fixed_delay=nan("");
  double   fixed_m=nan("");
  double   fixed_n=nan("");
  double   start_t=nan("");
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=svgfile[0]=(char)0;
  model=modelCodeIndex("feng2ms");
  /* 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(strncasecmp(cptr, "MODEL=", 6)==0) {
      cptr+=6;
      if(strcasecmp(cptr, "M2")==0) {model=modelCodeIndex("fengm2"); continue;}
      if(strcasecmp(cptr, "FENGM2")==0) {model=modelCodeIndex("fengm2"); continue;}
      if(strcasecmp(cptr, "M2S")==0) {model=modelCodeIndex("fengm2s"); continue;}
      if(strcasecmp(cptr, "FENG")==0) {model=modelCodeIndex("fengm2s"); continue;}
      if(strcasecmp(cptr, "LATEFDG")==0) {model=modelCodeIndex("surgefdgaif"); continue;}
      if(strcasecmp(cptr, "surgefdgaif")==0) {model=modelCodeIndex("surgefdgaif"); continue;}
      fprintf(stderr, "Error: invalid model selection.\n"); return(1);
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &fixed_delay) && fixed_delay>=0.0) continue;
    } else if(strncasecmp(cptr, "START=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &start_t) && start_t>=0.0) continue;
    } else if(strncasecmp(cptr, "DT=", 3)==0 && strlen(cptr)>3) {
      if(!atofCheck(cptr+3, &fixed_delay) && fixed_delay>=0.0) continue;
    } else if(strncasecmp(cptr, "M=", 2)==0 && strlen(cptr)>2) {
      if(!atofCheck(cptr+2, &fixed_m) && fixed_m>0.0) continue;
    } else if(strncasecmp(cptr, "N=", 2)==0 && strlen(cptr)>2) {
      if(!atofCheck(cptr+2, &fixed_n) && fixed_n>0.0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, 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);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    if(tacfile[0]) printf("tacfile := %s\n", tacfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(!isnan(fixed_m)) printf("fixed_m := %g\n", fixed_m);
    if(!isnan(fixed_n)) printf("fixed_n := %g\n", fixed_n);
    if(!isnan(fixed_delay)) printf("fixed_delay := %g\n", fixed_delay);
    if(!isnan(start_t)) printf("start_t := %g\n", start_t);
    printf("model := %s\n", modelCode(model));
    printf("weights := %d\n", weights);
  }
 
  /* Set parameter names and initial constraints */
  PAR plim; parInit(&plim); if(parAllocate(&plim, maxParNr, 1)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot initiate parameter list.\n"); return(1);}
  if(model==modelCodeIndex("surgefdgaif")) {
    plim.parNr=5;
    strcpy(plim.n[0].name, "dT");  plim.n[0].lim1=0.0;      plim.n[0].lim2=100.0;
    strcpy(plim.n[1].name, "A");   plim.n[1].lim1=1.0E-08;  plim.n[1].lim2=0.5;
    strcpy(plim.n[2].name, "B");   plim.n[2].lim1=1.0E-06;  plim.n[2].lim2=1.0E+06;
    strcpy(plim.n[3].name, "M");   plim.n[3].lim1=1.0E-06;  plim.n[3].lim2=1.0E+06;
    strcpy(plim.n[4].name, "N");   plim.n[4].lim1=1.0E-08;  plim.n[4].lim2=1.0E+01;
    if(!isnan(fixed_delay)) plim.n[0].lim1=plim.n[0].lim2=fixed_delay;
    if(!isnan(fixed_m)) plim.n[3].lim1=plim.n[3].lim2=fixed_m;
    if(!isnan(fixed_n)) plim.n[4].lim1=plim.n[4].lim2=fixed_n;
  } else if(model==modelCodeIndex("fengm2")) {
    plim.parNr=7;
    strcpy(plim.n[0].name, "dT");  plim.n[0].lim1=0.0;      plim.n[0].lim2=100.0;
    strcpy(plim.n[1].name, "A1");  plim.n[1].lim1=1.0E-06;  plim.n[1].lim2=1.0E+06;
    strcpy(plim.n[2].name, "L1");  plim.n[2].lim1=1.0E-02;  plim.n[2].lim2=1.0E+00;
    strcpy(plim.n[3].name, "A2");  plim.n[3].lim1=1.0E-06;  plim.n[3].lim2=1.0E+06;
    strcpy(plim.n[4].name, "L2");  plim.n[4].lim1=1.0E-03;  plim.n[4].lim2=1.0E-01;
    strcpy(plim.n[3].name, "A3");  plim.n[5].lim1=1.0E-06;  plim.n[5].lim2=1.0E+06;
    strcpy(plim.n[4].name, "L3");  plim.n[6].lim1=1.0E-08;  plim.n[6].lim2=1.0E-02;
  } else if(model==modelCodeIndex("fengm2s")) {
    plim.parNr=5;
    strcpy(plim.n[0].name, "dT");  plim.n[0].lim1=0.0;      plim.n[0].lim2=100.0;
    strcpy(plim.n[1].name, "A1");  plim.n[1].lim1=1.0E-06;  plim.n[1].lim2=1.0E+06;
    strcpy(plim.n[2].name, "L1");  plim.n[2].lim1=1.0E-02;  plim.n[2].lim2=1.0E+00;
    strcpy(plim.n[3].name, "A2");  plim.n[3].lim1=1.0E-06;  plim.n[3].lim2=1.0E+06;
    strcpy(plim.n[4].name, "L2");  plim.n[4].lim1=1.0E-06;  plim.n[4].lim2=1.0E-01;
  } else {
    fprintf(stderr, "Error: invalid model.\n");
    parFree(&plim); return(1);
  }
 
 
  /*
   *  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); parFree(&plim); 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); parFree(&plim); return(2);
  }
  /* Delete missing samples */
  tacDeleteMissingSamples(&tac);
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples for fit.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
  /* Refine dT limits */
  if(plim.n[0].lim2>plim.n[0].lim1 && plim.n[0].lim2>0.25*xmax) plim.n[0].lim2=0.25*xmax;
 
  /* If start time was given, then set the earlier sample values to zero */
  if(!isnan(start_t)) {
    for(int i=0; i<tac.sampleNr; i++) if(tac.x[i]<start_t) 
      for(int j=0; j<tac.tacNr; j++) tac.c[j].y[i]=0.0;
  }
 
  /* 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) {
    if(tacWByFreq(&tac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); return(2);
    }
  }
 
 
  /*
   *  Copy function parameters into PAR structure for printing and saving
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, maxParNr, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); parFree(&plim); return(3);
  }
  par.tacNr=tac.tacNr; par.parNr=modelParNr(model);
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<tac.tacNr; i++) {
    par.r[i].model=model;
    par.r[i].dataNr=tacWSampleNr(&tac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  for(int i=0; i<par.parNr; i++) strcpy(par.n[i].name, plim.n[i].name);
  /* set time and program name */
  {
    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);
  }
  /* set file names */
  iftPut(&par.h, "datafile", tacfile, 0, NULL);
 
 
  /*
   *  Fitting
   */
  if(verbose>1) printf("preparing for fitting\n");
 
  drandSeed(1);
 
  if(verbose>0) {printf("fitting...\n"); fflush(stdout);}
  for(int ci=0; ci<tac.tacNr; ci++) {
    if(verbose>1) printf("TAC %d: %s\n", 1+ci, tac.c[ci].name);
    /* Get y range */
    double ymax; int ymaxi;
    if(tacYRange(&tac, ci, NULL, &ymax, NULL, &ymaxi, NULL, NULL)!=0) {
      fprintf(stderr, "Error: invalid y (concentration) values.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); return(5);
    }
    if(verbose>3) printf("  ymax := %g\n", ymax);
    if(model==modelCodeIndex("surgefdgaif")) {
      /* Refine limits of b and m */
      if(plim.n[2].lim2>plim.n[2].lim1) plim.n[2].lim2=2.0*ymax;
      if(plim.n[3].lim2>plim.n[3].lim1) plim.n[3].lim2=2.0*ymax;
    } else if(model==modelCodeIndex("fengm2")) {
      /* Refine limits of a */
      if(plim.n[1].lim2>plim.n[1].lim1) plim.n[1].lim2=ymax/2.2;
      if(plim.n[3].lim2>plim.n[3].lim1) plim.n[3].lim2=0.2*ymax/2.2;
      if(plim.n[5].lim2>plim.n[5].lim1) plim.n[5].lim2=0.2*ymax/2.2;
    } else if(model==modelCodeIndex("fengm2s")) {
      /* Refine limits of a */
      if(plim.n[1].lim2>plim.n[1].lim1) plim.n[1].lim2=ymax/2.2;
      if(plim.n[3].lim2>plim.n[3].lim1) plim.n[3].lim2=0.2*ymax/2.2;
    }
    /* Set data pointers for the fit */
    double yfit[tac.sampleNr];
    FITDATA fitdata;
    fitdata.n=tac.sampleNr;
    fitdata.x=tac.x;
    fitdata.ymeas=tac.c[ci].y;
    fitdata.ysim=yfit;
    fitdata.w=tac.w;
    if(verbose>10) fitdata.verbose=verbose-10; else fitdata.verbose=0;
    /* Set NLLS options */
    NLOPT nlo; nloptInit(&nlo);
    if(nloptAllocate(&nlo, par.parNr)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); return(5);
    }
    if(model==modelCodeIndex("surgefdgaif")) nlo._fun=func_ipop;
    else if(model==modelCodeIndex("fengm2")) nlo._fun=func_ifengm2;
    else if(model==modelCodeIndex("fengm2s")) nlo._fun=func_ifengm2s;
    nlo.fundata=&fitdata;
    nlo.totalNr=par.parNr;
    for(int i=0; i<par.parNr; i++) {
      nlo.xlower[i]=plim.n[i].lim1;
      nlo.xupper[i]=plim.n[i].lim2;
    }
    for(int i=0; i<par.parNr; i++) {
      if(!(plim.n[i].lim2>plim.n[i].lim1)) nlo.xtol[i]=0.0;
      else nlo.xtol[i]=0.0001*(plim.n[i].lim2-plim.n[i].lim1);
    }
    /* Initial guess */
    if(model==modelCodeIndex("surgefdgaif")) {
      nlo.xfull[0]=0.8*nlo.xlower[0]+0.2*nlo.xupper[0];
      nlo.xfull[1]=0.7*nlo.xlower[1]+0.3*nlo.xupper[1];
      nlo.xfull[2]=0.6*nlo.xlower[2]+0.4*nlo.xupper[2];
      nlo.xfull[3]=0.6*nlo.xlower[3]+0.4*nlo.xupper[3];
      nlo.xfull[4]=0.8*nlo.xlower[4]+0.2*nlo.xupper[4];
    } else if(model==modelCodeIndex("fengm2")) {
      nlo.xfull[0]=0.8*nlo.xlower[0]+0.2*nlo.xupper[0];
      nlo.xfull[1]=0.5*ymax/6.6;
      nlo.xfull[2]=0.5;
      nlo.xfull[3]=0.5*ymax/6.6;
      nlo.xfull[4]=0.05;
      nlo.xfull[5]=0.02*ymax/6.6;
      nlo.xfull[6]=0.005;
    } else if(model==modelCodeIndex("fengm2s")) {
      nlo.xfull[0]=0.8*nlo.xlower[0]+0.2*nlo.xupper[0];
      nlo.xfull[1]=ymax/6.6;
      nlo.xfull[2]=0.5;
      nlo.xfull[3]=2.0*nlo.xfull[1]/15.0;
      nlo.xfull[4]=0.05;
    }
    nlo.maxFunCalls=50000;
    /* Fit */
    if(verbose>4) nloptWrite(&nlo, stdout);
    int orig_verbose=status.verbose; if(verbose>2) status.verbose=1; else status.verbose=0;
#if(0)
    if(nloptSimplexARRS(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); nloptFree(&nlo); return(6);
    }
#else
    if(nloptIATGO(&nlo, 1, 0, 0.25, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); nloptFree(&nlo); return(6);
    }
#endif
    status.verbose=orig_verbose;
    nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>2) nloptWrite(&nlo, stdout);
    if(verbose>3) {
      printf("measured and fitted TAC:\n");
      for(int i=0; i<tac.sampleNr; i++)
        printf("\t%g\t%g\t%g\n", tac.x[i], tac.c[ci].y[i], yfit[i]);
    }
    //if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
    /* Copy parameters */
    par.r[ci].wss=nlo.funval;
    if(model==modelCodeIndex("surgefdgaif")) {
      for(int i=0; i<par.parNr; i++) par.r[ci].p[i]=nlo.xfull[i];
    } else if(model==modelCodeIndex("fengm2")) {
      par.r[ci].p[0]=nlo.xfull[1];
      par.r[ci].p[1]=-nlo.xfull[2];
      par.r[ci].p[2]=nlo.xfull[3];
      par.r[ci].p[3]=-nlo.xfull[4];
      par.r[ci].p[4]=nlo.xfull[5];
      par.r[ci].p[5]=-nlo.xfull[6];
      par.r[ci].p[6]=nlo.xfull[0];
    } else if(model==modelCodeIndex("fengm2s")) {
      par.r[ci].p[0]=nlo.xfull[1];
      par.r[ci].p[1]=-nlo.xfull[2];
      par.r[ci].p[2]=nlo.xfull[3];
      par.r[ci].p[3]=-nlo.xfull[4];
      par.r[ci].p[4]=nlo.xfull[0];
    }
    nloptFree(&nlo);
  } // next TAC
 
  /* Print and save parameters */
  if(verbose>0 || !parfile[0]) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  if(parfile[0]) {
    /* 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); parFree(&plim); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_FIT, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); 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"); fflush(stdout);}
    /* Prepare place for function TAC */
    double sdist;
    tacGetSampleInterval(&tac, 1.0E-03, &sdist, NULL);
    int snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end-par.r[0].start, 2, NULL);
    TAC ftac; tacInit(&ftac);
    if(verbose>3) {printf("  allocating memory for %d samples\n", snr); fflush(stdout);}
    int ret=tacAllocate(&ftac, snr, tac.tacNr);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(ret));
      tacFree(&tac); parFree(&plim); parFree(&par); return(21);
    }
    ftac.tacNr=tac.tacNr; ftac.sampleNr=snr;
    if(verbose>3) {printf("  copying TAC headers\n"); fflush(stdout);}
    ret=tacCopyHdr(&tac, &ftac);
    if(ret==TPCERROR_OK) for(int ci=0; ci<tac.tacNr && ret==0; ci++)
      ret=tacCopyTacchdr(tac.c+ci, ftac.c+ci);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(ret));
      tacFree(&tac); parFree(&plim); 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; note that we are computing integral this time */
    ret=0;
    if(verbose>3) {printf("  computing function(s)\n"); fflush(stdout);}
    for(int ci=0; ci<tac.tacNr && ret==0; ci++)
      ret=mfEvalInt(modelCode(par.r[ci].model), par.parNr, par.r[ci].p, 
                  ftac.sampleNr, ftac.x, ftac.c[ci].y, verbose-5);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate function values.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(23);
    }
    if(verbose>10) tacWrite(&ftac, stdout, TAC_FORMAT_UNKNOWN, 1, NULL);
    /* Plot */
    if(verbose>3) {printf("  plotting function(s)\n"); fflush(stdout);}
//    ret=tacPlotFitSVG(&tac, &ftac, "", nan(""), nan(""), nan(""), nan(""), svgfile, NULL);
    ret=tacPlotFitSVG(&tac, &ftac, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(ret));
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(24);
    }
    tacFree(&ftac);
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
  }
 
 
  tacFree(&tac); 
  parFree(&plim);
  parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_ipop(int parNr, double *p, void *fdata)
{
  double v, wss=0.0;
  FITDATA *d=(FITDATA*)fdata;
  double a=p[1];
  double b=p[2];
  double m=p[3];
  double n=p[4];
 
  /* Calculate the curve values and weighted SS */
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=0.0;
    if(isnan(d->ymeas[i]) || isnan(d->x[i])) continue;
    double x=d->x[i]-p[0];
    if(x>0.0) {
      d->ysim[i] = m*(1.0 - (n*a*x + 1.0)*exp(-n*a*x))/(n*n*a*a);
      d->ysim[i] += (1.0 - exp(-a*x))/a;
      d->ysim[i] *= b;
    }
    v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
  if(d->verbose>0) {
    for(unsigned int i=0; i<(unsigned int)parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
  return(wss);
}
/*****************************************************************************/
double func_ifengm2(int parNr, double *p, void *fdata)
{
  double v, wss=0.0;
  FITDATA *d=(FITDATA*)fdata;
  double A1=p[1];
  double L1=p[2];
  double A2=p[3];
  double L2=p[4];
  double A3=p[5];
  double L3=p[6];
 
  /* Calculate the curve values and weighted SS */
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=0.0;
    if(isnan(d->ymeas[i]) || isnan(d->x[i])) continue;
    double x=d->x[i]-p[0];
    if(x>0.0) {
      v=0.0;
      if(L1!=0.0) {
        double a=exp(-L1*x);
        v=(A1-L1*(A2+A3))*(1.0-a)/(L1*L1); if(isfinite(v)) d->ysim[i]+=v;
        v=(A1/L1)*x*a; if(isfinite(v)) d->ysim[i]-=v;
      }
      if(L2!=0.0) v=(A2/L2)*(1.0-exp(-L2*x)); else v=A2*x;
      if(isfinite(v)) d->ysim[i]+=v;
      if(L3!=0.0) v=(A3/L3)*(1.0-exp(-L3*x)); else v=A3*x;
      if(isfinite(v)) d->ysim[i]+=v;
    }
    v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
  if(d->verbose>0) {
    for(unsigned int i=0; i<(unsigned int)parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
  return(wss);
}
/*****************************************************************************/
double func_ifengm2s(int parNr, double *p, void *fdata)
{
  double v, wss=0.0;
  FITDATA *d=(FITDATA*)fdata;
  double A1=p[1];
  double L1=p[2];
  double A2=p[3];
  double L2=p[4];
 
  /* Calculate the curve values and weighted SS */
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=0.0;
    if(isnan(d->ymeas[i]) || isnan(d->x[i])) continue;
    double x=d->x[i]-p[0];
    if(x>0.0) {
      v=0.0;
      if(L1!=0.0) {
        double a=exp(-L1*x);
        v=(A1-L1*A2)*(1.0-a)/(L1*L1); if(isfinite(v)) d->ysim[i]+=v;
        v=(A1/L1)*x*a; if(isfinite(v)) d->ysim[i]-=v;
      }
      if(L2!=0.0) v=(A2/L2)*(1.0-exp(-L2*x)); else v=A2*x;
      if(isfinite(v)) d->ysim[i]+=v;
    }
    v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
  if(d->verbose>0) {
    for(unsigned int i=0; i<(unsigned int)parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/