NoiseGenerators/test/CorrelatedNoisifierTest.cpp

0001 #include "CLHEP/Random/JamesRandom.h"
0002 #include "CommonTools/Statistics/interface/AutocorrelationAnalyzer.h"
0003 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0004 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.icc"
0005
0006 namespace CLHEP {
0007   class HepRandomEngine;
0008 }
0009
0010 typedef math::ErrorD<10>::type MyMat;
0011
0012 template class CorrelatedNoisifier<MyMat>;
0013
0014 template void CorrelatedNoisifier<MyMat>::noisify(std::vector<double> &,
0015                                                   CLHEP::HepRandomEngine *,
0016                                                   const std::vector<double> *) const;
0017
0018 int main() {
0019   math::ErrorD<10>::type input;
0020   for (int k = 0; k < 10; k++) {
0021     for (int kk = k; kk < 10; kk++) {
0022       input(k, kk) = kk == k       ? 1
0023                      : kk == k + 1 ? 0.67
0024                      : kk == k + 2 ? 0.53
0025                      : kk == k + 3 ? 0.44
0026                      : kk == k + 4 ? 0.39
0027                      : kk == k + 5 ? 0.36
0028                      : kk == k + 6 ? 0.38
0029                      : kk == k + 7 ? 0.35
0030                      : kk == k + 8 ? 0.36
0031                      : kk == k + 9 ? 0.32
0032                                    : 0.;
0033     }
0034   }
0035   CLHEP::HepJamesRandom engine;
0036
0037   typedef math::ErrorD<10>::type MatType;
0038   typedef CorrelatedNoisifier<MatType> Noisifier;
0039
0040   MatType chol;
0041
0042   for (unsigned int itry(0); itry != 2; ++itry) {
0043     //      Noisifier noisifier ( input, &engine ) ;
0044     Noisifier noisifier(0 == itry ? Noisifier(input) : Noisifier(nullptr, chol));
0045
0046     if (0 == itry)
0047       chol = noisifier.cholMat();
0048
0049     AutocorrelationAnalyzer analyzer(10);
0050
0051     const unsigned int nTotal = 1000000;
0052     for (unsigned int i = 0; i < nTotal; ++i) {
0053       std::vector<double> samples(10);
0054       noisifier.noisify(samples, &engine);
0055       analyzer.analyze(samples);
0056     }
0057
0058     double big(-999);
0059     math::ErrorD<10>::type ratdif;
0060     for (int k = 0; k < 10; k++) {
0061       for (int kk = k; kk < 10; kk++) {
0062         const double diff(input(k, kk) - analyzer.correlation(k, kk));
0063         ratdif(k, kk) = diff / input(k, kk);
0064         if (fabs(ratdif(k, kk)) > big)
0065           big = fabs(ratdif(k, kk));
0066       }
0067     }
0068
0069     std::cout << "In " << nTotal << " trials, the biggest fractional deviation\n"
0070               << "of observed correlations from input correlations =" << big << std::endl;
0071
0072     std::cout << std::endl
0073               << "Initial correlations:"
0074               << "\n"
0075               << input << std::endl;
0076     ;
0077
0078     std::cout << analyzer << std::endl;
0079
0080     std::cout << ratdif << std::endl;
0081     std::cout << std::endl << "\nSQUARE of initial matrix:\n" << input * input << std::endl;
0082   }
0083 }