tacslope.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 "tpcift.h"
#include "tpctac.h"
#include "tpcmodels.h"
#include "tpcpar.h"
#include "tpcfunc.h"
#include "tpctacmod.h"
#include "tpclinopt.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Find the highest slope in TAC file. The intercept with x axis",
  "can be used as an estimate of tracer appearance time.",
  " ",
  "Usage: @P [options] tacfile [parfile]",
  " ",
  "Options:",
  " -n=<fitted sample number>",
  "     Number of samples used to determine the slopes; by default 3.",
  " -positive",
  "     Slope is verified to be positive, and if not, then error is returned.",
  " -rmbkg",
  "     Sample values before estimated appearance time are set to zero.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: tactime, fitdelay, taccut, tacpeak, inpstart, tacrange",
  " ",
  "Keywords: TAC, input, time delay, time",
  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   tacfile[FILENAME_MAX], parfile[FILENAME_MAX];
  int    rmBkg=0; // 1=set values before appearance time to zero
  int    reqPos=0; // 1=error, if slope is not positive
  int    nFit=3; // Nr of samples used in line fits
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile[0]=parfile[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, "RMBKG", 2)==0) {
      rmBkg=1; reqPos=1; continue;
    } else if(strncasecmp(cptr, "POSITIVE", 3)==0) {
      reqPos=1; continue;
    } else if(strncasecmp(cptr, "N=", 2)==0) {
      if(atoiCheck(cptr+2, &nFit)==0 && nFit>1) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'\n", argv[ai]);
    return(1);
  } else break; // tac name argument may start with '-'
 
  TPCSTATUS status; statusInit(&status);
  statusSet(&status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  status.verbose=verbose-1;
  
  /* 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);}
 
 
  /* Arguments */
  if(ai<argc) strlcpy(tacfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(parfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: too many arguments.\n"); return(1);}
 
  /* Is something missing? */
  if(!tacfile[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
 
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    for(ai=0; ai<argc; ai++) printf("%s ", argv[ai]); 
    printf("\n");
    printf("tacfile := %s\n", tacfile);
    printf("nFit := %d\n", nFit);
    printf("reqPos := %d\n", reqPos);
    printf("rmBkg := %d\n", rmBkg);
  }
 
 
  /*
   *  Read the file
   */
  if(verbose>1) printf("reading %s\n", tacfile);
  TAC tac; tacInit(&tac);
  if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  if(verbose>2) {
    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));
  }
  /* Sort the data by sample times (x) */
  if(tacSortByTime(&tac, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(2);
  }
  /* Get x range */
  double xmin, xmax;
  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);
  }
  /* 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);
  }
  /* Check that fitting required nr of samples seems possible */
  if(nFit>=tac.sampleNr) {
    fprintf(stderr, "Error: too few samples for fitting.\n");
    tacFree(&tac); return(2);
  }
 
 
  /*
   *  Prepare space for results
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, 4, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  par.tacNr=tac.tacNr;
  par.parNr=4;
  par.format=PAR_FORMAT_TSV_UK;
  for(int i=0; i<tac.tacNr; i++) {
    par.r[i].model=0;
    par.r[i].dataNr=tacWSampleNr(&tac);
    par.r[i].start=xmin; 
    par.r[i].end=xmax;
  }
  /* Set parameter names */
  strcpy(par.n[0].name, "sMax");
  strcpy(par.n[1].name, "yIc"); par.n[1].unit=tac.cunit;
  strcpy(par.n[2].name, "Ta"); par.n[2].unit=tac.tunit;
  strcpy(par.n[3].name, "TsMax"); par.n[3].unit=tac.tunit;
 
 
  /*
   *  Calculate max slopes
   */
  if(verbose>1) {printf("searching for maximal slopes...\n"); fflush(stdout);}
  for(int i=0; i<tac.tacNr; i++) {
    double slope, yi, xi, xh;
    int ret=highestSlope(tac.x, tac.c[i].y, tac.sampleNr, nFit, nan(""), &slope, &yi, &xi, &xh);
    if(ret) {
      fprintf(stderr, "Error: cannot calculate max slope.\n");
      if(verbose>2) fprintf(stderr, "ret := %d\n", ret);
      tacFree(&tac); parFree(&par); return(4);
    }
    par.r[i].p[0]=slope;
    par.r[i].p[1]=yi;
    par.r[i].p[2]=xi;
    par.r[i].p[3]=xh;
    if(reqPos && !(slope>0.0)) {
      fprintf(stderr, "Error: slope is not positive.\n");
      tacFree(&tac); parFree(&par); return(5);
    }
  }
  if(verbose>1) {printf("... done.\n"); fflush(stdout);}
  /* Current, possibly edited, TAC data is not needed later */
  tacFree(&tac);
 
 
  /* 
   *  If requested, save parameters
   */
  if(parfile[0]) {
    /* Save file */
    if(verbose>1) printf("  saving %s\n", parfile);
    FILE *fp=fopen(parfile, "w");
    if(fp==NULL) {
      fprintf(stderr, "Error: cannot open file for writing.\n");
      parFree(&par); return(11);
    }
    int ret=parWrite(&par, fp, PAR_FORMAT_UNKNOWN, 1, &status);
    fclose(fp);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      parFree(&par); return(12);
    }
    if(verbose>0) printf("parameters saved in %s\n", parfile);
  }
 
  /*
   *  Print on screen the mean results, in a format identified by tactime.
   */
  for(int pi=0; pi<par.parNr; pi++) {
    double psum=0.0; int n=0;
    for(int i=0; i<par.tacNr; i++) if(isfinite(par.r[i].p[pi])) {psum+=par.r[i].p[pi]; n++;}
    if(n>0) printf("%s := %g\n", par.n[pi].name, psum/(double)n);
  }
 
 
  /*
   *  If requested, set to zero the values before Ta
   */
  if(rmBkg) {
    if(verbose>1) printf("setting values before Ta to zero...\n");
    /* Read TAC file again */
    if(tacRead(&tac, tacfile, &status)!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      tacFree(&tac); parFree(&par); return(6);
    }
    /* Set to zero the values before Ta */
    int nChanged=0;
    for(int ri=0; ri<tac.tacNr; ri++) {
      for(int fi=0; fi<tac.sampleNr; fi++)
        if(!(tac.x[fi]>=par.r[ri].p[2]) && !isnan(tac.c[ri].y[fi]) && tac.c[ri].y[fi]!=0.0) {
          tac.c[ri].y[fi]=0.0; nChanged++;
        }
    }
    /* If changes were made, then overwrite the original TAC file */
    if(nChanged==0) {
      if(verbose>0) printf("no changes made to %s\n", tacfile);
    } else {
      if(verbose>1) printf("writing %s\n", tacfile);
      FILE *fp; fp=fopen(tacfile, "w");
      if(fp==NULL) {
        fprintf(stderr, "Error: cannot open file for writing (%s)\n", tacfile);
        tacFree(&tac); parFree(&par); return(21);
      }
      if(tacWrite(&tac, fp, TAC_FORMAT_UNKNOWN, 1, &status)!=TPCERROR_OK) {
        fprintf(stderr, "Error: %s\n", errorMsg(status.error));
        fclose(fp); tacFree(&tac); parFree(&par); return(22);
      }
      fclose(fp);
    }
    tacFree(&tac);
    if(verbose>0) printf("edited TAC file saved.\n");
  }
  parFree(&par);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/