fit_av.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 "tpcpar.h"
#include "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_av(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int  in;
  double       *ix;
  double       *iy;
  double       *kernel;
  double       *cy;

  unsigned int n;
  double       *x;
  double       *y;
  double       *ysim;
  double       *w;
  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the parameters of compartmental model for",
  "A-V difference.",
  " ",
  "Usage: @P [Options] arterytac venatac [parfile]",
  " ",
  "Options:",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     output data file are used, if available.",
  " -wf",
  "     Weight by output TAC sampling interval.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "If TAC files contain more than one TAC, only the first is used.",
  " ",
  "See also: sim_av, fit_xexp, convexpf",
  " ",
  "Keywords: input, curve fitting",
  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 simC3vb(
  double *t,
  double *cab,
  const int nr,
  double f,
  double k1,
  double k2,
  double k3,
  double k4,
  double k5,
  double k6,
  double *cvb,
  double *ct
) {
  int i;
  double b, c, d, w, z, dt2;
  double cai, ca_last, t_last;
  double ct1, ct1_last, ct2, ct2_last, ct3, ct3_last;
  double ct1i, ct1i_last, ct2i, ct2i_last, ct3i, ct3i_last;
 
 
  /* Check for data */
  if(nr<2) return 1;
  if(cvb==NULL) return 2;
 
  /* Check actual parameter number */
  if(!(f>=0.0) || !(k1>=0.0) || k1>f) return 3;
  if(k3<=0.0) {k3=0.0; if(k2<=0.0) k2=0.0;}
  else if(k5<=0.0) {k5=0.0; if(k4<=0.0) k4=0.0;}
  else {if(k6<=0.0) k6=0.0;}
 
  /* Calculate curves */
  t_last=0.0; if(t[0]<t_last) t_last=t[0]; 
  cai=ca_last=0.0;
  ct1_last=ct2_last=ct3_last=ct1i_last=ct2i_last=ct3i_last=0.0;
  ct1=ct2=ct3=ct1i=ct2i=ct3i=0.0;
  for(i=0; i<nr; i++) {
    /* delta time / 2 */
    dt2=0.5*(t[i]-t_last);
    /* calculate values */
    if(dt2<0.0) {
      return 5;
    } else if(dt2>0.0) {
      /* arterial integral */
      cai+=(cab[i]+ca_last)*dt2;
      /* partial results */
      b=ct1i_last+dt2*ct1_last;
      c=ct2i_last+dt2*ct2_last;
      d=ct3i_last+dt2*ct3_last;
      w=k4 + k5 - (k5*k6*dt2)/(1.0+k6*dt2);
      z=1.0+w*dt2;
      /* 1st tissue compartment and its integral */
      ct1 = (
          + k1*z*cai + (k3*k4*dt2 - (k2+k3)*z)*b
          + k4*c + k4*k6*dt2*d/(1.0+k6*dt2)
        ) / ( z*(1.0 + dt2*(k2+k3)) - k3*k4*dt2*dt2 );
      ct1i = ct1i_last + dt2*(ct1_last+ct1);
      /* 2nd tissue compartment and its integral */
      ct2 = (k3*ct1i - w*c + k6*d/(1.0+k6*dt2)) / z;
      ct2i = ct2i_last + dt2*(ct2_last+ct2);
      /* 3rd tissue compartment and its integral */
      ct3 = (k5*ct2i - k6*d) / (1.0 + k6*dt2);
      ct3i = ct3i_last + dt2*(ct3_last+ct3);
    }
    if(f>0.0) {
      double dct = k1*cab[i] - k2*ct1; 
      cvb[i] = cab[i] - dct/f;
    } else
      cvb[i]=cab[i];
 
    /* copy values to argument arrays; set very small values to zero */
    if(ct!=NULL) {ct[i]=ct1+ct2+ct3; if(fabs(ct[i])<1.0e-12) ct[i]=0.0;}
 
    /* prepare to the next loop */
    t_last=t[i]; ca_last=cab[i];
    ct1_last=ct1; ct1i_last=ct1i;
    ct2_last=ct2; ct2i_last=ct2i;
    ct3_last=ct3; ct3i_last=ct3i;
  }
 
  return 0;
}

/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int     ai, help=0, version=0, verbose=1;
  char    afile[FILENAME_MAX], vfile[FILENAME_MAX], parfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  int     weights=0; // 0=default, 1=no weighting, 2=frequency
  int     model=0;
  int     ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  afile[0]=vfile[0]=parfile[0]=svgfile[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, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "SER2")==0) {
      model=1; continue;
    } else if(strcasecmp(cptr, "MET")==0) {
      model=2; 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(afile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(vfile, 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(!vfile[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("afile := %s\n", afile);
    printf("vfile := %s\n", vfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    printf("weights := %d\n", weights);
    printf("model := %d\n", model);
    fflush(stdout);
  }
 
 
  /*
   *  Read input and output data
   */
  if(verbose>1) printf("reading TAC data\n");
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  TAC atac, vtac; tacInit(&atac); tacInit(&vtac);
  int fitSampleNr;
  double fitdur=1.0E+12;
  ret=tacReadModelingData(afile, vfile, NULL, NULL, &fitdur, 0,
                          &fitSampleNr, &atac, &vtac, &status);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); tacFree(&vtac); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(atac.format));
    printf("a.sampleNr := %d\n", atac.sampleNr);
    printf("v.sampleNr := %d\n", vtac.sampleNr);
    printf("fitSampleNr := %d\n", fitSampleNr);
    printf("xunit := %s\n", unitName(atac.tunit));
    printf("yunit := %s\n", unitName(atac.cunit));
    printf("fitdur := %g s\n", fitdur);
  }
  if(fitSampleNr<4 || atac.sampleNr<4) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&atac); tacFree(&vtac); return(2);
  }
  /* Set places for simulated and fitted TACs */
  atac.tacNr=vtac.tacNr=1;
  if(tacAllocateMore(&atac, 1)!=TPCERROR_OK || tacAllocateMore(&vtac, 1)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); tacFree(&vtac); return(2);
  }
  /* Check and set weights */
  if(weights==0) {
    if(!tacIsWeighted(&vtac)) {
      vtac.weighting=WEIGHTING_OFF;
      for(int i=0; i<vtac.sampleNr; i++) vtac.w[i]=1.0;
    } 
  } else if(weights==1) {
    vtac.weighting=WEIGHTING_OFF;
    for(int i=0; i<vtac.sampleNr; i++) vtac.w[i]=1.0;
  } else if(weights==2) {
    ret=tacWByFreq(&vtac, ISOTOPE_UNKNOWN, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&atac); tacFree(&vtac); return(2);
    }
  }
 
 
  /*
   *  Fitting
   */
  if(verbose>1) printf("preparing for fitting\n");
 
  /* Set data pointers for the fit */
  FITDATA fitdata;
  fitdata.verbose=0;
  //fitdata.kernel=kernel;
  fitdata.in=atac.sampleNr;
  fitdata.ix=atac.x;
  fitdata.iy=atac.c[0].y;
  fitdata.cy=atac.c[1].y;
 
  fitdata.n=vtac.sampleNr;
  fitdata.x=vtac.x;
  fitdata.y=vtac.c[0].y;
  fitdata.ysim=vtac.c[1].y;
  fitdata.w=vtac.w;
 
 
 
  /* Set parameters for nonlinear optimization */
  if(verbose>2) printf("process initial parameter values\n");
  int parNr=4; if(model==1) parNr=6;
  NLOPT nlo; nloptInit(&nlo);
  ret=nloptAllocate(&nlo, parNr);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); tacFree(&vtac); return(3);
  }
  /* Set function */
  nlo._fun=func_av;
  nlo.fundata=&fitdata;
  nlo.totalNr=parNr; // delay, f, K1/f, K1/k2; k3, k3/k4
  { // delay [min]
    nlo.xlower[0]=-0.03*fitdur; 
    nlo.xupper[0]=0.1*fitdur; 
    nlo.xfull[0]=0.0;
    nlo.xdelta[0]=0.1;
    nlo.xtol[0]=0.005;
  }
  { // f [mL/(mL*min)]
    nlo.xlower[1]=0.01; 
    nlo.xupper[1]=3.0; 
    nlo.xfull[1]=0.10;
    nlo.xdelta[1]=0.05;
    nlo.xtol[1]=0.01;
  }
  { // K1/f
    nlo.xlower[2]=0.01; 
    nlo.xupper[2]=1.0; 
    nlo.xfull[2]=0.50;
    nlo.xdelta[2]=0.10;
    nlo.xtol[2]=0.01;
  }
  { // K1/k2
    nlo.xlower[3]=0.01; 
    nlo.xupper[3]=100.0; 
    nlo.xfull[3]=1.0;
    nlo.xdelta[3]=0.1;
    nlo.xtol[3]=0.01;
  }
  if(model==1) {
    // k3
    nlo.xlower[4]=0.0; 
    nlo.xupper[4]=1.0; 
    nlo.xfull[4]=0.1;
    nlo.xdelta[4]=0.05;
    nlo.xtol[4]=0.01;
    // k3/k4
    nlo.xlower[5]=1.0E-03; 
    nlo.xupper[5]=10.0; 
    nlo.xfull[5]=1.0;
    nlo.xdelta[5]=0.1;
    nlo.xtol[5]=0.02;
  }
  nlo.maxFunCalls=20000;
 
  /* Get nr of fixed parameters */
  unsigned int fixedNr=nloptFixedNr(&nlo);
  if(verbose>2) printf("fixedNr := %d\n", fixedNr); 
 
#if(0)
  // call function for testing
  fitdata.verbose=10;
  func_av(nlo.totalNr, nlo.xfull, &fitdata);
 
  nloptFree(&nlo);
  tacFree(&atac); tacFree(&vtac);
 
#else
 
  /* Fit */
  drandSeed(1);
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>0) printf("fitting\n");
  if(nloptITGO2(&nlo, 1, 0, 0, 0, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&atac); tacFree(&vtac); nloptFree(&nlo); return(5);
  }
  nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
  if(verbose>4) {
    printf("measured and fitted TAC:\n");
    for(int i=0; i<vtac.sampleNr; i++)
      printf("\t%g\t%g\t%g\n", vtac.x[i], vtac.c[0].y[i], vtac.c[1].y[i]);
  }
  if(verbose>2) nloptWrite(&nlo, stdout);
  if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
 
  /*
   *  Plot measured and fitted data, if requested
   */
  if(svgfile[0]) {
    if(verbose>1) printf("plotting measured and fitted data\n");
    for(int i=0; i<atac.sampleNr; i++) atac.c[0].y[i]=atac.c[1].y[i];
    /* Plot */
    ret=tacPlotFitSVG(&vtac, &atac, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&atac); tacFree(&vtac); nloptFree(&nlo); //parFree(&par);
      return(24);
    }
 
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
  }
 
 
 
  nloptFree(&nlo); //parFree(&par);
  tacFree(&atac); tacFree(&vtac);
 
#endif
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_av(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(d->verbose>20) printf("p[]: %g %g %g %g\n", p[0], p[1], p[2], p[3]);
 
  /* Process parameters */
  double delay, f, K1, k2, k3=0.0, k4=0.0;
  delay=p[0];
  f=p[1];
  K1=f*p[2];
  k2=K1/p[3];
  if(parNr==6) {
    k3=p[4];
    k4=k3/p[5];
  }
 
  /* Delay corrected sample times */
  double dx[d->in];
  for(unsigned int i=0; i<d->in; i++) dx[i]=d->ix[i]+delay;
 
  /* Simulate venous TAC */
  if(simC3vb(d->ix, d->iy, d->in, f, K1, k2, k3, k4, 0.0, 0.0, d->cy, NULL)) {
    fprintf(stderr, "Error: cannot simulate.\n");
    return(nan(""));
  }
 
  /* Interpolate to venous sample times */
  if(d->verbose>3) {fprintf(stdout, "interpolation...\n"); fflush(stdout);}
  if(liInterpolate(dx, d->cy, d->in, d->x, d->ysim, NULL, NULL, d->n, 0, 1, 0)!=0) {
    fprintf(stderr, "Error: cannot interpolate.\n");
    return(nan(""));
  }
  if(d->verbose>3) {
    printf("\nData x\ty\tSimulated\n");
    for(unsigned int i=0; i<d->n; i++)
      printf("%g\t%g\t%g\n", d->x[i], d->y[i], d->ysim[i]);
  }
 
  /* 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->n; i++) {
    double v=d->ysim[i] - d->y[i];
    wss+=d->w[i]*v*v;
  }
  if(d->verbose>1) {
    for(int i=0; i<parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/