taccat.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
/*****************************************************************************/
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Catenate the TACs from the second TAC file to the first one, or to a new",
  "datafile. Sample times or correction for physical decay are not changed.",
  "Each TAC file must contain the same number of TACs, and if more than one TAC,",
  "TAC names must match, too.",
  " ",
  "Usage: @P [options] tacfile1 tacfile2 [catenated_file]",
  " ",
  "Options:",
  " -both | -first | -second | -cut=<time>",
  "     In case of overlapping samples, either samples from both (-both),",
  "     first (-first), or second (-second, default) TAC are saved in",
  "     combined file, or specified cut time is used.",
  "     Cut time can also be given in file with keys 'x' or 'time'.",
  "     By default, overlap will lead to an error.",
  " --force",
  "     Program does not mind if the time or calibration units cannot be",
  "     converted to match, or if TAC names do not match.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "Example:",
  "     @P t455ap_pump.kbq t455ap_manual.kbq t455ap_combined.kbq",
  " ",
  "See also: tacblend, taccut, tacadd0, tacadd, tacunit, tactime, taccross",
  " ",
  "Keywords: TAC, tool, input, blood",
  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;
/*****************************************************************************/
 
/*****************************************************************************/
enum {OVL_ERROR, OVL_USE_BOTH, OVL_USE_FIRST, OVL_USE_SECOND, OVL_CUTTIME};

int main(int argc, char **argv)
{
  int ai, help=0, version=0, verbose=1;
  char tacfile1[FILENAME_MAX], tacfile2[FILENAME_MAX], outfile[FILENAME_MAX];
  int forceMode=0;
  int overlapMode=OVL_ERROR;
  double cutTime=nan("");
 
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  tacfile1[0]=tacfile2[0]=outfile[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(strcasecmp(cptr, "F")==0 || strcasecmp(cptr, "FORCE")==0) {
      forceMode=1; continue;
    } else if(strncasecmp(cptr, "BOTH", 1)==0) {
      overlapMode=OVL_USE_BOTH; continue;
    } else if(strncasecmp(cptr, "FIRST", 2)==0) {
      overlapMode=OVL_USE_FIRST; continue;
    } else if(strncasecmp(cptr, "SECOND", 2)==0) {
      overlapMode=OVL_USE_SECOND; continue;
    } else if(strncasecmp(cptr, "CUT=", 4)==0) {
      overlapMode=OVL_CUTTIME;
      cptr+=4;
      /* Try as a value first */
      cutTime=atofVerified(cptr);
      if(cutTime>0.0) continue;
      /* Then try as a IFT file */
      {
        FILE *fp=fopen(cptr, "r");
        if(fp==NULL) {fprintf(stderr, "Error: invalid cut time.\n"); return(1);}
        IFT ift; iftInit(&ift);
        int ret=iftRead(&ift, fp, 1, 0, NULL);
        fclose(fp);
        if(ret==0) {
          if(iftGetDoubleValue(&ift, "x", 0, &cutTime)>=0 && cutTime>0.0) {iftFree(&ift); continue;}
          if(iftGetDoubleValue(&ift, "time", 0, &cutTime)>=0 && cutTime>0.0) {iftFree(&ift); continue;}
          if(iftGetDoubleValue(&ift, "cut", 0, &cutTime)>=0 && cutTime>0.0) {iftFree(&ift); continue;}
          iftFree(&ift);
        }
      }
      fprintf(stderr, "Error: invalid cut time.\n"); return(1);
    }
    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-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(tacfile1, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(tacfile2, 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(!tacfile2[0]) {tpcPrintUsage(argv[0], info, stdout); return(1);}
  if(!outfile[0]) strcpy(outfile, tacfile1);
  if(overlapMode==OVL_CUTTIME && !(cutTime>0.0)) {fprintf(stderr, "Error: invalid cut time.\n"); 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("tacfile1 := %s\n", tacfile1);
    printf("tacfile2 := %s\n", tacfile2);
    printf("outfile := %s\n", outfile);
    printf("overlapMode := %d\n", overlapMode);
    if(overlapMode==OVL_CUTTIME) printf("cutTime := %g\n", cutTime);
    printf("forceMode := %d\n", forceMode);
    fflush(stdout);
  }
 
 
  /*
   *  Read the files
   */
  TAC tac1, tac2;
  tacInit(&tac1); tacInit(&tac2);
  if(verbose>1) printf("reading %s\n", tacfile1);
  if(tacRead(&tac1, tacfile1, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); tacFree(&tac2); return(2);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac1.format));
    printf("tacNr := %d\n", tac1.tacNr);
    printf("sampleNr := %d\n", tac1.sampleNr);
    printf("xunit := %s\n", unitName(tac1.tunit));
    printf("yunit := %s\n", unitName(tac1.cunit));
  }
  if(verbose>1) printf("reading %s\n", tacfile2);
  if(tacRead(&tac2, tacfile2, &status)!=TPCERROR_OK) {
    fprintf(stderr, "Error: %s\n", errorMsg(status.error));
    tacFree(&tac1); tacFree(&tac2); return(3);
  }
  if(verbose>2) {
    printf("fileformat := %s\n", tacFormattxt(tac2.format));
    printf("tacNr := %d\n", tac2.tacNr);
    printf("sampleNr := %d\n", tac2.sampleNr);
    printf("xunit := %s\n", unitName(tac2.tunit));
    printf("yunit := %s\n", unitName(tac2.cunit));
  }
 
  /* Check if the TAC number is different */
  if(verbose>1) printf("checking data...\n");
  if(tac1.tacNr!=tac2.tacNr) {
    fprintf(stderr, "Error: different number of TACs.\n");
    tacFree(&tac1); tacFree(&tac2); return(4);
  }
 
  /* Check that TAC names match, if more than one TAC */
  if(tac1.tacNr>1 && tacCompareNames(&tac1, &tac2, -1, &status)!=0) {
    if(forceMode) {
      fprintf(stderr, "Warning: TAC names do not match.\n");
    } else {
      /* Without force mode this is an error, unless sorting helps */
      tacSortByName(&tac1, NULL);
      tacSortByName(&tac2, NULL);
      if(tacCompareNames(&tac1, &tac2, -1, &status)!=0) {
        fprintf(stderr, "Error: TAC names do not match.\n");
        tacFree(&tac1); tacFree(&tac2); return(4);
      }
      fprintf(stderr, "Note: TACs sorted by name.\n");
    }
  }
 
  /* Delete missing samples */
  tacDeleteMissingSamples(&tac1);
  if(tacNotNaNs(&tac1, -1)<1) {
    fprintf(stderr, "Error: no valid samples in %s\n", tacfile1);
    tacFree(&tac1); tacFree(&tac2); return(2);
  }
  tacDeleteMissingSamples(&tac2);
  if(tacNotNaNs(&tac2, -1)<1) {
    fprintf(stderr, "Error: no valid samples in %s\n", tacfile2);
    tacFree(&tac1); tacFree(&tac2); return(3);
  }
 
  /* Sort by sample times */
  tacSortByTime(&tac1, NULL);
  tacSortByTime(&tac2, NULL);
 
  /* Check time units */
  if(tac1.tunit==UNIT_UNKNOWN || tac2.tunit==UNIT_UNKNOWN) {
    fprintf(stderr, "Warning: unknown time units.\n");
    if(forceMode) { // In force mode use what we got
      if(tac1.tunit!=UNIT_UNKNOWN) tac2.tunit=tac1.tunit;
      else if(tac2.tunit!=UNIT_UNKNOWN) tac1.tunit=tac2.tunit;
    }
  }
  if(tac2.tunit!=UNIT_UNKNOWN) {
    /* Try to convert time units to time units in the first file */
    int ret=tacXUnitConvert(&tac2, tac1.tunit, &status);
    if(ret!=TPCERROR_OK) {
      if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
      if(forceMode) {
        fprintf(stderr, "Warning: non-compatible TAC time units.\n");
      } else {
        fprintf(stderr, "Error: non-compatible TAC time units.\n");
        tacFree(&tac1); tacFree(&tac2); return(4);
      }
    }
  }
  /* If either TAC file only has frame mid times, then both have */
  if(tac1.isframe==0) tac2.isframe=0;
  else if(tac2.isframe==0) tac1.isframe=0;
 
  /* Check concentration units */
  if(tac1.cunit==UNIT_UNKNOWN || tac2.cunit==UNIT_UNKNOWN) {
    fprintf(stderr, "Warning: unknown concentration units.\n");
    if(forceMode) { // In force mode use what we got
      if(tac1.cunit!=UNIT_UNKNOWN) tac2.cunit=tac1.cunit;
      else if(tac2.cunit!=UNIT_UNKNOWN) tac1.cunit=tac2.cunit;
    }
  }
  if(tac2.cunit!=UNIT_UNKNOWN) {
    /* Try to convert concentration units to units in the first file */
    int ret=tacYUnitConvert(&tac2, tac1.cunit, &status);
    if(ret!=TPCERROR_OK) {
      if(verbose>2) fprintf(stderr, "Status: %s\n", errorMsg(status.error));
      if(forceMode) {
        fprintf(stderr, "Warning: non-compatible TAC concentration units.\n");
      } else {
        fprintf(stderr, "Error: non-compatible TAC concentration units.\n");
        tacFree(&tac1); tacFree(&tac2); return(4);
      }
    }
  }
 
  /* If one of TAC files only has frame mid times, then we only use mid times */
  if(tac1.isframe==0 || tac2.isframe==0) tac1.isframe=tac2.isframe=0;
  /* Remove weighting unless both have similar weighting */
  if(tac1.weighting!=tac2.weighting) tac1.weighting=tac2.weighting=WEIGHTING_OFF;
 
 
  /*
   *  Check if there is overlap in sample times of the two files
   */
  if(verbose>1) printf("checking for time overlap\n");
  int isOverlap=0;
  double xmin1, xmax1, xmin2, xmax2, xmin, xmax;
  if(tacXRange(&tac1, &xmin1, &xmax1) || tacXRange(&tac2, &xmin2, &xmax2)) {
    fprintf(stderr, "Error: invalid sample times.\n");
    tacFree(&tac1); tacFree(&tac2); return(5);
  }
  if(verbose>3) {
    printf("x range 1 := %g - %g\n", xmin1, xmax1);
    printf("x range 2 := %g - %g\n", xmin2, xmax2);
  }
  if(xmax1>xmin2) isOverlap=1;
  else if(!tac1.isframe && xmax1>=xmin2) isOverlap=1;
  if(isOverlap && overlapMode==OVL_ERROR) {
    fprintf(stderr, "Error: overlapping sample times.\n");
    tacFree(&tac1); tacFree(&tac2); return(5);
  }
  xmin=xmin1; if(xmin2<xmin) xmin=xmin2;
  xmax=xmax1; if(xmax2>xmax) xmax=xmax2;
  if(overlapMode==OVL_CUTTIME && (cutTime<xmin || cutTime>xmax)) {
    fprintf(stderr, "Error: cut time outside of data time range %g - %g\n", xmin, xmax);
    tacFree(&tac1); tacFree(&tac2); return(5);
  }
  if(isOverlap) {
    int fail=0;
    if(overlapMode==OVL_USE_FIRST && xmin1>xmin2) fail=1;
    if(overlapMode==OVL_USE_SECOND && xmin2<=xmin1) fail=2;
    if(overlapMode==OVL_CUTTIME && (xmin1>cutTime || xmax2<cutTime)) fail=3;
    if(fail) {
      fprintf(stderr, "Error: incompatible sample time ranges or wrong file order.\n");
      tacFree(&tac1); tacFree(&tac2); return(5);
    }
  }
  /* Do not care about overlap if all data is used anyway */
  if(isOverlap) {
    if(overlapMode==OVL_USE_BOTH) {
      isOverlap=0;
      fprintf(stderr, "Warning: overlap in sample times is ignored.\n");
    } else {
      fprintf(stderr, "Warning: overlap in sample times will be removed.\n");
    }
  }
 
  /* Set cut time based on TAC preference, if that is given */
//  if(overlapMode==OVL_USE_FIRST) cutTime=xmax1;
//  if(overlapMode==OVL_USE_SECOND) cutTime=xmin2;
  if(overlapMode==OVL_USE_FIRST) {cutTime=xmax1; if(xmax1==xmin2) cutTime+=1.0E-12;}
  if(overlapMode==OVL_USE_SECOND) {cutTime=xmin2; if(xmax1==xmin2) cutTime-=1.0E-12;}
  if(verbose>4 && isOverlap && (overlapMode==OVL_USE_FIRST || overlapMode==OVL_USE_SECOND))
    printf("cutTime := %g\n", cutTime);
 
 
  /* Delete overlapping part from first dataset, if necessary */
  if(isOverlap && overlapMode!=OVL_USE_FIRST && xmax1>=cutTime) {
    if(verbose>1) printf("deleting overlap from the first dataset\n");
    if(tac1.isframe) {
      for(int i=tac1.sampleNr-1; i>=0; i--) 
        if(tac1.x2[i]>cutTime) tacDeleteSample(&tac1, i);
    } else {
      for(int i=tac1.sampleNr-1; i>=0; i--) 
        if(tac1.x[i]>cutTime) tacDeleteSample(&tac1, i);
    }
    if(tac1.sampleNr<1 || tacXRange(&tac1, &xmin1, &xmax1)) {
      fprintf(stderr, "Error: no data from file 1 included.\n");
      tacFree(&tac1); tacFree(&tac2); return(5);
    }
    if(verbose>3) printf("final x range 1 := %g - %g\n", xmin1, xmax1);
  }
  /* Delete overlapping part from second dataset, if necessary */
  if(isOverlap && overlapMode!=OVL_USE_SECOND && xmin2<=cutTime) {
    if(verbose>1) printf("deleting overlap from the second dataset\n");
    if(tac1.isframe) {
      for(int i=tac2.sampleNr-1; i>=0; i--)
        if(tac2.x1[i]<cutTime) tacDeleteSample(&tac2, i);
    } else {
      for(int i=tac2.sampleNr-1; i>=0; i--)
        if(tac2.x[i]<cutTime) tacDeleteSample(&tac2, i);
    }
    if(tac2.sampleNr<1 || tacXRange(&tac2, &xmin2, &xmax2)) {
      fprintf(stderr, "Error: no data from file 2 included.\n");
      tacFree(&tac1); tacFree(&tac2); return(5);
    }
    if(verbose>3) printf("final x range 2 := %g - %g\n", xmin2, xmax2);
  }
 
 
  /*
   *  Add data from second file in to the first one
   */
  if(verbose>1) printf("combining data\n");
  if(tacAllocateMoreSamples(&tac1, tac2.sampleNr)) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    tacFree(&tac1); tacFree(&tac2); return(6);
  }
  for(int i=0; i<tac2.sampleNr; i++) {
    tac1.x[tac1.sampleNr]=tac2.x[i];
    tac1.x1[tac1.sampleNr]=tac2.x1[i];
    tac1.x2[tac1.sampleNr]=tac2.x2[i];
    for(int j=0; j<tac1.tacNr; j++) tac1.c[j].y[tac1.sampleNr]=tac2.c[j].y[i];
    tac1.w[tac1.sampleNr]=tac2.w[i];
    tac1.sampleNr++;
  }
  /* Sort by sample time */
  tacSortByTime(&tac1, NULL);
 
  /* Add interfile type header contents that are stored as comments in file 2 */
  iftCopyItems(&tac1.h, &tac2.h, 1, 1, 2, NULL);
  /* Remove duplicate header fields */
  iftDeleteDuplicateKeys(&tac1.h, NULL);
 
  /* No need for the second dataset */
  tacFree(&tac2);
 
  /*
   *  Save data 
   */
  if(verbose>1) printf("writing %s\n", outfile);
  {
  FILE *fp; fp=fopen(outfile, "w");
  if(fp==NULL) {
    fprintf(stderr, "Error: cannot open file for writing (%s)\n", outfile);
    tacFree(&tac1); return(11);
  }
  int ret=tacWrite(&tac1, fp, TAC_FORMAT_UNKNOWN, 1, &status);
  fclose(fp);
  if(ret!=TPCERROR_OK) {
    fprintf(stderr, "Error (%d): %s\n", ret, errorMsg(status.error));
    tacFree(&tac1); return(12);
  }
  if(verbose>=0) {
    printf("%d samples saved in %s.\n", tac1.sampleNr, outfile);
  }
  }
  tacFree(&tac1); 
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/