fit_ratf.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"
/*****************************************************************************/
double func_ratf(int parNr, double *p, void*);
/*****************************************************************************/
double *x, *ymeas, *yfit, *w;
int dataNr=0, parNr=0;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of 1st-3rd degree rational function to PET plasma and",
  "tissue time-activity curves (TACs).",
  " ",
  "Function:",
  " ",
  "         p1 + p3*x + p5*x^2 + p7*x^3  ",
  "  f(x) = ---------------------------  ",
  "         p2 + p4*x + p6*x^2 + p8*x^3  ",
  " ",
  ", where p1=0 and p2=1.",
  " ",
  "Usage: @P [Options] tacfile fitfile",
  " ",
  "Options:",
  " -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.",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Options for selecting the rational function:",
  " -p3p3  3rd degree polynomial/3rd degree polynomial (default)",
  " -p3p2  3rd degree polynomial/2nd degree polynomial",
  " -p2p2  2nd degree polynomial/2nd degree polynomial",
  " -p2p1  2nd degree polynomial/1st degree polynomial",
  " -p1p1  1st degree polynomial/1st degree polynomial",
  " ",
  "Program writes the fit start and end times, nr of points, WSS,",
  "and parameters (p1, p2, ...) of the fitted function to the fit file.",
  " ",
  "See also: fit2dat, fit_exp, fit_sigm, fit_ppf",
  " ",
  "Keywords: curve fitting, TAC, simulation, rational function",
  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, pi, ri, type=-1, ret;
  char        *cptr, datfile[FILENAME_MAX], fitfile[FILENAME_MAX],
               svgfile[FILENAME_MAX];
  int          weights=0; // 0=default, 1=no weighting, 2=frequency
  DFT          dft;
  FIT          fit;
  double       a, b, c, tstart, tstop, miny, maxy;
  enum         {MODEL_UNKNOWN, MODEL_P3P3, MODEL_P3P2, MODEL_P2P2, 
                MODEL_P2P1, MODEL_P1P1};
  int          model=MODEL_UNKNOWN;
  static char *model_name[] = {"Unknown", "pol3/pol3", "pol3/pol2", "pol2/pol2",
                               "pol2/pol1", "pol1/pol1",  0};
  int func_par_nr[6]={0,6,5,4,3,2};
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=fitfile[0]=svgfile[0]=(char)0;
  dftInit(&dft); fitInit(&fit);
  /* 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, "P3P3")==0) {
      model=MODEL_P3P3; type=233; continue;
    } else if(strcasecmp(cptr, "P3P2")==0) {
      model=MODEL_P3P2; type=232; continue;
    } else if(strcasecmp(cptr, "P2P2")==0) {
      model=MODEL_P2P2; type=222; continue;
    } else if(strcasecmp(cptr, "P2P1")==0) {
      model=MODEL_P2P1; type=221; continue;
    } else if(strcasecmp(cptr, "P1P1")==0) {
      model=MODEL_P1P1; type=211; 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;
    }
    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(!datfile[0]) {strcpy(datfile, argv[ai]); continue;}
    else if(!fitfile[0]) {strcpy(fitfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%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 function was not selected by options, use default */
  if(type<=0) type=233; 
  if(model==MODEL_UNKNOWN) model=MODEL_P3P3;
  parNr=func_par_nr[model];
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("fitfile := %s\n", fitfile);
    printf("svgfile := %s\n", svgfile);
    printf("type := %d\n", type);
    printf("parNr := %d\n", parNr);
    printf("weights := %d\n", weights);
  }
 
  /*
   *  Set parameter constraints
   *  note that these apply to the [0,1] scaled data!
   */
  pmin[0]=-1.0e3;  pmax[0]=5.0e3;
  pmin[1]=-5.0e2;  pmax[1]=1.0e3;
  pmin[2]=-5.0e2;  pmax[2]=1.0e4;
  pmin[3]=-0.10;   pmax[3]=1.0e4;
  pmin[4]=-5.00;   pmax[4]=5.0e3;
  pmin[5]=-1.0e2;  pmax[5]=1.0e4;
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", datfile);
  if(dftRead(datfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", datfile, dfterrmsg);
    return(2);
  }
  if(dft.frameNr<3) {
    fprintf(stderr, "Error: not enough samples for decent fitting.\n");
    dftEmpty(&dft); return(3);
  }
  dataNr=dft.frameNr;
 
  /* Check for NA's */
  if(dft_nr_of_NA(&dft) > 0) {
    fprintf(stderr, "Error: missing sample(s) in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
 
  /* Sort the samples by time in case data is catenated from several curves */
  (void)dftSortByFrame(&dft);
  if(verbose>30) dftPrint(&dft);
 
  /* 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", datfile);
    dftEmpty(&dft); return(2);
  }
  if(dataNr<3 || tstop<=0.0 || maxy<=0.0) {
    fprintf(stderr, "Error: invalid contents in %s\n", datfile);
    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); return(4);
  }
  /* Set necessary information in fit data structure */
  fit.voiNr=dft.voiNr;
  strncpy(fit.datafile, datfile, FILENAME_MAX);
  tpcProgramName(argv[0], 1, 1, fit.program, 1024);
  strcpy(fit.unit, dft.unit);
  fit.time=time(NULL);
  for(ri=0; ri<dft.voiNr; ri++) {
    fit.voi[ri].type=type; fit.voi[ri].parNr=2+parNr;
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=dataNr;
  }
 
  /* 
   *  Scale data between [0,1] during the fitting
   */
  /* Set the scale factors */
  double xscale, yscale[dft.voiNr];
  xscale=1.0/tstop;
  for(ri=0; ri<dft.voiNr; ri++) {
    a=dft.voi[ri].y[0];
    for(fi=1; fi<dft.frameNr; fi++)
      if(dft.voi[ri].y[fi]>a) a=dft.voi[ri].y[fi];
    yscale[ri]=1.0/a;
  }
  if(verbose>1) {
    printf("xscale=%g yscale=", xscale);
    for(ri=0; ri<dft.voiNr; ri++) printf("%g ", yscale[ri]);
    printf("\n");
  }
  /* Allocate memory for scaled data for fitting */
  double *scaled_data;
  scaled_data=malloc(sizeof(double)*3*dft.frameNr);
  if(scaled_data==NULL) {
    fprintf(stderr, "Error: out of memory.\n");
    dftEmpty(&dft); fitEmpty(&fit);
    return(5);
  }
  /* Set pointers */
  x=scaled_data; ymeas=scaled_data+dft.frameNr; yfit=ymeas+dft.frameNr;
  w=dft.w;
  /* Calculate scaled x values */
  for(fi=0; fi<dft.frameNr; fi++) x[fi]=xscale*dft.x[fi];
 
 
  /*
   *  Fit each TAC
   */
  if(verbose>1) {printf("fitting\n"); fflush(stdout);}
  for(ri=0; ri<dft.voiNr; ri++) {
    if(verbose>0 && dft.voiNr>1) printf("%s\n", dft.voi[ri].name);
 
    /* Scale y values for this TAC */
    for(fi=0; fi<dft.frameNr; fi++) 
      ymeas[fi]=yscale[ri]*dft.voi[ri].y[fi];
 
    /* Set the initial values for parameters */
    fit.voi[ri].p[0]=0.0; fit.voi[ri].p[1]=1.0;
    fit.voi[ri].wss=3.402823e+38;
 
    /* fit */
    TGO_LOCAL_INSIDE=0;
    TGO_SQUARED_TRANSF=1;
    ret=tgo(pmin, pmax, func_ratf, NULL, parNr, 15, &fit.voi[ri].wss, 
            fit.voi[ri].p+2, 1000, 3, verbose-8);
    if(ret) {
      fprintf(stderr, "Error %d in TGO.\n", ret);
      dftEmpty(&dft); fitEmpty(&fit); free(scaled_data);
      return(6);
    }
 
    /* Print measured and fitted curve (scaled down) */
    if(verbose>6) {
      printf("     Measured  Fitted:\n");
      for(fi=0; fi<dft.frameNr; fi++)
        printf("  %2d: %8.2e   %8.2e %8.2e\n", fi+1, x[fi], ymeas[fi], yfit[fi]);
      printf("\n");
    }
    if(verbose>2) {
      printf(" fitted parameters (scaled), with limits:\n");
      for(pi=0; pi<parNr; pi++)
        printf(" %g [%g, %g]\n", fit.voi[ri].p[pi+2], pmin[pi], pmax[pi]);
      printf("\n");
    }
 
    /* Force parameter values inside limits, as they are in the simulation */ 
    modelCheckParameters(parNr, pmin, pmax,
                         fit.voi[ri].p+2, fit.voi[ri].p+2, &a);
    fit.voi[ri].wss/=a; // remove penalty from reported WSS
    /* Calculate the function values again in case parameters were changed */
    func_ratf(parNr, fit.voi[ri].p+2, NULL);
 
    /* Scale the fitted y values back */
    for(fi=0; fi<dft.frameNr; fi++) dft.voi[ri].y2[fi]=yfit[fi]/yscale[ri];
 
    /* Calculate the WSS from data in original scale */
    a=c=0.0;
    for(fi=0; fi<dft.frameNr; fi++) {
      b=dft.voi[ri].y2[fi]-dft.voi[ri].y[fi];
      a+=w[fi]*b*b;
      if(fi==0 || fabs(b)>c) c=fabs(b);
    }
    fit.voi[ri].wss=a;
    if(verbose>1) printf("  max_dev := %g\n", c);
    if(verbose>4) printf("\n  last_y := %g\n\n", dft.voi[ri].y2[dft.frameNr-1]);
 
    /* Scale back the parameters */
    for(pi=0; pi<parNr; pi+=2) 
      fit.voi[ri].p[pi+2]/=yscale[ri];
    for(pi=0; pi<parNr; pi++) 
      fit.voi[ri].p[pi+2]*=(double)pow(xscale, floor(1+pi/2));
 
  } /* next TAC */
  free(scaled_data);
 
 
  /*
   *  Print fit results on screen
   */
  if(verbose>1) {
    printf("\n"); 
    fitWrite(&fit, "stdout");
    printf("\n"); 
  }
 
  /*
   *  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);
  }
 
  /*
   *  Plotting of fitted TACs, if requested
   */
  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 plot 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_ratf(int parNr, double *p, void *fdata)
{
  int i;
  double v, dd, dr, x2, x3, 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 */
  x2=x3=1.0;
  for(i=0, wss=0.0; i<dataNr; i++) {
    /* Calculate x^2 and x^3 */
    if(parNr>=3) {x2=x[i]*x[i]; if(parNr>=5) x3=x2*x[i];}
    /* Calculate the divider */
    dr=1.0;
    if(parNr>=2) dr+=pa[1]*x[i];
    if(parNr>=4) dr+=pa[3]*x2;
    if(parNr>=6) dr+=pa[5]*x3;
    /* Calculate the dividend */
    dd=0.0;
    if(parNr>=1) dd+=pa[0]*x[i];
    if(parNr>=3) dd+=pa[2]*x2;
    if(parNr>=5) dd+=pa[4]*x3;
    /* Calculate the curve value */
    yfit[i]=dd/dr; 
    /* Calculate weighted SS */
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  if(!isfinite(wss)) wss=nan(""); 
  wss*=penalty;
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/