fit_wcbf.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_wcbf(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;

  double        pGM, pWM;

  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the radiowater model with PVC to brain TTACs using BTAC",
  "as the input function (Iida et al., 2000).",
  " ",
  "Usage: @P [Options] btacfile ttacfile [parfile]",
  " ",
  "Options:",
  " -Va=<value>",
  "     Arterial blood volume is constrained to given value (mL/mL).",
  " -dt=<value>",
  "     Delay time (activity appearance time in tissue minus that in blood)",
  "     in seconds. Dispersion in BTAC is corrected in relation to delay",
  "     as tau=0.31*dt-0.30 (Iida et al., 1989), when tau>0.",
  " -pGM=<value>",
  "     Partition coefficient of water in grey matter; by default 1.03 mL/mL.",
  " -pWM=<value>",
  "     Partition coefficient of water in white matter; by default 0.86 mL/L.",
  " -fWM=<value> | fit",
  "     Blood flow in white matter; by default 0.21 mL/(min*mL).",
  "     Enter 'fit' to estimate fW as a free parameter.",
  " -w1 | -wf",
  "     By default, weights in data file are used, if available.",
  "     With these options all weights can be set to 1.0 (no weighting)",
  "     or based on frame durations.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -fit=<Filename>",
  "     Fitted TTACs are written in specified TAC 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.",
  " ",
  "Reference:",
  "1. Iida H et al., Quantitation of regional cerebral blood flow corrected",
  "   for partial volume effect using O-15 water and PET: I. Theory, error",
  "   analysis, and stereologic comparison.",
  "   J Cereb Blood Flow Metab. 2000;20(8):1237-1251.",
  "   https://doi.org/10.1097/00004647-200008000-00009",
  "2. Iida et al., A determination of the regional brain/blood partition",
  "   coefficient of water using dynamic positron emission tomography.",
  "   J Cereb Blood Flow Metab. 1989;9(6):874-885.",
  "   https://doi.org/10.1038/jcbfm.1989.121",
  " ",
  "See also: fit_h2o, bfmh2o, fitk2",
  " ",
  "Keywords: TAC, modelling, perfusion, brain, 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], fitfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  double  fixed_dt=nan("");
  double  fixed_va=nan("");
  double  fixed_pg=1.03;
  double  fixed_pw=0.86;
  double  fixed_fw=0.21;
 
  drandSeed(1);
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  btacfile[0]=ttacfile[0]=parfile[0]=svgfile[0]=fitfile[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, "VA=", 3)==0 && strlen(cptr)>3) {
      if(!atofCheck(cptr+3, &fixed_va) && fixed_va>=0.0 && fixed_va<1.0) continue;
    } else if(strncasecmp(cptr, "DT=", 3)==0 && strlen(cptr)>3) {
      if(!atofCheck(cptr+3, &fixed_dt) && fixed_dt>-50.0 && fixed_dt<50.0) continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &fixed_dt) && fixed_dt>-50.0 && fixed_dt<50.0) continue;
    } else if(strncasecmp(cptr, "PGM=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &fixed_pg) && fixed_pg>0.0 && fixed_pg<1.05) continue;
    } else if(strncasecmp(cptr, "PWM=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &fixed_pw) && fixed_pw>0.0 && fixed_pw<1.05) continue;
    } else if(strncasecmp(cptr, "FWM=FIT", 5)==0) {
      fixed_fw=nan(""); continue;
    } else if(strncasecmp(cptr, "FWM=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &fixed_fw) && fixed_fw>=0.0 && fixed_fw<1.0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "FIT=", 4)==0 && strlen(cptr)>4) {
      strlcpy(fitfile, 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(fitfile[0]) printf("fitfile := %s\n", fitfile);
    if(simfile[0]) printf("simfile := %s\n", simfile);
    printf("weights := %d\n", weights);
    printf("fixed_pg := %g\n", fixed_pg);
    printf("fixed_pw := %g\n", fixed_pw);
    if(!isnan(fixed_va)) printf("fixed_va := %g\n", fixed_va);
    if(!isnan(fixed_fw)) printf("fixed_fw := %g\n", fixed_fw);
    if(!isnan(fixed_dt)) printf("fixed_dt := %g\n", fixed_dt);
    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 two TACs.\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);
  /* Fix any gaps and overlap in data if possible and if not then quit with error */
  if(tacSetXContiguous(&ttac) || tacSetXContiguous(&btac)) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&ttac); tacFree(&ttac); return(2);
  }
  /* 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);
    printf("final_xunit := %s\n", unitName(ttac.tunit));
  }
 
 
  /* 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);
    }
  }
 
 
  /*
   *  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, 6, &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=6;
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<par.tacNr; i++) {
    par.r[i].model=modelCodeIndex("radiowater-brain");
    par.r[i].dataNr=tacWSampleNr(&ttac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "aGM"); par.n[0].unit=UNIT_UNITLESS;
  strcpy(par.n[1].name, "fGM"); par.n[1].unit=UNIT_ML_PER_ML_MIN;
  strcpy(par.n[2].name, "aWM"); par.n[2].unit=UNIT_UNITLESS;
  strcpy(par.n[3].name, "fWM"); par.n[3].unit=UNIT_ML_PER_ML_MIN;
  strcpy(par.n[4].name, "Va");  par.n[4].unit=UNIT_ML_PER_ML;
  strcpy(par.n[5].name, "dT");  par.n[5].unit=UNIT_SEC;
  /* set file names */
  iftPut(&par.h, "inputfile", btacfile, 0, NULL);
  iftPut(&par.h, "datafile", ttacfile, 0, NULL);
  /* set fixed parameters */
  iftPutDouble(&par.h, "pGM", fixed_pg, 0, NULL);
  iftPutDouble(&par.h, "pWM", fixed_pw, 0, NULL);
 
 
  /*
   *  Fit radiowater model with PVC for the brain 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;
    fitdata.pGM=fixed_pg;
    fitdata.pWM=fixed_pw;
    if(verbose>10) fitdata.verbose=verbose-10; else fitdata.verbose=0;
    /* Set NLLS options */
    NLOPT nlo; nloptInit(&nlo);
    if(nloptAllocate(&nlo, 6)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n"); fflush(stderr); failed++;
      nloptFree(&nlo); continue;
    }
    nlo._fun=func_wcbf; nlo.maxFunCalls=50000;
    nlo.fundata=&fitdata;
    nlo.totalNr=6;
    /* Set initial values and limits */
    // alphaGM
    nlo.xlower[0]=0.00; nlo.xupper[0]=1.0; nlo.xfull[0]=0.5; nlo.xtol[0]=0.002; 
    // fGM
    nlo.xlower[1]=0.40; nlo.xupper[1]=1.4; nlo.xfull[1]=0.8; nlo.xtol[1]=0.002; 
    // parameter[2] = alphaWM/(1-alphaGM-Va), i.e. alphaWM=parameter[2]*(1-alpgaGM-Va) 
    nlo.xlower[2]=0.00; nlo.xupper[2]=1.0; nlo.xfull[2]=0.6; nlo.xtol[2]=0.002; 
    // fWM
    if(isnan(fixed_fw)) {
      nlo.xlower[3]=0.0; nlo.xupper[3]=0.5; nlo.xfull[3]=0.2; nlo.xtol[3]=0.002; 
    } else {
      nlo.xlower[3]=nlo.xupper[3]=nlo.xfull[3]=fixed_fw; nlo.xtol[3]=0.0; 
    }
    // Va
    if(isnan(fixed_va)) {
      nlo.xlower[4]=0.0; nlo.xupper[4]=0.1; nlo.xfull[4]=0.02; nlo.xtol[4]=0.001; 
    } else {
      nlo.xlower[4]=nlo.xupper[4]=nlo.xfull[4]=fixed_va; nlo.xtol[4]=0.0; 
    }
    // Delay
    if(isnan(fixed_dt)) {
      nlo.xlower[5]=-1.0; nlo.xupper[5]=10.0; nlo.xfull[5]=5.0; nlo.xtol[5]=0.01; 
    } else {
      nlo.xlower[5]=nlo.xupper[5]=nlo.xfull[5]=fixed_dt; nlo.xtol[5]=0.0; 
    }
 
 
#if(0)
    // call function for testing
    fitdata.verbose=10;
    double wss=func_wcbf(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);}
 
#else
    /* Fit */
    if(verbose>2) {
      printf("initial guess\n");
      nlo.funval=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
      nloptWrite(&nlo, stdout);
      fflush(stdout);
    }
 
    //status.verbose=10;
//    if(nloptSimplexARRS(&nlo, 0, &status)!=TPCERROR_OK) {
//  if(nloptIATGO(&nlo, 1, 150, 0.10, &status)!=TPCERROR_OK) {
    if(nloptSimplexMS(&nlo, 500, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error)); fflush(stderr); failed++;
      nloptFree(&nlo); continue;
    }
    //status.verbose=1;
    double wss=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(" wss := %g\n", wss);
 
#endif
 
 
    /* Copy parameters */
    for(int i=0; i<par.parNr; i++) par.r[ri].p[i]=nlo.xfull[i];
    par.r[ri].p[2]*=1.0-par.r[ri].p[0]-par.r[ri].p[4]; // alphaWM=fraction*(1-alphaGM-Va)
    par.r[ri].wss=wss;
 
    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);
  }
 
 
  /*
   *  Save and/or plot fitted TACs, if requested
   */
  if(svgfile[0] || fitfile[0]) {
 
    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];
 
    /*
     *  SVG plot of fitted and original data
     */
    if(svgfile[0]) {
      if(verbose>1) printf("plotting measured and fitted data\n");
      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);
    }
 
    /*
     *  Save fitted TTACs
     */
    if(fitfile[0]) {
      if(verbose>1) printf("writing %s\n", fitfile);
      int ret=TPCERROR_OK;
      FILE *fp; fp=fopen(fitfile, "w");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file for writing fitted TTACs.\n");
        ret=TPCERROR_FAIL;
      } else {
        ret=tacWrite(&fit, fp, TAC_FORMAT_PMOD, 1, &status);
        fclose(fp);
        if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      }
      if(ret!=TPCERROR_OK) {
        tacFree(&ttac); tacFree(&btac); parFree(&par); tacFree(&fit);
        return(22);
      }
      if(verbose>0) printf("fitted TACs saved in %s.\n", fitfile);
    }
 
    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;
    fitdata.pGM=fixed_pg;
    fitdata.pWM=fixed_pw;
    fitdata.verbose=0;
 
    double pp[6];
 
    /* Simulate one TAC at a time */
    for(int ri=0; ri<sim.tacNr; ri++) {
      fitdata.yt=fitdata.syt=sim.c[ri].y;
      for(int j=0; j<6; j++) pp[j]=par.r[ri].p[j];
      pp[2]=par.r[ri].p[2]/(1.0-par.r[ri].p[0]-par.r[ri].p[4]); // fraction = alphaWM/(1-alphaGM-Va)
 
      func_wcbf(6, pp, &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_wcbf(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s() %d\n", __func__, d->verbose); fflush(stdout);}
  if(parNr!=6 || 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 */
  double aGM=p[0];
  double fGM=p[1]/60.0;
  double fWM=p[3]/60.0;
  double Va=p[4];
  double dT=p[5];
  double tau=0.31*dT-0.30; if(!(tau>=0.0)) tau=0.0;
  double aWM=p[2]*(1.0-aGM-Va); // p[2] = fraction of aWM from remaining vol 
  //printf("pGM=%g pWM=%g\n", d->pGM, d->pWM);
  
 
  /* Make input TAC with delay and dispersion */
  double x[d->ni], y[d->ni], sy[d->ni];
  for(unsigned int i=0; i<d->ni; i++) x[i]=d->xi[i]+dT;
  for(unsigned int i=0; i<d->ni; i++) y[i]=d->yi[i];
  if(simDispersion(x, y, d->ni, tau, 0.0, sy)!=0) return(nan(""));
 
  /* Simulate tissue and PET curves at input sample times */
  double cgm[d->ni];
  if(simC1_d(x, y, d->ni, fGM, fGM/d->pGM, cgm)!=0) return(nan(""));
  double cwm[d->ni];
  if(simC1_d(x, y, d->ni, fWM, fWM/d->pWM, cwm)!=0) return(nan(""));
  for(unsigned int i=0; i<d->ni; i++)
    sy[i] = Va*y[i] + aGM*cgm[i] + aWM*cwm[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);
}
/*****************************************************************************/
 
/*****************************************************************************/