NoiseGenerators/interface/CorrelatedNoisifier.icc

0001 #include "FWCore/Utilities/interface/Exception.h"
0002 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0003
0004 #include "CLHEP/Random/RandGaussQ.h"
0005
0006 template <class M>
0007 const double CorrelatedNoisifier<M>::k_precision(1.e-7);
0008
0009 template <class M>
0010 CorrelatedNoisifier<M>::~CorrelatedNoisifier() {}
0011
0012 template <class M>
0013 CorrelatedNoisifier<M>::CorrelatedNoisifier(const M &symCorMat) : m_H(ROOT::Math::SMatrixIdentity()) {
0014   init(symCorMat);
0015 }
0016
0017 template <class M>
0018 CorrelatedNoisifier<M>::CorrelatedNoisifier(int *dummy, const M &cholDecMat) : m_H(cholDecMat) {
0019   initChol();
0020 }
0021
0022 template <class M>
0023 void CorrelatedNoisifier<M>::initChol() {
0024   M symCorMat;
0025   for (unsigned int i = 0; i < m_H.kRows; ++i) {
0026     for (unsigned int j = 0; j < m_H.kRows; ++j) {
0027       double sum(0);
0028       for (unsigned int k(0); k != m_H.kRows; ++k) {
0029         const double hik(k > i ? 0 : m_H(i, k));
0030         const double htkj(k > j ? 0 : m_H(k, j));
0031         sum += hik * htkj;
0032       }
0033       symCorMat(i, j) = sum;
0034     }
0035   }
0036   checkOffDiagonal(symCorMat);
0037 }
0038
0039 template <class M>
0040 void CorrelatedNoisifier<M>::resetCorrelationMatrix(const M &symCorMat) {
0041   m_H = ROOT::Math::SMatrixIdentity();
0042   init(symCorMat);
0043 }
0044
0045 template <class M>
0046 void CorrelatedNoisifier<M>::resetCholDecompMatrix(const M &cholDecMat) {
0047   m_H = cholDecMat;
0048   initChol();
0049 }
0050
0051 template <class M>
0052 void CorrelatedNoisifier<M>::init(const M &symCorMat) {
0053   checkOffDiagonal(symCorMat);
0054
0055   const bool okDecomp(m_isIdentity || computeDecomposition(symCorMat));
0056
0057   if (!okDecomp)
0058     throw cms::Exception("CorrelatedNoisifier") << "Bad covariance matrix. " << symCorMat;
0059
0060   M HHtDiff;
0061   const bool check(checkDecomposition(symCorMat, HHtDiff));
0062   if (!check)
0063     throw cms::Exception("CorrelatedNoisifier") << "Decomposition failed, difference = " << HHtDiff;
0064 }
0065
0066 template <class M>
0067 void CorrelatedNoisifier<M>::checkOffDiagonal(const M &symCorMat) {
0068   m_isDiagonal = true;
0069   m_isIdentity = true;
0070
0071   for (unsigned int i = 0; i < m_H.kRows; i++) {
0072     if (symCorMat(i, i) < 0.)
0073       throw cms::Exception("CorrelatedNoisifier") << "Bad correlation matrix.  Negative diagonal";
0074     if (m_isIdentity && fabs(symCorMat(i, i) - 1.) > k_precision)
0075       m_isIdentity = false;
0076     for (unsigned int j = 0; j < i; ++j) {
0077       if (fabs(symCorMat(i, j)) > k_precision) {
0078         m_isDiagonal = false;
0079         m_isIdentity = false;
0080         return;
0081       }
0082     }
0083   }
0084 }
0085
0086 template <class M>
0087 bool CorrelatedNoisifier<M>::computeDecomposition(const M &symCorMat) {
0088   const double s00 = std::sqrt(symCorMat(0, 0));
0089
0090   for (unsigned int i = 0; i < m_H.kRows; ++i) {
0091     m_H(i, 0) = symCorMat(i, 0) / s00;
0092     if (0 < i) {
0093       if (!m_isDiagonal) {
0094         for (unsigned int j = 1; j < i; ++j) {
0095           double sum(symCorMat(i, j));
0096           for (unsigned int k = 0; k < j; ++k)
0097             sum -= m_H(i, k) * m_H(j, k);
0098           m_H(i, j) = sum / m_H(j, j);
0099         }
0100       }
0101
0102       double hii = symCorMat(i, i);
0103       for (unsigned int j = 0; j < i; ++j) {
0104         const double hij(m_H(i, j));
0105         hii -= hij * hij;
0106       }
0107
0108       if (hii <= k_precision)
0109         return false;  // protect sqrt in next line
0110
0111       m_H(i, i) = sqrt(hii);
0112     }
0113   }
0114   return true;
0115 }
0116
0117 template <class M>
0118 bool CorrelatedNoisifier<M>::checkDecomposition(const M &symCorMat, M &HHtDiff) const {
0119   bool equal(true);
0120   for (unsigned int i = 0; i < m_H.kRows; ++i) {
0121     for (unsigned int j = 0; j < m_H.kRows; ++j) {
0122       double sum(0);
0123       for (unsigned int k(0); k != m_H.kRows; ++k) {
0124         const double hik(k > i ? 0 : m_H(i, k));
0125         const double htkj(k > j ? 0 : m_H(k, j));
0126         sum += hik * htkj;
0127       }
0128       const double diff(symCorMat(i, j) - sum);
0129       HHtDiff(i, j) = diff;
0130       equal = equal && (fabs(diff) < k_precision);
0131     }
0132   }
0133   return equal;
0134 }
0135
0136 template <class M>
0137 template <class T>
0138 void CorrelatedNoisifier<M>::noisify(T &frame,
0139                                      CLHEP::HepRandomEngine *engine,
0140                                      const typename CorrelatedNoisifier<M>::VecDou *rangau) const {
0141   if ((unsigned int)frame.size() != (unsigned int)m_H.kRows)
0142     throw cms::Exception("Configuration") << "CorrelatedNoisifier::noisify(frame) was passed wrong size frame "
0143                                              "object";
0144
0145   if (nullptr != rangau) {
0146     assert(m_H.kRows == rangau->size());
0147     m_vecgau = *rangau;
0148   } else {
0149     if (m_H.kRows != m_vecgau.size())
0150       m_vecgau = VecDou(m_H.kRows, (double)0.0);
0151     CLHEP::RandGaussQ::shootArray(engine, m_H.kRows, &m_vecgau.front());
0152   }
0153
0154   for (unsigned int i(0); i < m_H.kRows; ++i) {
0155     frame[i] += (m_isIdentity ? m_vecgau[i] : m_H(i, i) * m_vecgau[i]);
0156     if (!m_isDiagonal) {
0157       for (unsigned int j = 0; j < i; ++j)
0158         frame[i] += m_H(j, i) * m_vecgau[j];
0159     }
0160   }
0161 }
0162
0163 template <class M>
0164 const M &CorrelatedNoisifier<M>::cholMat() const {
0165   return m_H;
0166 }
0167
0168 template <class M>
0169 const typename CorrelatedNoisifier<M>::VecDou &CorrelatedNoisifier<M>::vecgau() const {
0170   return m_vecgau;
0171 }