fit_wliv.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#ifdef HAVE_OMP_H
#include <omp.h>
#endif
/*****************************************************************************/
#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_wliv(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int  ni;
  double       *xi;
  double       *yi;

  unsigned int  nt;
  double       *xt1;
  double       *xt2;
  double       *xt;
  double       *yt;
  double       *w;
  double       *syt;

  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the liver radiowater model to TTACs using arterial",
  "BTAC as the input function. Small delay of arterial input function is",
  "fitted by default, as well as larger delay and dispersion representing",
  "the portal vein input from GI system.",
  " ",
  "Usage: @P [Options] btacfile ttacfile [parfile]",
  " ",
  "Options:",
  " -delay=<value>",
  "     Delay time (sec) between BTAC and liver is constrained to given value.",
  " -delayp=<value>",
  "     Delay time (sec) in GI system is constrained to given value.",
  " -kGI=<value>",
  "     Dispersion in GI system is constrained to given rate constant (1/min);",
  "     dispersion time constant tau=1/kGI.",
  " -Vb=<value>",
  "     Blood volume is constrained to given value (mL/mL of liver).",
  " -w1 | -wf",
  "     All weights are set to 1.0 (no weighting), or based on TTAC frame",
  "     lengths; by default, weights in TTAC file are used, if available.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -sim=<Filename>",
  "     Fitted TTACs at BTAC sample times are saved in specified TAC file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Sample times must be in seconds, unless units are specified in the file.",
  " ",
  "See also: fit_h2o, bfmh2o, fitk2, sim_wliv, simdisp",
  " ",
  "Keywords: TAC, liver, modelling, perfusion, 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    btacfile[FILENAME_MAX], ttacfile[FILENAME_MAX], parfile[FILENAME_MAX], 
          svgfile[FILENAME_MAX], simfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  double  fixed_kgi=nan("");
  double  fixed_vb=nan("");
  double  fixed_dt_liver=nan("");
  double  fixed_dt_gi=nan("");
 
  drandSeed(6789); //drandSeed(time(NULL));
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  btacfile[0]=ttacfile[0]=parfile[0]=svgfile[0]=simfile[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, "kGI=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &fixed_kgi) && fixed_kgi>=0.0) continue;
    } else if(strncasecmp(cptr, "VB=", 3)==0 && strlen(cptr)>3) {
      if(!atofCheck(cptr+3, &fixed_vb) && fixed_vb>=0.0 && fixed_vb<1.0) continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &fixed_dt_liver) && fixed_dt_liver>=0.0) continue;
    } else if(strncasecmp(cptr, "DELAYL=", 7)==0 && strlen(cptr)>7) {
      if(!atofCheck(cptr+7, &fixed_dt_liver) && fixed_dt_liver>=0.0) continue;
    } else if(strncasecmp(cptr, "DELAYP=", 7)==0 && strlen(cptr)>7) {
      if(!atofCheck(cptr+7, &fixed_dt_gi) && fixed_dt_gi>=0.0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "SIM=", 4)==0 && strlen(cptr)>4) {
      strlcpy(simfile, 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(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(ttacfile, 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(!ttacfile[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("btacfile := %s\n", btacfile);
    printf("ttacfile := %s\n", ttacfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(simfile[0]) printf("simfile := %s\n", simfile);
    printf("weights := %d\n", weights);
    if(!isnan(fixed_kgi)) printf("fixed_kGI := %g\n", fixed_kgi);
    if(!isnan(fixed_vb)) printf("fixed_vb := %g\n", fixed_vb);
    if(!isnan(fixed_dt_liver)) printf("fixed_dt_liver := %g\n", fixed_dt_liver);
    if(!isnan(fixed_dt_gi)) printf("fixed_dt_gi := %g\n", fixed_dt_gi);
    fflush(stdout);
  }
 
 
  /*
   *  Read the data
   */
  if(verbose>1) printf("reading TACs\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC ttac, btac; tacInit(&ttac); tacInit(&btac);
 
  if(tacRead(&ttac, ttacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  if(verbose>2) {
    printf("ttac.fileformat := %s\n", tacFormattxt(ttac.format));
    printf("ttacNr := %d\n", ttac.tacNr);
    printf("ttac.sampleNr := %d\n", ttac.sampleNr);
    fflush(stdout);
  }
 
  if(tacRead(&btac, btacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  if(verbose>2) {
    printf("btac.fileformat := %s\n", tacFormattxt(btac.format));
    printf("btacNr := %d\n", btac.tacNr);
    printf("btac.sampleNr := %d\n", btac.sampleNr);
    fflush(stdout);
  }
  if(btac.tacNr>1) {
    if(verbose>0) fprintf(stderr, "Warning: BTAC file contains more than one TAC.\n");
    btac.tacNr=1;
  }
 
  if(ttac.sampleNr<7 || btac.sampleNr<7) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&ttac)>0 || tacNaNs(&btac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&ttac, &status);
  tacSortByTime(&btac, &status);
  /* Convert sample times into seconds */
  if(tacXUnitConvert(&ttac, UNIT_SEC, &status)!=TPCERROR_OK ||
     tacXUnitConvert(&btac, UNIT_SEC, &status)!=TPCERROR_OK) 
  {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&ttac); return(2);
  }
  /* Convert BTAC concentrations into TTAC units */
  if(tacYUnitConvert(&btac, ttac.cunit, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: check and set the data units.\n");
    tacFree(&ttac); tacFree(&ttac); return(2);
  }
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&ttac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
 
 
  /* Set places for fitted TACs */
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(tacAllocateMore(&ttac, ttac.tacNr)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&ttac); return(2);
  }
 
  /* Check and set weights */
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(weights==0) {
    if(!tacIsWeighted(&ttac)) {
      ttac.weighting=WEIGHTING_OFF;
      for(int i=0; i<ttac.sampleNr; i++) ttac.w[i]=1.0;
    } 
  } else if(weights==1) {
    ttac.weighting=WEIGHTING_OFF;
    for(int i=0; i<ttac.sampleNr; i++) ttac.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&ttac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); return(2);
    }
  }
  /* Number of samples with positive weight is needed to calculate AIC */
  unsigned int wsampleNr=tacWSampleNr(&ttac);
  if(verbose>2) printf("wsampleNr := %u\n", wsampleNr);
  if(wsampleNr<5) {
    fprintf(stderr, "Error: too few samples for fitting.\n");
    tacFree(&ttac); tacFree(&btac); return(2);
  }
 
 
  /*
   *  Prepare PAR structure for printing and saving model parameters
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &ttac, 8, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&btac); parFree(&par); return(3);
  }
  iftFree(&par.h); // remove stupid header info
  /* 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);
  }
  par.tacNr=ttac.tacNr; par.parNr=8;
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<par.tacNr; i++) {
    par.r[i].model=modelCodeIndex("radiowater-liver");
    par.r[i].dataNr=tacWSampleNr(&ttac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "K1a"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
  strcpy(par.n[1].name, "K1p"); par.n[1].unit=UNIT_ML_PER_ML_MIN;
  strcpy(par.n[2].name, "k2"); par.n[2].unit=UNIT_PER_MIN;
  strcpy(par.n[3].name, "Vb"); par.n[3].unit=UNIT_ML_PER_ML;
  strcpy(par.n[4].name, "LdeltaT"); par.n[4].unit=UNIT_SEC;
  strcpy(par.n[5].name, "PdeltaT"); par.n[5].unit=UNIT_SEC;
  strcpy(par.n[6].name, "kGI"); par.n[6].unit=UNIT_PER_MIN;
  strcpy(par.n[7].name, "p"); par.n[7].unit=UNIT_ML_PER_ML;
  /* set file names */
  iftPut(&par.h, "inputfile", btacfile, 0, NULL);
  iftPut(&par.h, "datafile", ttacfile, 0, NULL);
 
 
  /*
   *  Fit radiowater model for the TTACs
   */
  if(verbose==1) {printf("\nfitting...\n"); fflush(stdout);}
  int failed=0;
//#pragma omp parallel for
  for(int ri=0; ri<ttac.tacNr; ri++) {
    if(verbose>1) {printf("\nfitting %s\n", ttac.c[ri].name); fflush(stdout);}
    /* Set data pointers for the fit */
    FITDATA fitdata;
    fitdata.ni=btac.sampleNr;
    fitdata.xi=btac.x;
    fitdata.yi=btac.c[0].y;
    fitdata.nt=ttac.sampleNr;
    if(ttac.isframe) {
      fitdata.xt1=ttac.x1;
      fitdata.xt2=ttac.x2;
      fitdata.xt=NULL;
    } else {
      fitdata.xt=ttac.x;
      fitdata.xt1=NULL;
      fitdata.xt2=NULL;
    }
    fitdata.yt=ttac.c[ri].y;
    fitdata.w=ttac.w;
    fitdata.syt=ttac.c[ttac.tacNr+ri].y;
    if(verbose>10) fitdata.verbose=verbose-10; else fitdata.verbose=0;
    /* Set NLLS options */
    NLOPT nlo; nloptInit(&nlo);
    if(nloptAllocate(&nlo, 7)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n"); fflush(stderr); failed++;
      nloptFree(&nlo); continue;
    }
    nlo._fun=func_wliv; nlo.maxFunCalls=50000;
    nlo.fundata=&fitdata;
    nlo.totalNr=7;
    /* Set initial values and limits */
    // Kudomi et al 2008: kGI 0.497+-0.153 1/min, PdeltaT 0.7 +- 5.1 s, fA+fP 1-2 mL/(min*mL)
    // Slimani et al 2008: fA 0.15+-0.07, fP 1.11+-0.34, p 0.67 +- 0.03, Va 0.085+-0.054
    // Rijzewijk et al 2010: fA+fP 0.5-1.5
    nlo.xlower[0]=0.0;   nlo.xupper[0]=3.0;  nlo.xfull[0]=0.75; nlo.xtol[0]=0.00002; 
    nlo.xlower[1]=0.0;   nlo.xupper[1]=9.0;  nlo.xfull[1]=4.0;  nlo.xtol[1]=0.00005; 
    nlo.xlower[2]=0.01;  nlo.xupper[2]=20.0; nlo.xfull[2]=6.0;  nlo.xtol[2]=0.0001; 
    nlo.xlower[3]=0.0;   nlo.xupper[3]=0.80; nlo.xfull[3]=0.01; nlo.xtol[3]=0.00001; 
    nlo.xlower[4]=-30.0; nlo.xupper[4]=30.0; nlo.xfull[4]=11.0; nlo.xtol[4]=0.01; 
    nlo.xlower[5]=2.0;   nlo.xupper[5]=30.0; nlo.xfull[5]=14.0; nlo.xtol[5]=0.01; 
    nlo.xlower[6]=0.05;  nlo.xupper[6]=2.0;  nlo.xfull[6]=1.0;  nlo.xtol[6]=0.0001; 
    if(!isnan(fixed_vb)) {
      nlo.xlower[3]=nlo.xupper[3]=nlo.xfull[3]=fixed_vb; nlo.xtol[3]=0.0; 
    }
    if(!isnan(fixed_dt_liver)) {
      nlo.xlower[4]=nlo.xupper[4]=nlo.xfull[4]=fixed_dt_liver; nlo.xtol[4]=0.0; 
    }
    if(!isnan(fixed_dt_gi)) {
      nlo.xlower[5]=nlo.xupper[5]=nlo.xfull[5]=fixed_dt_gi; nlo.xtol[5]=0.0; 
    }
    if(!isnan(fixed_kgi)) {
      nlo.xlower[6]=nlo.xupper[6]=nlo.xfull[6]=fixed_kgi; nlo.xtol[6]=0.0; 
    }
 
 
#if(0)
    // call function for testing
    fitdata.verbose=10;
    double wss=func_wliv(nlo.totalNr, nlo.xfull, &fitdata);
    if(verbose>1) {
      printf("\nTime\tTTAC\tsimTTAC\n");
      for(unsigned int i=0; i<fitdata.nt; i++)
          printf("%g\t%g\t%g\n", ttac.x[i], ttac.c[ri].y[i], ttac.c[ri+ttac.tacNr].y[i]);
      printf("\n"); fflush(stdout);
    }
    if(verbose>1) {printf("\twss := %g\n", wss); fflush(stdout);}
 
    /* Copy parameters */
    par.r[ri].p[0]=nlo.xfull[0];
    par.r[ri].p[1]=nlo.xfull[1];
    par.r[ri].p[2]=nlo.xfull[2];
    par.r[ri].p[3]=nlo.xfull[3];
    par.r[ri].p[4]=nlo.xfull[4];
    par.r[ri].p[5]=nlo.xfull[5];
    par.r[ri].p[6]=nlo.xfull[6];
    par.r[ri].p[7]=(nlo.xfull[0]+nlo.xfull[1])/nlo.xfull[2];
    par.r[ri].wss=wss;
 
#else
 
    /* Fit */
    if(verbose>2) {
      printf("initial guess\n");
      nlo.funval=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
      nloptWrite(&nlo, stdout);
      fflush(stdout);
    }
    if(nloptSimplexARRS(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr); failed++;
      nloptFree(&nlo); continue;
    }
    double wss1=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>6) {
      printf("\nTime\tTTAC\tsimTTAC\n");
      for(unsigned int i=0; i<fitdata.nt; i++)
          printf("%g\t%g\t%g\n", ttac.x[i], ttac.c[ri].y[i], ttac.c[ri+ttac.tacNr].y[i]);
      printf("\n"); fflush(stdout);
    }
    if(verbose>2) nloptWrite(&nlo, stdout);
    if(verbose>2) printf(" wss1 := %g\n", wss1);
    /* Calculate AIC */
    double aic1=aicSS(wss1, wsampleNr, parFreeNr(nlo.totalNr, nlo.xlower, nlo.xupper));
    if(verbose>2) printf(" AIC1 := %g\n", aic1);
    
    /* Copy parameters */
    par.r[ri].p[0]=nlo.xfull[0];
    par.r[ri].p[1]=nlo.xfull[1];
    par.r[ri].p[2]=nlo.xfull[2];
    par.r[ri].p[3]=nlo.xfull[3];
    par.r[ri].p[4]=nlo.xfull[4];
    par.r[ri].p[5]=nlo.xfull[5];
    par.r[ri].p[6]=nlo.xfull[6];
    par.r[ri].p[7]=(nlo.xfull[0]+nlo.xfull[1])/nlo.xfull[2];
    par.r[ri].wss=wss1;
 
    /*
     *  Another fit, assuming only arterial input, in case tumours are included
     */
    nlo.xlower[0]=0.0;   nlo.xupper[0]=8.0;  nlo.xfull[0]=0.75; nlo.xtol[0]=0.00005; 
    nlo.xlower[1]=0.0;   nlo.xupper[1]=0.0;  nlo.xfull[1]=0.0;  nlo.xtol[1]=0.0; 
    nlo.xlower[2]=0.01;  nlo.xupper[2]=10.0; nlo.xfull[2]=1.0;  nlo.xtol[2]=0.00002; 
    nlo.xlower[3]=0.0;   nlo.xupper[3]=0.80; nlo.xfull[3]=0.01; nlo.xtol[3]=0.00001; 
    nlo.xlower[4]=-30.0; nlo.xupper[4]=30.0; nlo.xfull[4]=11.0; nlo.xtol[4]=0.01; 
    nlo.xlower[5]=0.0;   nlo.xupper[5]=0.0;  nlo.xfull[5]=0.0;  nlo.xtol[5]=0.0; 
    nlo.xlower[6]=0.0;   nlo.xupper[6]=0.0;  nlo.xfull[6]=0.0;  nlo.xtol[6]=0.0; 
    if(!isnan(fixed_vb)) {
      nlo.xlower[3]=nlo.xupper[3]=nlo.xfull[3]=fixed_vb; nlo.xtol[3]=0.0; 
    }
    if(!isnan(fixed_dt_liver)) {
      nlo.xlower[4]=nlo.xupper[4]=nlo.xfull[4]=fixed_dt_liver; nlo.xtol[4]=0.0; 
    }
    if(verbose>2) {
      printf("initial guess\n");
      nlo.funval=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
      nloptWrite(&nlo, stdout);
      fflush(stdout);
    }
    if(nloptSimplexARRS(&nlo, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr); failed++;
      nloptFree(&nlo); continue;
    }
    double wss2=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>6) {
      printf("\nTime\tTTAC\tsimTTAC\n");
      for(unsigned int i=0; i<fitdata.nt; i++)
          printf("%g\t%g\t%g\n", ttac.x[i], ttac.c[ri].y[i], ttac.c[ri+ttac.tacNr].y[i]);
      printf("\n"); fflush(stdout);
    }
    if(verbose>2) nloptWrite(&nlo, stdout);
    if(verbose>2) printf(" wss2 := %g\n", wss2);
    /* Calculate AIC */
    double aic2=aicSS(wss2, wsampleNr, parFreeNr(nlo.totalNr, nlo.xlower, nlo.xupper));
    if(verbose>2) printf(" AIC2 := %g\n", aic2);
 
    /* Use the results with better AIC */
    if(!(aic1<aic2)) {
      /* Copy parameters */
      par.r[ri].p[0]=nlo.xfull[0];
      par.r[ri].p[1]=nlo.xfull[1];
      par.r[ri].p[2]=nlo.xfull[2];
      par.r[ri].p[3]=nlo.xfull[3];
      par.r[ri].p[4]=nlo.xfull[4];
      par.r[ri].p[5]=nlo.xfull[5];
      par.r[ri].p[6]=nlo.xfull[6];
      par.r[ri].p[7]=(nlo.xfull[0]+nlo.xfull[1])/nlo.xfull[2];
      par.r[ri].wss=wss2;
    } else {
      /* We need to compute the fitted TACs again with the results from the first fit */
      nlo.xfull[0]=par.r[ri].p[0];
      nlo.xfull[1]=par.r[ri].p[1];
      nlo.xfull[2]=par.r[ri].p[2];
      nlo.xfull[3]=par.r[ri].p[3];
      nlo.xfull[4]=par.r[ri].p[4];
      nlo.xfull[5]=par.r[ri].p[5];
      nlo.xfull[6]=par.r[ri].p[6];
      nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    }
#endif
 
    nloptFree(&nlo);
  }
  if(failed>0) {tacFree(&ttac); tacFree(&btac); parFree(&par); return(5);}
 
 
  /* Print and save the parameters */
  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(&ttac); tacFree(&btac); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); 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(&ttac, &fit);
    for(int r=0; r<ttac.tacNr; r++)
      for(int i=0; i<ttac.sampleNr; i++)
        fit.c[r].y[i]=ttac.c[r+ttac.tacNr].y[i];
    /* Plot */
    if(tacPlotFitSVG(&ttac, &fit, "", nan(""), nan(""), nan(""), nan(""), svgfile, NULL)!=TPCERROR_OK) 
    {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); parFree(&par); tacFree(&fit);
      return(21);
    }
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
    tacFree(&fit);
  }
 
 
  /*
   *  Compute and save fitted TTACs at BTAC sample times, if requested
   */
  if(simfile[0]) {
    if(verbose>1) printf("calculating simulated TTACs\n");
 
    /* Allocate memory for simulated TTACs */
    TAC sim; tacInit(&sim);
    if(tacDuplicate(&ttac, &sim) || tacAllocateMoreSamples(&sim, btac.sampleNr-ttac.sampleNr)) {
      fprintf(stderr, "Error: cannot allocate space for simulated TTACs\n");
      tacFree(&ttac); tacFree(&btac); parFree(&par); tacFree(&sim);
      return(31);
    }
    sim.sampleNr=btac.sampleNr;
    sim.isframe=btac.isframe;
    (void)tacXCopy(&btac, &sim, 0, sim.sampleNr-1);
    sim.weighting=WEIGHTING_OFF;
 
    /* Set data pointers for the function data structure */
    FITDATA fitdata;
    fitdata.ni=btac.sampleNr;
    fitdata.xi=btac.x;
    fitdata.yi=btac.c[0].y;
    fitdata.nt=btac.sampleNr;
    if(btac.isframe) {
      fitdata.xt1=btac.x1;
      fitdata.xt2=btac.x2;
      fitdata.xt=NULL;
    } else {
      fitdata.xt=btac.x;
      fitdata.xt1=NULL;
      fitdata.xt2=NULL;
    }
    fitdata.w=btac.w;
 
    /* Simulate one TAC at a time */
    for(int i=0; i<sim.tacNr; i++) {
      fitdata.yt=fitdata.syt=sim.c[i].y;
      func_wliv(7, par.r[i].p, &fitdata);
    }
 
    /* Save simulated TTACs */
    if(verbose>1) printf("writing %s\n", simfile);
    FILE *fp; fp=fopen(simfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing (%s)\n", simfile);
      tacFree(&ttac); tacFree(&btac); parFree(&par); tacFree(&sim); return(32);
    }
    int ret=tacWrite(&sim, fp, TAC_FORMAT_PMOD, 1, &status);
    fclose(fp); tacFree(&sim);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
      tacFree(&ttac); tacFree(&btac); parFree(&par); return(33);
    }
    if(verbose>=0) {printf("%s saved.\n", simfile); fflush(stdout);}
  }
 
 
  tacFree(&ttac); tacFree(&btac); parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_wliv(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=7 || p==NULL || fdata==NULL || d->ni<1 || d->nt<1) return(nan(""));
  if(d->verbose>9) {
    printf("p[]: %g", p[0]);
    for(int i=1; i<parNr; i++) printf(" %g", p[i]);
    printf("\n"); fflush(stdout);
  }
 
  /* Process parameters; data is in seconds, therefore rate constants are converted to per sec */
  double K1a=p[0]/60.0;
  double K1p=p[1]/60.0;
  double k2=p[2]/60.0;
  double Vb=p[3];
  double LdT=p[4];
  double PdT=p[5];
  double kGI=p[6]/60.0;
 
  /* Make arterial input TAC with delay */
  double x[d->ni], *ca=d->yi;
  for(unsigned int i=0; i<d->ni; i++) x[i]=d->xi[i]+LdT;
 
  /* Make portal vein input TAC with both delays and dispersion */
  double x2[d->ni], cp[d->ni], sy[d->ni];
  for(unsigned int i=0; i<d->ni; i++) x2[i]=x[i]+PdT;
  if(kGI<1.0E-20) { // nothing comes out of portal vein
    for(unsigned int i=0; i<d->ni; i++) cp[i]=0.0;
  } else {
    for(unsigned int i=0; i<d->ni; i++) sy[i]=d->yi[i];
    if(simDispersion(x2, sy, d->ni, 1.0/kGI, 0.0, cp)!=0) return(nan(""));
    /* interpolate to the same sample times as the arterial input */
    if(liInterpolate(x2, sy, d->ni, x, cp, NULL, NULL, d->ni, 3, 1, 0)) return(nan(""));
  }
 
  /* Simulate tissue curve at input sample times */
  if(simC1DI(x, ca, cp, d->ni, K1a, K1p, k2, sy)!=0) return(nan(""));
  /* Simulate vascular contribution */
  for(unsigned int i=0; i<d->ni; i++) sy[i]*=(1.0-Vb);
  if(K1a+K1p > 0.0)
    for(unsigned int i=0; i<d->ni; i++) sy[i]+= Vb*(K1a*ca[i] + K1p*cp[i])/(K1a+K1p);
  else
    for(unsigned int i=0; i<d->ni; i++) sy[i]+= Vb*ca[i];
 
  /* Interpolate simulated tissue curve to the sample times of measured tissue */
  if(d->xt1==NULL || d->xt2==NULL) {
    if(liInterpolate(x, sy, d->ni, d->xt, d->syt, NULL, NULL, d->nt, 3, 1, 0))
      return(nan(""));
  } else {
    if(liInterpolateForPET(x, sy, d->ni, d->xt1, d->xt2, d->syt, NULL, NULL, d->nt, 3, 1, 0))
      return(nan(""));
  }
 
  /* 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->nt; i++) {
    double v=d->syt[i] - d->yt[i];
    wss+=d->w[i]*v*v;
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/