fit_exp.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"
/*****************************************************************************/
 
/*****************************************************************************/
/* Local functions */
double func_exp(int parNr, double *p, void*);
/*****************************************************************************/
double *x, *ymeas, *yfit, *w;
int dataNr=0, parNr=0;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
int start_index, end_index;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Non-linear fitting of the sum of 1-3 exponential functions to PET plasma and",
  "tissue time-activity curves (TACs).",
  " ",
  "Function:",
  " ",
  "  f(x) = p1*exp(p2*x) + p3*exp(p4*x) + p5*exp(p6*x)",
  " ",
  "Usage: @P [Options] tacfile fitfile",
  " ",
  "Options:",
  " -3[e] | -2[e] | -1[e]",
  "     Fit to sum of three (default), two, or single exponential function.",
  " -nonneg",
  "     p1, p3 and p5 are constrained to values >= 0.",
  " -expneg",
  "     p2, p4 and p6 are constrained to values <= 0.",
  " -starttime=<<Time (min)>|Peak>",
  "     Set the sample time from where fitting is started (0 by default),",
  "     or 'peak' to start fitting from common peak time.",
  " -endtime=<Time (min)>",
  "     Set the sample time where fitting is stopped (last sample by default).",
  " -svg=<Filename>",
  "     Fitted and measured TACs are plotted in specified SVG file.",
  " -autoname",
  "     Names for fit file and fitted TAC file are set automatically",
  " -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.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Function parameters are written in fit (2) or result (3) format,",
  "depending on filename extension (.fit or .res).",
  " ",
  "References:",
  "1. https://www.turkupetcentre.net/petanalysis/doc/format_tpc_dft.html",
  "2. https://www.turkupetcentre.net/petanalysis/doc/format_tpc_fit.html",
  "3. https://www.turkupetcentre.net/petanalysis/doc/format_tpc_res.html",
  " ",
  "See also: fit_fexp, fit_dexp, llsqe3, fit_ratf, fit2dat, extrapol, tacln",
  " ",
  "Keywords: TAC, curve fitting, input, simulation, extrapolation",
  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      n, fi, pi, ri, type=303, ret;
  int      is_nonneg=0, is_expneg=0;
  // Weights: 0=left as it is in datafile, 1=set to 1, 2=sampling frequency
  int      weighting=0; 
  char     datfile[FILENAME_MAX], outfile[FILENAME_MAX], svgfile[FILENAME_MAX],
          *cptr, temp[512];
  DFT      dft;
  FIT      fit;
  double   a, b, c, tstart, tstop;
  double   starttime=-1.0, endtime=-1.0;
  int      peak_start=0; // 1=start fit from peak
  int      tgo_nr, neigh_nr, iter_nr;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=outfile[0]=svgfile[0]=(char)0;
  dftInit(&dft); fitInit(&fit);
  parNr=6;
  /* Options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; 
    if(strcasecmp(cptr, "3")==0 || strcasecmp(cptr, "3e")==0) {
      type=303; parNr=6; continue;
    } else if(strcasecmp(cptr, "2")==0 || strcasecmp(cptr, "2e")==0) {
      type=302; parNr=4; continue;
    } else if(strcasecmp(cptr, "1")==0 || strcasecmp(cptr, "1e")==0) {
      type=301; parNr=2; continue;
    } else if(strncasecmp(cptr, "NONNEG", 2)==0) {
      is_nonneg=1; continue;
    } else if(strncasecmp(cptr, "EXPNEG", 2)==0) {
      is_expneg=1; continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weighting=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weighting=2; continue;
    } else if(strncasecmp(cptr, "STARTTIME=", 10)==0) {
      if(strcasecmp(cptr+10, "PEAK")==0) {peak_start=1; continue;}
      starttime=atof_dpi(cptr+10); if(starttime>=0.0) continue;
    } else if(strncasecmp(cptr, "ENDTIME=", 8)==0) {
      endtime=atof_dpi(cptr+8); if(endtime>0.0) 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(!outfile[0]) {strcpy(outfile, argv[ai]); continue;}
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
  if(!outfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
 
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("datfile := %s\n", datfile);
    printf("outfile := %s\n", outfile);
    printf("svgfile := %s\n", svgfile);
    printf("type := %d\nparNr := %d\n", type, parNr);
    printf("weighting := %d\n", weighting);
    if(is_nonneg) printf("non-negativity constraint was set for p1, p3 and p5\n");
    if(is_expneg) printf("negativity constraint was set for p2, p4 and p6\n");
    if(peak_start==0)
      printf("starttime := %g min\n", starttime);
    else
      printf("starttime := peak\n");
    printf("endtime := %g min\n", endtime);
  }
  if(verbose>6) MATHFUNC_TEST=verbose-6; else MATHFUNC_TEST=0;
 
 
  /*
   *  Set parameter constraints;
   *  note that these apply to the [0,1] scaled data!
   */
  pmin[0]=-2.0e2;  pmax[0]=5.0e2;  if(is_nonneg) {pmin[0]=0.0; pmax[0]*=5.0;}
  pmin[1]=-5.0e2;  pmax[1]=5.0e2;  if(is_expneg) {pmax[1]=0.0; pmin[1]*=2.0;}
  pmin[2]=-1.0e2;  pmax[2]=5.0e1;  if(is_nonneg) {pmin[2]=0.0; pmax[2]*=5.0;}
  pmin[3]=-0.9;    pmax[3]=+0.9;   if(is_expneg) {pmin[3]=0.0;}
  pmin[4]=-2.0e2;  pmax[4]=5.0e1;  if(is_nonneg) {pmin[4]=0.0; pmax[4]*=5.0;}
  pmin[5]=-0.9;    pmax[5]=+0.9;   if(is_expneg) {pmin[5]=0.0;}
 
 
  /*
   *  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);
  }
 
  /*
   *  Delete the TACs that have name 'SD' (produced at least by avgbolus)
   *  unless all TACs would be deleted.
   */
  for(ri=0, n=0; ri<dft.voiNr; ri++) {
    dft.voi[ri].sw2=0;
    if(strcasecmp(dft.voi[ri].voiname, "SD")==0) {
      dft.voi[ri].sw2=1; n++; continue;}
    if(strcasecmp(dft.voi[ri].hemisphere, "SD")==0) {
      dft.voi[ri].sw2=1; n++; continue;}
  }
  if(n>0 && n<dft.voiNr) {
    if(verbose>0) printf("%d TAC(s) with SDs are not used.\n", n);
    ri=0;
    while(ri<dft.voiNr) {
      if(dft.voi[ri].sw2) dftDelete(&dft, ri); else ri++;
    }
  }
 
  /* 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 */
  dftSortByFrame(&dft);
 
  /* If fit is to be started at peak time, then find the common peak time */
  if(peak_start!=0) {
    ret=dftRobustMinMaxTAC(&dft, -1, NULL, &starttime, NULL, NULL, NULL, NULL,
                           NULL, NULL, verbose-10);
    if(ret!=0) {
      fprintf(stderr, "Error: can not determine peak time.\n");
      if(verbose>1) printf("ret := %d\n", ret);
      dftEmpty(&dft); return(3);
    }
    /* Convert it to minutes, if sample times are in seconds */
    if(dft.timeunit==TUNIT_SEC) starttime/=60.;
    if(verbose>1) printf("peak_time := %g min\n", starttime);
  }
 
  /* Set the fit time */
  if(endtime<=0.0) endtime=1.0E20;
  dataNr=fittime_from_dft(&dft, &starttime, &endtime, &start_index, &end_index, verbose-10);
  if(dataNr<0) {
    fprintf(stderr, "Error: check the fit time range.\n");
    if(verbose>1) printf("ret := %d\n", -dataNr);
    dftEmpty(&dft); return(3);
  }
  if(verbose>2) {
    printf("start_index := %d\n", start_index);
    printf("end_index := %d\n", end_index);
  }
  if(dataNr<3) {
    fprintf(stderr, "Error: too few samples in the fit range.\n");
    dftEmpty(&dft); return(3);
  }
  /* Get start and end times in the same units as in the data */
  if(dft.timetype==DFT_TIME_STARTEND) {
    tstart=dft.x1[start_index]; tstop=dft.x2[end_index];
  } else {
    tstart=dft.x[start_index]; tstop=dft.x[end_index];
  }
  /* Check times */
  if(dataNr<3 || tstop<=tstart) {
    fprintf(stderr, "Error: bad time scale in %s\n", datfile);
    dftEmpty(&dft); return(2);
  }
 
  /* Check and set weights */
  if(weighting==0) {
    if(dft.isweight==0) {
      dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;}
  } else if(weighting==1) {
    dft.isweight=1; for(fi=0; fi<dft.frameNr; fi++) dft.w[fi]=1.0;
  } else if(weighting==2) {
    dftWeightByFreq(&dft);
  }
  if(verbose>2) {
    printf("data_weights := %g", dft.w[0]);
    for(fi=1; fi<dft.frameNr; fi++) printf(", %g", dft.w[fi]);
    printf("\n");
  }
 
 
  /*
   *  Allocate memory for fits
   */
  if(verbose>1) printf("allocating memory for fits.\n");
  if((ret=fit_allocate_with_dft(&fit, &dft))!=0) {
    fprintf(stderr, "Error: cannot prepare memory for fit parameters (%d)", ret);
    dftEmpty(&dft); return(4);
  }
  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=parNr;
    fit.voi[ri].start=tstart; fit.voi[ri].end=tstop;
    fit.voi[ri].dataNr=dataNr;
  }
 
  /* 
   *  Use scaled [0,1] data 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 || (verbose>0 && dft.voiNr>1)) printf("fitting TACs\n");
  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];
 
    tgo_nr=100+200*parNr;
    neigh_nr=2*parNr;
    iter_nr=parNr;
    TGO_SQUARED_TRANSF=1;
    TGO_LOCAL_INSIDE=1;
 
    /* Set the initial values for parameters */
    fit.voi[ri].wss=3.402823e+38;
 
    /* fit */
    ret=tgo(pmin, pmax, func_exp, NULL, parNr, neigh_nr, &fit.voi[ri].wss,
            fit.voi[ri].p, tgo_nr, iter_nr, 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], pmin[pi], pmax[pi]);
      printf("\n"); fflush(stdout);
    }
 
    /* Correct fitted parameters to match constraints like inside the function */
    /* No need to remove penalty from WSS because WSS is calculated later again */
    (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) {fprintf(stdout, "ok\n");}
      else fprintf(stderr, "fit collided with the limits.\n");
    }
 
    /* Calculate the function values again in case parameters were changed */
    func_exp(parNr, fit.voi[ri].p, 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=1; pi<parNr; pi+=2) {
      fit.voi[ri].p[pi-1]/=yscale[ri];
    }
    fit.voi[ri].p[1]*=xscale;
    /* Convert exponentials to the original values from the relative values */
    if(parNr>3) fit.voi[ri].p[3]*=fit.voi[ri].p[1];
    if(parNr>5) fit.voi[ri].p[5]*=fit.voi[ri].p[3];
 
  } /* next TAC */
  free(scaled_data);
 
 
  /*
   *  Print fit results on screen
   */
  if(verbose>1 || strcasecmp(outfile, "stdout")==0) {
    printf("\n"); 
    fitWrite(&fit, "stdout");
    printf("\n"); 
  }
 
 
  /*
   *  Save fit results, if required
   */
  if(outfile[0] && strcasecmp(outfile, "stdout")!=0) {
    if(verbose>1) printf("writing fit results in %s\n", outfile);
    cptr=strrchr(outfile, '.'); if(cptr!=NULL) cptr++;
    if(strcasecmp(cptr, "FIT")==0) { // Save in FIT format
      if(fitWrite(&fit, outfile)) {
        fprintf(stderr, "Error in writing '%s': %s\n", outfile, fiterrmsg);
        dftEmpty(&dft); fitEmpty(&fit); return(11);
      }
    } else { // Save in RES format
      RES res; resInit(&res);
      ret=fitToResult(&fit, &res, temp);
      if(ret!=0) {
        fprintf(stderr, "Error in making results: %s\n", temp);
        dftEmpty(&dft); fitEmpty(&fit); resEmpty(&res); return(11);
      }
      if(parNr==6) strcpy(res.titleline, "Func p1 p2 p3 p4 p5 p6 WSS");
      else if(parNr==4) strcpy(res.titleline, "Func a1 c1 a2 c2 WSS");
      else strcpy(res.titleline, "Func a1 c1 WSS");
      strcpy(res.program, fit.program);
      sprintf(res.datarange, "%g - %g", tstart, tstop);
      res.isweight=dft.isweight;
      if(resWrite(&res, outfile, verbose-3)) {
        fprintf(stderr, "Error in writing '%s': %s\n", outfile, fiterrmsg);
        dftEmpty(&dft); fitEmpty(&fit); resEmpty(&res); return(11);
      }
      resEmpty(&res);
    }
  }
 
 
  /*
   *  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];
    int fj;
    DFT dft2; dftInit(&dft2);
    ret=dftAllocateWithHeader(&dft2, 1+end_index-start_index, dft.voiNr, &dft);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&dft); fitEmpty(&fit);
      return(21);
    }
    for(fi=start_index, fj=0; fi<=end_index; fi++) {
      dft2.x[fj]=dft.x[fi]; dft2.x1[fj]=dft.x1[fi]; dft2.x2[fj]=dft.x2[fi];
      for(ri=0; ri<dft.voiNr; ri++) dft2.voi[ri].y[fj]=dft.voi[ri].y2[fi];
      fj++;
    }
 
    /* Save SVG plot of fitted and original data */
    if(svgfile[0]) {
      tpcProgramName(argv[0], 0, 0, tmp, 64); strcat(tmp, " ");
      if(strlen(dft.studynr)>0) strcat(tmp, dft.studynr);
      ret=plot_fit_svg(&dft, &dft2, tmp, svgfile, verbose-10);
      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_exp(int parNr, double *p, void *fdata)
{
  int i;
  double v, wss;
  double pa[MAX_PARAMETERS], penalty=1.0;
 
  /* Check parameters against the constraints */
  (void)modelCheckParameters(parNr, pmin, pmax, p, pa, &penalty);
  if(fdata) {}
 
  /* Calculate the curve values and weighted SS */
  if(parNr==6) {
    for(i=start_index, wss=0.0; i<=end_index; i++) {
      if(isnan(ymeas[i]) || isnan(x[i])) continue;
      yfit[i]=pa[0]*exp(pa[1]*x[i]) + 
              pa[2]*exp(pa[1]*pa[3]*x[i]) + 
              pa[4]*exp(pa[1]*pa[3]*pa[5]*x[i]);
      v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
    }
  } else if(parNr==4) {
    for(i=start_index, wss=0.0; i<=end_index; i++) {
      if(isnan(ymeas[i]) || isnan(x[i])) continue;
      yfit[i]=pa[0]*exp(pa[1]*x[i]) + 
              pa[2]*exp(pa[1]*pa[3]*x[i]);
      v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
    }
  } else if(parNr==2) {
    for(i=start_index, wss=0.0; i<=end_index; i++) {
      if(isnan(ymeas[i]) || isnan(x[i])) continue;
      yfit[i]=pa[0]*exp(pa[1]*x[i]);
      v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
    }
  } else {printf("Program error: parNr=%d\n", parNr); exit(10);}
  wss*=penalty;
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/