fit_wrlv.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 "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
#include "tpcbfm.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_disp(int parNr, double *p, void*);
double func_wrlv(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int  n;
  double       *x;
  double       *x1;
  double       *x2;
  double       *yrv;
  double       *ylv;
  double       *w;

  double       *fyrv;
  double       *fylv;

  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of exponential input function with paired, repeated",
  "eigenvalues, incorporating injection schedule, to RV and LV cavity BTACs",
  "from a radiowater PET study, where LV BTAC is modelled with simple",
  "delay and dispersion of the RV BTAC.",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -Ti=<infusion time>",
  "     Duration of intravenous infusion in seconds; 0, if short bolus.",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     output data file are used, if available.",
  " -wf",
  "     Weight by output TAC sampling interval.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "TAC file must contain two TACs, for RV and LV, in this order.",
  "Sample times must be in seconds, unless units are specified in the file.",
  " ",
  "See also: fit_disp, fit_sinf, fit_wpul, disp4dft, convexpf",
  " ",
  "Keywords: input, curve fitting, myocardium, radiowater",
  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
  double  fixedTi=nan("");
  int     ret;
 
 
  drandSeed(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, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "Ti=", 3)==0) {
      fixedTi=atofVerified(cptr+3); if(fixedTi>=0.0) 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]) { // note that parameter file is optional
    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(!isnan(fixedTi)) printf("fixedTi := %g\n", fixedTi);
    printf("weights := %d\n", weights);
    fflush(stdout);
  }
 
 
  /*
   *  Read the data
   */
  if(verbose>1) printf("reading BTAC data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  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<2) {
    fprintf(stderr, "Error: file does not contain RV and LV TACs.\n");
    tacFree(&tac); return(2);
  } else if(tac.tacNr>2) {
    if(verbose>0) fprintf(stderr, "Warning: file contains more than two TACs.\n");
    tac.tacNr=2;
  }
  if(tac.sampleNr<8) {
    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);
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    printf("orig_xmin := %g\n", xmin);
    printf("orig_xmax := %g\n", xmax);
  }
  if(!(xmin>=0.0)) {
    fprintf(stderr, "Error: negative sample times.\n");
    tacFree(&tac); return(2);
  }
  /* Convert sample times into seconds */
  if(tac.tunit==UNIT_UNKNOWN) {
    if(xmax<10.0) tac.tunit=UNIT_MIN; else tac.tunit=UNIT_SEC;
  }
  if(tacXUnitConvert(&tac, UNIT_SEC, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
    printf("final_xunit := %s\n", unitName(tac.tunit));
  }
 
  /* Find the sample time and index of max RV concentration */
  double ymin, ymax; int imax;
  ret=tacYRange(&tac, 0, &ymin, &ymax, NULL, &imax, NULL, NULL);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
    tacFree(&tac); return(2);
  }
  double tmax=tac.x[imax];
  if(verbose>1) {
    printf("peak_x := %g\n", tmax);
    printf("peak_y := %g\n", ymax);
    printf("peak_index := %d\n", imax);
  }
  if((tac.sampleNr-imax)<3 || !(ymin<0.5*ymax)) {
    fprintf(stderr, "Error: invalid TAC shape in %s\n", tacfile);
    tacFree(&tac); return(2);
  }
 
  /* Find the sample time of max LV concentration */
  double lv_tmax;
  {
    int imax;
    ret=tacYRange(&tac, 1, NULL, NULL, NULL, &imax, NULL, NULL);
    if(ret!=0) {
      fprintf(stderr, "Error: invalid contents in %s\n", tacfile);
      tacFree(&tac); return(2);
    }
    lv_tmax=tac.x[imax];
    if(verbose>1) {
      printf("LV_peak_x := %g\n", lv_tmax);
    }
  }
 
  /* Set places for fitted TACs */
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(tacAllocateMore(&tac, 2)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
 
  /* Check and set weights */
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  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); return(2);
    }
  }
 
 
  /*
   *  Fit delay and dispersion between RV and LV TACs
   */
  double delay, tau;
  {
    if(verbose>1) printf("fitting dispersion and delay between RV and LV TACs\n");
    FITDATA fitdata;
    /* Set data pointers for the fit */
    fitdata.n=tac.sampleNr;
    fitdata.x=tac.x;
    fitdata.yrv=tac.c[0].y;
    fitdata.ylv=tac.c[1].y;
    fitdata.fyrv=tac.c[2].y;
    fitdata.fylv=tac.c[3].y;
    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, 2)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n");
      tacFree(&tac); return(3);
    }
    nlo._fun=func_disp;
    nlo.fundata=&fitdata;
    nlo.totalNr=2;
    /* Set initial values and limits */
    double peakdif=lv_tmax-tmax;
    if(!(peakdif>0.0)) peakdif=0.0;
    if(peakdif>50.0) peakdif=50.0;
    nlo.xlower[0]=0.0; nlo.xupper[0]=60.0; nlo.xfull[0]=10.0; nlo.xdelta[0]=1.0; nlo.xtol[0]=0.001; 
    nlo.xlower[1]=0.0; nlo.xupper[1]=60.0; nlo.xfull[1]=peakdif; nlo.xdelta[1]=0.5; nlo.xtol[1]=0.001; 
 
    if(verbose>2) printf("non-linear optimization\n");
    if(nloptSimplex(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); nloptFree(&nlo); return(3);
    }
    if(verbose>3) {
      nloptWrite(&nlo, stdout);
      (void)func_disp(nlo.totalNr, nlo.xfull, nlo.fundata);
      double wss=0.0;
      for(unsigned int i=0; i<fitdata.n; i++) {
        double v=fitdata.fylv[i]-fitdata.ylv[i];
        wss+=fitdata.w[i]*v*v;
      }
      printf(" wss := %g\n", wss);
    }
    tau=nlo.xfull[0];
    delay=nlo.xfull[1];
    if(verbose>1) {
      printf(" fitted_tau := %g\n", tau);
      printf(" fitted_delay := %g\n", delay);
    }
    nloptFree(&nlo);
  }
 
 
  /*
   *  Get initial estimates for radiowater appearance time based on highest RV TAC slope
   */
  if(verbose>1) printf("slope method for estimating appearance time\n");
  double Ta=nan("");
  {
    double slope;
    if(highestSlope(tac.x, tac.c[0].y, tac.sampleNr, 3, 0.0, &slope, NULL, &Ta, NULL) || !(slope>0.0))
    {
      fprintf(stderr, "Error: invalid TAC shapes.\n");
      tacFree(&tac); return(2);
    }
    if(verbose>1) printf("  initial_Ta := %g\n", Ta);
  }
 
 
  /*
   *  Get initial estimates for exponential function parameters from the descending part
   *  of the RV BTAC
   */
  if(verbose>1) printf("initial exponential fitting\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  int enr=2;
  double init_a[enr], init_k[enr];
  /* Get the shortest frame length or sampling interval */
  double sdist;
  if(tacGetSampleInterval(&tac, 1.0E-03, &sdist, NULL)!=0) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>1) printf("shortest frame length := %g\n", sdist);
  {
    /* Get the data after the peak */
    int sampleNr=tac.sampleNr-imax-1;
    double *y=tac.c[0].y+imax+1;
    double *w=tac.w+imax+1;
    double x[sampleNr]; for(int i=0; i<sampleNr; i++) x[i]=tac.x[i+imax+1]-tmax;
    if(verbose>3) {
      printf("\n\tRV data for decaying exponential fitting\n");
      for(int i=0; i<sampleNr; i++) printf("\t%g\t%g\t%g\n", x[i], y[i], w[i]);
      printf("\n"); fflush(stdout);
    }
    /* Spectral analysis to find better range of parameters */
    if(verbose>4) printf("estimating initial k value range\n");
    int bfNr=100;
    double kmin, kmax; // kmin must be >0
    double xmin, xmax; doubleRange(x, sampleNr, &xmin, &xmax);
    kmin=1.0E-06/xmax; kmax=3.0/(xmin+sdist); if(kmax<100.0*kmin) kmax=100.0*kmin;
/*
    kmin=0.000001/sdist; kmax=10.0/sdist; if(kmax>100.) kmax=100.;
*/
    if(verbose>2) printf(" kmin := %g\n kmax := %g\n", kmin, kmax);
    double pa[bfNr], pk[bfNr], yfit[sampleNr];
    ret=spectralDExp(x, y, w, 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>5) {
      printf("\n\tSolutions for each k:\n");
      for(int bi=0; bi<bfNr; bi++) printf("\t%g\t%g\n", pk[bi], pa[bi]);
    }
    /* Set better k value range, if necessary */
    if(!(pa[0]>0.0 && pa[bfNr-1]>0.0)) {
      if(verbose>4) 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);
      }
      ret=spectralDExp(x, y, w, 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>5) {
        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("\n\tMeasured and fitted TAC:\n");
      for(int i=0; i<sampleNr; i++) printf("\t%g\t%g\n", y[i], yfit[i]);
      printf("\n"); fflush(stdout);
    }
    /* Copy the parameter clusters, combining some of those if necessary */
    ret=spectralBFExtract(pk, pa, bfNr, init_k, init_a, enr);
    if(ret==0) {
      fprintf(stderr, "Error: cannot estimate initial parameters.\n");
      tacFree(&tac); return(3);
    }
    if(verbose>3) {
      printf("\n\tInitial a and k values\n");
      for(int i=0; i<enr; i++) printf("\t%g\t%g\n", init_a[i], init_k[i]);
      printf("\n"); fflush(stdout);
    }
  }
 
 
 
  /*
   *  Fitting
   */
  if(verbose>1) printf("preparing for fitting\n");
 
  /* Set data pointers for the fit */
  FITDATA fitdata;
  fitdata.verbose=verbose-10;
  fitdata.n=tac.sampleNr;
  fitdata.x=tac.x;
  if(tac.isframe) {fitdata.x1=tac.x1; fitdata.x2=tac.x2;} else {fitdata.x1=NULL; fitdata.x1=NULL;}
  fitdata.yrv=tac.c[0].y;
  fitdata.ylv=tac.c[1].y;
  fitdata.fyrv=tac.c[2].y;
  fitdata.fylv=tac.c[3].y;
  fitdata.w=tac.w;
 
 
  /* Set parameters for nonlinear optimization */
  if(verbose>2) printf("process initial parameter values\n");
  int parNr=8;
  NLOPT nlo; nloptInit(&nlo);
  ret=nloptAllocate(&nlo, parNr);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  /* Set function */
  nlo._fun=func_wrlv;
  nlo.fundata=&fitdata;
  nlo.totalNr=parNr; // Ta, Ti, A1, L1, A2, L2, dT, tau
  nlo.maxFunCalls=50000;
  { // Appearance time Ta
    nlo.xlower[0]=0.8*Ta; 
    nlo.xupper[0]=1.2*Ta;
    nlo.xfull[0]=Ta;
    nlo.xdelta[0]=0.01*(nlo.xupper[0]-nlo.xlower[0]);
    nlo.xtol[0]=0.02;
  }
  if(fixedTi>=0.0) { // Infusion duration Ti
    nlo.xlower[1]=nlo.xupper[1]=nlo.xfull[1]=fixedTi; 
    nlo.xdelta[1]=0.0;
    nlo.xtol[1]=0.0;
  } else {
    nlo.xlower[1]=1.0;
    nlo.xfull[1]=tmax-Ta; if(nlo.xfull[1]<2.0) nlo.xfull[1]=2.0;
    nlo.xupper[1]=1.5*nlo.xfull[1];
    nlo.xdelta[1]=0.01*(nlo.xupper[1]-nlo.xlower[1]);
    nlo.xtol[1]=0.1*nlo.xdelta[1];
  }
  { // Delay between RV and LV
    nlo.xlower[6]=1.0; 
    nlo.xupper[6]=10.0;
    nlo.xfull[6]=delay; 
    if(nlo.xfull[6]<=nlo.xlower[6] || nlo.xfull[6]>=nlo.xupper[6]) nlo.xfull[6]=0.5*(nlo.xlower[6]+nlo.xupper[6]);
    nlo.xdelta[6]=0.02;
    nlo.xtol[6]=0.001;
  }
  { // Dispersion between RV and LV
    nlo.xlower[7]=1.0; 
    nlo.xupper[7]=15.0;
    nlo.xfull[7]=tau;
    if(nlo.xfull[7]<=nlo.xlower[7] || nlo.xfull[7]>=nlo.xupper[7]) nlo.xfull[7]=0.5*(nlo.xlower[7]+nlo.xupper[7]);
    nlo.xdelta[7]=0.02;
    nlo.xtol[7]=0.001;
  }
  { // Exponent function parameters
    for(int i=0; i<2; i++) {
      int ai=2+2*i;
      int ki=ai+1;
 
      nlo.xfull[ai]=init_a[i];
      nlo.xlower[ai]=0.02*init_a[i];
      nlo.xupper[ai]=10.0*init_a[i];
      nlo.xdelta[ai]=1.0E-02*init_a[i];
      nlo.xtol[ai]=1.0E-04*init_a[i];
 
      nlo.xfull[ki]=init_k[i];
      nlo.xlower[ki]=0.011*init_k[i];
      nlo.xupper[ki]=1.5*init_k[i];
      nlo.xdelta[ki]=1.0E-02*init_k[i]; if(!(nlo.xdelta[ki]>1.0E-05)) nlo.xdelta[ki]=1.0E-05;
      nlo.xtol[ki]=1.0E-06*init_k[i]; if(!(nlo.xtol[ki]>1.0E-06)) nlo.xtol[ki]=1.0E-06;
    }
  }
 
  /* Get nr of fixed parameters */
  unsigned int fixedNr=nloptFixedNr(&nlo);
  if(verbose>2) printf("fixedNr := %u\n", fixedNr); 
 
#if(0)
  // call function for testing
  fitdata.verbose=10;
  double wss=func_wrlv(nlo.totalNr, nlo.xfull, &fitdata);
  if(verbose>1) {
    printf("\nTime\tRV\tLV\tsimRV\tsimLV\n");
    for(unsigned int i=0; i<fitdata.n; i++) printf("%g\t%g\t%g\t%g\t%g\n", fitdata.x[i], 
        fitdata.yrv[i], fitdata.ylv[i], fitdata.fyrv[i], fitdata.fylv[i]);
    printf("\n"); fflush(stdout);
  }
  if(verbose>1) {printf("\twss := %g\n", wss); fflush(stdout);}
 
  nloptFree(&nlo);
  tacFree(&tac);
 
#else
 
  /* Fit */
  nlo.funval=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>0) {printf("fitting\n"); fflush(stdout);}
 
  if(nloptSimplex(&nlo, 50000, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&nlo); return(5);
  }
  nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
  if(verbose>6) {
    printf("\nTime\tRV\tLV\tsimRV\tsimLV\n");
    for(int i=0; i<tac.sampleNr; i++) printf("%g\t%g\t%g\t%g\t%g\n", tac.x[i], 
        tac.c[0].y[i], tac.c[1].y[i], tac.c[2].y[i], tac.c[3].y[i]);
  }
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
  /* Simplex again */
  for(unsigned int i=0; i<nlo.totalNr; i++) {
    nlo.xdelta[i]*=8.97E-02;
    nlo.xtol[i]*=1.0E-02;
  }
  if(nloptSimplex(&nlo, 50000, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&nlo); return(5);
  }
  nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
  if(verbose>6) {
    printf("\nTime\tRV\tLV\tsimRV\tsimLV\n");
    for(int i=0; i<tac.sampleNr; i++) printf("%g\t%g\t%g\t%g\t%g\n", tac.x[i], 
        tac.c[0].y[i], tac.c[1].y[i], tac.c[2].y[i], tac.c[3].y[i]);
  }
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
  /* Simplex once more */
  for(unsigned int i=0; i<nlo.totalNr; i++) {
    nlo.xdelta[i]*=1.23E-01;
    nlo.xtol[i]*=1.5E-02;
  }
  if(nloptSimplex(&nlo, 50000, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&nlo); return(5);
  }
  nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
  if(verbose>6) {
    printf("\nTime\tRV\tLV\tsimRV\tsimLV\n");
    for(int i=0; i<tac.sampleNr; i++) printf("%g\t%g\t%g\t%g\t%g\n", tac.x[i], 
        tac.c[0].y[i], tac.c[1].y[i], tac.c[2].y[i], tac.c[3].y[i]);
  }
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
 
 
  /*
   *  Copy function parameters into PAR structure for printing and saving
   */
  PAR par; parInit(&par);
  ret=parAllocateWithTAC(&par, &tac, nlo.totalNr, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); nloptFree(&nlo); return(10);
  }
  par.tacNr=2; par.parNr=nlo.totalNr;
  par.format=PAR_FORMAT_TSV_UK;
  {
    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);
  }
  iftPut(&par.h, "datafile", tacfile, 0, NULL);
  par.r[0].model=modelCodeIndex("ebolinf");
  par.r[0].dataNr=tacWSampleNr(&tac);
  par.r[0].start=xmin; 
  par.r[0].end=xmax;
  par.r[0].wss=nlo.funval;
  par.r[0].fitNr=nlo.totalNr-fixedNr;
  par.r[1].model=modelCodeIndex("ebolinfdd");
  par.r[1].dataNr=tacWSampleNr(&tac);
  par.r[1].start=xmin; 
  par.r[1].end=xmax;
  par.r[1].wss=nlo.funval;
  par.r[1].fitNr=nlo.totalNr-fixedNr;
  /* Set parameter names and units */
  strcpy(par.n[0].name, "Ta"); par.n[0].unit=tac.tunit;
  strcpy(par.n[1].name, "Ti"); par.n[1].unit=tac.tunit;
  strcpy(par.n[2].name, "A1"); par.n[2].unit=tac.cunit;
  strcpy(par.n[3].name, "L1"); par.n[3].unit=unitInverse(tac.tunit);
  strcpy(par.n[4].name, "A2"); par.n[4].unit=tac.cunit;
  strcpy(par.n[5].name, "L2"); par.n[5].unit=unitInverse(tac.tunit);
  strcpy(par.n[6].name, "dT"); par.n[6].unit=tac.tunit;
  strcpy(par.n[7].name, "tau"); par.n[7].unit=tac.tunit;
  /* Set parameter values in correct order, and write lambda's as positive values */
  par.r[0].p[0]=par.r[1].p[0]=nlo.xfull[0];
  par.r[0].p[1]=par.r[1].p[1]=nlo.xfull[1];
  par.r[0].p[2]=par.r[1].p[2]=nlo.xfull[2];
  par.r[0].p[3]=par.r[1].p[3]=+nlo.xfull[3];
  par.r[0].p[4]=par.r[1].p[4]=nlo.xfull[4];
  par.r[0].p[5]=par.r[1].p[5]=+nlo.xfull[5];
  par.r[0].p[6]=0.0; par.r[1].p[6]=nlo.xfull[6];
  par.r[0].p[7]=0.0; par.r[1].p[7]=nlo.xfull[7];
  /* Print and save */
  if(verbose>0) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  if(parfile[0]) {
    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;
    /* 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(&nlo); 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(&nlo); 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");
    TAC fit; tacInit(&fit);
    (void)tacDuplicate(&tac, &fit);
    for(int i=0; i<tac.sampleNr; i++) {fit.c[0].y[i]=tac.c[2].y[i]; fit.c[1].y[i]=tac.c[3].y[i];}
    /* Plot */
    ret=tacPlotFitSVG(&tac, &fit, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&fit); nloptFree(&nlo); parFree(&par);
      return(21);
    }
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
    tacFree(&fit);
  }
 
  nloptFree(&nlo); parFree(&par);
  tacFree(&tac);
 
#endif
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_disp(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=2) return(nan(""));
  if(d->verbose>20) printf("p[]: %g %g\n", p[0], p[1]);
 
  /* Add dispersion to the input TAC */
  double dtac[2*d->n];
  double *buf=dtac+d->n;
  for(unsigned int i=0; i<d->n; i++) dtac[i]=d->yrv[i];
  if(simDispersion(d->x, dtac, d->n, p[0], 0.0, buf)) return(nan(""));
 
  /* Add delay and interpolate to output sample times */
  for(unsigned int i=0; i<d->n; i++) buf[i]=d->x[i]+p[1];
  if(liInterpolate(buf, dtac, d->n, d->x, d->fylv, NULL, NULL, d->n, 0, 1, 0)!=0) return(nan(""));
 
  /* Calculate the weighted SS */
  double wss=0.0;
  for(unsigned int i=0; i<d->n; i++) {
    double v=d->fylv[i]-d->ylv[i];
    wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
 
/*****************************************************************************/
 
 
double func_wrlv(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=8 || p==NULL || fdata==NULL || d->n<1) return(nan(""));
  if(d->verbose>20) {
    printf("p[]: %g", p[0]);
    for(int i=1; i<parNr; i++) printf(" %g", p[i]);
    printf("\n"); fflush(stdout);
  }
 
  /* Process parameters */
  double Ta=p[0];
  double Ti=p[1];
  double A[3], L[3], ApL[3]; 
  A[0]=p[2]; L[0]=p[3]; A[1]=p[4]; L[1]=p[5]; A[2]=0.0; L[2]=0.0;
  double deltaT=p[6];
  double tau=p[7];
  /* derived parameters */
  double one_per_Ti; if(Ti>0.0) one_per_Ti=1.0/Ti; else one_per_Ti=1.0;
  for(int i=0; i<3; i++) if(L[i]>1.0E-12) ApL[i]=A[i]/L[i]; else ApL[i]=A[i];
 
  /* Compute function value at each sample time */
  for(unsigned int fi=0; fi<d->n; fi++) {
    double t=d->x[fi];
    if(t<=Ta) {
      d->fyrv[fi]=0.0;
    } else if(t<Ta+Ti) {
      double f=0.0;
      for(int i=0; i<3; i++) {
        double e2=-L[i]*(t-Ta);
        if(e2>-0.000001) e2=1.0; else if(e2<-50) e2=0.0; else e2=exp(e2);
        f+=ApL[i]*(1.0-e2);
      }
      if(Ti>0.0) f*=one_per_Ti;
      d->fyrv[fi]=f;
    } else {
      double f=0.0;
      for(int i=0; i<3; i++) {
        double e1=-L[i]*(t-Ta-Ti);
        if(e1>-0.000001) e1=1.0; else if(e1<-50) e1=0.0; else e1=exp(e1);
        double e2=-L[i]*(t-Ta);
        if(e2>-0.000001) e2=1.0; else if(e2<-50) e2=0.0; else e2=exp(e2);
        f+=ApL[i]*(e1-e2);
      }
      if(Ti>0.0) f*=one_per_Ti;
      d->fyrv[fi]=f;
    }
  }
 
  /* Copy RV TAC to LV TAC with delay */
  if(deltaT>0.0) {
    double xt[d->n];
    for(unsigned int fi=0; fi<d->n; fi++) xt[fi]=d->x[fi]+deltaT;
    if(liInterpolate(xt, d->fyrv, d->n, d->x, d->fylv, NULL, NULL, d->n, 0, 1, 0)!=0) return(nan(""));
  } else {
    for(unsigned int fi=0; fi<d->n; fi++) d->fylv[fi]=d->fyrv[fi];
  }
 
  /* Add dispersion to LV TAC if required */
  if(tau>0.0) (void)simDispersion(d->x, d->fylv, d->n, tau, 0.0, NULL);
 
  /* Calculate the weighted SS */
  if(d->verbose>2) {fprintf(stdout, "computing WSS...\n"); fflush(stdout);}
  double wss=0.0;
  for(unsigned i=0; i<d->n; i++) {
    double v1=d->fyrv[i] - d->yrv[i];
    double v2=d->fylv[i] - d->ylv[i];
    wss+=d->w[i]*(v1*v1 + v2*v2);
  }
  if(d->verbose>1) {
    for(int i=0; i<parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss); fflush(stdout);
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/