tacx.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include "tpcift.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpctac.h"
/*****************************************************************************/
 
/*****************************************************************************/
int tacXCopy(
  TAC *tac1,
  TAC *tac2,
  int i1,
  int i2
) {
  if(tac1==NULL || tac2==NULL) return(TPCERROR_FAIL);
  if(i1<0 || i1>i2) return(TPCERROR_FAIL);
  if(i2>=tac1->_sampleNr || i2>=tac2->_sampleNr) return(TPCERROR_FAIL);
  
  for(int i=i1; i<=i2; i++) {
    tac2->x[i]=tac1->x[i];
    tac2->x1[i]=tac1->x1[i];
    tac2->x2[i]=tac1->x2[i];
  }
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacCorrectFrameOverlap(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(d==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  /* If no frame start and end times, then just check the order */
  if(d->isframe==0) {
    if(tacVerifyTimeOrder(d, status)) {
      if(verbose>0) fprintf(stderr, "Error: sample times are not in order.\n");
      return(status->error);
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
  /* Check each frame */  
  double overlap, flen1, flen2, overlap_limit=0.0;
  for(int fi=0; fi<d->sampleNr-1; fi++) {
    overlap=d->x2[fi]-d->x1[fi+1];
    if(overlap==0.0) continue; // no gap or overlap
    /* Calculate the frame length of current frame and the next frame */
    flen1=d->x2[fi]-d->x1[fi]; flen2=d->x2[fi+1]-d->x1[fi+1];
    if(flen1<0.0 || flen2<0.0) return(1);
    /* Set the limit */
    if(flen1<flen2) overlap_limit=0.2*flen1; else overlap_limit=0.2*flen2;
    /* Check if gap or overlap is too large to be fixed automatically */
    if(overlap<-overlap_limit) continue; // gap is too large, then do nothing
    if(overlap>overlap_limit) { // overlap is too large: error
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OVERLAPPING_DATA);
      return(TPCERROR_OVERLAPPING_DATA);
    }
    /* Correct the small gap/overlap by making frame durations more similar */
    if(overlap>0.0) { // overlap
      if(flen1>flen2) d->x2[fi]=d->x1[fi+1]; else d->x1[fi+1]=d->x2[fi];
    } else { // gap
      if(flen1>flen2) d->x1[fi+1]=d->x2[fi]; else d->x2[fi]=d->x1[fi+1];
    }
  }
  
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacXRange(
  TAC *d,
  double *xmin,
  double *xmax
) {
  if(xmin!=NULL) *xmin=nan("");
  if(xmax!=NULL) *xmax=nan("");
  /* Check the data */
  if(d==NULL || d->sampleNr<1) return(1);
  /* Find the min and max time */
  double mi, ma, *xi, *xa;
  mi=ma=nan("");
  if(d->isframe) {xi=d->x1; xa=d->x2;} else xi=xa=d->x;
  for(int i=0; i<d->sampleNr; i++) {
    if(!isfinite(xi[i]) || !isfinite(xa[i])) continue;
    if(isnan(mi) || isnan(ma)) {mi=xi[i]; ma=xa[i]; continue;}
    if(xi[i]<mi) mi=xi[i]; else if(xa[i]>ma) ma=xa[i];
  }
  if(xmin!=NULL) *xmin=mi;
  if(xmax!=NULL) *xmax=ma;
  if(!isfinite(mi) || !isfinite(ma)) return(2);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSampleXRange(
  TAC *d,
  double *xmin,
  double *xmax
) {
  if(xmin!=NULL) *xmin=nan("");
  if(xmax!=NULL) *xmax=nan("");
  /* Check the data */
  if(d==NULL || d->sampleNr<1) return(1);
  /* Find the min and max time */
  double mi, ma, *xi, *xa;
  mi=ma=nan("");
  if(d->isframe) {xi=d->x1; xa=d->x2;} else xi=xa=d->x;
  int weighted=tacIsWeighted(d);
  int i, j;
  for(i=0; i<d->sampleNr; i++) {
    if(!isfinite(xi[i]) || !isfinite(xa[i])) continue;
    if(weighted && !(d->w[i]>0.0)) continue;
    for(j=0; j<d->tacNr; j++) if(isfinite(d->c[j].y[i])) break;
    if(j==d->tacNr) continue;
    if(isnan(mi) || isnan(ma)) {mi=xi[i]; ma=xa[i]; continue;}
    if(xi[i]<mi) mi=xi[i]; else if(xa[i]>ma) ma=xa[i];
  }
  if(xmin!=NULL) *xmin=mi;
  if(xmax!=NULL) *xmax=ma;
  if(!isfinite(mi) || !isfinite(ma)) return(2);
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacMinX(
  TAC *d
) {
  if(d==NULL || !(d->sampleNr>0)) return(-1);
  int imin=-1; double xmin=nan("");
  for(int i=0; i<d->sampleNr; i++) {
    double t1; if(d->isframe) t1=0.5*(d->x1[i]+d->x2[i]); else t1=d->x[i];
    if(!isfinite(t1)) continue;
    if(!isfinite(xmin) || imin<0 || t1<xmin) {xmin=t1; imin=i; continue;}
  }
  return(imin);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacIsX(
  TAC *d
) {
  if(d==NULL || d->sampleNr<1) return(0);
  double a, b;
  if(tacXRange(d, &a, &b)) return(0);
  if(fabs(a)<1.0E-30 && fabs(b)<1.0E-30) return(0);
  if(fabs(a)<1.0E-30 && fabs(a-b)<1.0E-30) return(0);
  return(1);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacXMatch(
  TAC *d1,
  TAC *d2,
  const int verbose
) {
  if(verbose>0) printf("\n%s()\n", __func__);
 
  /* Check data */
  if(d1==NULL || d2==NULL) return 0;
  double xmin[2], xmax[2];
  if(tacXRange(d1, xmin, xmax) || tacXRange(d2, xmin+1, xmax+1)) return(0);
 
  /* Which has less samples? */
  int smaller_sampleNr;
  if(d1->sampleNr<d2->sampleNr) smaller_sampleNr=d1->sampleNr;
  else smaller_sampleNr=d2->sampleNr;
  if(verbose>1) printf("smaller_sampleNr := %d\n", smaller_sampleNr);
  if(smaller_sampleNr<=0) return 0;
  
  /* Get conversion factors for sample times to get those in sec */
  double ts1, ts2;
  ts1=unitConversionFactor(d1->tunit, UNIT_SEC);
  if(isnan(ts1)) {
    if(verbose>1) printf("cannot convert tac1 sample times to sec.\n");
    ts1=1.0; if(d1->tunit==UNIT_UNKNOWN) ts1=60.0; // assuming min
  }
  ts2=unitConversionFactor(d2->tunit, UNIT_SEC);
  if(isnan(ts2)) {
    if(verbose>1) printf("cannot convert tac2 sample times to sec.\n");
    ts2=1.0; if(d2->tunit==UNIT_UNKNOWN) ts2=60.0; // assuming min
  }
  /* Convert time ranges to sec */
  xmin[0]*=ts1; xmax[0]*=ts1;
  xmin[1]*=ts2; xmax[1]*=ts2;
  if(verbose>1) {
    printf("tac1->isframe := %d\n", d1->isframe);
    if(verbose>2) printf("time_range1 := %g - %g\n", xmin[0], xmax[0]);
    printf("tac2->isframe := %d\n", d2->isframe);
    if(verbose>2) printf("time_range2 := %g - %g\n", xmin[1], xmax[1]);
  }
 
  /* Set the accepted time difference */
  double f, accepted_timedif=2.2; // sec
  /* With short study the accepted time difference must be shorter */
  f=xmax[0]-xmin[0]; if((xmax[1]-xmin[1])>f) f=xmax[1]-xmin[1];
  f*=0.01; if(accepted_timedif>f) accepted_timedif=f;
  if(verbose>1) {
    printf("accepted_timedif := %g [s]\n", accepted_timedif);
    fflush(stdout);
  }
 
  /* Compare sample times frame-by-frame */
  int fi, n=0;
  if(d1->isframe==0 && d2->isframe==0) {
    if(verbose>2) {printf("frame mid times\n"); fflush(stdout);}
    for(fi=0; fi<smaller_sampleNr; fi++) {
      if(d1->x[fi]*ts1==d2->x[fi]*ts2) continue;
      if(isnan(d1->x[fi]) && isnan(d2->x[fi])) continue; // match if both NaN
      if(isnan(d1->x[fi]) || isnan(d2->x[fi])) {n++; continue;} // no match
      f=fabs(d1->x[fi]*ts1-d2->x[fi]*ts2);
      if(verbose>3) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>4) printf("  %g vs %g\n", ts1*d1->x[fi], ts2*d2->x[fi]);
      if(f>accepted_timedif) n++;
    }
  } else if(d1->isframe!=0 && d2->isframe!=0) {
    if(verbose>2) {printf("frame start and end times\n"); fflush(stdout);}
    for(fi=0; fi<smaller_sampleNr; fi++) {
      if(d1->x1[fi]*ts1==d2->x1[fi]*ts2) continue;
      if(isnan(d1->x1[fi]) && isnan(d2->x1[fi])) continue; // match if both NaN
      if(isnan(d1->x1[fi]) || isnan(d2->x1[fi])) {n++; continue;} // no match
      f=fabs(d1->x1[fi]*ts1-d2->x1[fi]*ts2);
      if(verbose>3) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>4) printf("  %g vs %g\n", ts1*d1->x1[fi], ts2*d2->x1[fi]);
      if(f>accepted_timedif) n++;
    }
    for(fi=0; fi<smaller_sampleNr; fi++) {
      if(d1->x2[fi]*ts1==d2->x2[fi]*ts2) continue;
      if(isnan(d1->x2[fi]) && isnan(d2->x2[fi])) continue; // match if both NaN
      if(isnan(d1->x2[fi]) || isnan(d2->x2[fi])) {n++; continue;} // no match
      f=fabs(d1->x2[fi]*ts1-d2->x2[fi]*ts2);
      if(verbose>3) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>4) printf("  %g vs %g\n", ts1*d1->x2[fi], ts2*d2->x2[fi]);
      if(f>accepted_timedif) n++;
    }
  } else {
    if(verbose>2) {printf("mixed frame times\n"); fflush(stdout);}
    double x1, x2;
    for(fi=0; fi<smaller_sampleNr; fi++) {
      if(d1->isframe) x1=0.5*(d1->x1[fi]+d1->x2[fi]); else x1=d1->x[fi];
      if(d2->isframe) x2=0.5*(d2->x1[fi]+d2->x2[fi]); else x2=d2->x[fi];
      if(x1*ts1==x2*ts2) continue;
      if(isnan(x1) && isnan(x2)) continue; // match if both NaN
      if(isnan(x1) || isnan(x2)) {n++; continue;} // no match
      f=fabs(x1*ts1-x2*ts2);
      if(verbose>3) printf("timedif[%d] := %g\n", fi, f);
      if(verbose>4) printf("  %g vs %g\n", x1, x2);
      if(f>accepted_timedif) n++;
    }
  }
  if(verbose>2) {printf("nr of different frame times := %d\n", n); fflush(stdout);}
  if(n==0) return 1; else return 0; 
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacAddZeroSample(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(d==NULL || d->_sampleNr<1 || d->_tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
 
  /* Check if we already have the zero sample */
  if(d->sampleNr>0) {
    if((d->isframe && d->x1[0]<=0.0) || (!d->isframe && d->x[0]<=0.0)) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
      return TPCERROR_OK;
    }
  }
  
  /* If we have no samples at all, but space for samples, then just make one */
  if(d->sampleNr==0 && d->_sampleNr>0) {
    d->x[0]=d->x1[0]=d->x2[0]=0.0; d->sampleNr=1;
    for(int j=0; j<d->tacNr; j++) d->c[j].y[0]=0.0;
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
  
  /* Add space for one more sample */
  int ret;
  ret=tacAllocateMoreSamples(d, 1);
  if(ret!=TPCERROR_OK) {statusSet(status, __func__, __FILE__, __LINE__, ret); return ret;}
  
  /* Move all data to make room for the zero sample */
  for(int fi=d->sampleNr; fi>0; fi--) {
    d->x[fi]=d->x[fi-1]; d->x1[fi]=d->x1[fi-1]; d->x2[fi]=d->x2[fi-1];
    for(int ri=0; ri<d->tacNr; ri++) d->c[ri].y[fi]=d->c[ri].y[fi-1];
    d->w[fi]=d->w[fi-1];
  } // previous sample
  d->sampleNr++;
  /* ... and then add the zero sample */
  for(int ri=0; ri<d->tacNr; ri++) d->c[ri].y[0]=0.0;
  if(!d->isframe) {
    d->x[0]=d->x1[0]=d->x2[0]=0.0;
  } else {
    d->x1[0]=0.0; d->x2[0]=d->x1[1]; d->x[0]=0.5*(d->x1[0]+d->x2[0]);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacDeleteSample(
  TAC *d,
  int i
) {
  if(d==NULL || i<0 || i>=d->sampleNr) return TPCERROR_FAIL;
  int fi, ri;
  for(fi=i+1; fi<d->sampleNr; fi++) {
    d->x[fi-1]=d->x[fi]; d->x1[fi-1]=d->x1[fi]; d->x2[fi-1]=d->x2[fi];
    for(ri=0; ri<d->tacNr; ri++) d->c[ri].y[fi-1]=d->c[ri].y[fi];
    d->w[fi-1]=d->w[fi];
  }
  d->sampleNr--;
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacDeleteMissingSamples(
  TAC *d
) {
  if(d==NULL) return TPCERROR_FAIL;
  int fi, ri, n, delthat=0;
  fi=d->sampleNr-1;
  while(fi>=0) {
    delthat=0;
    /* check the x value(s) for this sample */
    if(d->isframe) {
      if(!isfinite(d->x1[fi]) || !isfinite(d->x2[fi])) delthat=1;
    } else {
      if(!isfinite(d->x[fi])) delthat=1;
    }
    /* check the y value(s) for this sample */
    for(ri=n=0; ri<d->tacNr; ri++) if(isfinite(d->c[ri].y[fi])) n++;
    if(n==0) delthat=1;
    if(delthat!=0) {
      /* delete sample */
      tacDeleteSample(d, fi);
    }
    /* go to previous sample */
    fi--;
  }
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacExtractRange(
  TAC *d1,
  TAC *d2,
  double startT,
  double endT
) {
  if(d1==NULL || d2==NULL) return TPCERROR_FAIL;
  if(d1->sampleNr<1) return TPCERROR_NO_DATA;
  if(endT<startT) return TPCERROR_INVALID_VALUE;
 
  /* Duplicate source TAC, in case target points to the same data */
  TAC tac; tacInit(&tac); if(tacDuplicate(d1, &tac)!=TPCERROR_OK) return TPCERROR_FAIL;
 
  /* Remove any old contents in target structure */
  tacFree(d2);
 
  /* Get the range of TAC frame middle times */
  double beg, end;
  if(tac.isframe) {
    beg=0.5*(tac.x1[0]+tac.x2[0]); 
    end=0.5*(tac.x1[tac.sampleNr-1]+tac.x2[tac.sampleNr-1]); 
  } else {
    beg=tac.x[0]; end=tac.x[tac.sampleNr-1];
  }
  
  /* Check whether required range is outside the TAC range */
  if(startT>=end || endT<beg) {tacFree(&tac); return TPCERROR_NO_DATA;}
 
  /* If the whole TAC is inside required range, then just copy it */
  if(startT<=beg && endT>=end) {
    int ret=tacDuplicate(&tac, d2);
    tacFree(&tac); return(ret);
  }
  
  /* Count how many samples are inside the required range */
  /* And get the first and last sample indices */
  int i, n=0, i1=-1, i2=-1;
  double x;
  for(i=0; i<tac.sampleNr; i++) {
    if(tac.isframe) x=0.5*(tac.x1[i]+tac.x2[i]); else x=tac.x[i];
    if(x<startT) continue;
    if(x>endT) break;
    i2=i; if(i1<0) i1=i;
    n++;
  }
  if(n==0 || i1<0 || i2<i1) {tacFree(&tac); return TPCERROR_NO_DATA;}
  
  /* Allocate memory for the data to be extracted */
  int ret;
  ret=tacAllocate(d2, n, tac.tacNr); if(ret!=TPCERROR_OK) {tacFree(&tac); return(ret);}
  d2->sampleNr=n;
  d2->tacNr=tac.tacNr;
 
  /* Copy the contents */
  ret=tacCopyHdr(&tac, d2); if(ret!=TPCERROR_OK) {tacFree(&tac); return(ret);}
  for(int j=0, ret=0; j<tac.tacNr && ret==0; j++)
    ret=tacCopyTacchdr(&tac.c[j], &d2->c[j]);
  if(ret!=TPCERROR_OK) {tacFree(&tac); return(ret);}
  d2->isframe=tac.isframe;
  d2->weighting=tac.weighting;
  for(i=i1, n=0; i<=i2 && n<d2->sampleNr; i++, n++) {
    d2->x[n]=tac.x[i]; d2->x1[n]=tac.x1[i]; d2->x2[n]=tac.x2[i];
    d2->w[n]=tac.w[i];
    for(int j=0; j<tac.tacNr; j++) d2->c[j].y[n]=tac.c[j].y[i];
  }
  tacFree(&tac);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacExtractSamples(
  TAC *d1,
  TAC *d2,
  int si,
  int ei
) {
  if(d1==NULL || d2==NULL) return TPCERROR_FAIL;
  if(d1->sampleNr<1) return TPCERROR_NO_DATA;
  if(si<0 || ei<si || si>d1->sampleNr-1) return TPCERROR_INVALID_VALUE;
  if(ei>d1->sampleNr-1) ei=d1->sampleNr-1;
 
  /* Remove any old contents in target struct */
  tacFree(d2);
 
  /* If the whole TAC is inside required range, then just copy it */
  if(si==0 && ei==d1->sampleNr-1) return(tacDuplicate(d1, d2));
  
  /* Count how many samples are inside the required range */
  int n=1+ei-si; if(n==0) return TPCERROR_NO_DATA;
  
  /* Allocate memory for the data to be extracted */
  int ret=tacAllocate(d2, n, d1->tacNr); if(ret!=TPCERROR_OK) return(ret);
  d2->sampleNr=n;
  d2->tacNr=d1->tacNr;
 
  /* Copy the contents */
  ret=tacCopyHdr(d1, d2); if(ret!=TPCERROR_OK) return(ret);
  for(int j=0, ret=0; j<d1->tacNr && ret==0; j++) ret=tacCopyTacchdr(&d1->c[j], &d2->c[j]);
  if(ret!=TPCERROR_OK) return(ret);
  d2->isframe=d1->isframe;
  d2->weighting=d1->weighting;
  for(int i=si, n=0; i<=ei && n<d2->sampleNr; i++, n++) {
    d2->x[n]=d1->x[i]; d2->x1[n]=d1->x1[i]; d2->x2[n]=d1->x2[i];
    d2->w[n]=d1->w[i];
    for(int j=0; j<d1->tacNr; j++) d2->c[j].y[n]=d1->c[j].y[i];
  }
 
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacCheckX1X2X(
  TAC *d
) {
  if(d==NULL || d->sampleNr<1) return(0);
  int i, n;
  double x1, x2, x, dx;
 
  for(i=n=0; i<d->sampleNr; i++) {
    x1=d->x1[i]; x2=d->x2[i]; dx=x2-x1; x=0.5*(x1+x2);
    /* if not similar, or both missing (NaN), that is considered as non-match */
    if(!doubleMatch(x, d->x[i], 0.15*dx)) {n++; break;}
    /* check that there is not much overlap with previous sample */
    if(i>0 && (x1+0.15*dx)<d->x2[i-1]) {n++; break;}
    /* check that there is not much gap on both sides */
    if(i>0 && i<d->sampleNr-1) {
      if(x1-0.15*dx>d->x2[i-1] && x2+0.15*dx<d->x1[i+1]) {n++; break;}
    }
  }
  if(n>0) return(0); else return(1);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSetX(
  TAC *d,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
 
  if(d==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  if(d->sampleNr==0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
  /* Check that some x values are there to begin with */ 
  if(!tacIsX(d)) { 
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_X);
    return TPCERROR_NO_X;
  }
 
  /* If x1 and x2 are set, then calculate x for each sample */
  if(d->isframe) {
    if(verbose>1) printf("calculating x[] based on x1[] and x2[].\n");
    int i, n, ret;
    for(i=n=0; i<d->sampleNr; i++) {
      d->x[i]=0.5*(d->x1[i]+d->x2[i]);
      if(isfinite(d->x[i])) n++;
    }
    if(n>0) ret=TPCERROR_OK; else ret=TPCERROR_NO_X; 
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
 
  /* We should start guessing x1 and x2 values based on their averages in x */
  if(verbose>1) printf("estimating x1[] and x2[] based on x[].\n");
 
  /* If reasonable x1 and x2 values already are set, then quit */
  if(tacCheckX1X2X(d)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
 
  /* If only one sample, then set x1=0 and x2=2*x, and quit */
  if(d->sampleNr==1) {
    if(verbose>2) printf("just one sample.\n");
    if(d->x[0]<=0.0) {d->x1[0]=2.0*d->x[0]; d->x2[0]=0.0;}
    else {d->x1[0]=0.0; d->x2[0]=2.0*d->x[0];}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
 
  /* If only two samples */
  if(d->sampleNr==2) {
    if(verbose>2) printf("just two samples.\n");
    if(d->x[0]>0.0 && 2.0*d->x[0]<0.95*d->x[1]) {
      d->x1[0]=0.0; d->x2[0]=2.0*d->x[0];
      d->x1[1]=d->x2[0]; d->x2[1]=d->x2[0]+2.0*(d->x[1]-d->x2[0]);
    } else {
      double f; f=d->x[1]-d->x[0];
      d->x1[0]=d->x[0]-0.5*f; d->x2[0]=d->x[0]+0.5*f;
      d->x1[1]=d->x2[0]; d->x2[1]=d->x2[0]+2.0*f;
    }
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return TPCERROR_OK;
  }
 
  int i, n, ret;
 
  /* Allocate playground as TAC struct t */
  TAC t; tacInit(&t);
  n=d->tacNr; d->tacNr=1; ret=tacDuplicate(d, &t); d->tacNr=n;
  if(ret==0) ret=tacAllocateMore(&t, 1);
  if(ret) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  double *x, *x1, *x2, *fa, *fb;
  x=t.x; x1=t.x1; x2=t.x2; fa=t.c[0].y; fb=t.c[1].y;
  n=t.sampleNr;
 
  /* Convert playground times to sec if possible */
  tacXUnitConvert(&t, UNIT_SEC, NULL);
 
  /* Set each x1 and x1 to x */
  for(i=0; i<n; i++) x1[i]=x2[i]=x[i];
 
  /* Calculate the min gap, excluding zeroes */
  double gap, v;
  gap=nan("");
  for(i=0; i<n-1; i++) {
    v=x[i+1]-x[i]; if(!isfinite(v) || v<0.001) continue;
    if(isnan(gap)) {gap=v; continue;}
    if(v<gap) gap=v;
  }
 
  if(!isfinite(gap)) {
    tacFree(&t);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_X);
    return(TPCERROR_NO_X);
  }
  if(verbose>2) printf("smallest gap := %g\n", gap);
  /* Set gap to smallest step value for frame lengths */
  if(t.tunit!=UNIT_SEC) {
    gap*=0.01;
  } else {
    /* Set to 0.5 sec, unless calculated gap is smaller */
    if(gap>0.5) gap=0.5;
  }
  if(verbose>2) printf("stepsize := %g\n", gap);
 
  /* Initiate frame lengths to 0 */
  for(i=0; i<n; i++) fa[i]=fb[i]=0.0;
  int fixed[n]; for(i=0; i<n; i++) fixed[i]=0;
 
  /* Iterate until gaps are as small as possible */
  int anr=0;
  do {
    anr=0;
    /* Set tentative frames to the last accepted frame lengths */
    for(i=0; i<n; i++) fa[i]=fb[i];
    /* Add fraction of gap to frames that have not yet been fixed */
    for(i=0; i<n; i++) if(!fixed[i]) fa[i]+=gap;
    /* Calculate x1 and x2 based on the tentative frames */
    for(i=0; i<n; i++) {x1[i]=x[i]-0.5*fa[i]; x2[i]=x[i]+0.5*fa[i];}
    /* If no overlap, then copy frame */
    i=0; 
    if(x2[i]<=x1[i+1]+0.001*gap && (x[i]<=0.0 || x1[i]>-0.001*gap)) {
      fb[i]=fa[i]; if(!fixed[i]) anr++;
    } else 
      fixed[i]=1;
    for(i=1; i<n-1; i++) {
      if(x2[i]<=x1[i+1]+0.001*gap && x1[i]>=x2[i-1]-0.001*gap) {
        fb[i]=fa[i]; if(!fixed[i]) anr++;} else fixed[i]=1;
    }
    if(x1[i]>=x2[i-1]-0.001*gap) {
      fb[i]=fa[i]; if(!fixed[i]) anr++;} else fixed[i]=1;
    if(verbose>5) {
      printf("frames:");
      for(i=0; i<n; i++) printf(" %g", fb[i]);
      printf("\n");
    }
  } while(anr>0);
 
  if(verbose>6) {
    printf("temp frames:");
    for(i=0; i<n; i++) printf(" %g", fa[i]);
    printf("\n");
  }
 
  /* Calculate final frames */
  for(i=0; i<n; i++) {x1[i]=x[i]-0.5*fb[i]; x2[i]=x[i]+0.5*fb[i];}
 
  /* Convert playground to original time units */
  tacXUnitConvert(&t, d->tunit, NULL);
 
  /* Copy the guessed x1 and x2 values from playground */
  t.isframe=1; tacXCopy(&t, d, 0, d->sampleNr-1);
 
  /* Delete playground TAC */
  tacFree(&t); 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return TPCERROR_OK;
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacGetSampleInterval(
  TAC *d,
  double ilimit,
  double *minfdur,
  double *maxfdur
) {
  if(minfdur!=NULL) *minfdur=nan("");
  if(maxfdur!=NULL) *maxfdur=nan("");
  /* Check the data */
  if(d==NULL || d->sampleNr<1) return(1);
  double mi=nan(""), ma=nan("");
  if(d->isframe) {
    for(int i=0; i<d->sampleNr; i++) {
      double f=d->x2[i]-d->x1[i];
      if(!isnan(ilimit) && !(f>=ilimit)) continue;
      if(isnan(mi) || f<mi) mi=f;
      if(isnan(ma) || f>ma) ma=f;
    }
  } else {
    for(int i=1; i<d->sampleNr; i++) {
      double f=d->x[i]-d->x[i-1];
      if(!isnan(ilimit) && !(f>=ilimit)) continue;
      if(isnan(mi) || f<mi) mi=f;
      if(isnan(ma) || f>ma) ma=f;
    }
    if(isnan(mi) || isnan(ma)) { // as the last resort use the time of first sample
      double f=d->x[0]; 
      if(f>0.0 && !isnan(ilimit) && f>=ilimit) {
        if(isnan(mi)) mi=f;
        if(isnan(ma)) ma=f;
      }
    }
  }
  if(minfdur!=NULL) *minfdur=mi;
  if(maxfdur!=NULL) *maxfdur=ma;
  
  if(minfdur!=NULL && isnan(mi)) return(2);
  if(maxfdur!=NULL && isnan(ma)) return(2);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacToBars(
  TAC *tac1,
  TAC *tac2
) {
  if(tac1==NULL || tac2==NULL || tac1->sampleNr<1 || tac1->tacNr<1) return(1);
  /* Make sure that frame start and end times are present */
  if(tac1->isframe==0) return(2);
  /* Remove any old contents in target struct */
  tacFree(tac2);
 
  /* Make a copy of the TAC data */
  if(tacDuplicate(tac1, tac2)!=0) return(3);
  /* Add room for more frames; we plot each frame as its own bar
     and therefore we may need 4x the original space */
  if(tacAllocateMoreSamples(tac2, 3*tac1->sampleNr)!=0) return(4);
 
  /* Make separate 'samples' for each frame start and end time */
  tac2->isframe=0;
  int fj=0;
  for(int fi=0; fi<tac1->sampleNr; fi++) {
    /* Do nothing if missing values */
    if(isnan(tac1->x1[fi]) || isnan(tac1->x2[fi])) continue;
    if(tac1->tacNr==1 && isnan(tac1->c[0].y[fi])) continue;
    /* Otherwise make bar of each frame */
    tac2->x[fj]=tac1->x1[fi];
    for(int ri=0; ri<tac2->tacNr; ri++) tac2->c[ri].y[fj]=0.0;
    fj++;
    tac2->x[fj]=tac1->x1[fi];
    for(int ri=0; ri<tac2->tacNr; ri++) tac2->c[ri].y[fj]=tac1->c[ri].y[fi];
    fj++;
    tac2->x[fj]=tac1->x2[fi];
    for(int ri=0; ri<tac2->tacNr; ri++) tac2->c[ri].y[fj]=tac1->c[ri].y[fi];
    fj++;
    tac2->x[fj]=tac1->x2[fi];
    for(int ri=0; ri<tac2->tacNr; ri++) tac2->c[ri].y[fj]=0.0;
    fj++;
  }
  tac2->sampleNr=fj;
  if(fj<1) {tacFree(tac2); return(5);}
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
#if(1)
int tacFramesToSteps(
  TAC *inp,
  TAC *out,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || out==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  /* Check the function input */
  if(inp->sampleNr<1 || inp->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check that input does not have any missing values */
  if(tacNaNs(inp)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  }
 
  /* Check that there are no gaps or overlap between frames */
  {
    //printf("  checking that there are no gaps or overlap\n"); fflush(stdout);
    int ret=tacIsXContiguous(inp);
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    if(ret!=TPCERROR_OK) return(ret);
  }
 
  /* First, duplicate the TAC data */
  //printf("  copying data\n"); fflush(stdout);
  if(tacDuplicate(inp, out)!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
 
  /* If frame start and end times are NOT present, then consider that work is done */
  if(inp->isframe==0) return(TPCERROR_OK);
 
  /* Add room for more frames */
  //printf("  allocating memory for output\n"); fflush(stdout);
  if(tacAllocateMoreSamples(out, 2+2*out->sampleNr)!=0) {
    tacFree(out);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
 
  /* Make separate 'samples' for each frame start and end time */
  out->isframe=0;
  out->sampleNr=0;
  /* Add zero sample at zero time if data starts later */
  if(inp->x1[0]>0.01*(inp->x2[0]-inp->x1[0])) {
    if(verbose>10) printf("  add zero sample\n");
    out->x[0]=0.0;
    for(int ri=0; ri<out->tacNr; ri++) out->c[ri].y[0]=0.0;
    out->sampleNr++;
  }
  /* Add zero sample at the start of the first frame */
  if(verbose>10) printf("  add first frame zero sample\n");
  out->x[out->sampleNr]=inp->x1[0];
  for(int ri=0; ri<out->tacNr; ri++) out->c[ri].y[out->sampleNr]=0.0;
  out->sampleNr++;
 
  /* Add frame start and end times */
  for(int fi=0; fi<inp->sampleNr; fi++) {
    if(verbose>10) printf("  frame %d: %g - %g\n", 1+fi, inp->x1[fi], inp->x2[fi]);
    /* add start sample, if any of y values differs from previous one */
    out->x[out->sampleNr]=out->x[out->sampleNr-1];
//    if(fi==0) out->x[out->sampleNr]=inp->x1[fi];
//    else out->x[out->sampleNr]=out->x[out->sampleNr-1];
    int dn=0;
    for(int ri=0; ri<out->tacNr; ri++) {
      out->c[ri].y[out->sampleNr]=inp->c[ri].y[fi];
      if(out->c[ri].y[out->sampleNr]!=out->c[ri].y[out->sampleNr-1]) dn++;
    }
    if(dn>0) out->sampleNr++;
    /* add end sample */
    out->x[out->sampleNr]=inp->x2[fi];
    for(int ri=0; ri<out->tacNr; ri++)
      out->c[ri].y[out->sampleNr]=inp->c[ri].y[fi];
    out->sampleNr++;
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
#else
int tacFramesToSteps(
  TAC *inp,
  TAC *out,
  TPCSTATUS *status
) {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) printf("%s()\n", __func__);
  if(inp==NULL || out==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
  /* Check the function input */
  if(inp->sampleNr<1 || inp->tacNr<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return TPCERROR_NO_DATA;
  }
 
  /* Check that input does not have any missing values */
  if(tacNaNs(inp)>0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_MISSING_DATA);
    return TPCERROR_MISSING_DATA;
  } 
 
  /* First, duplicate the TAC data */
  if(tacDuplicate(inp, out)!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return TPCERROR_FAIL;
  }
 
  /* Make sure that frame start and end times are present */
  if(inp->isframe==0) return(TPCERROR_OK); // if not, consider that work is done
 
  /* Add room for more frames */
  if(tacAllocateMoreSamples(out, 1+2*out->sampleNr)!=0) {
    tacFree(out);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return TPCERROR_OUT_OF_MEMORY;
  }
 
  /* Make separate 'samples' for each frame start and end time */
  out->isframe=0;
  int fi=0, fj=0;
  if(fabs(inp->x1[fi])>1.0E-20) { // add zero sample, if data does not already have it
    out->x[fj]=inp->x1[fi];
    for(int ri=0; ri<out->tacNr; ri++) out->c[ri].y[fj]=0.0;
    fj++;
  }
  for(fi=0; fi<inp->sampleNr; fi++) {
    if(fi>0 && inp->x1[fi]<inp->x2[fi-1]) { // overlap
      double x=0.5*(inp->x1[fi]+inp->x2[fi-1]);
      if(x<out->x[fj-2]) { // too large overlap
        tacFree(out);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OVERLAPPING_DATA);
        return TPCERROR_OVERLAPPING_DATA;
      }
      out->x[fj-1]=out->x[fj]=x;
    } else { // no overlap
      out->x[fj]=inp->x1[fi];
    }
    for(int ri=0; ri<out->tacNr; ri++) out->c[ri].y[fj]=inp->c[ri].y[fi];
    fj++;
    out->x[fj]=inp->x2[fi];
    for(int ri=0; ri<out->tacNr; ri++) out->c[ri].y[fj]=inp->c[ri].y[fi];
    fj++;
  }
  out->sampleNr=fj;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}
#endif
/*****************************************************************************/
 
/*****************************************************************************/
int tacIsXContiguous(
  TAC *d
) {
  if(d==NULL) return TPCERROR_FAIL;
  if(tacXNaNs(d)>0) return(TPCERROR_INVALID_X);
 
  /* If no frame start and end times, then just check the order */
  if(d->isframe==0) {
    if(tacVerifyTimeOrder(d, NULL)) return(TPCERROR_INVALID_XRANGE);
    return TPCERROR_OK;
  }
 
  /* Check each frame */
  for(int fi=0; fi<d->sampleNr-1; fi++) {
    double e=d->x2[fi]-d->x1[fi+1];
    if(e<0.0) return(TPCERROR_LARGE_GAP);
    if(e>0.0) return(TPCERROR_OVERLAPPING_DATA);
  }
 
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int tacSetXContiguous(
  TAC *d
) {
  if(d==NULL) return TPCERROR_FAIL;
  if(tacXNaNs(d)>0) return(TPCERROR_INVALID_X);
 
  /* If no frame start and end times, then just make the times increasing */
  if(d->isframe==0) return(tacSortByTime(d, NULL));
 
  /* If there are no gaps or overlap, then we're done */
  if(tacIsXContiguous(d)==TPCERROR_OK) return(TPCERROR_OK);
 
  /* Make a duplicate for editing */
  TAC tac; tacInit(&tac);
  int ret=tacDuplicate(d, &tac); if(ret!=TPCERROR_OK) return(ret);
 
  /* Check for small frame overlaps or gaps; fix small, and return error if large overlap */
  int fixed=0, large=0;
  for(int fi=0; fi<d->sampleNr-1; fi++) {
    double dif=d->x1[fi+1]-d->x2[fi]; if(dif==0.0) continue; // not gap or overlap
    double fdur1=d->x2[fi]-d->x1[fi];
    double fdur2=d->x2[fi+1]-d->x1[fi+1];
    if(dif>0) { // gap
      if(dif<0.2*fdur1 && dif<0.2*fdur2) { // small gap; increase the shorter frame duration
        if(fdur1<fdur2) tac.x2[fi]=d->x1[fi+1]; else tac.x1[fi+1]=d->x2[fi];
        fixed++;
      } else {
        large++; // large gap not fixed yet
      }
      continue;
    }
    /* Overlap */
    dif=-dif;
    if(dif>0.2*fdur1 || dif>0.2*fdur2) { // too large overlap
      tacFree(&tac); return(TPCERROR_OVERLAPPING_DATA);
    }
    /* Reduce the longer frame length */
    if(fdur1>fdur2) tac.x2[fi]=d->x1[fi+1]; else tac.x1[fi+1]=d->x2[fi];
    fixed++;
  }
 
  /* If we fixed something, then copy the fixed frame times */
  if(fixed>0) {
    for(int i=0; i<d->sampleNr; i++) {
      d->x1[i]=tac.x1[i]; d->x2[i]=tac.x2[i]; d->x[i]=0.5*(d->x1[i]+d->x2[i]);
    }
  }
 
  /* If there were no large gaps, we're done */
  if(large==0) {tacFree(&tac); return(TPCERROR_OK);}
 
  /* Fill large gaps with new frames */
  d->weighting=WEIGHTING_OFF;
  if(tacAllocateMoreSamples(d, large)!=0) {
    tacFree(&tac); return(TPCERROR_OUT_OF_MEMORY);
  }
  d->sampleNr=1;
  for(int fi=1; fi<tac.sampleNr; fi++) {
    double dif=tac.x1[fi]-tac.x2[fi-1];
    if(dif>0.0) {
      /* Add new frame times */
      d->x1[d->sampleNr]=tac.x2[fi-1]; d->x2[d->sampleNr]=tac.x1[fi];
      d->x[fi]=0.5*(d->x1[d->sampleNr]+d->x2[d->sampleNr]);
      /* Calculate y values for the new frame */
      double x1=0.5*(tac.x1[fi-1]+tac.x2[fi-1]);
      double x2=0.5*(tac.x1[fi]+tac.x2[fi]);
      for(int ri=0; ri<tac.tacNr; ri++) {
        double s=(tac.c[ri].y[fi]-tac.c[ri].y[fi-1])/(x2-x1);
        d->c[ri].y[d->sampleNr]=tac.c[ri].y[fi-1] + s*(d->x[fi]-x1);
      }
      d->sampleNr++;
    }
    /* copy the existing frame */
    d->x1[d->sampleNr]=tac.x1[fi]; d->x2[d->sampleNr]=tac.x2[fi];
    for(int ri=0; ri<tac.tacNr; ri++) d->c[ri].y[d->sampleNr]=tac.c[ri].y[fi];
    d->sampleNr++;
  }
 
  tacFree(&tac);
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/