RandomMultiGauss.cc

Go to the documentation of this file.

//#include "CommonDet/DetUtilities/interface/DetExceptions.h"
#include "CommonTools/Statistics/interface/RandomMultiGauss.h"
#include "CLHEP/Random/RandGauss.h"
 
#include <cfloat>
//
// constructor with means and covariance
//
RandomMultiGauss::RandomMultiGauss(const AlgebraicVector& aVector, const AlgebraicSymMatrix& aMatrix)
    : theSize(aMatrix.num_row()), theMeans(aVector), theTriangle(theSize, theSize, 0) {
  //
  // Check consistency
  //
  if (theMeans.num_row() == theSize) {
    initialise(aMatrix);
  } else {
    //    throw DetLogicError("RandomMultiGauss: size of vector and matrix do not match");
    theMeans = AlgebraicVector(theSize, 0);
  }
}
//
// constructor with covariance (mean = 0)
//
RandomMultiGauss::RandomMultiGauss(const AlgebraicSymMatrix& aMatrix)
    : theSize(aMatrix.num_row()), theMeans(theSize, 0), theTriangle(theSize, theSize, 0) {
  //
  initialise(aMatrix);
}
//
// construct triangular matrix (Cholesky decomposition)
//
void RandomMultiGauss::initialise(const AlgebraicSymMatrix& aMatrix) {
  //
  // Cholesky decomposition with protection against empty rows/columns
  //
  for (int i1 = 0; i1 < theSize; i1++) {
    if (fabs(aMatrix[i1][i1]) < FLT_MIN)
      continue;
 
    for (int i2 = i1; i2 < theSize; i2++) {
      if (fabs(aMatrix[i2][i2]) < FLT_MIN)
        continue;
 
      double sum = aMatrix[i2][i1];
      for (int i3 = i1 - 1; i3 >= 0; i3--) {
        if (fabs(aMatrix[i3][i3]) < FLT_MIN)
          continue;
        sum -= theTriangle[i1][i3] * theTriangle[i2][i3];
      }
 
      if (i1 == i2) {
        //
        // check for positive definite input matrix, but allow for effects
        // due to finite precision
        //
        if (sum <= 0) {
          //      if ( sum<-FLT_MIN )  throw DetLogicError("RandomMultiGauss: input matrix is not positive definite");
          sum = FLT_MIN;
        }
        theTriangle[i1][i1] = sqrt(sum);
      } else {
        theTriangle[i1][i2] = 0.;
        theTriangle[i2][i1] = sum / theTriangle[i1][i1];
      }
    }
  }
}
//
// generate vector of random numbers
//
AlgebraicVector RandomMultiGauss::fire() {
  AlgebraicVector vRandom(theSize, 0);
  for (int i = 0; i < theSize; i++) {
    if (theTriangle[i][i] != 0)
      vRandom[i] = CLHEP::RandGauss::shoot();
  }
  return theTriangle * vRandom + theMeans;
}