Belle II Software light-2406-ragdoll
MultivariateNormalGenerator.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <framework/utilities/MultivariateNormalGenerator.h>
10#include <framework/logging/Logger.h>
11
12using namespace Belle2;
13
15{
16 m_mean.resize(0);
17 m_transform.resize(0, 0);
18}
19
20bool MultivariateNormalGenerator::setMeanCov(int n, const double* mean, const double* cov)
21{
22 Eigen::VectorXd emean(n);
23 Eigen::MatrixXd ecov(n, n);
24 //Fill the vector and matrix from the buffers in row major layout
25 for (int i = 0; i < n; ++i) {
26 emean[i] = mean[i];
27 for (int j = 0; j < n; ++j) {
28 ecov(i, j) = cov[i * n + j];
29 }
30 }
31 //And delegate to the correct function
32 return setMeanCov(emean, ecov);
33}
34
35
36bool MultivariateNormalGenerator::setMeanCov(const Eigen::VectorXd& mean, const Eigen::MatrixXd& cov)
37{
38 reset();
39 if (mean.rows() != cov.rows()) {
40 B2ERROR("Mean values and covariance matrix need to be of the same dimension");
41 return false;
42 }
43 if (cov.rows() != cov.cols()) {
44 B2ERROR("Covariance matrix needs to be a square matrix");
45 return false;
46 }
47
48 // The usual way to calulate multivariate normal distributed values is to use
49 // the cholesky decomposition C = ll' and then calculate y = mean + l * x
50 // where x is a vector of standard normal distributed values. However this
51 // only works for positive definite matrices so we use the LDLT composition
52 // which gives us C = P'LDL'P and we compute L=P'Lsqrt(D) to get the same
53 // result but in a more robust way which works also for semi-definite
54 // matrices
55 auto ldlt = cov.ldlt();
56 if (ldlt.info() != Eigen::Success) {
57 B2ERROR("Cannot compute LDLT decomposition of covariance "
58 "matrix, maybe not positive semi-definite?");
59 return false;
60 }
61 Eigen::MatrixXd L = ldlt.matrixL();
62 Eigen::MatrixXd D = ldlt.vectorD().asDiagonal();
63 if (D.minCoeff() < 0) {
64 B2ERROR("MultivariateNormalGenerator: Negative values when computing LDL^T "
65 "decomposition, cannot compute M=AA^T, resulting random numbers "
66 "will not be correct");
67 return false;
68 }
69 auto P = ldlt.transpositionsP().transpose();
70 m_transform = P * L * D.cwiseSqrt();
71 //And save the mean values
72 m_mean = mean;
73 return true;
74}
75
Eigen::VectorXd m_mean
Member to store the mean values of the distribution.
bool setMeanCov(int n, const double *mean, const double *cov)
set the mean and covariance for the distribution with array interface: mean and covariance are passed...
void reset()
reset the generator setting the size to 0.
Eigen::MatrixXd m_transform
Member to store the transformation matrix for standard normal distributed random values.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24