integrate.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 "tpcextensions.h"
/*****************************************************************************/
#include "tpcli.h"
/*****************************************************************************/
 
/*****************************************************************************/
int liIntegrate(
  double *x,
  double *y,
  const int nr,
  double *yi,
  const int se,
  const int verbose
) {
  if(verbose>0) printf("liIntegrate(x, y, %d, yi, %d)\n", nr, se);
  /* Check for data */
  if(nr<1) return 1;
  if(x==NULL || y==NULL || yi==NULL) return 1;
 
  /* Initiate the first integral to 0 */
  yi[0]=0.0;
  
  /* Extrapolate the initial phase */
  if(se==0) {
    /* No extrapolation */
  } else if(se==1) {
    /* Extrapolation with a rectangle */
    if(x[0]>0.0) yi[0]+=x[0]*y[0];
  } else if(se==2) {
    /* Extrapolation with a triangle starting from zero */
    if(x[0]>0.0 && y[0]>0.0) yi[0]+=0.5*x[0]*y[0];
  } else if(se==3) {
    /* Extrapolation with a triangle, if small gap */
    if(x[0]>0.0 && y[0]>0.0) {
      if(nr<2 || (2.0*x[0]-x[1])<1.0E-10)
        yi[0]+=0.5*x[0]*y[0];
    }
  } else if(se==4) {
    /* Extrapolation with a triangle */
    if(x[0]>0.0 && y[0]>0.0) {
      if(nr<2 || (2.0*x[0]-x[1])<1.0E-10)
        yi[0]+=0.5*x[0]*y[0];
      else {
        double nx;
        nx=2.0*x[0]-x[1];
        yi[0]+=0.5*(x[0]-nx)*y[0];        
      }
    }
  }
 
  /* Dot-to-dot for the rest of data */
  for(int j=1; j<nr; j++) yi[j]=yi[j-1]+0.5*(y[j]+y[j-1])*(x[j]-x[j-1]);
 
  if(verbose>3) printf("liIntegrate() ready\n");
  return 0;
}
/*****************************************************************************/
 
/*****************************************************************************/
int liIntegratePET(
  double *x1,
  double *x2,
  double *y,
  int nr,
  double *ie,
  double *iie,
  const int verbose
) {
  if(verbose>0) printf("liIntegratePET(x1, x2, y, %d, i, ii)\n", nr);
 
  /* Check for data */
  if(nr<1) return(1);
  if(x1==NULL || x2==NULL || y==NULL) return(1);
  /* No use to do anything if both output arrays are NULL */
  if(ie==NULL && iie==NULL) return(0);
 
  int i;
  double x, xs, fdur, half_fdur, gap, gap_y, midi, midii, fulli;
  double box_integral=0.0, box_integral2=0.0;
  double gap_integral=0.0, gap_integral2=0.0;
  double prev_x2=x1[0], prev_fdur=0.0;
 
  for(i=0; i<nr; i++) {
    /* Duration of this time frame */
    fdur=x2[i]-x1[i]; if(!(fdur>=0.0)) return(2);
    half_fdur=0.5*fdur;
    x=0.5*(x1[i]+x2[i]);
    if(verbose>3) {
      printf("  x=%g\n", x);
      printf("  fdur=%g\n", fdur);
    }
    /* Integral of a possible gap between frames */
    gap=x1[i]-prev_x2;
    if(i>0 && gap>1.0E-10) {
      /* Estimate y during gap as weighted mean of previous and this frame */
      gap_y=fdur*y[i]+prev_fdur*(y[i-1]);
      xs=fdur+prev_fdur; if(xs>0.0) gap_y/=xs;
      gap_integral=gap*gap_y;
      gap_integral2=gap*0.5*(box_integral+box_integral+gap_integral);
      if(verbose>3) {
        printf("  gap=%g\n", gap);
        printf("  gap_y=%g\n", gap_y);
        printf("  gap_integral=%g\n", gap_integral);
        printf("  gap_integral2=%g\n", gap_integral2);
      }
      /* Add gap integrals */
      box_integral+=gap_integral;
      box_integral2+=gap_integral2;
    }
    /* Calculate integrals at frame middle time */
    midi=y[i]*half_fdur;
    midii=0.5*(box_integral+box_integral+midi)*half_fdur;
    if(verbose>4) {
      printf("  midframe_integral=%g\n", midi);
      printf("  midframe_integral2=%g\n", midii);
    }
    if(ie!=NULL) ie[i]=box_integral+midi; 
    if(iie!=NULL) iie[i]=box_integral2+midii;
    /* Add full frame integrals */
    fulli=fdur*y[i];
    box_integral2+=0.5*(box_integral+box_integral+fulli)*fdur;
    box_integral+=fulli;
    if(verbose>3) {
      printf("  box_integral=%g\n", box_integral);
      printf("  box_integral2=%g\n", box_integral2);
    }
    /* Prepare for the next frame */
    prev_x2=x2[i]; prev_fdur=fdur;
  } // next time frame
 
  if(verbose>3) printf("liIntegratePET() ready\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int liIntegrateFE(
  double *x1,
  double *x2,
  double *y,
  int nr,
  double *ie,
  double *iie,
  const int verbose
) {
  if(verbose>0) {printf("liIntegrateFE(x1, x2, y, %d, i, ii, %d)\n", nr, verbose); fflush(stdout);}
 
  /* Check for data */
  if(nr<1 || x1==NULL || x2==NULL || y==NULL) return(TPCERROR_NO_DATA);
  /* No use to do anything if both output arrays are NULL */
  if(ie==NULL && iie==NULL) return(TPCERROR_OK);
 
  double prev_fdur=0.0;
  double prev_iy=0.0, prev_iiy=0.0;
  double iy=0.0, iiy=0.0;
  for(int i=0; i<nr; i++) {
    /* Duration of this time frame */
    double fdur=x2[i]-x1[i]; 
    if(!isfinite(fdur)) return(TPCERROR_INVALID_X);
    if(fabs(fdur)<1.0E-10) { // frame length is about zero
      if(ie!=NULL) ie[i]=prev_iy;
      if(iie!=NULL) iie[i]=prev_iiy;
      continue;
    }
    if(fdur<0.0) return(TPCERROR_INVALID_X);
    /* Check for gap */
    if(i>0) {
      double gap=x1[i]-x2[i-1];
      if(gap<0.0) return(TPCERROR_INVALID_X); // overlap not accepted
      if(gap>1.0E-10) {
        /* Estimate the integral of during the gap;
           estimate y during gap as weighted mean of previous and this frame */
        double gap_y=fdur*y[i]+prev_fdur*(y[i-1]);
        double xs=fdur+prev_fdur; if(xs>0.0) gap_y/=xs;
        double gap_integral=gap*gap_y;
        double gap_integral2=gap*0.5*(iy+iy+gap_integral);
        /* Add gap integrals */
        iy+=gap_integral;
        iiy+=gap_integral2;
      }
    }
    /* Calculate integrals at frame end */
    double box_integral=fdur*y[i];
    double box_integral2=fdur*0.5*(iy+iy+box_integral);
    /* Add frame integrals */
    iy+=box_integral;
    iiy+=box_integral2;
    if(ie!=NULL) ie[i]=iy; 
    if(iie!=NULL) iie[i]=iiy;
    /* Prepare for the next frame */
    prev_fdur=fdur; prev_iy=iy; prev_iiy=iiy;
  } // next time frame
 
  return(TPCERROR_OK);
}
/*****************************************************************************/
 
/*****************************************************************************/
int liIntegrateHalfFrame(
  double *x1,
  double *x2,
  double *y,
  const size_t nr,
  double *fhi,
  const int verbose
) {
  if(verbose>0) 
    printf("liIntegrateHalfFrame(x1, x2, y, %zu, fhi, %d)\n", nr, verbose);
 
  /* Check for data */
  if(nr<1) return 1;
  if(x1==NULL || x2==NULL || y==NULL || fhi==NULL) return 1;
 
  /* If just one frame */
  if(nr==1) {
    if(verbose>2) printf("just one frame.\n");
    fhi[0]=y[0]*(x2[0]-x1[0])/2.0;
    return 0;
  }
 
  /* Calculate frame middle times and durations / 2.
     Check that frames durations are >=0.
     Check that frame times are in ascending order.
  */
  double hfdur[nr], x[nr];
  for(size_t i=0; i<nr; i++) {
    x[i]=0.5*(x1[i]+x2[i]);
    hfdur[i]=0.5*(x2[i]-x1[i]); 
    if(hfdur[i]<0.0) {
      if(verbose>0) fprintf(stderr, "negative frame length.\n");
      return 2;
    }
    if(i>0 && x[i]<x[i-1]) {
      if(verbose>0) fprintf(stderr, "frame times not ascending.\n");
      return 3;
    }
  }
 
  /* If just two frames */
  if(nr<3) {
    if(verbose>2) printf("just two frames.\n");
    double k, dx;
    dx=x[1]-x[0]; 
    if(dx<1.0E-100) {fhi[0]=y[0]*hfdur[0]; fhi[1]=y[1]*hfdur[1]; return 0;}
    k=(y[1]-y[0])/dx;
    fhi[0]=hfdur[0]*(y[0]-k*0.5*hfdur[0]);
    fhi[1]=hfdur[1]*(y[1]-k*0.5*hfdur[1]);
    return 0; 
  }
 
  /* Compute data matrix for solving (very simply) the y values at frame middle
     times */
  double A[2][nr], B[nr], xdif;
  A[0][0]=1.0; A[1][0]=0.0; B[0]=y[0];
  for(size_t i=1; i<nr-1; i++) {
    xdif=x[i+1]-x[i];
    if(!(hfdur[i]>1.0E-100) || !(xdif>1.0E-100)) { // catches both zero and NaN
      A[0][i]=0.0; A[1][i]=0.0; B[i]=y[i];
    } else {
      A[0][i]=4.0*(x[i]-x[i-1])/hfdur[i] - (x[i+1]-x[i-1])/xdif;
      A[1][i]=(x[i]-x[i-1])/xdif;
      B[i]=4.0*y[i]*(x[i]-x[i-1])/hfdur[i];
    }
  }
  A[0][nr-1]=1.0; A[1][nr-1]=0.0; B[nr-1]=y[nr-1]; 
  if(verbose>4) {
    printf("A[1]\tA[2]\t\tB\n");
    for(size_t i=0; i<nr; i++) printf("%g\t%g\t\t%g\n", A[0][i], A[1][i], B[i]);
  }
 
  /* Solve */
  for(size_t i=1; i<nr-1; i++) {
    if(fabs(A[0][i])<1.0E-100) continue;
    A[0][i]-=A[1][i-1]; B[i]-=B[i-1];
    A[1][i]/=A[0][i]; B[i]/=A[0][i]; A[0][i]=1.0; 
  }
  if(verbose>5) {
    printf("A[1]\tA[2]\t\tB\n");
    for(size_t i=0; i<nr; i++) printf("%g\t%g\t\t%g\n", A[0][i], A[1][i], B[i]);
  }
  fhi[0]=y[0]; fhi[nr-1]=y[nr-1];
  for(size_t i=nr-2; i>0; i--) {
    fhi[i]=B[i]-fhi[i+1]*A[1][i];
    if(verbose>5) printf("X[%zu]=%g\n", i+1, fhi[i]);
  }
  /* Set any NaNs to original values (which of course may be NaN as well) */
  for(size_t i=1; i<nr-1; i++) if(isnan(fhi[i])) fhi[i]=y[i];
  if(verbose>10) {
    printf("estimated frame mid values:\n");
    for(size_t i=0; i<nr; i++) printf("  %g\n", fhi[i]);
  }
 
  /* Compute the first half integrals */
  double newi, lasty=fhi[0];
  double dx=x[1]-x[0];
  if(!(dx>1.0E-100)) fhi[0]*=hfdur[0];
  else fhi[0]=hfdur[0]*(fhi[0]-0.5*hfdur[0]*(fhi[1]-fhi[0])/dx);
  for(size_t i=1; i<nr; i++) {
    dx=x[i]-x[i-1];
    if(!(dx>1.0E-100)) newi=fhi[i]*hfdur[0];
    else newi=hfdur[i]*(fhi[i]-0.5*hfdur[i]*(fhi[i]-lasty)/dx);
    lasty=fhi[i]; fhi[i]=newi;    
  }
 
  if(verbose>1) printf("liIntegrateHalfFrame() ready\n");
  return 0;                     
}
/*****************************************************************************/
 
/*****************************************************************************/
int liDerivate(
  double *x,
  double *y,
  const int nr,
  double *d,
  double *dd,
  const int verbose
) {
  if(verbose>0) printf("liDerivate(x, y, %d, *d, *dd, %d)\n", nr, verbose);
  /* Check for data */
  if(nr<2) return(1);
  if(x==NULL || y==NULL || d==NULL) return(1);
 
  double dx;
 
  /* Estimate the first data point */
  d[0]=0.0; if(dd!=NULL) dd[0]=0.0;
  if(fabs(y[0])>1.0E-20) {
    dx=x[1]-x[0]; if(!(dx>0.0)) return(2);
    d[0]=y[0]/dx; if(!isfinite(d[0])) return(3);
    if(dd!=NULL) {dd[0]=d[0]/dx; if(!isfinite(dd[0])) return(4);}
  }
 
  /* Following data points */
  for(int i=1; i<nr; i++) {
    dx=x[i]-x[i-1]; if(!(dx>0.0)) return(2);
    d[i]=(y[i]-y[i-1])/dx; if(!isfinite(d[i])) return(3);
    if(dd!=NULL) {dd[i]=(d[i]-d[i-1])/dx; if(!isfinite(dd[i])) return(4);}
  }
 
  if(verbose>3) printf("liDerivate() ready\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int liDerivate3(
  double *x,
  double *y,
  const int nr,
  double *d,
  double *dd,
  const int verbose
) {
  if(verbose>0) printf("liDerivate3(x, y, %d, *d, *dd, %d)\n", nr, verbose);
  /* Check for data */
  if(nr<3) return(1);
  if(x==NULL || y==NULL || d==NULL) return(1);
 
  int i=0;
  /* Estimate the first data point */
  if(linRegr3p(x[i]+x[i]-x[i+1], 0.0, x[i], y[i], x[i+1], y[i+1], &d[i], NULL)) return(2);
  /* Following data points except the last */
  for(++i; i<nr-1; i++)
    if(linRegr3p(x[i-1], y[i-1], x[i], y[i], x[i+1], y[i+1], &d[i], NULL)) return(3);
  /* The last point */
  d[i]=(y[i]-y[i-1])/(x[i]-x[i-1]);
  if(!isfinite(d[i])) return(4);
 
  /* Second derivative */
  if(dd!=NULL) {
    return(liDerivate3(x, d, nr, dd, NULL, verbose)); 
  }
 
  if(verbose>3) printf("liDerivate3() ready\n");
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/