tacpeak.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"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Find peak time and value in TAC file.",
  " ",
  "Usage: @P [options] tacfile [outputfile]",
  " ",
  "Options:",
  " -decr=<Percentage>",
  "     Estimate time when TAC decreases to specified percentage of peak value.",
  " -wt",
  "     Peak time is reported as weighted average of subsequent samples.",
  " -swa",
  "     TAC mean is subtracted from peak value.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "See also: imgpeak, dftmax, inpstart, tacrange, taccut, tacslope",
  " ",
  "Keywords: TAC, peak, max, 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];
  double decr=nan("");
  int    timeAsWAvg=0; // 1=report peak time as activity weighted mean.
  int    subtractWAvg=0; // 1=subtract TAC mean from the peak value.
 
 
  /*
   *  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, "DECR=", 5)==0) {
      if(atofCheck(cptr+5, &decr)==0) {
        if(decr>1.0 && decr<100.0) continue;
      }
    } else if(strcasecmp(cptr, "WT")==0) {
      timeAsWAvg=1; continue;
    } else if(strcasecmp(cptr, "SWA")==0) {
      subtractWAvg=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);
    if(parfile[0]) printf("parfile := %s\n", parfile);
    if(!isnan(decr)) printf("decr_percent := %g\n", decr);
    printf("timeAsWAvg := %d\n", timeAsWAvg);
    printf("subtractWAvg := %d\n", subtractWAvg);
  }
 
 
  /*
   *  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);
  }
 
  /*
   *  Prepare space for results
   */
  if(verbose>1) printf("preparing space for parameters\n");
  PAR par; parInit(&par);
  if(parAllocateWithTAC(&par, &tac, 3, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac); return(3);
  }
  par.tacNr=tac.tacNr;
  par.parNr=2; if(!isnan(decr)) par.parNr++;
  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, "Peak");
  strcpy(par.n[1].name, "PeakT");
  sprintf(par.n[2].name, "T%g%%", roundf(decr));
 
  /*
   *  Search peak values
   */
  if(verbose>1) printf("searching peaks\n");
  for(int ci=0; ci<tac.tacNr; ci++) {
 
    if(verbose>2) printf("%s\n", tac.c[ci].name);
    for(int i=0; i<par.parNr; i++) par.r[ci].p[i]=nan("");
 
    /* Find peak value and its sample index */
    double ymax; int imax;
    if(tacYRange(&tac, ci, NULL, &ymax, NULL, &imax, NULL, NULL)!=0) {
      fprintf(stderr, "Error: no maximum found for TAC '%s'.\n", tac.c[ci].name);
      continue;
    }
    if(verbose>2) {
      printf("  maxv := %g\n", ymax);
      printf("  maxx := %g\n", tac.x[imax]);
    }
    par.r[ci].p[0]=ymax;
    par.r[ci].p[1]=tac.x[imax];
 
    /* If requested, refine peak time to activity weighted avg */
    if(timeAsWAvg && imax>0 && imax<tac.sampleNr-1) {
      double s=0.0, ws=0.0;
      int i1, i2; i1=imax-1; i2=imax+1; if(i1>0 && i2<tac.sampleNr-1) {i1--; i2++;}
      for(int i=i1; i<=i2; i++) {s+=tac.c[ci].y[i]*tac.x[i]; ws+=tac.c[ci].y[i];}
      s/=ws; if(verbose>2) printf("  weighted_maxx := %g\n", s);
      if(!isnan(s)) par.r[ci].p[1]=s;
    }
 
    /* If requested, search time after peak when TAC drops to specified level */
    if(!isnan(decr)) {
      double ylim=(decr/100.0)*ymax;
      double x1=tac.x[imax];
      double y1=tac.c[ci].y[imax];
      double x2, y2; x2=y2=nan("");
      for(int i=imax+1; i<tac.sampleNr; i++) {
        if(isnan(tac.x[i])) continue;
        if(isnan(tac.c[ci].y[i])) continue;
        if(tac.c[ci].y[i]>ylim) {x1=tac.x[i]; y1=tac.c[ci].y[i]; continue;}
        x2=tac.x[i]; y2=tac.c[ci].y[i];
        break;
      }
      if(isnan(x2) || isnan(y2)) {
        fprintf(stderr, "Error: TAC '%s' does not drop to specified level.\n", tac.c[ci].name);
      } else {
        if(verbose>2) {
          printf("  prev_sample := (%g, %g)\n", x1, y1);
          printf("  next_sample := (%g, %g)\n", x2, y2);
        }
        if(fabs(y2-y1)>1.0E-20) par.r[ci].p[2]=x1+(ylim-y1)*(x2-x1)/(y2-y1);
        else par.r[ci].p[2]=0.5*(x1+x2);
      }
    }
 
    /* If requested, compute TAC mean and subtract it from the peak value */
    if(subtractWAvg) {
      double auc=tacAUC(&tac, ci, xmin, xmax, NULL);
      if(isnan(auc)) {
        fprintf(stderr, "Error: invalid mean of TAC '%s'.\n", tac.c[ci].name);
      } else {
        double avg=auc/(xmax-xmin);
        par.r[ci].p[0]-=avg;
      }
    }
 
  } // next TAC
  tacFree(&tac);
 
 
  /* 
   *  Print and save parameters
   */
  if(verbose>0 || !parfile[0]) parWrite(&par, stdout, PAR_FORMAT_TSV_UK, 1, NULL);
  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); parFree(&par);
    if(ret!=TPCERROR_OK) {
      fprintf(stderr, "Error: %s\n", errorMsg(status.error));
      return(12);
    }
    if(verbose>0) printf("parameters saved in %s\n", parfile);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/