fit_wpul.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_wpul(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 radiowater model to pulmonary TTACs using BTAC",
  "from RV cavity as the input function. Small delay and dispersion of input",
  "function prior to lungs is allowed and fitted by default. Regional",
  "radioactivity concentration is modelled as Cpet(t) = Vb*Cb(t) + Ct(t),",
  "and thus the estimated K1 = (1-Va-Vb)*f. Air volume fraction (Va) could",
  "be estimated from transmission or CT scan.",
  " ",
  "Usage: @P [Options] btacfile ttacfile [parfile]",
  " ",
  "Options:",
  " -tau=<value>",
  "     Dispersion time constant is constrained to given value (sec).",
  " -Vb=<value>",
  "     Blood volume is constrained to given value (mL/mL of lung).",
  " -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, fit_wrlv, fit_disp",
  " ",
  "Keywords: TAC, modelling, perfusion, lungs, 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_tau=nan("");
  double  fixed_vb=nan("");
 
  drandSeed(1);
 
 
  /*
   *  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, "TAU=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &fixed_tau) && fixed_tau>=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, "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_tau)) printf("fixed_tau := %g\n", fixed_tau);
    if(!isnan(fixed_vb)) printf("fixed_vb := %g\n", fixed_vb);
    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==2) {
    if(verbose>0) fprintf(stdout, "Note: BTAC file assumed to contain RV and LV BTACs.\n");
  } else if(btac.tacNr>2) {
    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);
  }
 
 
  /* 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);
    }
  }
 
 
#if(0)
  if(btac.isframe) {
    TAC tmp; tacInit(&tmp);
    if(tacFramesToSteps(&btac, &tmp, 0.01)) {
      fprintf(stderr, "Error: invalid input sample times.\n");
      tacFree(&ttac); tacFree(&btac); return(2);
    }
    tacFree(&btac);
    if(tacDuplicate(&tmp, &btac)) {
      fprintf(stderr, "Error: invalid input sample times.\n");
      tacFree(&ttac); tacFree(&btac); tacFree(&tmp); return(2);
    }
    tacFree(&tmp);
  }
#endif
 
 
  /*
   *  If LV cavity BTAC is available, then add 5% of it to RV input, with delay.
   *  LV input should probably have its own compartment, but for now like this.
   */
  if(btac.tacNr==2) {
    double x[btac.sampleNr], y[btac.sampleNr];
    for(int i=0; i<btac.sampleNr; i++) x[i]=btac.x[i]+4.0;
    if(liInterpolate(x, btac.c[1].y, btac.sampleNr, btac.x, y, NULL, NULL, btac.sampleNr, 3, 1, 0))
    {
      fprintf(stderr, "Error: invalid input sample times.\n");
      tacFree(&ttac); tacFree(&btac); return(2);
    }
    for(int i=0; i<btac.sampleNr; i++) btac.c[0].y[i]=0.95*btac.c[0].y[i]+0.05*y[i];
  }
 
 
 
  /*
   *  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-lung");
    par.r[i].dataNr=tacWSampleNr(&ttac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "tau"); par.n[0].unit=UNIT_SEC;
  strcpy(par.n[1].name, "deltaT"); par.n[1].unit=UNIT_SEC;
  strcpy(par.n[2].name, "Vb"); par.n[2].unit=UNIT_ML_PER_ML;
  strcpy(par.n[3].name, "K1"); par.n[3].unit=UNIT_ML_PER_ML_MIN;
  strcpy(par.n[4].name, "k2"); par.n[4].unit=UNIT_PER_MIN;
  strcpy(par.n[5].name, "K1/k2"); par.n[5].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 lung 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, 5)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n"); fflush(stderr); failed++;
      nloptFree(&nlo); continue;
    }
    nlo._fun=func_wpul; nlo.maxFunCalls=10000;
    nlo.fundata=&fitdata;
    nlo.totalNr=5;
    /* Set initial values and limits */
    // dispersion time tau
    if(isnan(fixed_tau)) {
      nlo.xlower[0]=0.0; nlo.xupper[0]=8.0; nlo.xfull[0]=2.0; nlo.xtol[0]=0.01; 
    } else {
      nlo.xlower[0]=nlo.xupper[0]=nlo.xfull[0]=fixed_tau; nlo.xtol[0]=0.0; 
    }
    // delay time deltaT
    nlo.xlower[1]=0.0; nlo.xupper[1]=4.0; nlo.xfull[1]=2.0; nlo.xtol[1]=0.001; 
    // Vb (of total ROI volume including air)
    if(isnan(fixed_vb)) {
      nlo.xlower[2]=0.00; nlo.xupper[2]=0.20; nlo.xfull[2]=0.01; nlo.xtol[2]=0.00001; 
    } else {
      nlo.xlower[2]=nlo.xupper[2]=nlo.xfull[2]=fixed_vb; nlo.xtol[2]=0.0; 
    }
    // K1
    nlo.xlower[3]=0.00; nlo.xupper[3]=5.0; nlo.xfull[3]=1.4; nlo.xtol[3]=0.00001; 
    // k2
    nlo.xlower[4]=0.00; nlo.xupper[4]=20.0; nlo.xfull[4]=5.0; nlo.xtol[4]=0.00001; 
 
 
#if(0)
    // call function for testing
    fitdata.verbose=10;
    double wss=func_wpul(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, 30, &status)!=TPCERROR_OK) {
//    if(nloptIATGO(&nlo, 1, 50, 0.05, &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 */
    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]=par.r[ri].p[3]/par.r[ri].p[4];
    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);
  }
 
 
  /*
   *  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_wpul(5, 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_wpul(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=5 || 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 tau=p[0]; if(tau<0.0) tau=0.0;
  double deltaT=p[1];
  double Vb=p[2]; if(Vb<0.0) Vb=0.0;
  double K1=p[3]/60.0; if(K1<0.0) K1=0.0;
  double k2=p[4]/60.0; if(k2<0.0) k2=0.0;
 
  /* 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]+deltaT;
  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 curve at input sample times */
  if(simC1(x, y, d->ni, K1, k2, sy)!=0) return(nan(""));
  for(unsigned int i=0; i<d->ni; i++) sy[i]+=Vb*y[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);
}
/*****************************************************************************/
 
/*****************************************************************************/