Belle II Software  release-05-02-19
DecorrelationMatrixHelper.h
1 /*************************************************************************
2 * BASF2 (Belle Analysis Framework 2) *
3 * Copyright(C) 2015 - Belle II Collaboration *
4 * *
5 * Author: The Belle II Collaboration *
6 * Contributors: Thomas Madlener *
7 * *
8 * This software is provided "as is" without any warranty. *
9 **************************************************************************/
10 
11 #pragma once
12 
13 #include <Eigen/Core>
14 #include <framework/logging/Logger.h>
15 
16 #include <array>
17 #include <vector>
18 #include <numeric>
19 
20 namespace Belle2 {
38  template<size_t Ndims>
39  const Eigen::Matrix<double, Ndims, Ndims, Eigen::RowMajor> calculateCovMatrix(std::array<std::vector<double>, Ndims> inputData)
40  {
41  // typedefs used in function
42  using RVector = Eigen::Matrix<double, 1, Ndims, Eigen::RowMajor>; // typedef for row vector type of one measurement
43  using DMatrix = Eigen::Matrix<double, Eigen::Dynamic, Ndims, Eigen::ColMajor>; // typedef for dataMatrix type
44  using DVector = Eigen::Matrix<double, Eigen::Dynamic, 1>; // typedef for column vector in dataMatrix
45 
46  // check if the inputData represent a "full" data matrix
47  size_t nSamples = inputData[0].size();
48  for (size_t i = 1; i < Ndims; ++i) {
49  if (inputData[i].size() != nSamples) {
50  B2ERROR("The input data is no MxN matrix, cannot calculate covariance matrix! Returning identity");
51  return Eigen::Matrix<double, Ndims, Ndims, Eigen::RowMajor>::Identity();
52  }
53  }
54 
55  // calculate the mean of every column and store them in a row vector
56  RVector meanVector{};
57  for (size_t i = 0; i < Ndims; ++i) {
58  meanVector(i) = std::accumulate(inputData[i].begin(), inputData[i].end(), 0.0) / inputData[i].size();
59  }
60 
61  // map the data to a more easily accessible Eigen::Matrix where each row is one measurement (consisting of Ndims values)
62  DMatrix dataMatrix(nSamples, Ndims);
63  for (size_t i = 0; i < Ndims; ++i) {
64  dataMatrix.col(i) = Eigen::Map<const DVector>(inputData[i].data(), inputData[i].size());
65  }
66 
67  // define a matrix where each row holds the mean vector to be able to substract the mean value from each column directly
68  DMatrix meanMatrix(nSamples, Ndims);
69  for (size_t i = 0; i < Ndims; ++i) {
70  meanMatrix.col(i) = DVector::Ones(nSamples) * meanVector(i);
71  }
72 
73  // calculate the cov matrix as product of dataMatrix reduced by meanValues in order to (hopefully) utilize Eigens optimizations
74  // COULDDO: test if this can be accelerated by splitting the process in smaller matrices
75  return (dataMatrix - meanMatrix).transpose() * (dataMatrix - meanMatrix) / (nSamples);
76  }
77 
79 } // end namespace
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::calculateCovMatrix
const Eigen::Matrix< double, Ndims, Ndims, Eigen::RowMajor > calculateCovMatrix(std::array< std::vector< double >, Ndims > inputData)
calculates the empirical covariance matrix from the inputData.
Definition: DecorrelationMatrixHelper.h:47