fit_1tcm.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"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_1tcm(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int  ni;
  double       *xi;
  double       *yi;
  double       *yb;

  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 1TCM to regional TTACs using PTAC as input function",
  "and BTAC for vascular volume correction. Delay time between TTACs and",
  "PTAC/BTAC is fitted by default. Regional radioactivity concentration is",
  "modelled as Cpet(t) = Vb*Cb(t) + Ct(t),",
  "and thus the estimated parameters are reported per total ROI volume.",
  " ",
  "Usage: @P [Options] ptacfile btacfile ttacfile [parfile]",
  " ",
  "Options:",
  " -end=<value>",
  "     Tissue data length, used in the fitting, is restricted to",
  "     given time (min).",
  " -Vb=<value>",
  "     Blood volume is constrained to given value (mL/mL).",
  "     If Vb is set to zero, then btacfile is not read, but a place filler",
  "     must nonetheless be given.",
  " -Delay=<value>",
  "     Delay time is constrained to given value (sec).",
  " -wf | -wfd | -w1",
  "     With -wf Weights are based on frame length or sampling interval;",
  "     with -wfd the late frames are given less weight by using formula",
  "     weight=(frame duration)*exp(-t*ln(2)/halflife) (Thiele et al, 2008);",
  "     with -w1 all weights are set to 1.0 (no weighting);",
  "     by default, weights in TTAC file are used, if available.",
  " -i=<Isotope code>",
  "     Isotope, for example C-11, in case it is not found inside TTAC, but",
  "     is needed with option -wfd.",
  " -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
  " ",
  "If time units are not specified in files, minutes are assumed.",
  " ",
  "See also: fitk2, fit_h2o, bfmh2o, fit_wrlv",
  " ",
  "Keywords: TAC, modelling, compartmental model",
  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 tacReadModelingData2(
  const char *tissuefile,
  const char *inputfile1,
  const char *inputfile2,
  const char *inputfile3,
  double *fitdur,
  int cutInput,
  int *fitSampleNr,
  TAC *tis,
  TAC *inp,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(tis==NULL || inp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(tissuefile==NULL || inputfile1==NULL || strnlen(inputfile1, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check the function input */
  int input_nr=1;
  if(inputfile2!=NULL && strnlen(inputfile2, 1)>0) input_nr++;
  if(inputfile3!=NULL && strnlen(inputfile3, 1)>0) {
    if(input_nr<2) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      return TPCERROR_NO_DATA;
    }
    input_nr++;
  }
  if(verbose>2) printf("input_nr := %d\n", input_nr);
 
  /* Delete any previous data */
  tacFree(tis); tacFree(inp);
  if(fitSampleNr!=NULL) *fitSampleNr=0;
 
  int ret;
  /*
   *  Read tissue data
   */
  if(verbose>1) printf("reading tissue data in %s\n", tissuefile);
  ret=tacRead(tis, tissuefile, status);
  if(ret!=TPCERROR_OK) return(ret);
 
  /* Do not check frame number; static scan is fine here */
  
  /* Check for NaN's */
  if(tacNaNs(tis)>0) {
    tacFree(tis);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
 
  /* Sort the tissue data by increasing sample times */
  ret=tacSortByTime(tis, status);
  if(ret!=TPCERROR_OK) {
    tacFree(tis);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  
  /* Make sure that there is no overlap in sample frame times; samples must
     be sorted before this */
  ret=tacCorrectFrameOverlap(tis, status);
  if(ret!=TPCERROR_OK) {tacFree(tis); return ret;}
 
  /*
   *  Read first input data
   */
  if(verbose>1) printf("reading input data 1 in %s\n", inputfile1);
  ret=tacRead(inp, inputfile1, status);
  if(ret!=TPCERROR_OK) {tacFree(tis); return ret;}
  /* Check and correct the sample time unit */
  if(tis->tunit==UNIT_UNKNOWN) tis->tunit=inp->tunit;
  else if(inp->tunit==UNIT_UNKNOWN) inp->tunit=tis->tunit;
  if(inp->tunit==UNIT_UNKNOWN && verbose>0) {
    fprintf(stderr, "Warning: input sample time units not known.\n");}
  if(tis->tunit!=inp->tunit && tacXUnitConvert(inp, tis->tunit, status)) {
    tacFree(tis); tacFree(inp); return ret;}
  /* Check TAC nr */
  if(inp->tacNr>1) {
    if(verbose>0) fprintf(stderr, "Warning: using only first TAC in %s\n", inputfile1);
    inp->tacNr=1;
  }
  /* Check sample nr */
  if(inp->sampleNr<4) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    tacFree(tis); tacFree(inp); return TPCERROR_TOO_FEW;
  }
  /* Sort the data by increasing sample times */
  tacSortByTime(inp, NULL);
  /* If inp contains isotope and tis does not, then copy it */
  if(tacGetIsotope(tis)!=ISOTOPE_UNKNOWN) {
    int isotope=tacGetIsotope(inp);
    if(isotope!=ISOTOPE_UNKNOWN) tacSetIsotope(tis, isotope);
  }
 
  /*
   *  Read following input files, if available
   */
  char *fname;
  TAC tmptac; tacInit(&tmptac);
  for(int ii=2; ii<=input_nr; ii++) {
    if(ii==2) fname=(char*)inputfile2; else fname=(char*)inputfile3;
    if(verbose>1) printf("reading input data %d in %s\n", ii, fname);
    ret=tacRead(&tmptac, fname, status);
    if(ret!=TPCERROR_OK) {tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;}
    /* Check TAC nr */
    if(tmptac.tacNr>1) {
      if(verbose>0) fprintf(stderr, "Warning: using only first TAC in %s\n", fname);
      tmptac.tacNr=1;
    }
    /* Check sample nr */
    if(tmptac.sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return TPCERROR_TOO_FEW;
    }
    /* Sort the data by increasing sample times */
    tacSortByTime(&tmptac, NULL);
 
    /* Check and correct the sample time unit */
    if(tis->tunit==UNIT_UNKNOWN) tis->tunit=tmptac.tunit;
    else if(tmptac.tunit==UNIT_UNKNOWN) tmptac.tunit=tis->tunit;
    if(inp->tunit!=tmptac.tunit && tacXUnitConvert(&tmptac, inp->tunit, status)) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
 
    /* Check and correct the sample concentration unit */
    if(inp->cunit==UNIT_UNKNOWN) inp->cunit=tmptac.cunit;
    else if(tmptac.cunit==UNIT_UNKNOWN) tmptac.cunit=inp->cunit;
    if(inp->cunit!=tmptac.cunit && tacYUnitConvert(&tmptac, inp->cunit, status)) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
 
    /* Copy to input data */
    if(tacInterpolateInto(&tmptac, inp, NULL, NULL, status)!=TPCERROR_OK) {
      tacFree(tis); tacFree(inp); tacFree(&tmptac); return ret;
    }
    
    tacFree(&tmptac);
  } // next input file
 
  /* Set time unit to min */
  ret=tacXUnitConvert(tis, UNIT_MIN, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check that regional data times are in minutes.\n");
  }
  ret=tacXUnitConvert(inp, UNIT_MIN, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check that input data times are in minutes.\n");
  }
  /* Check that input and tissue time ranges are about the same */
  double iend, tend;
  {
    ret=tacXRange(inp, NULL, &iend); if(ret==0) ret=tacXRange(tis, NULL, &tend);
    if(ret || iend<=0.0 || tend<=0.0) {
      tacFree(tis); tacFree(inp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
      return TPCERROR_INVALID_XRANGE;
    }
    if(tend>10.0*iend || tend<0.10*iend) {
      if(verbose>0) fprintf(stderr, "Warning: check the sample time units.\n");
    }
  }
  ret=tacYUnitConvert(inp, tis->cunit, status);
  if(ret && verbose>0) {
    fprintf(stderr, "Warning: check the calibration units.\n");
  }
 
  /*
   *  Check and set fit time length
   */
  if(verbose>1) printf("checking and setting fit time length\n");
  /* Set fit duration */
  double starttime=0, endtime=*fitdur;
  int fnr=tacFittime(tis, &starttime, &endtime, NULL, NULL, status);
  if(verbose>3) {
    fprintf(stdout, "tis.sampleNr := %d\n", tis->sampleNr);
    fprintf(stdout, "starttime := %g\n", starttime);
    fprintf(stdout, "endtime := %g\n", endtime);
    //fprintf(stdout, "first := %d\n", first);
    //fprintf(stdout, "last := %d\n", last);
    fprintf(stdout, "fitSampleNr := %d\n", fnr);
  }
  *fitdur=endtime;
  if(fitSampleNr!=NULL) *fitSampleNr=fnr;
 
  /* Check that input data does not end much before fitdur */
  if(*fitdur>1.2*iend) {
    tacFree(tis); tacFree(inp);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
    return TPCERROR_TOO_FEW;
  }
 
  /* Cut off too many input samples to make calculation faster */
  if(cutInput && iend>*fitdur) {
    if(verbose>1) printf("cut off too many input samples\n");
    int i; double f;
    for(i=0; i<inp->sampleNr; i++) {
      if(inp->isframe) f=0.5*(inp->x1[i]+inp->x2[i]); else f=inp->x[i];
      if(f>(*fitdur)) break;
    }
    if(i<inp->sampleNr) i++; 
    inp->sampleNr=i;
    if(inp->sampleNr<4) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_TOO_FEW);
      tacFree(tis); tacFree(inp); return TPCERROR_TOO_FEW;
    }
  }
  if(verbose>2) fprintf(stdout, "inp.sampleNr := %d\n", inp->sampleNr);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
 
 
/*****************************************************************************/
 
/*****************************************************************************/
int nnlsRTCM(
  TAC *ttac,
  TAC *ptac,
  const int tcn,
  const int vbc,
  PAR *par,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(ttac, ptac, %d, %d, par)\n", __func__, tcn, vbc);
  if(ttac==NULL || ptac==NULL || par==NULL || tcn<1 || tcn>3) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(ttac->tacNr<1 || ttac->sampleNr<3 || ptac->tacNr<1 || ptac->sampleNr<3) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
 
 
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    ptacfile[FILENAME_MAX], btacfile[FILENAME_MAX], ttacfile[FILENAME_MAX],
          parfile[FILENAME_MAX], svgfile[FILENAME_MAX], simfile[FILENAME_MAX];
  int     weightMethod=WEIGHTING_UNKNOWN;
  double  fixed_dt=nan("");
  double  fixed_vb=nan("");
  double  tstop=nan(""); // fit end time
  isotope isot=ISOTOPE_UNKNOWN;
 
  drandSeed(1);
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  ptacfile[0]=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) {
      weightMethod=WEIGHTING_OFF; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weightMethod=WEIGHTING_ON_F; continue;
    } else if(strcasecmp(cptr, "WFD")==0) {
      weightMethod=WEIGHTING_ON_FD; continue;
    } else if(strncasecmp(cptr, "I=", 2)==0) {
      if((isot=isotopeIdentify(cptr+2))==ISOTOPE_UNKNOWN) {
        fprintf(stderr, "Error: invalid isotope '%s'.\n", cptr+2); return(1);}
      continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &fixed_dt) && fixed_dt>=0.0) continue;
    } else if(strncasecmp(cptr, "DT=", 3)==0 && strlen(cptr)>3) {
      if(!atofCheck(cptr+3, &fixed_dt) && fixed_dt>=0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0 && strlen(cptr)>4) {
      if(!atofCheck(cptr+4, &tstop) && tstop>0.0) continue;
    } else if(strncasecmp(cptr, "STOP=", 5)==0 && strlen(cptr)>5) {
      if(!atofCheck(cptr+5, &tstop) && tstop>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(ptacfile, argv[ai++], FILENAME_MAX);
  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);
  }
  /* If Vb is set to zero, then btacfile is not needed */
  if(fixed_vb==0.0) btacfile[0]=(char)0;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ptacfile := %s\n", ptacfile);
    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("weightMethod := %d\n", weightMethod);
    if(!isnan(tstop)) printf("tstop := %g\n", tstop);
    if(!isnan(fixed_dt)) printf("fixed_dt := %g\n", fixed_dt);
    if(!isnan(fixed_vb)) printf("fixed_vb := %g\n", fixed_vb);
    fflush(stdout);
  }
 
 
  /*
   *  Read tissue and input data
   */
  if(verbose>0) {printf("reading tissue and input data\n"); fflush(stdout);}
  TAC ttac, input; tacInit(&ttac); tacInit(&input);
  int fitSampleNr=0;
  double fitdur=1.0E+10; if(tstop>0.01) fitdur=tstop;
 
  if(tacReadModelingData2(ttacfile, ptacfile, btacfile, NULL, &fitdur, 0,
                          &fitSampleNr, &ttac, &input, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&input); return(2);
  }
  if(isot!=ISOTOPE_UNKNOWN) tacSetIsotope(&ttac, isot);
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(ttac.format));
    printf("tacNr := %d\n", ttac.tacNr);
    printf("tac.sampleNr := %d\n", ttac.sampleNr);
    printf("input.sampleNr := %d\n", input.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(ttac.tunit));
    printf("yunit := %s\n", unitName(ttac.cunit));
    printf("fitdur := %g s\n", fitdur);
    printf("isotope :=  %s\n", isotopeName(tacGetIsotope(&ttac)));
    fflush(stdout);
  }
 
 
  /* Check and set weights */
  if(verbose>0) {printf("setting weights\n"); fflush(stdout);}
  if(tacSetWeights(&ttac, weightMethod, fitSampleNr, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&input); return(3);
  }
  /* 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(&input); return(2);
  }
 
  /* Get x range, weights considered */
  double xmin, xmax;
  if(tacSampleXRange(&ttac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&ttac); tacFree(&input); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
 
 
  /*
   *  Prepare PAR structure for printing and saving model parameters
   */
  if(verbose>1) {printf("preparing space for parameters\n"); fflush(stdout);}
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &ttac, 5, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&ttac); tacFree(&input); 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=5;
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<par.tacNr; i++) {
    par.r[i].model=modelCodeIndex("1TCM");
    par.r[i].dataNr=tacWSampleNr(&ttac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "K1"); par.n[0].unit=UNIT_ML_PER_ML_MIN;
  strcpy(par.n[1].name, "k2"); par.n[1].unit=UNIT_PER_MIN;
  strcpy(par.n[2].name, "Vb"); par.n[2].unit=UNIT_ML_PER_ML;
  strcpy(par.n[3].name, "dT"); par.n[3].unit=UNIT_SEC;
  strcpy(par.n[4].name, "K1/k2"); par.n[4].unit=UNIT_ML_PER_ML;
  /* set file names */
  iftPut(&par.h, "datafile", ttacfile, 0, NULL);
  iftPut(&par.h, "plasma", ptacfile, 0, NULL);
  iftPut(&par.h, "blood", btacfile, 0, NULL);
 
 
  /*
   *  Delay correction, if fixed delay time was given
   */
  if(fixed_dt<0.0 || fixed_dt>0.0) {
    if(verbose>1) {printf("correcting data for user-defined delay time\n"); fflush(stdout);}
    if(tacDelay(&input, -fixed_dt/60.0, -1, &status) != TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&input); parFree(&par); return(5);
    }
    if(verbose>1) {printf("  input data moved %g s\n", -fixed_dt); fflush(stdout);}
  }
 
 
  /*
   *  Vascular volume correction, if fixed Vb was given
   */
  if(fixed_vb>0.0) {
    if(verbose>1) {printf("correcting data for user-defined Vb\n"); fflush(stdout);}
    /* Separate BTAC */
    TAC btac; tacInit(&btac);
    int ret=tacExtract(&input, &btac, 1);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(ret));
      tacFree(&ttac); tacFree(&input); parFree(&par); tacFree(&btac); return(6);
    }
    /* Interpolate BTAC to TTAC sample times */
    TAC bitac; tacInit(&bitac);
    ret=tacInterpolate(&btac, &ttac, &bitac, NULL, NULL, &status);
    tacFree(&btac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s.\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&input); parFree(&par); tacFree(&bitac); return(6);
    }
    /* Vb correction */
    ret=tacVb(&ttac, -1, &bitac, fixed_vb, 0, 1, &status);
    tacFree(&bitac);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&input); parFree(&par); return(6);
    }
    if(verbose>1) {printf("  corrected for Vb %g\n", fixed_vb); fflush(stdout);}
    tacFree(&btac);
  }
 
 
  /*
   *  Get good initial guess by calculating multi-linear solution
   */
  {
    if(verbose>2) {printf("multi-linear solution to get good initial guess\n"); fflush(stdout);}
    PAR mlpar; parInit(&mlpar);
    int ret;
    if(fixed_vb) ret=nnlsRTCM(&ttac, &input, 1, 0, &mlpar, &status);
    else ret=nnlsRTCM(&ttac, &input, 1, 1, &mlpar, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&ttac); tacFree(&input); parFree(&par); parFree(&mlpar); return(7);
    }
 
 
    parFree(&mlpar);
  }
 
 
  /*
   *  Fit model to the TTACs
   */
  if(verbose==1) {printf("\n fitting...\n"); fflush(stdout);}
  for(int ri=0; ri<ttac.tacNr; ri++) {
    if(verbose>1) {printf("\n fitting %s\n", ttac.c[ri].name); fflush(stdout);}
 
    // to be written
 
  } // next TTAC
 
 
 
 
  tacFree(&ttac); tacFree(&input); parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_1tcm(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);
}
/*****************************************************************************/
 
/*****************************************************************************/