fit_suri.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_suri(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *ymeas;
  double       *ysim;
  double       *w;
} FITDATA;
const int maxParNr=4;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the Surge function plus integral (recirculation)",
  "to PET time-activity curves (TACs).",
  " ",
  "Function:",
  "  when x<=dT, f(x) = 0 ",
  "  when x>dT,  f(x) = A*(x-dT)*exp(-B*(x-dT)) + ",
  "                     (C*A/B^2)*(1-exp(-B*(x-dT))-B*(x-dT)*exp(-B*(x-dT)))",
  ", where A>=0, B>=0, C>=0, and dT>=0.",
  " ",
  "Usage: @P [Options] tacfile [parfile]",
  " ",
  "Options:",
  " -w1 | -wf",
  "     All weights are set to 1.0 (no weighting), or weights are based on",
  "     sampling interval. By default, weights in tacfile are used, if available.",
  " -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.",
  "     Setting feasible limits are highly recommended.",
  " -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_ratf, fit2dat",
  " ",
  "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;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[0]=svgfile[0]=limfile[0]=(char)0;
  model=modelCodeIndex("surgerecircd");
  /* 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(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);
    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, "A");  plim.n[0].lim1=0.0;      plim.n[0].lim2=1.0E+04;
  strcpy(plim.n[1].name, "B");  plim.n[1].lim1=1.0E-06;  plim.n[1].lim2=1.0E+01;
  strcpy(plim.n[2].name, "C");  plim.n[2].lim1=1.0E-06;  plim.n[2].lim2=1.0E+03;
  strcpy(plim.n[3].name, "D");  plim.n[3].lim1=0.0;      plim.n[3].lim2=1.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);
 
 
 
  /*
   *  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 structure 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);
    }
    nlo._fun=func_suri;
    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(1)
    if(nloptSimplexARRS(&nlo, 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_suri(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[3];
    if(x>0.0) {
      double a=p[0];
      double b=p[1];
      double c=p[2];
      double e=exp(-b*x);
      d->ysim[i]=x*e + (c/(b*b))*(1.0-(b*x+1.0)*e);
      d->ysim[i]*=a;
    }
    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);
}
/*****************************************************************************/
 
/*****************************************************************************/