RandMultiGauss.cc

Go to the documentation of this file.

// $Id: RandMultiGauss.cc,v 1.3 2003/08/13 20:00:13 garren Exp $
// -*- C++ -*-
//
// -----------------------------------------------------------------------
//                             HEP Random
//                          --- RandMultiGauss ---
//                      class implementation file
// -----------------------------------------------------------------------
 
// =======================================================================
// Mark Fischler  - Created: 17th September 1998
// =======================================================================
 
// Some theory about how to get the Multivariate Gaussian from a bunch 
// of Gaussian deviates.  For the purpose of this discussion, take mu = 0.
//
// We want an n-vector x with distribution (See table 28.1 of Review of PP)
//
//            exp ( .5 * x.T() * S.inverse() * x ) 
//   f(x;S) = ------------------------------------
//                 sqrt ( (2*pi)^n * S.det() )
//
// Suppose S = U * D * U.T() with U orthogonal ( U*U.T() = 1 ) and D diagonal.
// Consider a random n-vector y such that each element y(i)is distributed as a
// Gaussian with sigma = sqrt(D(i,i)).  Then the distribution of y is the
// product of n Gaussains which can be written as 
//
//            exp ( .5 * y.T() * D.inverse() * y ) 
//   f(y;D) = ------------------------------------
//                 sqrt ( (2*pi)^n * D.det() )
//
// Now take an n-vector x = U * y (or y = U.T() * x ).  Then substituting,
//
//            exp ( .5 * x * U * D.inverse() U.T() * x ) 
// f(x;D,U) = ------------------------------------------
//                 sqrt ( (2*pi)^n * D.det() )
//
// and this simplifies to the desired f(x;S) when we note that 
//  D.det() = S.det()  and U * D.inverse() * U.T() = S.inverse()
//
// So the strategy is to diagonalize S (finding U and D), form y with each 
// element a Gaussian random with sigma of sqrt(D(i,i)), and form x = U*y.
// (It turns out that the D information needed is the sigmas.)
// Since for moderate or large n the work needed to diagonalize S can be much
// greater than the work to generate n Gaussian variates, we save the U and
// sigmas for both the default S and the latest S value provided.  
 
#include "CLHEP/RandomObjects/RandMultiGauss.h"
#include "CLHEP/RandomObjects/defs.h"
#include <cmath>    // for log()
#include <iostream>
 
namespace CLHEP {
 
// ------------
// Constructors
// ------------
 
RandMultiGauss::RandMultiGauss( HepRandomEngine & anEngine, 
                const HepVector & mu,
                const HepSymMatrix & S ) 
      : localEngine(&anEngine), 
    deleteEngine(false), 
    set(false),
    nextGaussian(0.0)
{ 
  if (S.num_row() != mu.num_row()) {
    std::cerr << "In constructor of RandMultiGauss distribution: \n" <<
            "      Dimension of mu (" << mu.num_row() << 
        ") does not match dimension of S (" << S.num_row() << ")\n";
    std::cerr << "---Exiting to System\n";
    exit(1);
  }
  defaultMu = mu;
  defaultSigmas = HepVector(S.num_row());
  prepareUsigmas (S,  defaultU, defaultSigmas);
}
  
RandMultiGauss::RandMultiGauss( HepRandomEngine * anEngine, 
                const HepVector & mu,
                const HepSymMatrix & S ) 
      : localEngine(anEngine), 
    deleteEngine(true), 
    set(false),
    nextGaussian(0.0)
{ 
  if (S.num_row() != mu.num_row()) {
    std::cerr << "In constructor of RandMultiGauss distribution: \n" <<
            "      Dimension of mu (" << mu.num_row() << 
        ") does not match dimension of S (" << S.num_row() << ")\n";
    std::cerr << "---Exiting to System\n";
    exit(1);
  }
  defaultMu = mu;
  defaultSigmas = HepVector(S.num_row());
  prepareUsigmas (S,  defaultU, defaultSigmas);
}
  
RandMultiGauss::RandMultiGauss( HepRandomEngine & anEngine ) 
      : localEngine(&anEngine), 
    deleteEngine(false), 
    set(false),
    nextGaussian(0.0)
{ 
  defaultMu = HepVector(2,0);
  defaultU  = HepMatrix(2,1);
  defaultSigmas = HepVector(2);
  defaultSigmas(1) = 1.;
  defaultSigmas(2) = 1.;
}
 
RandMultiGauss::RandMultiGauss( HepRandomEngine * anEngine ) 
      : localEngine(anEngine), 
    deleteEngine(true), 
    set(false),
    nextGaussian(0.0)
{ 
  defaultMu = HepVector(2,0);
  defaultU  = HepMatrix(2,1);
  defaultSigmas = HepVector(2);
  defaultSigmas(1) = 1.;
  defaultSigmas(2) = 1.;
}
 
RandMultiGauss::~RandMultiGauss() {
  if ( deleteEngine ) delete localEngine;
}
 
// ----------------------------
// prepareUsigmas()
// ----------------------------
 
void RandMultiGauss::prepareUsigmas( const HepSymMatrix & S,
                            HepMatrix & U, 
                            HepVector & sigmas ) { 
 
  HepSymMatrix tempS ( S ); // Since diagonalize does not take a const s
                // we have to copy S.
 
  U = diagonalize ( &tempS );           // S = U Sdiag U.T()
  HepSymMatrix D = S.similarityT(U);        // D = U.T() S U = Sdiag
  for (int i = 1; i <= S.num_row(); i++) {
    double s2 = D(i,i);
    if ( s2 > 0 ) {
    sigmas(i) = sqrt ( s2 );
    } else {
      std::cerr << "In RandMultiGauss distribution: \n" <<
            "      Matrix S is not positive definite.  Eigenvalues are:\n";
      for (int ixx = 1; ixx <= S.num_row(); ixx++) {
    std::cerr << "      " << D(ixx,ixx) << std::endl;
      }
      std::cerr << "---Exiting to System\n";
      exit(1);
    }
  }
} // prepareUsigmas
 
// -----------
// deviates()
// -----------
 
HepVector RandMultiGauss::deviates (    const HepMatrix & U,
                    const HepVector & sigmas,
                    HepRandomEngine * engine,
                                bool& available,
                    double& next)
{
  // Returns vector of gaussian randoms based on sigmas, rotated by U,
  // with means of 0.
 
  int n = sigmas.num_row(); 
  HepVector v(n);  // The vector to be returned
 
  double r,v1,v2,fac;
  
  int i = 1;
  if (available) {
    v(1) = next;
    i = 2;
    available = false;
  }
    
  while ( i <= n ) {
    do {
      v1 = 2.0 * engine->flat() - 1.0;
      v2 = 2.0 * engine->flat() - 1.0;
      r = v1*v1 + v2*v2;
    } while ( r > 1.0 );
    fac = sqrt(-2.0*log(r)/r);
    v(i++) = v1*fac;
    if ( i <= n ) {
      v(i++) = v2*fac;
    } else {
      next = v2*fac;
      available = true;
    }
  } 
 
  for ( i = 1; i <= n; i++ ) {
    v(i) *= sigmas(i);
  }
 
  return U*v;
 
} // deviates() 
 
// ---------------
// fire signatures
// ---------------
 
HepVector RandMultiGauss::fire() {
  // Returns a pair of unit normals, using the S and mu set in constructor,
  // utilizing the engine belonging to this instance of RandMultiGauss.
 
  return defaultMu + deviates ( defaultU, defaultSigmas, 
                    localEngine, set, nextGaussian );
 
} // fire();
 
 
HepVector RandMultiGauss::fire( const HepVector& mu, const HepSymMatrix& S ) {
 
  HepMatrix U;  
  HepVector sigmas(mu.num_row());
 
  if (mu.num_row() == S.num_row()) {
    prepareUsigmas ( S, U, sigmas );
    return mu + deviates ( U, sigmas, localEngine, set, nextGaussian );
  } else {
    std::cerr << "In firing RandMultiGauss distribution with explicit mu and S: \n"
         << "      Dimension of mu (" << mu.num_row() << 
        ") does not match dimension of S (" << S.num_row() << ")\n";
    std::cerr << "---Exiting to System\n";
    exit(1);
  }
  return mu;    // This line cannot be reached.  But without returning 
        // some HepVector here, KCC 3.3 complains.
 
} // fire(mu, S);
 
 
// --------------------
// fireArray signatures
// --------------------
 
void RandMultiGauss::fireArray( const int size, HepVector* array ) {
 
  int i;
  for (i = 0; i < size; ++i) {
    array[i] = defaultMu + deviates ( defaultU, defaultSigmas, 
                    localEngine, set, nextGaussian );
  } 
 
} // fireArray ( size, vect )
 
 
void RandMultiGauss::fireArray( const int size, HepVector* array,
                                 const HepVector& mu, const HepSymMatrix& S ) {
 
  // For efficiency, we diagonalize S once and generate all the vectors based
  // on that U and sigmas.
 
  HepMatrix U;  
  HepVector sigmas(mu.num_row());
  HepVector mu_ (mu);
 
  if (mu.num_row() == S.num_row()) {
    prepareUsigmas ( S, U, sigmas );
  } else {
    std::cerr << 
    "In fireArray for RandMultiGauss distribution with explicit mu and S: \n"
         << "      Dimension of mu (" << mu.num_row() << 
        ") does not match dimension of S (" << S.num_row() << ")\n";
    std::cerr << "---Exiting to System\n";
    exit(1);
  }
 
  int i;
  for (i=0; i<size; ++i) {
    array[i] = mu_ + deviates(U, sigmas, localEngine, set, nextGaussian);
  }
 
} // fireArray ( size, vect, mu, S )
 
// ----------
// operator()
// ----------
 
HepVector RandMultiGauss::operator()() {
  return fire();
} 
 
HepVector RandMultiGauss::operator()
            ( const HepVector& mu, const HepSymMatrix& S ) {
  return fire(mu,S);
} 
 
 
}  // namespace CLHEP