LinearPredictor.cc

Go to the documentation of this file.

#include "LinearPredictor.h"

#include <cmath>
#include <vector>

#ifndef __python_test__
#include <utl/AugerException.h>
#endif

#include <Eigen/Dense> 
using namespace Eigen;

#ifndef __python_test__
using namespace utl;
#endif
using namespace std;

namespace RdChannelLinearPredictorRFISuppressor {
  LinearPredictor::LinearPredictor(unsigned int numCoeffPerChannel, unsigned int delayLine, unsigned int numChannels, unsigned int optimization)
  :
    numCoeffPerChannel(numCoeffPerChannel),
    numChannels(numChannels),
    delayLine(delayLine),
    //channels(numChannels), // unused. LN.
    dimCovMatr((numCoeffPerChannel + 1)*numChannels),
    dimPredictorCoefficients(numCoeffPerChannel*numChannels),
    totalPredictorOffset(numCoeffPerChannel + delayLine),
    optimization(optimization),
    covMatr(dimCovMatr, dimCovMatr),
    predictorCoefficients()
  {
    if (!numChannels) {
      #ifndef __python_test__
      throw InvalidConfigurationException("Need at least one channel to work on.");
      #else
      cout << "Need at least one channel to work on." << endl;
      exit(1);
      #endif
    };
    covMatr.setZero();
  }

  void LinearPredictor::train(std::vector< std::vector<double> > const &channels, double fudgeFactor) {
    if (channels.size() != numChannels) {
      #ifndef __python_test__
      throw InvalidConfigurationException("The number of provided (active) channels is not equal to the number of configured channels.");
      #else
      cout << "The number of provided (active) channels is not equal to the number of configured channels." << endl;
      exit(1);
      #endif
    };
    unsigned int channelLength = channels[0].size();
    unsigned int maxTraceStartIndex = channelLength - totalPredictorOffset;
    VectorXd samplingVector(dimCovMatr);

    // covMatr contains all relevant (co)variances.
    // The complete matrix looks like:
    // X ... X Y ... Y
    // . .   . Y .   Y
    // .  .  . Y  .  Y
    // .   . . Y   . Y
    // X ... X Y ... Y
    // Z ... Z A ... A
    // . .   . . .   .
    // .  .  . .  .  .
    // .   . . .   . .
    // Z ... Z A ... A
    // XX is a block-Toeplitz matrix
    // Where X is a square matrix with as dimension the total number of coefficients dimX
    #define dimX dimPredictorCoefficients
    // Where A is a square matrix with as dimension the total number of channels dimA
    #define dimA numChannels
    // X and Y are used but Z and A are not used in the rest of this calculation
    
    for (unsigned int traceStartIndex = 0; traceStartIndex < maxTraceStartIndex; ++traceStartIndex) {
      // fill the first part of the sampling vector before the delay line into the sampling vector
      unsigned int samplingIndex = 0;
      unsigned int traceIndex = 0;
      for (
        traceIndex = traceStartIndex; 
        traceIndex < traceStartIndex + numCoeffPerChannel; 
        ++traceIndex
      ) {
        for (unsigned int channelNum = 0; channelNum < numChannels; ++channelNum) {
          samplingVector(samplingIndex++) = channels[channelNum][traceIndex];
        }
      }
      // fill the last samples after the delay line into the sampling vector
      traceIndex += delayLine;
      for (unsigned int channelNum = 0; channelNum < numChannels; ++channelNum, ++samplingIndex) {
        samplingVector(samplingIndex) = channels[channelNum][traceIndex]; 
      }
      if (optimization == 0) {
        covMatr += samplingVector * samplingVector.transpose();
      } else {
        // Optimize exploiting the block-band-diagonal structure of the matrix
        #define S samplingVector
        covMatr.block(0, 0, dimA, dimX) += S.block(0, 0, dimA, 1) * S.block(0, 0, dimX, 1).transpose();
        covMatr.block(dimA, 0, dimX, dimA) += S.block(dimA, 0, dimX, 1) * S.block(0, 0, dimA, 1).transpose();
        covMatr.block(dimX, dimA, dimA, dimX) += S.block(dimX, 0, dimA, 1) * S.block(dimA, 0, dimX, 1).transpose();
        covMatr.block(0, dimX, dimX, dimA) += S.block(0, 0, dimX, 1) * S.block(dimX, 0, dimA, 1).transpose();
        #undef S
      }
    }
    if (optimization > 0) {
      for (unsigned int i = 1; i < numCoeffPerChannel; ++i) {
        for (unsigned int j = 1; j < numCoeffPerChannel; ++j) {
          if (i > j) {
            covMatr.block(i*dimA, j*dimA, dimA, dimA) = covMatr.block((i-j)*dimA, 0, dimA, dimA);
          } else {
            covMatr.block(i*dimA, j*dimA, dimA, dimA) = covMatr.block(0, (j-i)*dimA, dimA, dimA);          
          }
        }
      }
    }
    // We take the sub-matrix X from the covariance matrix
    MatrixXd subMatrixX = covMatr.block(0, 0, dimX, dimX);
    // We multiply the block diagonal with the fudge factor
    for (unsigned int i = 0; i < numCoeffPerChannel; ++i) { 
      for (unsigned int a = 0; a < numChannels; ++a) { 
        for (unsigned int b = 0; b < numChannels; ++b) {
          subMatrixX(i*numChannels+a, i*numChannels+b) *= (1.0+fudgeFactor);
        }
      }
    }
    // // This is the old rather inefficient way of finding the solution to the igenvalue problem
    // // We calculate the inverse of the sub-matrix X
    // MatrixXd inverseSubMatrixX(subMatrixX.inverse());
    // // We take the sub-matrix Y from the covariance matrix
    // MatrixXd subMatrixY = covMatr.block(0, dimX, dimX, dimA);
    // // The coefficients are the inner product of these two
    // predictorCoefficients = inverseSubMatrixX * subMatrixY;
    // The matrix is positive semi-definite so we choose LDLT.
    MatrixXd subMatrixY = covMatr.block(0, dimX, dimX, dimA);
    predictorCoefficients = subMatrixX.ldlt().solve(subMatrixY);
  }    

  std::vector< std::vector<double> > LinearPredictor::predict(std::vector< std::vector<double> > const &channels) const {
    if (channels.size() != numChannels) {
      #ifndef __python_test__
      throw InvalidConfigurationException("The number of provided (active) channels is not equal to the number of configured channels.");
      #else
      cout << "The number of provided (active) channels is not equal to the number of configured channels." << endl;
      exit(1);
      #endif
    };
    unsigned int channelLength = channels[0].size();
    unsigned int maxTraceStartIndex = channelLength - totalPredictorOffset;
    VectorXd samplingVector(dimCovMatr);

    std::vector< std::vector<double> > predicted(numChannels, vector<double>(channelLength, 0.0));

    for (unsigned int traceStartIndex = 0; traceStartIndex < maxTraceStartIndex; ++traceStartIndex) {
      unsigned int traceIndexPredict = traceStartIndex + totalPredictorOffset;
      for (unsigned int channelNum1 = 0; channelNum1 < numChannels; ++channelNum1) {
        unsigned int samplingIndex = 0;
        for (
          unsigned int traceIndex = traceStartIndex; 
          traceIndex < traceStartIndex + numCoeffPerChannel; 
          ++traceIndex
        ) {
          for (unsigned int channelNum2 = 0; channelNum2 < numChannels; ++channelNum2) {
            predicted[channelNum1][traceIndexPredict] += predictorCoefficients(samplingIndex++, channelNum1) * channels[channelNum2][traceIndex];
          }
        }
      }
    }
    return predicted;
  }
}