tacorder.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacVerifyTimeOrder(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check that required data exists */
  if(d==NULL || d->sampleNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  /* All sample times must be available */
  if(tacXNaNs(d)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
  /* Less than two samples are always in order */
  if(d->sampleNr<2) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
 
  int fi;
  double p, s;
 
  fi=0; if(d->isframe==0) p=d->x[fi]; else p=0.5*(d->x1[fi]+d->x2[fi]);
  for(fi=1; fi<d->sampleNr; fi++) {
    if(d->isframe==0) s=d->x[fi]; else s=0.5*(d->x1[fi]+d->x2[fi]);
    if(s<p) {
      if(verbose>2) printf("x[%d]=%g\nx[%d]=%g\n", fi-1, p, fi, s);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OVERLAPPING_DATA);
      return(TPCERROR_OVERLAPPING_DATA);
    }
    p=s;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSortByTime(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check that required data exists */
  if(d==NULL || d->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  /* All sample times must be available */
  if(tacXNaNs(d)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(d->sampleNr<1) return TPCERROR_OK;
  
  TACC *c;
  int ri, fi, fj;
  double s, a1, a2;
  for(fi=0; fi<d->sampleNr-1; fi++) for(fj=fi+1; fj<d->sampleNr; fj++) {
    if(d->isframe==0) {a1=d->x[fi]; a2=d->x[fj];}
    else {a1=0.5*(d->x1[fi]+d->x2[fi]); a2=0.5*(d->x1[fj]+d->x2[fj]);}
    if(a2>=a1) continue;
    s=d->x[fi];  d->x[fi]=d->x[fj];   d->x[fj]=s;
    s=d->x1[fi]; d->x1[fi]=d->x1[fj]; d->x1[fj]=s;
    s=d->x2[fi]; d->x2[fi]=d->x2[fj]; d->x2[fj]=s;
    s=d->w[fi];  d->w[fi]=d->w[fj];   d->w[fj]=s;
    for(ri=0; ri<d->tacNr; ri++) {
      c=d->c+ri;
      s=c->y[fi];  c->y[fi]=c->y[fj];   c->y[fj]=s;
    }
  }  
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSortByConc(
  TAC *d,
  const int i,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check that required data exists */
  if(d==NULL || d->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  if(i<0 || i>=d->_tacNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(d->sampleNr<1) return TPCERROR_OK;
 
  TACC *c;
  int ri, fi, fj;
  double s, a1, a2;
  for(fi=0; fi<d->sampleNr-1; fi++) for(fj=fi+1; fj<d->sampleNr; fj++) {
    a1=d->c[i].y[fi]; a2=d->c[i].y[fj];
    if(a2<=a1) continue;
    s=d->x[fi];  d->x[fi]=d->x[fj];   d->x[fj]=s;
    s=d->x1[fi]; d->x1[fi]=d->x1[fj]; d->x1[fj]=s;
    s=d->x2[fi]; d->x2[fi]=d->x2[fj]; d->x2[fj]=s;
    s=d->w[fi];  d->w[fi]=d->w[fj];   d->w[fj]=s;
    for(ri=0; ri<d->_tacNr; ri++) {
      c=d->c+ri;
      s=c->y[fi];  c->y[fi]=c->y[fj];   c->y[fj]=s;
    }
  }
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
 
static int tacQSortName(const void *c1, const void *c2)
{
  int res;
  res=strcasecmp( ((TACC*)c1)->name, ((TACC*)c2)->name );
  return(res);
}
 
int tacSortByName(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check that required data exists */
  if(d==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(d->tacNr<2) return TPCERROR_OK;
  qsort(d->c, d->tacNr, sizeof(TACC), tacQSortName);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSortByAUC(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  /* Check that required data exists */
  if(d==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  if(d->tacNr<2) return TPCERROR_OK;
  /* Calculate AUCs for each TAC */
  double auc[d->tacNr], fdur, v;
  int i, j;
  for(i=0; i<d->tacNr; i++) {
    auc[i]=0.0;
    for(j=0; j<d->sampleNr; j++) {
      fdur=d->x2[0]-d->x1[0]; v=d->c[i].y[j]; 
      if(isnan(fdur) || isnan(v)) continue;
      auc[i]+=fdur*v;
    }
  }
  /* Sort */
  double a;
  for(i=0; i<d->tacNr-1; i++)
    for(j=i+1; j<d->tacNr; j++)
      if(auc[j]>auc[i]) {
        tacSwapTACCs(d, j, i); 
        a=auc[i]; auc[i]=auc[j]; auc[j]=a;
      }
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacMoveTACC(
  TAC *d,
  int from,
  int to
) {
  int i;
  TACC tacc;
  size_t taccSize;
 
  if(d==NULL || from<0 || to<0) return(TPCERROR_FAIL);
  if(from>=d->_tacNr || to>=d->_tacNr) return(TPCERROR_FAIL);
  if(from==to) return(TPCERROR_OK); // nothing to do
  taccSize=sizeof(TACC);
  memcpy(&tacc, d->c+from, taccSize);
  if(from>to) for(i=from; i>to; i--)
    memcpy(d->c+i, d->c+(i-1), taccSize);
  else for(i=from; i<to; i++)
    memcpy(d->c+i, d->c+(i+1), taccSize);
  memcpy(d->c+i, &tacc, taccSize);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSwapTACCs(
  TAC *d,
  int i1,
  int i2
) {
  TACC tacc;
  size_t taccSize;
 
  if(d==NULL || i1<0 || i2<0) return(TPCERROR_FAIL);
  if(i1>=d->_tacNr || i2>=d->_tacNr) return(TPCERROR_FAIL);
  if(i1==i2) return(TPCERROR_OK); // nothing to do
  taccSize=sizeof(TACC);
  memcpy(&tacc, d->c+i1, taccSize);
  memcpy(d->c+i1, d->c+i2, taccSize);
  memcpy(d->c+i2, &tacc, taccSize);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacDeleteTACC(
  TAC *d,
  int i
) {
  /* Check that region exists */
  if(d==NULL || i<0 || i>=d->tacNr) return(TPCERROR_FAIL);
  /* If the last TACC is to be deleted, then just decrease the tacNr */
  if(i==d->tacNr-1) {d->tacNr--; return(TPCERROR_OK);}
  /* Otherwise move it to the last position, and then decrease tacNr */
  if(tacMoveTACC(d, i, d->tacNr-1)!=TPCERROR_OK) return(TPCERROR_FAIL);
  d->tacNr--;
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacMultipleSamples(
  TAC *d1,
  const int fixMode,
  TAC *d2,
  const int verbose
) {
  if(verbose>0) {printf("\n%s(%d)\n", __func__, fixMode); fflush(stdout);}
 
  /* Check data */
  if(fixMode!=0 && (d1==NULL || d2==NULL || d1->sampleNr<1)) return(2);
  if(d1==NULL || d1->sampleNr<2 || !tacIsX(d1)) return(0);
 
  /* Set limit for similarity of sample times; the same sample times can have been saved in 
     files with variable precision. */
  const double simlim=1.0E-03; if(verbose>1) printf("  limit := %g\n", simlim);
 
  /* Check if there are any multiple samples */
  int isMultiple=0;
  for(int i=0; i<d1->sampleNr-1 && !isMultiple; i++) {
    double t1; if(d1->isframe) t1=0.5*(d1->x1[i]+d1->x2[i]); else t1=d1->x[i];
    for(int j=i+1; j<d1->sampleNr; j++) {
      double t2; if(d1->isframe) t2=0.5*(d1->x1[j]+d1->x2[j]); else t2=d1->x[j];
      double dt=fabs(t2-t1); 
      if(verbose>4) printf("  |%g - %g| = %g\n", t1, t2, dt);
      if(dt<simlim) {isMultiple=1; break;}
    }
  }
  if(verbose>1) {printf("  isMultiple := %d\n", isMultiple); fflush(stdout);}
  if(fixMode==0) return(isMultiple); // done, if fixing not asked for
 
  /* If there were no multiples, copy the input to output and quit */
  if(!isMultiple) {
    TAC itac; tacInit(&itac); 
    if(tacDuplicate(d1, &itac)!=TPCERROR_OK) return(3);
    if(tacDuplicate(&itac, d2)!=TPCERROR_OK) {tacFree(&itac); return(3);}
    tacFree(&itac); return(0);
  }
 
  /* Ok there are multiples that need to be fixed. */
  /* We'll start by making a copy of the input data, and removing any existing data in 
     the output; copy is needed, because output structure can be the same as input */
  TAC itac; tacInit(&itac); if(tacDuplicate(d1, &itac)!=TPCERROR_OK) return(4);
  tacFree(d2);
  /* Delete missing samples from the input */
  tacDeleteMissingSamples(&itac); if(verbose>2) printf("  itac.sampleNr := %d\n", itac.sampleNr);
  /* Allocate memory for the output data */
  if(tacAllocate(d2, itac.sampleNr, itac.tacNr)!=TPCERROR_OK) {tacFree(&itac); return(5);}
  if(tacCopyHdr(&itac, d2)!=TPCERROR_OK) {tacFree(&itac); return(6);}
  for(int r=0; r<itac.tacNr; r++)
    if(tacCopyTacchdr(&itac.c[r], &d2->c[r])!=TPCERROR_OK) {tacFree(&itac); return(7);}
  d2->weighting=WEIGHTING_OFF;
  d2->tacNr=itac.tacNr;
 
  /* Go through the input data */
  while(itac.sampleNr>0) {
    /* Find the smallest sample time */
    int i=tacMinX(&itac); if(i<0) break;
    if(verbose>3) printf("  smallest x[%d]=%g\n", i, itac.x[i]);
    /* Copy values to the first free output sample */
    for(int r=0; r<itac.tacNr; r++) d2->c[r].y[d2->sampleNr]=0.0;
    int rn[itac.tacNr]; for(int r=0; r<itac.tacNr; r++) rn[r]=0;
    for(int r=0; r<itac.tacNr; r++)
      if(isfinite(itac.c[r].y[i])) {
        d2->c[r].y[d2->sampleNr]=itac.c[r].y[i];
        rn[r]++;
      }
    d2->x[d2->sampleNr]=itac.x[i]; d2->x1[d2->sampleNr]=itac.x1[i]; d2->x2[d2->sampleNr]=itac.x2[i];
    double t1; if(itac.isframe) t1=0.5*(itac.x1[i]+itac.x2[i]); else t1=itac.x[i];
    /* Delete the sample from input */
    tacDeleteSample(&itac, i);
    /* Search the input for matching sample times */
    while(itac.sampleNr>0) {
      /* Find the smallest sample time left */
      int i=tacMinX(&itac); if(i<0) break;
      if(verbose>3) printf("    smallest x[%d]=%g\n", i, itac.x[i]);
      double t2; if(itac.isframe) t2=0.5*(itac.x1[i]+itac.x2[i]); else t2=itac.x[i];
      double dt=fabs(t2-t1); if(verbose>4) printf("  |%g - %g| = %g\n", t1, t2, dt);
      if(!(dt<simlim)) break; // if this is not close enough, then none will be
      /* Add values to the first free output sample */
      for(int r=0; r<itac.tacNr; r++)
        if(isfinite(itac.c[r].y[i])) {
          d2->c[r].y[d2->sampleNr]+=itac.c[r].y[i];
          rn[r]++;
        }
      /* Delete also this sample from input */
      tacDeleteSample(&itac, i);
    }
    /* Divide values by the summed number of samples */
    for(int r=0; r<itac.tacNr; r++) d2->c[r].y[d2->sampleNr]/=(double)rn[r];
    /* Move to next output sample */
    d2->sampleNr++;
  }
  tacFree(&itac);
 
  /* Check that we got something left */
  if(d2->sampleNr<1) return(10);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/