KfMatrixSym.h

Go to the documentation of this file.

/* Copyright (C) 2007-2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Sergey Gorbunov [committer] */
 

 
#ifndef KF_UTILS_KfMatrixSym_h
#define KF_UTILS_KfMatrixSym_h 1
 
#include "KfDefs.h"
#include "KfMath.h"
#include "KfUtils.h"
 
#include <boost/serialization/access.hpp>
 
#include <iomanip>
#include <sstream>
#include <string>
 
#if !defined(NO_ROOT) && !XPU_IS_HIP_CUDA
#include <Rtypes.h>  // for ClassDef
#endif
 
namespace cbm::algo::kf
{

  template<typename T, int N>
  class MatrixSym : public std::array<T, N*(N + 1) / 2> {
   public:
    static constexpr int kDimension{N};                    
    static constexpr int kNofElements{(N) * (N + 1) / 2};  

    MatrixSym() { this->fill(0.); }

    ~MatrixSym() = default;

    T operator()(int i, int j) const { return (*this)[math::SymIndex(i, j)]; }

    template<int i, int j>
    T operator()(Tag<i, j>) const
    {
      return (*this)[math::SymIndex(i, j)];
    }

    T& operator()(int i, int j) { return (*this)[math::SymIndex(i, j)]; }
 

    template<int i, int j>
    T& operator()(Tag<i, j>)
    {
      return (*this)[math::SymIndex(i, j)];
    }


    template<typename... Args>
    void Reset(Args... args)
    {
      static_assert(sizeof...(args) == kDimension,
                    "MatrixSym::Reset: number of arguments does not match matrix dimension ");
      this->fill(0.);
      SetDiagonal<0>(args...);
    }

    std::string ToString(int iv = -1) const;

    std::string ToStringCorrelations(int iv = -1) const;

    bool IsFinite() const { return IsFinite<DoPrintDebug::N>(nullptr); }

    bool IsSimdEntryConsistent(int iv) const { return IsSimdEntryConsistent<DoPrintDebug::N>(nullptr, iv); }

    bool IsConsistent(int nSimdFilled = -1) const { return IsConsistent<DoPrintDebug::N>(nullptr, nSimdFilled); }

    template<DoPrintDebug FlagPrintDebug = DoPrintDebug::Y>
    bool IsFinite(std::stringstream* ss) const;

    template<DoPrintDebug FlagPrintDebug = DoPrintDebug::Y>
    bool IsSimdEntryConsistent(std::stringstream* ss, int iv) const;

    template<DoPrintDebug FlagPrintDebug = DoPrintDebug::Y>
    bool IsConsistent(std::stringstream* ss, int nSimdFilled = -1) const;

    template<class Archive>
    void serialize(Archive& ar, const unsigned int /*version*/)
    {
      const std::array<T, kNofElements>& self = *this;
      ar& self;
    }
 
   private:
    template<int index, typename... Args>
    void SetDiagonal(T val, Args... args)
    {
      if constexpr (index < kDimension) {
        (*this)(Tag<index, index>{}) = val;
      }
      if constexpr (index < kDimension - 1) {
        SetDiagonal<index + 1>(args...);
      }
    }
 
#if !defined(NO_ROOT) && !XPU_IS_HIP_CUDA
    //   ClassDefNV(MatrixSym, 1);
#endif
  };  // class MatrixSym
 
 
  template<int N>
  using MatrixSymS = MatrixSym<fscal, N>;
 
  template<int N>
  using MatrixSymF = MatrixSym<float, N>;
 
  template<int N>
  using MatrixSymD = MatrixSym<double, N>;
 

 
  template<typename T, int N>
  inline std::string MatrixSym<T, N>::ToString(int iv) const
  {
    std::stringstream s;
    s.setf(std::ios::scientific, std::ios::floatfield);
    s << std::setprecision(6);
 
    if constexpr (std::is_same_v<T, fvec>) {
      // SIMD vector type
      if (iv >= 0) {
        // print only one component of the SIMD vector
        for (int i = 0; i < kDimension; i++) {
          for (int j = 0; j <= i; j++) {
            s << " " << operator()(i, j)[iv];
          }
          s << std::endl;
        }
      }
      else {
        // print all components of the SIMD vector
        for (unsigned int ivec = 0; ivec < fvec::Size; ++ivec) {
          s << "Vector component " << ivec << ":" << std::endl;
          s << ToString(ivec);
        }
      }
    }
    else {
      // scalar data type
      for (int i = 0; i < kDimension; i++) {
        for (int j = 0; j <= i; j++) {
          s << " " << operator()(i, j);
        }
        s << std::endl;
      }
    }
    return s.str();
  }
 
 
  template<typename T, int N>
  inline std::string MatrixSym<T, N>::ToStringCorrelations(int iv) const
  {
    std::stringstream s;
    s.setf(std::ios::scientific, std::ios::floatfield);
    s << std::setprecision(6);
 
    if constexpr (std::is_same_v<T, fvec>) {
      // SIMD vector type
      if (iv >= 0) {
        // print only one component of the SIMD vector
        for (int i = 0; i < kDimension; i++) {
          double sigi = sqrt(operator()(i, i)[iv]);
          for (int j = 0; j < i; j++) {
            double sigj = sqrt(operator()(j, j)[iv]);
            double cij  = operator()(i, j)[iv];
            double corr = cij / (sigi * sigj);
            s << " " << corr;
          }
          s << " " << sigi << std::endl;
        }
      }
      else {
        // print all components of the SIMD vector
        for (unsigned int ivec = 0; ivec < fvec::Size; ++ivec) {
          s << "Vector component " << ivec << ":" << std::endl;
          s << ToStringCorrelations(ivec);
        }
      }
    }
    else {
      // scalar data type
      for (int i = 0; i < kDimension; i++) {
        double sigi = sqrt(operator()(i, i));
        for (int j = 0; j < i; j++) {
          double sigj = sqrt(operator()(j, j));
          double cij  = operator()(i, j);
          double corr = cij / (sigi * sigj);
          s << " " << corr;
        }
        s << " " << sigi << std::endl;
      }
    }
    return s.str();
  }
 
 
  template<typename T, int N>
  template<DoPrintDebug FlagPrintDebug>
  inline bool MatrixSym<T, N>::IsFinite(std::stringstream* ss) const
  {
    // verify that all the numbers in the object are valid floats
    bool ret = true;
 
    for (int i = 0; i < kDimension; i++) {
      for (int j = 0; j <= i; j++) {
        T val = this->operator()(i, j);
        if (!utils::IsFinite(val)) {
          ret = false;
          if constexpr (FlagPrintDebug == DoPrintDebug::Y) {
            if (ss) {
              (*ss) << " Cov matrix element (" << i << "," << j << ")  is undefined: " << val << std::endl;
            }
          }
        }
      }
    }
    return ret;
  }
 
 
  template<typename T, int N>
  template<DoPrintDebug FlagPrintDebug>
  inline bool MatrixSym<T, N>::IsSimdEntryConsistent(std::stringstream* ss, int iv) const
  {
    // verify that all the numbers in the object are valid floats
 
    bool ok = IsFinite<FlagPrintDebug>(ss);
    // verify diagonal elements.
    // Cii is a squared dispersion of i-th parameter, it must be positive
 
    for (int i = 0; i < kDimension; i++) {
      double val = utils::simd::GetSimdEntry(this->operator()(i, i), iv);
      if (val <= 0.) {
        ok = false;
        if constexpr (FlagPrintDebug == DoPrintDebug::Y) {
          if (ss) {
            (*ss) << "  C[" << i << "," << i << "] is not positive: " << val << std::endl;
          }
        }
      }
    }
 
    // verify non-diagonal elements.
    // Cij/sqrt(Cii*Cjj) is a correlation between i-th and j-th parameter,
    // it must belong to [-1,1]
 
    for (int i = 1; i < kDimension; i++) {
      for (int j = 0; j < i; j++) {
        double tolerance = 1.0;
        double cij       = utils::simd::GetSimdEntry(this->operator()(i, j), iv);
        double cii       = utils::simd::GetSimdEntry(this->operator()(i, i), iv);
        double cjj       = utils::simd::GetSimdEntry(this->operator()(j, j), iv);
        if (cij * cij > tolerance * (cii * cjj)) {
          ok = false;
          if constexpr (FlagPrintDebug == DoPrintDebug::Y) {
            if (ss) {
              (*ss) << " Correlation [" << i << "," << j << "] is too large: " << cij / sqrt(cii * cjj) << std::endl;
            }
          }
        }
      }
    }
 
    // verify triplets of correlations
    // Kxy * Kxy + Kxz * Kxz + Kyz * Kyz <= 1 + 2 * Kxy * Kxz * Kyz
 
    for (int i = 2; i < kDimension; i++) {
      for (int j = 1; j < i; j++) {
        for (int m = 0; m < j; m++) {
          double tolerance = 1.0;
          double Cxx       = utils::simd::GetSimdEntry(this->operator()(i, i), iv);
          double Cyy       = utils::simd::GetSimdEntry(this->operator()(j, j), iv);
          double Czz       = utils::simd::GetSimdEntry(this->operator()(m, m), iv);
          double Cxy       = utils::simd::GetSimdEntry(this->operator()(i, j), iv);
          double Cxz       = utils::simd::GetSimdEntry(this->operator()(i, m), iv);
          double Cyz       = utils::simd::GetSimdEntry(this->operator()(j, m), iv);
          if (Cxx * Cyz * Cyz + Cyy * Cxz * Cxz + Czz * Cxy * Cxy
              > tolerance * (Cxx * Cyy * Czz + 2. * Cxy * Cyz * Cxz)) {
            ok = false;
            if constexpr (FlagPrintDebug == DoPrintDebug::Y) {
              double Kxy = Cxy / sqrt(Cxx * Cyy);
              double Kxz = Cxz / sqrt(Cxx * Czz);
              double Kyz = Cyz / sqrt(Cyy * Czz);
              if (ss) {
                (*ss) << " Correlations between parameters " << i << ", " << j << ", " << m << " are wrong: " << Kxy
                      << " " << Kxz << " " << Kyz << std::endl
                      << " inequation: " << Kxy * Kxy + Kxz * Kxz + Kyz * Kyz << " > " << 1 + 2 * Kxy * Kxz * Kyz
                      << std::endl;
              }
            }
          }
        }
      }
    }
 
    if constexpr (FlagPrintDebug == DoPrintDebug::Y) {
      if (!ok) {
        if (ss) {
          (*ss) << "Matrix elements are not consistent: " << std::endl << std::endl;
          (*ss) << ToString() << std::endl;
        }
      }
    }
    return ok;
  }  // IsSimdEntryConsistent
 
 
  template<typename T, int N>
  template<DoPrintDebug FlagPrintDebug>
  inline bool MatrixSym<T, N>::IsConsistent(std::stringstream* ss, int nSimdFilled) const
  {
    int size = 1;
 
    if constexpr (std::is_same_v<T, fvec>) {
      size = fvec::size();
    }
 
    assert(nSimdFilled <= size);
    if (nSimdFilled < 0) {
      nSimdFilled = size;
    }
 
    bool ok = true;
    for (int i = 0; i < nSimdFilled; ++i) {
      ok = ok && IsSimdEntryConsistent<FlagPrintDebug>(ss, i);
    }
 
    if constexpr (FlagPrintDebug == DoPrintDebug::Y) {
      if (!ok) {
        if (ss) {
          (*ss) << "MatrixSym parameters are not consistent: " << std::endl << std::endl;
          if (nSimdFilled == size) {
            (*ss) << "  All vector elements are filled " << std::endl << std::endl;
          }
          else {
            (*ss) << "  Only first " << nSimdFilled << " vector elements are filled " << std::endl << std::endl;
          }
          (*ss) << ToString(-1) << std::endl;
        }
      }
    }
    return ok;
  }
 
}  // namespace cbm::algo::kf
 
 
#endif