fit_bpr.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <time.h>
/*****************************************************************************/
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpccurveio.h"
#include "libtpcsvg.h"
#include "libtpcmodext.h"
/*****************************************************************************/
 
/*****************************************************************************/
#define MAX_FUNC_NR 4
enum {MODEL_UNKNOWN, MODEL_IGV, MODEL_HILLB, MODEL_IHILLB};
int model=MODEL_UNKNOWN;
int func_par_nr[MAX_FUNC_NR]={0,5,5,5};
static char *model_name[] = {"Unknown", "IGV", "HILLB", "IHILLB", 0};
/*****************************************************************************/
#define MAX_MINMEAS_NR 3
enum {MINMEAS_UNKNOWN, MINMEAS_SS, MINMEAS_ABS};
int min_meas=MINMEAS_UNKNOWN;
/*****************************************************************************/
double *x, *ymeas, *yfit, *w; /* Local */
int dataNr=0, parNr=0;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double lastWSS;
/*****************************************************************************/
/* Local functions */
double (*func)(int, double*, void*);
double func_igv(int parNr, double *p, void*);
double func_hillb(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of a function to sampled blood-to-plasma ratio data.",
  " ",
  "Usage: @P [Options] datafile fitfile",
  " ",
  "Options:",
  " -model=<igv|hillb|ihillb>",
  "     Select the function; default is IGV.",
  " -delay[=<time>]",
  "     Fit also delay time, or, set delay time (0 by default).",
  " -bl=<baseline>",
  "     Constrain baseline level to specified value; fitted by default.",
  " -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.",
  " -fast or -safe",
  "     Speed up the fitting but increase the chance of failure, or",
  "     increase the reliability at the cost of computing time.",
  " -min=<ss|abs>",
  "     Sum-of-squares (SS) is minimized by default, but optionally",
  "     sum of absolute deviations can be selected.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Datafile must contain at least 2 columns, sample times and ratios;",
  "DFT file (1) can include also weights.",
  "Program writes the fit start and end times, nr of points, WSS/n,",
  "and parameters of the fitted function to the FIT file (2).",
  " ",
  "Example:",
  "     @P -model=igv -svg=ia765bprat.svg ia765bp.rat ia765bprat.fit",
  " ",
  "References:",
  "1. https://www.turkupetcentre.net/petanalysis/format_tpc_dft.html",
  "2. https://www.turkupetcentre.net/petanalysis/format_tpc_fit.html",
  " ",
  "See also: fit2dat, avgfract, b2plasma, bpr2cpr, tacinv, taccat, tacblend",
  " ",
  "Keywords: curve fitting, input, blood, plasma, modelling, simulation",
  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;
  int      fi, ri, pi, type, ret;
  int      weights=0; // 0=default, 1=no weighting, 2=frequency
  double   fixedDelay=0.0; // Fitted, if <0
  double   fixedBaseline=-1.0; // Fitted, if <0
  char     dftfile[FILENAME_MAX], fitfile[FILENAME_MAX], svgfile[FILENAME_MAX];
  char    *cptr;
  double   tstart, tstop, miny, maxy, par_tmp[MAX_PARAMETERS];
  int      tgo_nr, neighNr=0, reliability_level=2; /* 1=fast, 2=normal, 3=slow */
  DFT      dft;
  FIT      fit;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  dftInit(&dft); fitInit(&fit);
  dftfile[0]=fitfile[0]=svgfile[0]=(char)0;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    if(strcasecmp(cptr, "SAFE")==0) {
      reliability_level=3; continue;
    } else if(strcasecmp(cptr, "NORMAL")==0) {
      reliability_level=2; continue;
    } else if(strcasecmp(cptr, "FAST")==0) {
      reliability_level=1; continue;
    } else 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) {
      strcpy(svgfile, cptr+4); continue;
    } else if(strncasecmp(cptr, "MODEL=", 6)==0) {
      cptr+=6;
      if(strncasecmp(cptr, "IGV", 3)==0) {model=MODEL_IGV; continue;}
      if(strncasecmp(cptr, "HILLB", 3)==0) {model=MODEL_HILLB; continue;}
      if(strncasecmp(cptr, "IHILLB", 3)==0) {model=MODEL_IHILLB; continue;}
      fprintf(stderr, "Error: invalid ratio model.\n"); return(1);
    } else if(strncasecmp(cptr, "MIN=", 4)==0) {
      cptr+=4;
      if(strcasecmp(cptr, "SS")==0 || strcasecmp(cptr, "WSS")==0) {
        min_meas=MINMEAS_SS; continue;}
      if(strcasecmp(cptr, "ABS")==0 || strcasecmp(cptr, "WABS")==0) {
        min_meas=MINMEAS_ABS; continue;}
      fprintf(stderr, "Error: invalid minimization measure.\n"); return(1);
    } else if(strncasecmp(cptr, "DELAY=", 6)==0) {
      fixedDelay=atof_dpi(cptr+6); if(fixedDelay>=0.0) continue;
    } else if(strcasecmp(cptr, "DELAY")==0) {
      fixedDelay=-1.0; continue;
    } else if(strncasecmp(cptr, "BL=", 3)==0) {
      fixedBaseline=atof_dpi(cptr+3); if(fixedBaseline>=0.0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
  
  /* 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 */
  for(; ai<argc; ai++) {
    if(!dftfile[0]) {strcpy(dftfile, argv[ai]); continue;}
    if(!fitfile[0]) {strcpy(fitfile, argv[ai]); continue;}
    fprintf(stderr, "Error: too many arguments: '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Is something missing? */
  if(!fitfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(model==MODEL_UNKNOWN) model=MODEL_IGV;
  if(min_meas==MINMEAS_UNKNOWN) min_meas=MINMEAS_SS;
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("dftfile := %s\n", dftfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("model := %d\n", model);
    printf("min_meas := %d\n", min_meas);
    printf("reliability_level := %d\n", reliability_level);
    printf("weights := %d\n", weights);
    if(fixedDelay>=0) printf("fixedDelay := %g\n", fixedDelay);
    if(fixedBaseline>=0) printf("fixedBaseline := %g\n", fixedBaseline);
  }
  if(verbose>4) MATHFUNC_TEST=verbose-4; else MATHFUNC_TEST=0;
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading '%s'\n", dftfile);
  if(dftRead(dftfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", dftfile, dfterrmsg);
    return(2);
  }
  if(dft.frameNr<3) {
    fprintf(stderr, "Error: not enough samples for decent fitting.\n");
    dftEmpty(&dft); return(3);
  }
  dataNr=dft.frameNr;
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
  if(verbose>30) dftPrint(&dft);
 
  /* Check for NA's */
  if(dft_nr_of_NA(&dft) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", dftfile);
    dftEmpty(&dft); return(2);
  }
 
  /* Get min and max X and Y */
  ret=dftMinMax(&dft, &tstart, &tstop, &miny, &maxy);
  if(ret!=0) {
    fprintf(stderr, "Error: invalid contents in %s\n", dftfile);
    dftEmpty(&dft); return(2);
  }
  if(tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", dftfile);
    dftEmpty(&dft); return(2);
  }
  if(tstart<0.0) fprintf(stderr, "Warning: negative x value(s).\n");
  if(miny<0.0) fprintf(stderr, "Warning: negative y value(s).\n");
 
  /* Check and set weights */
  if(weights==0) {
    if(dft.isweight==0) {
      dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;}
  } else if(weights==1) {
    dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;
  } else if(weights==2) {
    dftWeightByFreq(&dft);
  }
  if(verbose==10) dftPrint(&dft);
 
 
  /*
   *  Allocate memory for fits
   */
  if(verbose>1) printf("allocating memory for fits.\n");
  ret=fit_allocate_with_dft(&fit, &dft); if(ret) {
    fprintf(stderr, "Error: cannot allocate space for fit results (%d).\n", ret);
    dftEmpty(&dft); fitEmpty(&fit); return(4);
  }
  /* Set necessary information in fit data structure */
  fit.voiNr=dft.voiNr;
  strncpy(fit.datafile, dftfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 256);
  fit.time=time(NULL);
  parNr=func_par_nr[model];
  if(model==MODEL_IGV) type=1402;
  else if(model==MODEL_HILLB) type=844;
  else if(model==MODEL_IHILLB) type=844;
  else type=0;
  for(ri=0; ri<fit.voiNr; ri++) {
    fit.voi[ri].type=type;
    fit.voi[ri].parNr=parNr;
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=dataNr;
  }
  if(verbose>8) fitPrint(&fit);
 
  /*
   *  Set fitting parameters
   */
  switch(reliability_level) {
    case 1: tgo_nr=150; neighNr=4; break;
    case 3: tgo_nr=1200; neighNr=6; break;
    case 2:
    default: tgo_nr=300; neighNr=5; break;
  }
  /* Set constraints and function */
  if(fixedDelay>=0.0) pmin[0]=pmax[0]=fixedDelay;
  else {pmin[0]=0.0; pmax[0]=0.10*tstop;}
  switch(model) {
    case MODEL_IGV:
      func=func_igv;
      /* baseline, BL/PL water contents ratio */
      if(fixedBaseline>=0.0) {pmin[1]=pmax[1]=fixedBaseline;}
      else {pmin[1]=0.80; pmax[1]=1.20;}
      /* scale */ 
      pmin[2]=0.001; pmax[2]=1.0;
      /* exponent */
      pmin[3]=0.0; pmax[3]=1.0;
      /* exp div */
      pmin[4]=1.0e-05; pmax[4]=1.0e+03;
      break;
    case MODEL_HILLB:
      func=func_hillb;
      /* scale */ 
      pmin[1]=0.0; pmax[1]=5.0;
      /* exponent */
      pmin[2]=1.0; pmax[2]=20.0;
      /* steepness */
      pmin[3]=1.0E-6; pmax[3]=1.0E+8;
      /* backround */
      if(fixedBaseline>=0.0) {pmin[4]=pmax[4]=fixedBaseline;}
      else {pmin[4]=0.5; pmax[4]=1.0;}
      break;
    case MODEL_IHILLB:
      func=func_hillb;
      /* scale */ 
      pmin[1]=-5.0; pmax[1]=0.0;
      /* exponent */
      pmin[2]=1.0; pmax[2]=20.0;
      /* steepness */
      pmin[3]=1.0E-6; pmax[3]=1.0E+8;
      /* backround */
      if(fixedBaseline>=0.0) {pmin[4]=pmax[4]=fixedBaseline;}
      else {pmin[4]=0.5; pmax[4]=5.0;}
      break;
  }
 
  /*
   *  Fit each TAC
   */
  if(verbose>0) {printf("fitting\n"); fflush(stdout);}
  for(ri=0; ri<dft.voiNr; ri++) {
 
    if(verbose>2) {
      fprintf(stdout, " fitting %s: ", dft.voi[ri].name);
      fflush(stdout);
    } else if(dft.voiNr>1 && verbose>0) {
      fprintf(stdout, "."); fflush(stdout);
    }
 
    /* Set data pointers */
    x=dft.x; ymeas=dft.voi[ri].y; yfit=dft.voi[ri].y2;
    w=dft.voi[ri].y3; for(fi=0; fi<dataNr; fi++) w[fi]=dft.w[fi];
 
    /* Set the initial values for parameters */
    fit.voi[ri].wss=3.402823e+38;
 
    /* Fit */
    ret=tgo(pmin, pmax, func, NULL, parNr, neighNr,
            &fit.voi[ri].wss, fit.voi[ri].p, tgo_nr, 0, verbose-5);
    if(ret) {
      fprintf(stderr, "Error %d in TGO.\n", ret);
      dftEmpty(&dft); fitEmpty(&fit); return(5);
    }
    fit.voi[ri].wss=lastWSS; // non-penalized WSS
    /* Correct fitted parameters to match constraints like inside the function */
    (void)modelCheckParameters(parNr, pmin, pmax, fit.voi[ri].p, fit.voi[ri].p, 
                               NULL);
 
    /* Warn user about parameters that collide with limits */
    if(verbose>1) {
      ret=modelCheckLimits(parNr, pmin, pmax, fit.voi[ri].p);
      if(ret==0) {if(verbose>2) fprintf(stdout, "ok\n");}
      else fprintf(stderr, "warning, fit collided with the limits.\n");
    }
    /*
     *  Print additional information in verbose mode
     */
    if(verbose>1) {
      double wss, a, aic;
      int n, k;
      char tmp[64];
      if(dft.voiNr==1) strcpy(tmp, ""); else sprintf(tmp, "[%d]", ri+1);
      /* Calculate WSS */
      wss=0.0;
      for(fi=0; fi<dataNr; fi++) {
        a=ymeas[fi]-yfit[fi];
        wss+=w[fi]*a*a;
      }
      printf("WSS%s := %g\n", tmp, wss);
      /* Nr of fitted parameters */
      for(pi=k=0; pi<parNr; pi++) if(pmax[pi]>pmin[pi]) k++;
      printf("Nr_of_fitted_parameters%s := %d\n", tmp, k);
      /* Nr of data points */
      for(fi=n=0; fi<dataNr; fi++) if(w[fi]>0.0) n++;
      printf("Nr_of_fitted_samples%s := %d\n", tmp, n);
      /* Calculate AIC */
      aic=aicSS(wss, n, k);  printf("AIC%s := %g\n", tmp, aic);
    }
 
    /* Print measured and fitted curve */
    if(verbose>5) {
      printf("\n  Frame   Time        Value       Fitted     Weight \n");
      for(fi=0; fi<dataNr; fi++)
        printf("  %02d   %8.3f  %12.4f  %12.4f  %8.4f\n",
          fi+1, x[fi], ymeas[fi], yfit[fi], w[fi]);
    }
    if(verbose>3) {
      printf(" fitted parameters:\n");
      for(pi=0; pi<fit.voi[ri].parNr; pi++)
        printf("  p[%d] := %g  [%g, %g]\n",
          pi+1, fit.voi[ri].p[pi], pmin[pi], pmax[pi]);
    }
  } /* next TAC */
  if(verbose>0) fprintf(stdout, "\n");
 
 
  /*
   *  Set fitted parameters as required by FIT format
   */
  for(ri=0; ri<dft.voiNr; ri++) {
    for(pi=0; pi<parNr; pi++) par_tmp[pi]=fit.voi[ri].p[pi];
    switch(model) {
      case MODEL_IGV:
        fit.voi[ri].p[0]=-par_tmp[2];
        fit.voi[ri].p[1]= par_tmp[3];
        fit.voi[ri].p[2]= par_tmp[4];
        fit.voi[ri].p[3]= par_tmp[0];
        fit.voi[ri].p[4]= par_tmp[1];
        break;
      case MODEL_HILLB:
      case MODEL_IHILLB:
        fit.voi[ri].p[0]=par_tmp[1];
        fit.voi[ri].p[1]=par_tmp[2];
        fit.voi[ri].p[2]=par_tmp[3];
        fit.voi[ri].p[3]=par_tmp[4];
        fit.voi[ri].p[4]=par_tmp[0];
        break;
    }
  }
 
 
  /*
   *  Print fit results on screen
   */
  if(verbose>0 && fit.voiNr<=50) fitWrite(&fit, "stdout");
 
  /*
   *  Save fit results
   */
  if(verbose>1) printf("saving results in %s\n", fitfile);
  if(fitWrite(&fit, fitfile)) {
    fprintf(stderr, "Error in writing '%s': %s\n", fitfile, fiterrmsg);
    dftEmpty(&dft); fitEmpty(&fit); return(11);
  }
 
  /*
   *  Saving and/or plotting of fitted TACs
   */
  if(svgfile[0]) {
    if(verbose>1) printf("saving SVG plot\n");
 
    /* Create a DFT containing fitted TACs */
    char tmp[128];
    DFT dft2;
    dftInit(&dft2); ret=dftdup(&dft, &dft2);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&dft); fitEmpty(&fit);
      return(21);
    }
    for(ri=0; ri<dft.voiNr; ri++) for(fi=0; fi<dft.frameNr; fi++)
      dft2.voi[ri].y[fi]=dft2.voi[ri].y2[fi];
 
    /* Save SVG plot of fitted and original data */
    tpcProgramName(argv[0], 0, 0, tmp, 64); strcat(tmp, " ");
    strcat(tmp, model_name[model]); strcat(tmp, " ");
    if(strlen(dft2.studynr)>1) strcat(tmp, dft.studynr);
    ret=plot_fit_svg(&dft, &dft2, tmp, svgfile, verbose-8);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
      dftEmpty(&dft2); dftEmpty(&dft); fitEmpty(&fit);
      return(30+ret);
    }
    if(verbose>0) printf("Plots written in %s\n", svgfile);
 
    dftEmpty(&dft2);
  }
 
  /* Free memory */
  dftEmpty(&dft); fitEmpty(&fit);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_igv(int parNr, double *p, void *fdata)
{
  int i;
  double xt, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
 
  /* Calculate the curve and WSS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-pa[0];
    if(xt<0.0) {
      yfit[i]=0.0;
    } else {
      yfit[i]=pa[2]*pow(xt, pa[3])*exp(-xt/pa[4]);
    }
    yfit[i]=pa[1]-yfit[i];
    /* Calculate weighted SS or ABS */
    v=yfit[i]-ymeas[i];
    if(min_meas==MINMEAS_SS) v*=v; else v=fabs(v);
    wss+=w[i]*v;
  }
  lastWSS=wss;
  wss*=penalty;
  return(wss);
}
double func_hillb(int parNr, double *p, void *fdata)
{
  int i;
  double xt, v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
  if(0) {
    for(i=0; i<parNr; i++) printf(" p[%d]=%g", i, p[i]);
    printf("\n");
  }
 
  /* Calculate the curve and WSS */
  for(i=0, wss=0.0; i<dataNr; i++) {
    xt=x[i]-pa[0];
    if(xt<=0.0) {
      yfit[i]=0.0;
    } else {
      v=pow(xt, pa[2]);
      yfit[i]=pa[1]*v/(v+pa[3]);
    }
    yfit[i]+=pa[4];
    /* Calculate weighted SS or ABS */
    v=yfit[i]-ymeas[i];
    if(min_meas==MINMEAS_SS) v*=v; else v=fabs(v);
    wss+=w[i]*v;
  }
  lastWSS=wss;
  wss*=penalty;
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/