EventPrimitivesCovarianceHelpers.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/

// EventPrimitivesHelpers.h, (c) ATLAS Detector software
 
#ifndef EVENTPRIMITIVES_EVENTPRIMITIVESCOVARIANCEHELPERS_H
#define EVENTPRIMITIVES_EVENTPRIMITIVESCOVARIANCEHELPERS_H
 
#include "EventPrimitives/EventPrimitives.h"
//
#include <Eigen/Cholesky>
//
#include <limits>
 
namespace Amg {

inline bool saneCovarianceElement(double ele) {
  constexpr double upper_covariance_cutoff = std::numeric_limits<float>::max();
  return !(std::isnan(ele) || std::isinf(ele) ||
           std::abs(ele) > upper_covariance_cutoff);
}

template <int N>
inline bool hasPositiveOrZeroDiagElems(const AmgSymMatrix(N) & mat) {
  constexpr int dim = N;
  for (int i = 0; i < dim; ++i) {
    if (mat(i, i) < 0.0 || !saneCovarianceElement(mat(i, i)))
      return false;
  }
  return true;
}
 
inline bool hasPositiveOrZeroDiagElems(const Amg::MatrixX& mat) {
  int dim = mat.rows();
  for (int i = 0; i < dim; ++i) {
    if (mat(i, i) < 0.0 || !saneCovarianceElement(mat(i, i)))
      return false;
  }
  return true;
}

template <int N>
inline bool hasPositiveDiagElems(const AmgSymMatrix(N) & mat) {
  constexpr double MIN_COV_EPSILON = std::numeric_limits<float>::min();
  constexpr int dim = N;
  for (int i = 0; i < dim; ++i) {
    if (mat(i, i) < MIN_COV_EPSILON || !saneCovarianceElement(mat(i, i)))
      return false;
  }
  return true;
}
inline bool hasPositiveDiagElems(const Amg::MatrixX& mat) {
  constexpr double MIN_COV_EPSILON = std::numeric_limits<float>::min();
  int dim = mat.rows();
  for (int i = 0; i < dim; ++i) {
    if (mat(i, i) < MIN_COV_EPSILON || !saneCovarianceElement(mat(i, i)))
      return false;
  }
  return true;
}

template <int N>
inline bool isPositiveSemiDefinite(const AmgSymMatrix(N) & mat) {
  if (!hasPositiveOrZeroDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::LDLT<AmgSymMatrix(N)> ldltCov(mat);
  return (ldltCov.info() == Eigen::Success && ldltCov.isPositive());
}
inline bool isPositiveSemiDefinite(const Amg::MatrixX& mat) {
  if (!hasPositiveOrZeroDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::LDLT<Amg::MatrixX> ldltCov(mat);
  return (ldltCov.info() == Eigen::Success && ldltCov.isPositive());
}

template <int N>
inline bool isPositiveDefinite(const AmgSymMatrix(N) & mat) {
  if (!hasPositiveDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::LLT<AmgSymMatrix(N)> lltCov(mat);
  return (lltCov.info() == Eigen::Success);
}
inline bool isPositiveDefinite(const Amg::MatrixX& mat) {
  if (!hasPositiveDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::LLT<Amg::MatrixX> lltCov(mat);
  return (lltCov.info() == Eigen::Success);
}

template <int N>
inline bool isPositiveSemiDefiniteSlow(const AmgSymMatrix(N) & mat) {
  if (!hasPositiveOrZeroDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::SelfAdjointEigenSolver<AmgSymMatrix(5)> eigensolver(mat);
  auto res = eigensolver.eigenvalues();
  for (size_t i = 0; i < N; ++i) {
    if (res[i] < 0) {
      return false;
    }
  }
  return true;
}
inline bool isPositiveSemiDefiniteSlow(const Amg::MatrixX& mat) {
  if (!hasPositiveOrZeroDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::SelfAdjointEigenSolver<Amg::MatrixX> eigensolver(mat);
  auto res = eigensolver.eigenvalues();
  int dim = mat.rows();
  for (int i = 0; i < dim; ++i) {
    if (res[i] < 0) {
      return false;
    }
  }
  return true;
}
template <int N>
inline bool isPositiveDefiniteSlow(const AmgSymMatrix(N) & mat) {
  if (!hasPositiveDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::SelfAdjointEigenSolver<AmgSymMatrix(5)> eigensolver(mat);
  auto res = eigensolver.eigenvalues();
  for (size_t i = 0; i < N; ++i) {
    if (res[i] <= 0) {
      return false;
    }
  }
  return true;
}
inline bool isPositiveDefiniteSlow(const Amg::MatrixX& mat) {
  if (!hasPositiveDiagElems(mat)) {
    // fast check for necessary condition
    return false;
  }
  Eigen::SelfAdjointEigenSolver<Amg::MatrixX> eigensolver(mat);
  auto res = eigensolver.eigenvalues();
  int dim = mat.rows();
  for (int i = 0; i < dim; ++i) {
    if (res[i] <= 0) {
      return false;
    }
  }
  return true;
}
// Chi squared statistic of a linear regression as defined in Frühwirth, 
// section 3.2.1, equation 3.19. Precision is defined as the inverse of 
// of the covariance.
template <typename T, typename U>
inline double chi2(const T& precision, const U& residual, const int sign = 1) {
  return sign * residual.transpose() * precision * residual;
}
 
}  // namespace Amg
 
#endif