source: Sophya/trunk/ArchTOIPipe/Processors/correl.cc@ 3689

#include <iostream.h>
#include <stdlib.h>
#include "correl.h"

CorrelEstimator::CorrelEstimator(int n, int wsize) {
  ncor = n;
  winsize = wsize;

  mT  = new int_4[winsize];
  mX  = new   r_8[winsize];
  mY  = new   r_8[winsize];
  
  mC  = new   r_8[ncor];
  mNC = new int_4[ncor];

  for (int i=0; i<ncor; i++) {
    mC[i] = mNC[i] = 0;
  }

  mNDat = 0;
  mINext = 0;
  mIFirst = 0;
}

CorrelEstimator::~CorrelEstimator() {
  delete[] mNC;
  delete[] mC;
  delete[] mY;
  delete[] mX;
  delete[] mT;
}

void CorrelEstimator::reset() {
  for (int i=0; i<ncor; i++) {
    mC[i] = mNC[i] = 0;
  }

  mNDat = 0;
  mINext = 0;
  mIFirst = 0;

}

void CorrelEstimator::push(int_4 t, r_8 x, r_8 y) {
  if (mNDat == winsize) {
    pop();
  }
  mT[mINext] = t;
  mX[mINext] = x;
  mY[mINext] = y;
  int iCur = mINext;
  mINext++; if (mINext == winsize) mINext=0;
  mNDat++;
  if (mINext > winsize) {cerr << "CorrelEstimator::push mINext > winsize"<< endl; abort();}

  
  for (int i=0; i<ncor; i++) {
    int ii=iCur;
    while (mT[ii] > t-i) {
      ii--;
      if (ii<0) ii = winsize-1;
      if ((mINext <= mIFirst  &&  (ii<mIFirst && ii>=mINext)) ||
          (mINext >  mIFirst  &&  (ii<mIFirst || ii>=mINext))) {  // not found
        ii = -1;
        break;
      }
    }
    if (ii>0) {
      mNC[i]++;
      mC[i] += mX[iCur] * mY[ii];
    }
  }
}

void CorrelEstimator::push(int_4 t, r_8 x) {
  push(t,x,x);
}

void CorrelEstimator::pop() {
  int t = mT[mIFirst];
  for (int i=0; i<ncor; i++) {
    int ii=mIFirst;
    while (mT[ii] < t+i) {
      ii++;
      if (ii>=winsize) ii = 0;
      if ((mINext <= mIFirst  &&  ii<mIFirst  &&  ii>=mINext) ||
          (mINext >  mIFirst  &&  ii>=mINext)) {  // not found
        ii = -1;
        break;
      }
    }
    if (ii>0) {
      mNC[i]--;
      mC[i] -= mX[ii] * mY[mIFirst];
    }
 
    mIFirst++;
    if (mIFirst >= winsize) {
      mIFirst = 0;
    }
  }
}

void CorrelEstimator::recompute() {
}

r_8 CorrelEstimator::correl(int k) {
  if (k<0 || k>= ncor) {
    cerr << "invalid parameter " << k << " for correl "
         << "(valid : [0;" << ncor-1<<"])" << endl;
    abort();
  }

  if (mNC[k] == 0) {
    return 0;
  } else {
    return mC[k]/mNC[k];
  }
}

int_4 CorrelEstimator::ns4correl(int k) {
  if (k<0 || k>= ncor) {
    cerr << "invalid parameter " << k << " for correl "
         << "(valid : [0;" << ncor-1<<"])" << endl;
    abort();
  }

  return mNC[k];
}