EigenUtils.cxx

Go to the documentation of this file.

#include "CAFAna/Core/EigenUtils.h"

#include "CAFAna/Core/Utilities.h" // For LogLikelihood(double, double)

#include <iostream>

namespace ana
{
  //----------------------------------------------------------------------
  Eigen::MatrixXd EigenMatrixXdFromTMatrixD(const TMatrixD* mat)
  {
    Eigen::MatrixXd ret(mat->GetNrows(), mat->GetNcols());
    // TMatrixD doesn't appear to have a GetArray()
    for(int i = 0; i < mat->GetNrows(); ++i){
      for(int j = 0; j < mat->GetNcols(); ++j){
        ret.coeffRef(i, j) = (*mat)(i, j);
      }
    }
    return ret;
  }

  //----------------------------------------------------------------------
  TMatrixD TMatrixDFromEigenMatrixXd(const Eigen::MatrixXd& mat)
  {
    TMatrixD ret(mat.rows(), mat.cols());
    // TMatrixD doesn't appear to have a GetArray()
    for(int i = 0; i < mat.rows(); ++i){
      for(int j = 0; j < mat.cols(); ++j){
        ret(i, j) = mat.coeffRef(i, j);
      }
    }
    return ret;
  }

  //----------------------------------------------------------------------
  double LogLikelihood(const Eigen::ArrayXd& ea, const Eigen::ArrayXd& oa, bool useOverflow)
  {
    assert(ea.size() == oa.size());

    double chi = 0;

    const int bufferBins = useOverflow ? 0 : -1;

    for(int i = 0; i < ea.size()+bufferBins; ++i){
      chi += LogLikelihood(ea[i], oa[i]);
    }

    return chi;
  }

  //----------------------------------------------------------------------
  double Chi2CovMx(const Eigen::ArrayXd& e, const Eigen::ArrayXd& o, const Eigen::MatrixXd& covmxinv)
  {
    assert(e.size() == covmxinv.rows()+2);

    Eigen::ArrayXd diff = e-o;
    // Drop underflow and overflow bins
    const Eigen::ArrayXd diffSub(Eigen::Map<Eigen::ArrayXd>(diff.data()+1, diff.size()-2));
    // dot collapses things down to a single number
    return diffSub.matrix().dot(covmxinv*diffSub.matrix());
  }

  //----------------------------------------------------------------------
  Eigen::MatrixXd SymmMxInverse(const Eigen::MatrixXd& mx)
  {
    // check if there are any null rows/columns.
    // if there are, they make the matrix singular.
    // we will remove them temporarily,
    // invert the matrix, then put them back afterwards.
    std::set<int> nullRows;
    for (int row = 0; row < mx.rows(); row++)
    {
      bool rowIsNull = true;
      for (int col = 0; col < mx.cols(); col++)
      {
        if (mx(row, col) != 0.)
        {
          rowIsNull = false;
          break;
        }
      }

      if (rowIsNull)
        nullRows.insert(row);
    }

    std::cerr << " Notice: covariance matrix has " << nullRows.size() << " null rows.\n"
        << "They will be removed before inverting and added back afterwards." << std::endl;

    // create a new matrix for inverting, skipping the null rows
    Eigen::MatrixXd invMx(mx.rows() - nullRows.size(),
                          mx.cols() - nullRows.size());
    unsigned int skippedRows = 0;
    for (int row = 0; row < mx.rows(); row++)
    {
      if (nullRows.find(row) != nullRows.end())
      {
        skippedRows++;
        continue;
      }
      unsigned int skippedCols = 0;
      for (int col = 0; col < mx.cols(); col++)
      {
        // since we assumed the matrix is symmetric,
        // we can just use the null rows list here
        if (nullRows.find(col) != nullRows.end())
        {
          skippedCols++;
          continue;
        }

        invMx(col-skippedCols, row-skippedRows) = invMx(row-skippedRows, col-skippedCols) = mx(row, col);
      }
    }

    invMx = invMx.inverse();

    // put back the empty rows if there were any
    if (!nullRows.empty())
    {
      skippedRows = 0;
      Eigen::MatrixXd retMx(mx.rows(), mx.cols());
      for (int row = 0; row < mx.rows(); row++)
      {
        if (nullRows.find(row) != nullRows.end())
        {
          skippedRows++;
          continue;
        }

        unsigned int skippedCols = skippedRows;
        for (int col = 0; col < mx.cols(); col++)
        {
          if (nullRows.find(col) != nullRows.end())
          {
            skippedCols++;
            continue;
          }

          retMx(col, row) = retMx(row, col) = invMx(row-skippedRows, col-skippedCols);
        }
      }

      return retMx;
    }

    return invMx;
  }
}