fit_wcdf.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 "tpctacmod.h"
#include "tpcrand.h"
#include "tpcnlopt.h"
#include "tpcfileutil.h"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_wcdf(int parNr, double *p, void*);
double func_wcdfd(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *ymeas;
  double       *ysim;
  double       *w;
} FITDATA;
const int maxParNr=5;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the Weibull CDF and Weibull PDF sum function to PET",
  "time-activity curves (TACs).",
  " ",
  "Function:",
  "  when x<=dT, f(x) = 0 ",
  "  when x>dT,  f(x) = A * (n/B) * ((x-dT)/B)^(n-1) * exp(-((x-dT)/B)^n) +",
  "                     k * A * (1-exp(-((x-dT)/B)^n))",
  ", where A>=0, B>0, n>0 (usually n>1), k>=0, and dT>=0.",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -CDF",
  "     WCDF function without WPDF (derivative of WCDF) is fitted:",
  "     when x>dT, f(x) = A * (1-exp(-((x-dT)/B)^n))",
  " -w1",
  "     All weights are set to 1.0 (no weighting); by default, weights in",
  "     data file are used, if available.",
  " -wf",
  "     Weight by sampling interval.",
  " -lim[=<filename>]",
  "     Specify the constraints for function parameters;",
  "     This file with default values can be created by giving this option",
  "     as the only command-line argument to this program.",
  "     Without file name the default values are printed on screen.",
  " -delay=<time>",
  "     Delay time (dT) is constrained to specified value (>=0).",
  " -k=<value>",
  "     Parameter k is constrained to given value; setting k to zero causes",
  "     the function to approach zero.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: fit_sinf, fit_feng, fit_suri, fit_ratf, fit2dat, extrapol",
  " ",
  "Keywords: TAC, 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 main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char     tacfile[FILENAME_MAX], parfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  char     limfile[FILENAME_MAX];
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  unsigned int model=0;
  double   fixed_delay=nan("");
  double   fixed_k=nan("");
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=svgfile[0]=limfile[0]=(char)0;
  model=modelCodeIndex("weibullcdfdd");
  /* 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(strcasecmp(cptr, "CDF")==0 || strcasecmp(cptr, "WCDF")==0) {
      model=modelCodeIndex("weibullcdfd"); continue;
    } else if(strcasecmp(cptr, "CDFD")==0 || strcasecmp(cptr, "WCDFD")==0) {
      model=modelCodeIndex("weibullcdfdd"); continue;
    } else if(strncasecmp(cptr, "DELAY=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &fixed_delay) && fixed_delay>=0.0) continue;
    } else if(strncasecmp(cptr, "DT=", 3)==0 && strlen(cptr)>3) {
      if(!atofCheck(cptr+3, &fixed_delay) && fixed_delay>=0.0) continue;
    } else if(strncasecmp(cptr, "K=", 2)==0 && strlen(cptr)>2) {
      if(!atofCheck(cptr+2, &fixed_k) && fixed_k>=0.0) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "LIM=", 4)==0 && strlen(cptr)>4) {
      strlcpy(limfile, cptr+4, FILENAME_MAX); continue;
    } else if(strcasecmp(cptr, "LIM")==0) {
      strcpy(limfile, "stdout"); 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(tacfile, 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);
  }
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    if(tacfile[0]) printf("tacfile := %s\n", tacfile);
    if(limfile[0]) printf("limfile := %s\n", limfile);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(svgfile[0]) printf("svgfile := %s\n", svgfile);
    if(!isnan(fixed_k)) printf("fixed_k := %g\n", fixed_k);
    if(!isnan(fixed_delay)) printf("fixed_delay := %g\n", fixed_delay);
    printf("model := %s\n", modelCode(model));
    printf("weights := %d\n", weights);
  }
 
  /* Set parameter names and initial constraints */
  PAR plim; parInit(&plim); if(parAllocate(&plim, maxParNr, 1)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot initiate parameter list.\n"); return(1);}
  plim.parNr=maxParNr;
  strcpy(plim.n[0].name, "dT");  plim.n[0].lim1=0.0;      plim.n[0].lim2=100.0;
  strcpy(plim.n[1].name, "A");   plim.n[1].lim1=1.0E-06;  plim.n[1].lim2=1.0E+03;
  strcpy(plim.n[2].name, "B");   plim.n[2].lim1=1.0E-06;  plim.n[2].lim2=20.0;
  if(model==modelCodeIndex("weibullcdfdd")) {
  strcpy(plim.n[3].name, "n");   plim.n[3].lim1=0.5;      plim.n[3].lim2=4.0;
  strcpy(plim.n[4].name, "k");   plim.n[4].lim1=0.0;      plim.n[4].lim2=5.0;
  } else {
  strcpy(plim.n[3].name, "n");   plim.n[3].lim1=1.0E-03;  plim.n[3].lim2=4.0;
  strcpy(plim.n[4].name, "k");   plim.n[4].lim1=0.0;      plim.n[4].lim2=0.0;
  }
  /* If only file name for parameter constraints was given, then write one and exit */
  if(limfile[0] && !tacfile[0]) {
    /* Check that initial value file does not exist */
    if(strcasecmp(limfile, "stdout")!=0 && fileExist(limfile)) {
      fprintf(stderr, "Error: parameter constraint file %s exists.\n", limfile);
      parFree(&plim); return(1);
    }
    /* Create parameter file */
    plim.parNr=modelParNr(model);
    if(verbose>1 && strcasecmp(limfile, "stdout")!=0) 
      printf("writing parameter constraints file\n");
    if(parWriteLimits(&plim, limfile, verbose-2)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot write constraints file.\n");
      parFree(&plim); return(1);
    }
    if(strcasecmp(limfile, "stdout")!=0)
      fprintf(stdout, "Parameter constraints written in %s\n", limfile);
    parFree(&plim); return(0);
  }
 
  /* Did we get all the information that we need? */
  if(!tacfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    parFree(&plim); return(1);
  }
  if(strcasecmp(limfile, "stdout")==0) limfile[0]=(char)0;
 
 
  /*
   *  Read parameter constraints if file for that was given
   */
  if(limfile[0]) {
    if(parReadLimits(&plim, limfile, verbose-2)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot read constraints file.\n");
      parFree(&plim); return(1);
    }
    if(verbose>1) parListLimits(&plim, stdout);
  }
  /* Set optional user-defined constraints */
  if(!isnan(fixed_delay)) plim.n[0].lim1=plim.n[0].lim2=fixed_delay;
  if(!isnan(fixed_k)) plim.n[4].lim1=plim.n[4].lim2=fixed_k;
  if(verbose>1) parListLimits(&plim, stdout);
 
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); parFree(&plim); return(2);
  }
  if(verbose>1) {
    printf("fileformat := %s\n", tacFormattxt(tac.format));
    printf("tacNr := %d\n", tac.tacNr);
    printf("sampleNr := %d\n", tac.sampleNr);
    printf("xunit := %s\n", unitName(tac.tunit));
    printf("yunit := %s\n", unitName(tac.cunit));
    printf("isframe := %d\n", tac.isframe);
  }
  if(tac.tacNr<1) {
    fprintf(stderr, "Error: file contains no data.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Get x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); parFree(&plim); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
  /* Refine dT limits */
  if(!limfile[0] && plim.n[0].lim2>plim.n[0].lim1 && xmax<plim.n[0].lim2/0.75) 
    plim.n[0].lim2=0.75*xmax;
 
  /* Check and set weights */
  if(weights==0) {
    if(!tacIsWeighted(&tac)) {
      tac.weighting=WEIGHTING_OFF;
      for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
    } 
  } else if(weights==1) {
    tac.weighting=WEIGHTING_OFF;
    for(int i=0; i<tac.sampleNr; i++) tac.w[i]=1.0;
  } else if(weights==2) {
    if(tacWByFreq(&tac, ISOTOPE_UNKNOWN, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); return(2);
    }
  }
 
 
  /*
   *  Copy function parameters into PAR struct for printing and saving
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, maxParNr, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); parFree(&plim); return(3);
  }
  par.tacNr=tac.tacNr; par.parNr=modelParNr(model);
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<tac.tacNr; i++) {
    par.r[i].model=model;
    par.r[i].dataNr=tacWSampleNr(&tac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  for(int i=0; i<par.parNr; i++) strcpy(par.n[i].name, plim.n[i].name);
  /* 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);
  }
  /* set file names */
  iftPut(&par.h, "datafile", tacfile, 0, NULL);
 
 
  /*
   *  Fitting
   */
  if(verbose>1) printf("preparing for fitting\n");
 
  drandSeed(1);
 
  if(verbose>0) printf("fitting...\n");
  for(int ci=0; ci<tac.tacNr; ci++) {
    if(verbose>1) printf("TAC %d: %s\n", 1+ci, tac.c[ci].name);
    /* Get y range */
    double ymax; int ymaxi;
    if(tacYRange(&tac, ci, NULL, &ymax, NULL, &ymaxi, NULL, NULL)!=0) {
      fprintf(stderr, "Error: invalid y (concentration) values.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); return(5);
    }
    if(verbose>3) printf("  ymax := %g\n", ymax);
    /* Set data pointers for the fit */
    double yfit[tac.sampleNr];
    FITDATA fitdata;
    fitdata.n=tac.sampleNr;
    fitdata.x=tac.x;
    fitdata.ymeas=tac.c[ci].y;
    fitdata.ysim=yfit;
    fitdata.w=tac.w;
    /* Set NLLS options */
    NLOPT nlo; nloptInit(&nlo);
    if(nloptAllocate(&nlo, par.parNr)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); return(5);
    }
    if(model==modelCodeIndex("weibullcdfdd")) nlo._fun=func_wcdfd; else nlo._fun=func_wcdf;
    nlo.fundata=&fitdata;
    nlo.totalNr=par.parNr;
    for(int i=0; i<par.parNr; i++) {
      nlo.xlower[i]=plim.n[i].lim1;
      nlo.xupper[i]=plim.n[i].lim2;
    }
    for(int i=0; i<par.parNr; i++) {
      if(!(plim.n[i].lim2>plim.n[i].lim1)) nlo.xtol[i]=0.0;
      else nlo.xtol[i]=0.0001*(plim.n[i].lim2-plim.n[i].lim1);
    }
    nlo.maxFunCalls=50000;
    /* Fit */
    if(verbose>4) nloptWrite(&nlo, stdout);
#if(0)
    if(nloptITGO2(&nlo, 1, 0, 0, 0, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); nloptFree(&nlo); return(6);
    }
#else
    if(nloptIATGO(&nlo, 1, 0, 0.25, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); nloptFree(&nlo); return(6);
    }
#endif
    nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>2) nloptWrite(&nlo, stdout);
    if(verbose>6) {
      printf("measured and fitted TAC:\n");
      for(int i=0; i<tac.sampleNr; i++)
        printf("\t%g\t%g\t%g\n", tac.x[i], tac.c[ci].y[i], yfit[i]);
    }
    //if(verbose>2) printf(" wss := %g\n", nlo.funval);
 
    /* Copy parameters */
    for(int i=0; i<par.parNr; i++) par.r[ci].p[i]=nlo.xfull[i];
    par.r[ci].wss=nlo.funval;
 
    nloptFree(&nlo);
  } // next TAC
 
  /* Print and save parameters */
  if(verbose>0 || !parfile[0]) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  if(parfile[0]) {
    /* 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(&tac); parFree(&plim); parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_FIT, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); 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");
    /* Prepare place for function TAC */
    double sdist;
    tacGetSampleInterval(&tac, 1.0E-03, &sdist, NULL);
    int snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end-par.r[0].start, 2, NULL);
    TAC ftac; tacInit(&ftac);
    int ret=tacAllocate(&ftac, snr, tac.tacNr);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); return(21);
    }
    ftac.tacNr=tac.tacNr; ftac.sampleNr=snr;
    ret=tacCopyHdr(&tac, &ftac);
    for(int ci=0; ci<tac.tacNr && ret==0; ci++)
      ret=tacCopyTacchdr(tac.c+ci, ftac.c+ci);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(22);
    }
    snr=(int)simSamples(0.2*sdist, 0.0, par.r[0].end-par.r[0].start, 2, ftac.x);
    ftac.isframe=0;
    for(int i=0; i<ftac.sampleNr; i++) ftac.x[i]+=par.r[0].start;
    /* Compute function values at sample times */
    ret=0;
    for(int ci=0; ci<tac.tacNr && ret==0; ci++)
      ret=mfEvalY(modelCode(par.r[ci].model), par.parNr, par.r[ci].p, 
                  ftac.sampleNr, ftac.x, ftac.c[ci].y, verbose-5);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot calculate function values.\n");
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(23);
    }
    if(verbose>10) tacWrite(&ftac, stdout, TAC_FORMAT_UNKNOWN, 1, NULL);
    /* Plot */
    ret=tacPlotFitSVG(&tac, &ftac, "", nan(""), nan(""), nan(""), nan(""), svgfile, &status);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&plim); parFree(&par); tacFree(&ftac);
      return(24);
    }
    tacFree(&ftac);
    if(verbose>0) printf("Measured and fitted data plotted in %s\n", svgfile);
  }
 
 
  tacFree(&tac); 
  parFree(&plim);
  parFree(&par);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_wcdfd(int parNr, double *p, void *fdata)
{
  double v, wss=0.0;
  FITDATA *d=(FITDATA*)fdata;
 
  /* Calculate the curve values and weighted SS */
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=0.0;
    if(isnan(d->ymeas[i]) || isnan(d->x[i])) continue;
    double x=d->x[i]-p[0];
    if(x>0.0) {
      double et=exp(-pow(x/p[2], p[3]));
      double wpdf=(p[3]/p[2]) * pow(x/p[2], p[3]-1.0) * et;
      d->ysim[i]=p[1]*(wpdf + p[4]*(1.0-et));
    }
    v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
  if(0) {
    for(unsigned int i=0; i<(unsigned int)parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
  return(wss);
}
 
double func_wcdf(int parNr, double *p, void *fdata)
{
  double v, wss=0.0;
  FITDATA *d=(FITDATA*)fdata;
 
  /* Calculate the curve values and weighted SS */
  for(unsigned int i=0; i<d->n; i++) {
    d->ysim[i]=0.0;
    if(isnan(d->ymeas[i]) || isnan(d->x[i])) continue;
    double x=d->x[i]-p[0];
    if(x>0.0) {
      d->ysim[i]=p[1]*(1.0 - exp(-pow(x/p[2], p[3])));
    }
    v=d->ysim[i]-d->ymeas[i];
    wss+=d->w[i]*v*v;
  }
  if(0) {
    for(unsigned int i=0; i<(unsigned int)parNr; i++) printf(" %g", p[i]);
    printf(" -> %g\n", wss);
  }
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/