b2ptrap.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 "tpcfileutil.h"
#include "tpcift.h"
#include "tpctac.h"
#include "tpcpar.h"
#include "tpcli.h"
#include "tpccm.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
#include "tpcrand.h"
#include "tpcmodels.h"
#include "tpcnlopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
const int maxParNr=1;
/* Local functions */
double func_b2p(int parNr, double *p, void*);
/*****************************************************************************/
typedef struct FITDATA {
  unsigned int n;
  double       *x;
  double       *ymeas;
  double       *ysim;
  double       *w;
  double       start;
  optimality_criterion optcrit;
  int           verbose;
} FITDATA;
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Calculate or fit PTAC based on BTAC in cell trapping model.",
  "Give parameter k in units 1/min, or 'fit' in order to try to fit it.",
  " ",
  "Usage: @P [Options] btacfile k ptacfile",
  " ",
  "Options:",
  " -ztime=<Time>",
  "     Time (in minutes) after the PTAC is assumed to be 0 and BTAC steady;",
  "     by default 3 min; only used when k is fitted.",
  " -svg=<Filename>",
  "     BTAC and (1-HCT*)PTAC are plotted in specified SVG file.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: fitmtrap, imgmtrap, b2rbc",
  " ",
  "Keywords: input, blood, plasma",
  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    btacfile[FILENAME_MAX], ptacfile[FILENAME_MAX];
  char    svgfile[FILENAME_MAX];
  double  fixed_k=nan("");
  double  zt=3.0;
  int     weights=1; // 0=file, 1=no weighting, 2=frequency
  optimality_criterion optcrit=OPTCRIT_OLS;
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  btacfile[0]=ptacfile[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) continue;
    if(strncasecmp(cptr, "SVG=", 4)==0 && strlen(cptr)>4) {
      strlcpy(svgfile, cptr+4, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "ZTIME=", 6)==0 && strlen(cptr)>6) {
      if(!atofCheck(cptr+6, &zt) && zt>0.0) continue;
    } else if(strcasecmp(cptr, "W1")==0) {
      weights=1; continue;
    } else if(strcasecmp(cptr, "WF")==0) {
      weights=2; 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(btacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    if(strcasecmp(argv[ai], "FIT")!=0) {
      if(atofCheck(argv[ai], &fixed_k) || !(fixed_k>0.0)) {
        fprintf(stderr, "Error: invalid k value '%s'.\n", argv[ai]);
        return(1);
      }
    }
    ai++;
  }
  if(ai<argc) strlcpy(ptacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {
    fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]);
    return(1);
  }
 
  /* Did we get all the information that we need? */
  if(!ptacfile[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("btacfile := %s\n", btacfile);
    printf("ptacfile := %s\n", ptacfile);
    if(!isnan(fixed_k)) printf("fixed_k := %g\n", fixed_k);
    else printf("ztime := %g min\n", zt);
    printf("optcrit := %s\n", optcritCode(optcrit));
    fflush(stdout);
  }
 
 
  /*
   *  Read TAC data
   */
  if(verbose>1) printf("reading %s\n", btacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, btacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); 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); return(2);
  } else if(tac.tacNr>1) {
    if(verbose>0) fprintf(stderr, "Warning: file contains more than one curve.\n");
    tac.tacNr=1;
  }
  /* Check NaNs */
  if(tacNaNs(&tac)>0) {
    fprintf(stderr, "Error: data contains missing values.\n");
    tacFree(&tac); return(2);
  }
  /* Sort data by sample time */
  tacSortByTime(&tac, &status);
  /* Take average of any duplicate samples */
  if(tacMultipleSamples(&tac, 1, &tac, verbose-2)!=0) {
    fprintf(stderr, "Error: cannot process duplicate samples.\n");
    tacFree(&tac); return(2);
  }
  /* Get initial x range */
  double xmin, xmax;
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  /* Convert sample times into minutes */
  if(tac.tunit==UNIT_UNKNOWN) {
    if(xmax<10.0) tac.tunit=UNIT_MIN; else tac.tunit=UNIT_SEC;
  }
  if(tacXUnitConvert(&tac, UNIT_MIN, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  /* Get final x range */
  if(tacXRange(&tac, &xmin, &xmax)!=0) {
    fprintf(stderr, "Error: invalid data sample times.\n");
    tacFree(&tac); return(2);
  }
  if(verbose>1) {
    printf("xmin := %g\n", xmin);
    printf("xmax := %g\n", xmax);
  }
  if(tac.sampleNr<3) {
    fprintf(stderr, "Error: too few samples.\n");
    tacFree(&tac); return(2);
  }
  /* Check that ztime setting is feasible with the data */
  if(isnan(fixed_k) && (zt>=xmax || zt<=xmin)) {
    fprintf(stderr, "Error: ztime outside of data range.\n");
    tacFree(&tac); return(2);
  }
 
 
  /* Add place for PTAC */
  if(verbose>1) fprintf(stdout, "allocating memory for PTAC\n");
  if(tacAllocateMore(&tac, 1)!=TPCERROR_OK) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    tacFree(&tac); return(3);
  }
  strcpy(tac.c[1].name, "PTAC");
  strcpy(tac.c[0].name, "BTAC");
  tac.tacNr=2;
 
 
  /* Add data weights, if requested */
  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); return(3);
    }
  }
 
 
 
  /* Set data pointers for the fit */
  FITDATA fitdata;
  fitdata.n=tac.sampleNr;
  fitdata.x=tac.x;
  fitdata.ymeas=tac.c[0].y;
  fitdata.ysim=tac.c[1].y;
  fitdata.w=tac.w;
  fitdata.start=zt;
  fitdata.optcrit=optcrit;
  if(verbose>10) fitdata.verbose=verbose-10; else fitdata.verbose=0;
 
 
  /*
   *  If k is given, then simulate PTAC
   */
  if(fixed_k>0.0) {
    func_b2p(1, &fixed_k, &fitdata);
  }
 
 
  /*
   *  Fitting, if k value was not given
   */
  if(isnan(fixed_k)) {
 
    if(verbose>1) printf("preparing for fitting\n");
    drandSeed(1);
 
    /* Set NLLS options */
    NLOPT nlo; nloptInit(&nlo);
    if(nloptAllocate(&nlo, 1)!=TPCERROR_OK) {
      fprintf(stderr, "Error: cannot initiate NLLS.\n");
      tacFree(&tac); return(5);
    }
    nlo._fun=func_b2p;
    nlo.fundata=&fitdata;
    nlo.totalNr=1;
    nlo.xlower[0]=0.0;
    nlo.xupper[0]=1.0;
    nlo.xtol[0]=0.00001;
    nlo.xfull[0]=0.25; // initial guess
    nlo.maxFunCalls=500000;
    /* Fit */
    if(verbose>4) nloptWrite(&nlo, stdout);
    if(nloptITGO1(&nlo, 1, 0, 10000, 20, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); nloptFree(&nlo); return(6);
    }
    nlo._fun(nlo.totalNr, nlo.xfull, nlo.fundata);
    if(verbose>2) nloptWrite(&nlo, stdout);
    if(verbose>6) {
      printf("measured and simulated TAC:\n");
      for(int i=0; i<tac.sampleNr; i++)
        printf("\t%g\t%g\t%g\n", tac.x[i], tac.c[0].y[i], tac.c[1].y[i]);
    }
    if(verbose>0) printf("WSS := %g\n", nlo.funval);
    if(verbose>0) printf("k := %g\n", nlo.xfull[0]);
 
    nloptFree(&nlo);
  }
 
  /*
   *  Write PTAC
   */
  {
    if(verbose>1) printf("saving (1-HCT)*Cp TAC\n");
    tacSwapTACCs(&tac, 0, 1);
    tac.tacNr=1;
    int ret=TPCERROR_OK;
    FILE *fp; fp=fopen(ptacfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing '%s'.\n", ptacfile);
      ret=TPCERROR_FAIL;
    } else {
      ret=tacWrite(&tac, fp, TAC_FORMAT_PMOD, 1, &status);
      fclose(fp);
      if(ret!=TPCERROR_OK) fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    }
    if(ret!=TPCERROR_OK) {
      tacFree(&tac);
      return(11);
    }
    if(verbose>0) printf("(1-HCT)*PTAC written in %s.\n", ptacfile);
  }
 
 
  tacFree(&tac);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************
 *
 *  Functions to be minimized
 *
 *****************************************************************************/
double func_b2p(int parNr, double *p, void *fdata)
{
  FITDATA *d=(FITDATA*)fdata;
 
  if(d->verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  if(parNr!=1) return(nan(""));
  if(d->verbose>20) printf("p[]: %g\n", p[0]);
 
  /* Calculate the curve values and weighted SS */
 
  double cpi=0.0;
  double t_last=0.0, cp_last=0.0, cpi_last=0.0;
 
  double wss=0.0;
  for(unsigned int i=0; i<d->n; i++) {
    /* delta time / 2 */
    double dt2=0.5*(d->x[i]-t_last);
    if(!isnan(d->ymeas[i]) && !isnan(d->x[i]) && dt2>0.0) {
      /* plasma */
      d->ysim[i] = (d->ymeas[i] - p[0]*(cpi_last + dt2*cp_last)) / (1.0 + dt2*p[0]);
      /* integral */
      cpi+=(d->ysim[i]+cp_last)*dt2;
    } else {
      d->ysim[i]=cp_last;
    }
    t_last=d->x[i];
    cp_last=d->ysim[i];
    cpi_last=cpi;
    double v=d->ysim[i]; // PTAC compared to zero
    if(d->x[i]>=d->start && !isnan(v)) {
      if(d->optcrit==OPTCRIT_LAD)
        wss+=d->w[i]*fabs(v);
      else
        wss+=d->w[i]*v*v;
    }
  }
 
  if(d->verbose>20) printf("p[]: %g   WSS=%g   start=%g\n", p[0], wss, d->start);
  return(wss);
}
/*****************************************************************************/
 
/*****************************************************************************/