spillinp.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcfileutil.h"
#include "tpcli.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "General geometrical model for cardiac spillover correction, based on",
  "the following equations describing the regional TACs of myocardial muscle",
  "and cavity with blood (Cb) and muscle tissue (Ct) curves:",
  " ",
  "Croi,myo(t) = fBV*Cb(t) + (1-fBV)*Ct(t) ",
  "Croi,cav(t) = beta*Cb(t) + (1-beta)*Ct(t) ",
  " ",
  ", where fBV<beta. The spillover corrected BTAC can be calculated as",
  " ",
  "        (1-fBV)*Croi,cav(t) + (1-beta)*Croi,myo(t)",
  "Cb(t) = ------------------------------------------",
  "                beta*(1-fBV) - fBV*(1-beta)",
  " ",
  "and, optionally, the corrected tissue TAC can be calculated as",
  " ",
  "        beta*Croi,myo(t) - fBV*Croi,cav(t)",
  "Ct(t) = ----------------------------------",
  "            beta*(1-fBV) - fBV*(1-beta)",
  " ",
  "Radioactivity concentrations in regions-of-interest is given in TAC file,",
  "including at least the myocardial muscle and cavity ROIs.",
  "The names or numbers of these two ROIs are given next.",
  " ",
  "Usage: @P [Options] tacfile MYO CAV fBV beta cbfile [ctfile]",
  " ",
  "Options:",
  " -sim[ulate]",
  "     Instead of PVC, PVE is simulated from Ct(t) and Cb(t).",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "  @P aim25.tac myoc cavity 0.2 0.9 aim25ab.tac",
  " ",
  "References:",
  "1. Feng et al. (1996), https://doi.org/10.1109/10.486290",
  "2. Li et al. (1998), https://doi.org/10.1007/BF02522867",
  "3. Fang et al. (2008), https://doi.org/10.2967/jnumed.107.047613",
  " ",
  "See also: heartcor, tacadd, b2plasma, tacformat",
  " ",
  "Keywords: input, blood, PVC, myocardium, cavity, spill-over",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_ggm(int parNr, double *p, void*);
double func_ggm2(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int  nr;
  double       *xr;
  double       *ymyo;
  double       *ycav;

  unsigned int  nb;
  double       *xb;
  double       *yb;

  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
int ggmEstim(
  TAC *rtac,
  const int iMyo,
  const int iCav,
  TAC *stac,
  double start,
  int method,
  double *fBV,
  double *beta,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s(*rtac, %d, %d, *stac, %g, %d)\n", __func__, iMyo, iCav, start, method);
  if(rtac==NULL || stac==NULL || fBV==NULL || beta==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  if(rtac->tacNr<2 || rtac->sampleNr<2 || stac->tacNr<1 || stac->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
  if(iMyo<0 || iMyo>=rtac->tacNr || iCav<0 || iCav>=rtac->tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
 
  /* Check/set start time */
  double rx1, rx2, sx1, sx2;
  if(tacXRange(rtac, &rx1, &rx2) || tacXRange(stac, &sx1, &sx2)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
    return(TPCERROR_INVALID_XRANGE);
  }
  if(!(start>0.0)) {
    if(verbose>1) printf("setting start time\n");
    if(method==0) {
      start=0.5*(rx1+rx2);
      if(sx1<start) start=sx1;
    } else {
      if(verbose>1) printf("finding peak\n");
      int i;
      double ymax;
      if(tacYRange(rtac, iCav, NULL, &ymax, NULL, &i, NULL, NULL)) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_XRANGE);
        return(TPCERROR_INVALID_XRANGE);
      }
      start=rtac->x[i];
      if(verbose>1) printf("Cavity peak %g at %g (frame %d)\n", ymax, start, i);
    }
  }
  if(verbose>1) {
    printf("ROI data range: %g - %g\n", rx1, rx2);
    printf("Blood sample data range: %g - %g\n", sx1, sx2);
    printf("Start time: %g\n", start);
  }
  /* Extract the data before start time */
  TAC ftac; tacInit(&ftac);
  int ret=tacExtractRange(rtac, &ftac, start, rx2);
  if(ret==TPCERROR_OK && ftac.sampleNr<2) ret=TPCERROR_INVALID_XRANGE;
  if(ret!=TPCERROR_OK) {
    tacFree(&ftac);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
 
  /* Set data for the fit */
  FITDATA fitdata;
  fitdata.verbose=0;
  fitdata.nr=ftac.sampleNr;
  fitdata.xr=ftac.x;
  fitdata.ymyo=ftac.c[iMyo].y;
  fitdata.ycav=ftac.c[iCav].y;
  fitdata.nb=stac->sampleNr;
  fitdata.xb=stac->x;
  fitdata.yb=stac->c[0].y;
 
  /* Set parameters for non-linear optimisation */
  int parNr;
  if(method==0) {
    parNr=2; // beta and fb; fBV=fb*beta
  } else {
    parNr=8; // beta, fb (fBV=fb*beta), and Feng function parameters
  }
  NLOPT nlo; nloptInit(&nlo);
  ret=nloptAllocate(&nlo, parNr);
  if(ret!=TPCERROR_OK) {
    tacFree(&ftac);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  /* Set function */
  nlo.fundata=&fitdata;
  nlo.totalNr=parNr;
  {
    nlo.xlower[0]=0.1; //0.1; 
    nlo.xupper[0]=1.0; //1.0; 
    nlo.xfull[0]=0.5; //0.45;
    nlo.xdelta[0]=0.001;
    nlo.xtol[0]=0.0001;
  }
  {
    nlo.xlower[1]=0.0; //0.0; 
    nlo.xupper[1]=0.9; //0.9
    nlo.xfull[1]=0.5; //0.69;
    nlo.xdelta[1]=0.001;
    nlo.xtol[1]=0.0001;
  }
  if(method==0) {
    nlo._fun=func_ggm;
    nlo.maxFunCalls=20000;
  } else {
    nlo._fun=func_ggm2;
    nlo.maxFunCalls=50000;
    int i=2; // A1
    nlo.xlower[i]=0.0;
    nlo.xupper[i]=5000.0; 
    nlo.xfull[i]=1300.0;
    nlo.xdelta[i]=50.0;
    nlo.xtol[i]=0.01;
    i=4; // A2
    nlo.xlower[i]=0.0;
    nlo.xupper[i]=100.0; 
    nlo.xfull[i]=7.5;
    nlo.xdelta[i]=0.5;
    nlo.xtol[i]=0.01;
    i=6; // A3
    nlo.xlower[i]=0.0;
    nlo.xupper[i]=100.0; 
    nlo.xfull[i]=6.8;
    nlo.xdelta[i]=0.4;
    nlo.xtol[i]=0.01;
    i=3; // L1
    nlo.xlower[i]=-8.;
    nlo.xupper[i]=0.0; 
    nlo.xfull[i]=-6.7;
    nlo.xdelta[i]=0.1;
    nlo.xtol[i]=0.001;
    i=5; // f*L1
    nlo.xlower[i]=0.0;
    nlo.xupper[i]=0.9; 
    nlo.xfull[i]=0.025;
    nlo.xdelta[i]=0.002;
    nlo.xtol[i]=0.0001;
    i=7; // g*f*L1
    nlo.xlower[i]=0.0;
    nlo.xupper[i]=0.9; 
    nlo.xfull[i]=0.1;
    nlo.xdelta[i]=0.02;
    nlo.xtol[i]=0.001;
  }
 
#if(0)
  // call function for testing
  fitdata.verbose=10;
  double ss=func_ggm(nlo.totalNr, nlo.xfull, &fitdata);
  printf("ss=%g\n", ss);
  nloptFree(&nlo);
  *fBV=0.25; *beta=0.50; // placeholder
#else
  /* Fit */
  if(verbose>2) {
    printf("initial guess\n");
    nlo.funval=nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    nloptWrite(&nlo, stdout);
    fflush(stdout);
  }
 
  //fitdata.verbose=10;
  //ret=nloptSimplexMS(&nlo, 500, status);
  //ret=nloptIATGO(&nlo, 1, 1000, 0.15, status);
  ret=nloptIATGO(&nlo, 1, 0, 0.15, status);
  if(ret!=TPCERROR_OK) {
    tacFree(&ftac); nloptFree(&nlo);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  *beta=nlo.xfull[0];
  *fBV=nlo.xfull[1]*nlo.xfull[0];
  if(verbose>1) {printf("  SS := %g\n", nlo.funval); fflush(stdout);}
#endif
 
  tacFree(&ftac); nloptFree(&nlo);
 
  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   tacfile[FILENAME_MAX], bfile[FILENAME_MAX], tfile[FILENAME_MAX];
  char   sfile[FILENAME_MAX];
  char   myocname[MAX_TACNAME_LEN+1], cavname[MAX_TACNAME_LEN+1];
  double fBV=nan(""), beta=nan("");
  int    mode=0; // 0=correct, 1=simulate PVE, 2=fit
  int    method=0; // fit method
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=bfile[0]=tfile[0]=sfile[0]=(char)0;
  myocname[0]=cavname[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(strncasecmp(cptr, "SIMULATE", 3)==0) {
      mode=1; continue;
    } else if(strncasecmp(cptr, "FFIT", 2)==0) {
      method=1; continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-2;
 
  /* 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);}
 
  /* The (non-option) arguments */
  if(ai<argc) {strlcpy(tacfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(myocname, argv[ai], MAX_TACNAME_LEN); ai++;}
  if(ai<argc) {strlcpy(cavname, argv[ai], MAX_TACNAME_LEN); ai++;}
  if(ai>argc) {fprintf(stderr, "Error: missing command-line argument(s).\n"); return(1);}
 
  if(ai<argc) { // next one or two arguments: either sample file, or fBV and beta
    if(!atofCheck(argv[ai], &fBV) && !atofCheck(argv[ai+1], &beta)) {
      if(!(fBV>=0.0 && fBV<=1.0)) {
        fprintf(stderr, "Error: invalid Fbv command-line argument.\n"); return(1);}
      if(!(beta>=0.0 && beta<=1.0)) {
        fprintf(stderr, "Error: invalid Beta command-line argument.\n"); return(1);}
      /* Check that fBV<beta */
      if(fBV>=beta) {
        fprintf(stderr, "Error: fBV must be lower than beta.\n"); return(1);}
      ai+=2;
    } else {
      if(!fileExist(argv[ai])) {
        fprintf(stderr, "Error: missing sample file %s\n", argv[ai]); return(1);}
      strlcpy(sfile, argv[ai], FILENAME_MAX);
      mode=3;
      ai++;
    }
  }
  if(ai<argc) {strlcpy(bfile, argv[ai], FILENAME_MAX); ai++;}
  if(ai<argc) {strlcpy(tfile, argv[ai], FILENAME_MAX); ai++;}
  /* check that there are no extra arguments */
  if(ai<argc) {fprintf(stderr, "Error: extra command-line argument.\n"); return(1);}
  /* Is something missing? */
  if(!bfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]); 
    printf("\n");
    printf("tacfile := %s\n", tacfile);
    if(sfile[0]) printf("sfile := %s\n", sfile);
    if(isfinite(fBV)) printf("fBV := %g\n", fBV);
    if(isfinite(beta)) printf("beta := %g\n", beta);
    printf("mode := %d\n", mode);
    if(mode==3) printf("method := %d\n", method);
    printf("bfile := %s\n", bfile);
    if(tfile[0]) printf("tfile := %s\n", tfile);
  }
 
  /* Check that fBV<beta */
  if(fBV>=beta) {
    fprintf(stderr, "Error: fBV must be lower than beta.\n");
    return(1);
  }
 
 
  /*
   *  Read the data
   */
  if(verbose>1) printf("reading TACs\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(tac.tacNr<2) {
    fprintf(stderr, "Error: file does not contain myocardial and cavity ROI.\n");
    tacFree(&tac); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
 
  /*
   *  Find the myocardial and cavity ROIs among the TACs
   */
  int iMyoc=-1;
  {
    int n=tacSelectTACs(&tac, myocname, 1, &status);
    if(n>0) iMyoc=tacSelectBestReference(&tac);
    if(iMyoc<0) {
      fprintf(stderr, "Error: cannot identify myocardial ROI.\n");
      tacFree(&tac); return(3);
    }
  }
  int iCav=-1;
  {
    int n=tacSelectTACs(&tac, cavname, 1, &status);
    if(n>0) iCav=tacSelectBestReference(&tac);
    if(iCav<0) {
      fprintf(stderr, "Error: cannot identify cavity ROI.\n");
      tacFree(&tac); return(3);
    }
  }
  if(verbose>1) {
    printf("myocardial ROI: %s\n", tac.c[iMyoc].name);
    printf("cavity ROI: %s\n", tac.c[iCav].name);
  }
 
 
  /*
   *  If blood sample file is available, then estimate fBV and beta
   */
  if(mode==3) {
    if(verbose>1) printf("estimating PVC parameters\n");
    /* Read the blood sample(s) */
    TAC bs; tacInit(&bs);
    if(tacRead(&bs, sfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&bs); return(2);
    }
    if(bs.tacNr>1) {
      fprintf(stderr, "Note: using only the first concentration column in sample file.\n");
      bs.tacNr=1;
    }
    /* Check NaNs */
    if(tacNaNs(&bs)>0) {
      fprintf(stderr, "Error: blood sample file contains missing values.\n");
      tacFree(&tac); tacFree(&bs); return(2);
    }
    /* Sort data by sample time */
    tacSortByTime(&bs, &status);
    /* Make sure that blood sample data are in same units as ROI TAC data */
    if((bs.tunit!=tac.tunit && tacXUnitConvert(&bs, tac.tunit, &status)) ||
       (bs.cunit!=tac.cunit && tacYUnitConvert(&bs, tac.cunit, &status))) {
      fprintf(stderr, "Error: missing or incompatible data units.\n");
    }
 
    /* Estimation */
    if(ggmEstim(&tac, iMyoc, iCav, &bs, 0.0, method, &fBV, &beta, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&bs); return(9);
    }
 
 
    tacFree(&bs);
    if(verbose>1) {
      printf("fBV := %g\n", fBV);
      printf("beta := %g\n", beta);
    }
    mode=0; // So that BTAC will be calculated next
  }
 
 
  /*
   *  Calculate BTAC or cavity TAC
   */
  TAC btac; tacInit(&btac);
  {
    /* Setup place for BTAC */
    int cn=tac.tacNr;
    tac.tacNr=1;
    if(tacDuplicate(&tac, &btac)) {
      fprintf(stderr, "Error: cannot set up BTAC data.\n");
      tacFree(&tac); tacFree(&btac); return(4);
    }
    tac.tacNr=cn;
    if(mode==0) {
      /* Calculate BTAC */
      strcpy(btac.c[0].name, "Blood");
      double d=beta*(1.0-fBV)-fBV*(1.0-beta); if(verbose>2) printf("GGM-divider := %g\n", d);
      if(!(fabs(d)>1.0E-010)) {
        fprintf(stderr, "Error: invalid fBV and beta.\n");
        tacFree(&tac); tacFree(&btac); return(4);
      }
      for(int i=0; i<tac.sampleNr; i++)
        btac.c[0].y[i] = ((1.0-fBV)*tac.c[iCav].y[i] - (1.0-beta)*tac.c[iMyoc].y[i]) / d;
    } else {
      /* Calculate cavity TAC */
      strcpy(btac.c[0].name, "Cav");
      for(int i=0; i<tac.sampleNr; i++)
        btac.c[0].y[i] = beta*tac.c[iCav].y[i] + (1.0-beta)*tac.c[iMyoc].y[i];
    }
    /* Write file */
    FILE *fp; fp=fopen(bfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot write %s\n", bfile);
      tacFree(&tac); tacFree(&btac); 
      return(11);
    }
    int ret=tacWrite(&btac, fp, TAC_FORMAT_UNKNOWN, 0, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&btac); 
      return(11);
    }
    if(verbose>0) printf("saved %s\n", bfile);
    tacFree(&btac);
  }
 
  /*
   *  If requested, calculate also TTAC or myocardial TAC
   */
  if(tfile[0]) {
    /* Setup place for TTAC */
    TAC ttac; tacInit(&ttac);
    int cn=tac.tacNr;
    tac.tacNr=1;
    if(tacDuplicate(&tac, &ttac)) {
      fprintf(stderr, "Error: cannot set up TTAC data.\n");
      tacFree(&tac); tacFree(&ttac); return(5);
    }
    tac.tacNr=cn;
    if(mode==0) {
      /* Calculate TTAC */
      strcpy(ttac.c[0].name, "Tissue");
      double d=beta*(1.0-fBV)-fBV*(1.0-beta); //if(verbose>2) printf("GGM-divider := %g\n", d);
      if(!(fabs(d)>1.0E-010)) {
        fprintf(stderr, "Error: invalid fBV and beta.\n");
        tacFree(&tac); tacFree(&ttac); return(5);
      }
      for(int i=0; i<tac.sampleNr; i++)
        ttac.c[0].y[i] = (beta*tac.c[iMyoc].y[i] - fBV*tac.c[iCav].y[i]) / d;
    } else {
      /* Calculate Myo */
      strcpy(ttac.c[0].name, "Myo");
      for(int i=0; i<tac.sampleNr; i++)
        ttac.c[0].y[i] = fBV*tac.c[iCav].y[i] + (1.0-fBV)*tac.c[iMyoc].y[i];
    }
    /* Write file */
    FILE *fp; fp=fopen(tfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot write %s\n", tfile);
      tacFree(&tac); tacFree(&ttac); 
      return(12);
    }
    int ret=tacWrite(&ttac, fp, TAC_FORMAT_UNKNOWN, 0, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); tacFree(&ttac); 
      return(12);
    }
    if(verbose>0) printf("saved %s\n", tfile);
    tacFree(&ttac);
  }
 
  tacFree(&tac);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_ggm(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(d->verbose>2) printf("p[]: %g %g\n", p[0], p[1]);
 
  /* Process parameters */
  double beta, fBV;
  beta=p[0];
  fBV=p[1]*beta;
 
  /* GGM PVC */
  double blood[d->nr], logblood[d->nr];
  double v=beta*(1.0-fBV)-fBV*(1.0-beta); if(d->verbose>3) printf("GGM-divider := %g\n", v);
  if(!(fabs(v)>1.0E-010)) return(nan(""));
  for(unsigned int i=0; i<d->nr; i++) 
    blood[i] = ((1.0-fBV)*d->ycav[i] - (1.0-beta)*d->ymyo[i]) / v;
  /* Line fit to log-transformed BTAC */
  for(unsigned int i=0; i<d->nr; i++) logblood[i]=log(blood[i]);
  double m, c;
  int n=fitLine(d->xr, logblood, d->nr, &m, &c); 
  if(n<2) {
    if(d->verbose>4) {
      printf("  fitLine() -> %d\n", n);
      for(unsigned int i=0; i<d->nr; i++)
        printf("  %g %g %g\n", d->xr[i], blood[i], logblood[i]);
    }
    return(nan(""));
  }
  if(d->verbose>3) printf("  fitted line f(x) = %g*x + %g, n=%d\n", m, c, n);
  /* Calculate SS between blood samples and fitted BTAC */
  double ss=0.0;
  n=0;
  for(unsigned int i=0; i<d->nb; i++) {
    double b = c + m*d->xb[i];
    b = exp(b) - d->yb[i];
    if(isfinite(b)) {ss+=b*b; n++;}
  }
  if(n==0) return(nan(""));
  if(d->verbose>1) {
    for(int i=0; i<parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", ss);
  }
 
  return(ss);
}
/*****************************************************************************/
 
/*****************************************************************************/
double func_ggm2(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
 
  /* Process parameters */
  double beta, fBV;
  beta=p[0];
  fBV=p[1]*beta;
  double A1, A2, A3, L1, L2, L3;
  A1=p[2]; L1=p[3];
  A2=p[4]; L2=p[5]*L1;
  A3=p[6]; L3=p[7]*L2;
 
  /* GGM PVC */
  double blood[d->nr];
  double v=beta*(1.0-fBV)-fBV*(1.0-beta); if(d->verbose>3) printf("GGM-divider := %g\n", v);
  if(!(fabs(v)>1.0E-010)) return(nan(""));
  for(unsigned int i=0; i<d->nr; i++) 
    blood[i] = ((1.0-fBV)*d->ycav[i] - (1.0-beta)*d->ymyo[i]) / v;
  /* SS between Feng function and BTAC */
  double ss=0.0;
  for(unsigned int i=0; i<d->nr; i++) {
    double x=d->xr[i];
    double f=(A1*x-A2-A3)*exp(L1*x)+A2*exp(L2*x)+A3*exp(L3*x);
    f-=blood[i];
    ss+=f*f;
  }
  /* Add SS between Feng function and blood sample(s) */
  for(unsigned int i=0; i<d->nb; i++) {
    double x=d->xb[i];
    double f=(A1*x-A2-A3)*exp(L1*x)+A2*exp(L2*x)+A3*exp(L3*x);
    f-=d->yb[i];
    ss+=10.0*f*f;
  }
  if(d->verbose>1) {
    for(int i=0; i<parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", ss);
  }
 
  return(ss);
}
/*****************************************************************************/
 
/*****************************************************************************/