paucinf.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(int parNr, double *p, void*);
/*****************************************************************************/
 
/*****************************************************************************/
double *x, *ymeas, *yfit, *w;
int dataNr=0, parNr=0;
double pmin[MAX_PARAMETERS], pmax[MAX_PARAMETERS];
double start_time, end_time;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Estimates the half-life (t1/2) and elimination rate constant (kEL) of",
  "a PET tracer in plasma, and the AUC of plasma time-activity curve (PTAC)",
  "from 0 to infinite time. AUC can be used e.g. to estimate the total",
  "clearance (ClT) of PET tracer after a single intravenous dose:",
  " ",
  "          Dose     ",
  "   CL = --------   ",
  "     T  AUC        ",
  "           0-Inf   ",
  " ",
  "The natural logarithm of tracer concentration is plotted against time",
  "from bolus infusion. The plot becomes linear in the end phase, as the",
  "tracer is eliminated according to the laws of first-order reaction",
  "kinetics. The slope of the linear part of the plot equals -kEL.",
  " ",
  "Usage: @P [Options] ptacfile [resultfile]",
  " ",
  "Options:",
  " -end=<sample time>",
  "     Use data from 0 to given time in the line fit; by default",
  "     search for the best line fit extends to the last PTAC sample.",
  " -start=<sample time>",
  "     Use data from given time to end in the line fit; by default,",
  "     the search for the best line fit extends to the first sample.",
  " -winnonlin",
  "     Best line fit is searched using algorithm that gives similar results",
  "     to WinNonlin with the same simple model.",
//" -nlfit",
//"     Result based on line fitting is refined with exp function fitting.",
  " -minnr=<sample number>",
  "     Enter the minimum number of samples to use in the fit; by default 3.",
  " -svg=<Filename>",
  "     Plots of log transformed TACs and fitted lines are written in specified",
  "     SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Plasmafile must contain a time column, and one or more concentration columns",
  "separated by space(s) or tabulator(s). The result is printed on screen.",
  " ",
  "See also: extrapol, interpol, fit_exp, fit_dexp, fit2dat, fit2auc",
  " ",
  "Keywords: input, modelling, pharmacokinetics, clearance, elimination rate",
  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     datfile[FILENAME_MAX], outfile[FILENAME_MAX], svgfile[FILENAME_MAX], temp[512];
  double   fittime=-1.0, starttime=0.0;
  int      check_for_improvement=0; // 1=WinNonlin compatibility mode
  int      nonlinfit=0; // 1=best line fit range is re-fitted (non-linear)
  int      forced_minNr=0;
  FILE    *logfp; // file pointer for extra output
  int      ret;
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  datfile[0]=outfile[0]=svgfile[0]=(char)0;
  /* 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==NULL) continue;
    if(strncasecmp(cptr, "E=", 2)==0 && strlen(cptr)>2) {
      fittime=atof_dpi(cptr+2); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "END=", 4)==0 && strlen(cptr)>4) {
      fittime=atof_dpi(cptr+4); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "ENDTIME=", 8)==0) {
      fittime=atof_dpi(cptr+8); if(fittime>0.0) continue;
    } else if(strncasecmp(cptr, "START=", 6)==0 && strlen(cptr)>6) {
      starttime=atof_dpi(cptr+6); if(starttime>=0.0) continue;
    } else if(strncasecmp(cptr, "STARTTIME=", 10)==0) {
      starttime=atof_dpi(cptr+10); if(starttime>=0.0) continue;
    } else if(strncasecmp(cptr, "WINNONLIN", 3)==0) {
      check_for_improvement=1; continue;
    } else if(strcasecmp(cptr, "NLFIT")==0) {
      nonlinfit=1; continue;
    } else if(strncasecmp(cptr, "MINNR=", 6)==0) {
      if(atoi_with_check(cptr+6, &forced_minNr)==0 && forced_minNr>=3) continue;
    } else if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); 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 */
  if(ai<argc) strlcpy(datfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(outfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
 
  /* Is something missing? */
  if(!datfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  /* WinNonlin does not do nonlinear fitting */
  if(check_for_improvement==1) nonlinfit=0;
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("ptacfile := %s\n", datfile);
    printf("svgfile := %s\n", svgfile);
    if(outfile[0]) printf("resfile := %s\n", outfile);
    if(fittime>=0.0) printf("endtime := %g\n", fittime);
    if(starttime>=0.0) printf("starttime := %g\n", starttime);
    if(forced_minNr>0) printf("forced_minNr := %d\n", forced_minNr);
    printf("starttime := %g\n", starttime);
    printf("check_for_improvement := %d\n", check_for_improvement);
    printf("nonlinfit := %d\n", nonlinfit);
  }
  if(verbose>1) logfp=stdout; else logfp=NULL;
 
 
  /*
   *  Read data
   */
  if(verbose>1) printf("reading %s\n", datfile);
  if(verbose>10) CSV_TEST=verbose-10; else CSV_TEST=0;
  DFT dft; dftInit(&dft);
  if(dftRead(datfile, &dft)) {
    fprintf(stderr, "Error in reading '%s': %s\n", datfile, dfterrmsg);
    return(2);
  }
  /* 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);
 
  /* Set the last sample to be fitted */
  int fitNr=dft.frameNr;
  if(fittime>0.0) {
    while(dft.x[fitNr-1]>fittime && fitNr>0) fitNr--;
  }
  fittime=dft.x[fitNr-1];
  if(verbose>1) printf("fittime := %g\n", dft.x[fitNr-1]);
 
  /* Set Max nr of fitted points */
  int max_nr=fitNr;
  int min_nr=3; if(forced_minNr>min_nr) min_nr=forced_minNr;
  if(starttime>0.0) {
    int i;
    for(i=fitNr-1; i>=0; i--) if(dft.x[i]<starttime) break;
    max_nr=fitNr-(i+1);
  }
  if(verbose>1) {
    printf("min_nr := %d\n", min_nr);
    printf("max_nr := %d\n", max_nr);
  }
  /* Check that there are at least 3 samples */
  if(max_nr<3 || max_nr<min_nr) {
    fflush(stdout);
    fprintf(stderr, "Error: too few samples for exponential fitting\n");
    fflush(stderr); dftEmpty(&dft); return(2);
  }
  /* or give a warning if less than 5 samples */
  if(max_nr<5) {
    fprintf(stderr, "Warning: few samples for exponential fitting\n");
  }
 
 
  /*
   *  Allocate memory for results
   */
  if(verbose>1) printf("allocating memory for results.\n");
  RES res; resInit(&res);
  if(res_allocate_with_dft(&res, &dft)!=0) {
    fprintf(stderr, "Error: cannot setup memory for results.\n");
    dftEmpty(&dft); resEmpty(&res); return(3);
  }
  /* Copy titles & filenames */
  tpcProgramName(argv[0], 1, 1, res.program, 128);
  strncpy(res.plasmafile, datfile, FILENAME_MAX);
  sprintf(res.datarange, "%g %s", fittime, petTunit(dft.timeunit));
  res.datanr=max_nr;
  res.time=time(NULL);
  res.isweight=0;
  /* Set the parameter names */
  res.parNr=5;
  {
  int pi;
  pi=0; strcpy(res.parname[pi], "kEL"); strcpy(res.parunit[pi], "");
  pi++; strcpy(res.parname[pi], "t1/2"); strcpy(res.parunit[pi], "");
  pi++; strcpy(res.parname[pi], "AUC"); strcpy(res.parunit[pi], "");
  pi++; strcpy(res.parname[pi], "Start"); strcpy(res.parunit[pi], "");
  pi++; strcpy(res.parname[pi], "C0"); strcpy(res.parunit[pi], "");
  }
 
  /*
   *  Make ln transformation for TACs
   */
  if(verbose>1) printf("allocating memory for ln transformed data\n");
  DFT dftln; dftInit(&dftln); 
  if(dftdup(&dft, &dftln)!=0) {
    fprintf(stderr, "Error: out of memory.\n");
    dftEmpty(&dft); resEmpty(&res); return(4);
  }
  if(verbose>1) printf("ln transformation\n");
  if(dft_ln(&dftln, NULL)!=0) {
    fprintf(stderr, "Error: cannot make ln transformation.\n");
    dftEmpty(&dft); resEmpty(&res); return(4);
  }
  if(verbose>10) {dftPrint(&dftln); fflush(stdout);}
 
 
  /*
   *  Compute best linear fit to the end of ln transformed TACs.
   */
  if(verbose>1) printf("linear fitting\n");
  FIT fit; fitInit(&fit); 
  {
  double f;
  int n;
  f=fittime; n=min_nr;
  ret=dft_end_line(&dftln, &f, &n, max_nr, -1.0, check_for_improvement,
                   &fit, logfp, temp);
  }
  if(ret!=0) {
    fprintf(stderr, "Error in linear fit: %s.\n", temp);
    dftEmpty(&dft); dftEmpty(&dftln); resEmpty(&res); fitEmpty(&fit);
    return(6);
  }
  if(verbose>1) {
    printf("Results of linear fitting:\n");
    fitPrint(&fit);
  }
  if(verbose>1) {
    if(fit.voiNr==1) printf("Adj_R^2 := %g (n=%d)\n", fit.voi[0].p[2], fit.voi[0].dataNr);
    else for(int ri=0; ri<fit.voiNr; ri++)
      printf("Adj_R^2(%s) := %g (n=%d)\n", fit.voi[ri].name, fit.voi[ri].p[2], fit.voi[ri].dataNr);
  }
 
  /* Non-linear fit to the previously found fit time range, if required.
   * Note that this is NOT done in WinNonlin
   */
  if(nonlinfit==1) {
    if(verbose>0) printf("non-linear fitting\n");
    for(int ri=0; ri<dft.voiNr; ri++) {
      int pi;
      double pfit[MAX_PARAMETERS];
      if(verbose>1) {fprintf(stdout, "%s: \n", dft.voi[ri].name); fflush(stdout);}
      /* Set data pointers */
      dataNr=dft.frameNr;
      x=dft.x; ymeas=dft.voi[ri].y; yfit=dft.voi[ri].y2; w=dft.w;
      start_time=fit.voi[ri].start;
      end_time=fit.voi[ri].end;
      /* Set the initial values for parameters */
      fit.voi[ri].wss=3.402823e+38;
      for(pi=0; pi<parNr; pi++) pfit[pi]=0.0;
      /* Set parameter constraints for 1-exp fit */
      pmin[0]=0.0; pmax[0]=2.0*exp(fit.voi[ri].p[0]);
      pmin[1]=2.0*fit.voi[ri].p[1]; pmax[1]=0.0;
      /* fit */
      ret=tgo(pmin, pmax, func, NULL, parNr, 4, &fit.voi[ri].wss, pfit, 0, 0, verbose-19);
      if(ret) {
        fprintf(stderr, "Error %d in TGO.\n", ret);
        dftEmpty(&dft); dftEmpty(&dftln); resEmpty(&res); fitEmpty(&fit);
        return(7);
      }
      /* Print measured and fitted curve */
      if(verbose>9) {
        printf("     Measured  Fitted:\n");
        for(int i=0; i<dft.frameNr; i++)
          printf("  %2d: %8.2e   %8.2e %8.2e\n", i+1, x[i], ymeas[i], yfit[i]);
      }
      if(verbose>2) {
        printf(" fitted parameters:");
        for(int pi=0; pi<parNr; pi++) printf(" %g", pfit[pi]); 
        printf("\n");
      }
      if(verbose>3) {
        printf(" parameter constraints:");
        for(int pi=0; pi<parNr; pi++) printf(" (%g - %g)", pmin[pi], pmax[pi]); 
        printf("\n");
      }
      /* Write fit results into FIT struct */
      fit.voi[ri].p[0]=log(pfit[0]);
      fit.voi[ri].p[1]=pfit[1];
    } /* next TAC */
  }
 
  /*  Calculate AUC using trapezoidal rule from time 0 to the time of last
   *  measured sample, where the extrapolation starts.
   *  Write fit parameters into results
   */
  for(int ri=0; ri<dftln.voiNr; ri++) {
    double auc_0_t, kel, c0, t12, auc_t_inf, auc_0_inf;
    auc_0_t=kel=c0=t12=auc_t_inf=auc_0_inf=0.0;
    if(verbose>0 && dft.voiNr>1) printf("%s :\n", dftln.voi[ri].name);
    /* Trapezoidal AUC calculation */
    integrate(dft.x, dft.voi[ri].y, dft.frameNr, dft.voi[ri].y3);
    auc_0_t=dft.voi[ri].y3[fitNr-1];
    if(verbose>2) printf("auc_0_t := %g\n", auc_0_t);
    /* Get parameters from the fit */
    kel=-fit.voi[ri].p[1];
    c0=exp(fit.voi[ri].p[0]);
    /* Calculate t(1/2) */
    t12=M_LN2/kel; if(verbose>1) printf("t(1/2) := %g\n", t12);
    /* Exp function integral from last (fitted) time point to infinity */
    auc_t_inf=c0*exp(-kel*fit.voi[ri].end)/kel;
    if(verbose>2) printf("auc_t_inf := %g\n", auc_t_inf);
    /* Fill the results */
    res.voi[ri].parameter[0]=kel;
    res.voi[ri].parameter[1]=t12;
    res.voi[ri].parameter[2]=auc_0_inf=auc_0_t+auc_t_inf;
    res.voi[ri].parameter[3]=fit.voi[ri].start; // selected fit start time
    res.voi[ri].parameter[4]=c0;
    /* Print results */
    if(verbose>0 || !outfile[0]) {
      fprintf(stdout, "  AUC(0-Inf)=%g [AUC(0-%g)=%g AUC(%g-Inf)=%g]\n",
        auc_0_inf, fit.voi[ri].end, auc_0_t, fit.voi[ri].end, auc_t_inf);
      fprintf(stdout, "  k(el)=%g t(1/2)=%g estimated between %g-%g\n",
        kel, t12, fit.voi[ri].start, fit.voi[ri].end);
    }
  }
  fitEmpty(&fit);
 
 
  /*
   *  Plotting of fitted TACs
   */
  if(svgfile[0]) {
 
    /* Create a DFT containing fitted TACs */
    double k, b;
    char tmp[64];
    DFT dft3;
    dftInit(&dft3); ret=dftdup(&dftln, &dft3);
    if(ret) {
      fprintf(stderr, "Error: cannot save fitted curves.\n");
      dftEmpty(&dft); dftEmpty(&dftln); resEmpty(&res);
      return(21);
    }
    for(int ri=0; ri<dft3.voiNr; ri++) {
      k=-res.voi[ri].parameter[0];
      b=log(res.voi[ri].parameter[4]);
      for(int fi=0; fi<dft3.frameNr; fi++) {
        dft3.voi[ri].y[fi]=nan("");
        if(dft3.x[fi]<res.voi[ri].parameter[3]) continue;
        if(dft3.x[fi]>fittime) continue;
        dft3.voi[ri].y[fi]=b+k*dft3.x[fi];
      }
    }
 
    /* Save SVG plot of fitted and original data */
    if(verbose>1) printf("saving SVG plot\n");
    sprintf(tmp, "Line fit to log-transformed data ");
    if(strlen(res.studynr)>0 && strcmp(res.studynr, ".")!=0) strcat(tmp, res.studynr);
    ret=plot_fit_svg(&dftln, &dft3, tmp, svgfile, verbose-30);
    if(ret) {
      fprintf(stderr, "Error (%d) in writing '%s'.\n", ret, svgfile);
      dftEmpty(&dftln); dftEmpty(&dft3); dftEmpty(&dft); resEmpty(&res);
      return(30+ret);
    }
    if(verbose>0) printf("Plots written in %s\n", svgfile);
 
    dftEmpty(&dft3);
  }
  dftEmpty(&dftln);
 
 
  /*
   *  Save the results, if required
   */
  if(outfile[0]) {
    if(verbose>1) printf("saving results in %s\n", outfile);
    ret=resWrite(&res, outfile, verbose-10);
    if(ret) {
      fprintf(stderr, "Error in writing '%s': %s\n", outfile, reserrmsg);
      dftEmpty(&dft); resEmpty(&res); return(11);
    }
  }
 
  /* Free memory */
  dftEmpty(&dft); resEmpty(&res);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func(int parNr, double *p, void *fdata)
{
  double v, wss=0.0;
  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 */
  for(int i=0; i<dataNr; i++) {
    if(x[i]<start_time || x[i]>end_time || w[i]<=0.0) continue;
    yfit[i]=pa[0]*exp(pa[1]*x[i]);
    v=yfit[i]-ymeas[i]; wss+=w[i]*v*v;
  }
  wss*=penalty;
 
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/