genie::utils::math Namespace Reference

Simple mathematical utilities not found in ROOT's TMath. More...

Functions
TMatrixD	CholeskyDecomposition (const TMatrixD &cov)
TVectorD	CholeskyGenerateCorrelatedParams (const TMatrixD &Lch, TVectorD &mean)
TVectorD	CholeskyGenerateCorrelatedParams (const TMatrixD &Lch, TVectorD &mean, TVectorD &g_uncorrelated)
TVectorD	CholeskyGenerateCorrelatedParamVariations (const TMatrixD &Lch)
TVectorD	CholeskyCalculateCorrelatedParamVariations (const TMatrixD &Lch, TVectorD &g_uncorrelated)
double	KahanSummation (double x[], unsigned int n)
double	KahanSummation (const vector< double > &x)
bool	AreEqual (double x1, double x2)
bool	AreEqual (float x1, float x2)
bool	IsWithinLimits (double x, Range1D_t range)
bool	IsWithinLimits (float x, Range1F_t range)
bool	IsWithinLimits (int i, Range1I_t range)
double	NonNegative (double x)
double	NonNegative (float x)

Detailed Description

Simple mathematical utilities not found in ROOT's TMath.

Author

Costas Andreopoulos <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

May 06, 2004

Copyright (c) 2003-2019, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Function Documentation

bool genie::utils::math::AreEqual	(	double	x1,
		double	x2
	)

Definition at line 243 of file MathUtils.cxx.

References dx, LOG, pINFO, and submit_syst::x2.

Referenced by genie::PathLengthList::AreAllZero(), genie::Spline::ClosestKnotValueIsZero(), genie::AxialFormFactor::Compare(), genie::ELFormFactors::Compare(), genie::QELFormFactors::Compare(), and genie::DISStructureFunc::Compare().

{
  double err = 0.001*DBL_EPSILON;
  double dx  = TMath::Abs(x1-x2);
  if(dx<err) {
    LOG("Math", pINFO) << x1 << " := " << x2;
    return true;
  }
  return false;;
}

bool genie::utils::math::AreEqual	(	float	x1,
		float	x2
	)

Definition at line 254 of file MathUtils.cxx.

References dx, LOG, pINFO, and submit_syst::x2.

{
  float err = FLT_EPSILON;
  float dx  = TMath::Abs(x1-x2);
  if(dx<err) {
    LOG("Math", pINFO) << x1 << " := " << x2;
    return true;
  }
  return false;;
}

TVectorD genie::utils::math::CholeskyCalculateCorrelatedParamVariations	(	const TMatrixD &	Lch,
		TVectorD &	g_uncorrelated
	)

Definition at line 193 of file MathUtils.cxx.

References ana::assert(), MECModelEnuComparisons::g, ncols, and fillBadChanDBTables::nrows.

{
  int ncols = cholesky_triangular.GetNcols();
  int nrows = cholesky_triangular.GetNrows();
  int npars = g_uncorrelated.GetNrows();

  assert(ncols==nrows);
  assert(npars==nrows);

  // create a vector of unit Gaussian variables
  // and multiply by Lt to introduce the appropriate correlations
  TVectorD g(g_uncorrelated);
  g *= cholesky_triangular;

  return g;
}

TMatrixD genie::utils::math::CholeskyDecomposition

(

const TMatrixD &

cov

)

Definition at line 27 of file MathUtils.cxx.

References ana::assert(), E, epsilon, exit(), stan::math::fabs(), MECModelEnuComparisons::i, calib::j, CLHEP::L, LOG, getGoodRuns4SAM::n, ncols, fillBadChanDBTables::nrows, pERROR, pINFO, ana::Sqrt(), and tmp.

{
// Perform a Cholesky decomposition of the input covariance matrix and
// return the lower triangular matrix
//
   const double epsilon = 1E-12;

   int ncols = cov_matrix.GetNcols();
   int nrows = cov_matrix.GetNrows();

   assert(ncols==nrows);

   int n = nrows;

   TMatrixD L(n, n);

   for (int i = 0; i < n; ++i) {

     // calculate the diagonal term first
     L(i,i) = cov_matrix(i,i);
     for (int k = 0; k < i; ++k) {
        double tmp = L(k,i);
        L(i,i) -= tmp*tmp;
     }//k

     if(L(i,i) <= 0) {
       if(fabs(L(i,i)) < epsilon){
         L(i,i)=epsilon;
         LOG("Cholesky", pINFO) 
           << "Changed element (" << i << ", " << i << ") to " << L(i,i);
       }
       else{
         LOG("Cholesky", pERROR) 
            << "Decomposed covariance matrix not positive-definite";
         LOG("Cholesky", pERROR) 
            << "L(" << i << "," << i << ") = " << L(i,i);
         exit(1);
       }
     }
     L(i,i) = TMath::Sqrt(L(i,i));
     // then the off-diagonal terms
     for (int j = i+1; j < n; ++j) {
        L(i,j) = cov_matrix(i,j);
        for (int k = 0; k < i; ++k) {
           L(i,j) -= L(k,i)*L(k,j);
        }
        L(i,j) /= L(i,i);
     }//j
  }//i

  // create the transpose of L
  TMatrixD LT(TMatrixD::kTransposed,L);

  return LT;
}

TVectorD genie::utils::math::CholeskyGenerateCorrelatedParams	(	const TMatrixD &	Lch,
		TVectorD &	mean
	)

Definition at line 83 of file MathUtils.cxx.

References exit(), MECModelEnuComparisons::g, MECModelEnuComparisons::i, LOG, getGoodRuns4SAM::n, ncols, fillBadChanDBTables::nrows, pERROR, and registry_explorer::v.

{
// Generate a vector of correlated params

  int ncols = cholesky_triangular.GetNcols();
  int nrows = cholesky_triangular.GetNrows();
  int npars = mean_params.GetNrows();

  if(ncols != nrows) {
    LOG("Cholesky", pERROR) 
        << "Mismatch between number of columns (" << ncols 
        << ") & rows (" << nrows << ")";
    exit(1);
  }
  if(npars != nrows) {
    LOG("Cholesky", pERROR) 
        << "Mismatch between number of parameters (" << npars
        << ") & array size (" << nrows << ")";
    exit(1);
  }

  int n = nrows;

  // create a vector of unit Gaussian variables
  // and multiply by Lt to introduce the appropriate correlations
  TVectorD g(n);
  for (int k = 0; k < n; ++k) {
    g(k) = RandomGen::Instance()->RndNum().Gaus();
  }
  g *= cholesky_triangular;

  // add the mean value offsets and store the results
  TVectorD correlated_params(n);
  for (int i = 0; i < n; ++i) {
     double v = mean_params[i];
     v += g(i);
     correlated_params[i] = v;
  }

  return correlated_params;
}

TVectorD genie::utils::math::CholeskyGenerateCorrelatedParams	(	const TMatrixD &	Lch,
		TVectorD &	mean,
		TVectorD &	g_uncorrelated
	)

Definition at line 126 of file MathUtils.cxx.

References exit(), MECModelEnuComparisons::i, LOG, getGoodRuns4SAM::n, ncols, fillBadChanDBTables::nrows, pERROR, and registry_explorer::v.

{
// Generate a vector of correlated params

  int ncols = cholesky_triangular.GetNcols();
  int nrows = cholesky_triangular.GetNrows();
  int npars = mean_params.GetNrows();
  int nunco = g_uncorrelated.GetNrows();

  if(ncols != nrows) {
    LOG("Cholesky", pERROR) 
        << "Mismatch between number of columns (" << ncols 
        << ") & rows (" << nrows << ")";
    exit(1);
  }
  if(npars != nrows) {
    LOG("Cholesky", pERROR) 
        << "Mismatch between number of parameters (" << npars
        << ") & array size (" << nrows << ")";
    exit(1);
  }
  if(nunco != nrows) {
    LOG("Cholesky", pERROR) 
        << "Mismatch between size of uncorrelated parameter vector (" << nunco
        << ") & array size (" << nrows << ")";
    exit(1);
  }

  int n = nrows;

  // create a vector of unit Gaussian variables
  // and multiply by Lt to introduce the appropriate correlations
  g_uncorrelated *= cholesky_triangular;

  // add the mean value offsets and store the results
  TVectorD correlated_params(n);
  for (int i = 0; i < n; ++i) {
     double v = mean_params[i];
     v += g_uncorrelated(i);
     correlated_params[i] = v;
  }

  return correlated_params;
}

TVectorD genie::utils::math::CholeskyGenerateCorrelatedParamVariations

(

const TMatrixD &

Lch

)

Definition at line 172 of file MathUtils.cxx.

References ana::assert(), MECModelEnuComparisons::g, getGoodRuns4SAM::n, ncols, and fillBadChanDBTables::nrows.

{
  int ncols = cholesky_triangular.GetNcols();
  int nrows = cholesky_triangular.GetNrows();

  assert(ncols==nrows);

  int n = nrows;

  // create a vector of unit Gaussian variables
  // and multiply by Lt to introduce the appropriate correlations
  TVectorD g(n);
  for (int k = 0; k < n; ++k) {
    g(k) = RandomGen::Instance()->RndNum().Gaus();
  }
  g *= cholesky_triangular;

  return g;
}

bool genie::utils::math::IsWithinLimits	(	double	x,
		Range1D_t	range
	)

Definition at line 265 of file MathUtils.cxx.

References genie::Range1D_t::max, and genie::Range1D_t::min.

Referenced by genie::utils::kinematics::ApplyCutsToKineLimits(), genie::DFRKinematicsGenerator::ComputeMaxXSec(), genie::KPhaseSpace::IsAllowed(), and genie::utils::kinematics::PhaseSpaceVolume().

{
  return ( x >= range.min && x <= range.max );
}

bool genie::utils::math::IsWithinLimits	(	float	x,
		Range1F_t	range
	)

Definition at line 270 of file MathUtils.cxx.

References genie::Range1F_t::max, and genie::Range1F_t::min.

{
  return ( x >= range.min && x <= range.max );
}

bool genie::utils::math::IsWithinLimits	(	int	i,
		Range1I_t	range
	)

Definition at line 275 of file MathUtils.cxx.

References genie::Range1I_t::max, and genie::Range1I_t::min.

{
  return ( i >= range.min && i <= range.max );
}

double genie::utils::math::KahanSummation	(	double	x[],
		unsigned int	n
	)

Definition at line 211 of file MathUtils.cxx.

References make_syst_table_plots::c, MECModelEnuComparisons::i, getGoodRuns4SAM::n, sum, confusionMatrixTree::t, and submit_syst::y.

{
// the Kahan summation algorithm - minimizes the error when adding a sequence
// of finite precision floating point numbers (compensated summation)

  double sum = x[0];
  double c   = 0.0;
  for(unsigned int i=1; i<n; i++) {
    double y = x[i]-c;
    double t = sum+y;
    c   = (t-sum) - y;
    sum = t;
  }
  return sum;
}

double genie::utils::math::KahanSummation

(

const vector< double > &

x

)

Definition at line 227 of file MathUtils.cxx.

References make_syst_table_plots::c, MECModelEnuComparisons::i, sum, confusionMatrixTree::t, and submit_syst::y.

{
// the Kahan summation algorithm - minimizes the error when adding a sequence
// of finite precision floating point numbers (compensated summation)

  double sum = x[0];
  double c   = 0.0;
  for(unsigned int i=1; i<x.size(); i++) {
    double y = x[i]-c;
    double t = sum+y;
    c   = (t-sum) - y;
    sum = t;
  }
  return sum;
}

double genie::utils::math::NonNegative

(

double

x

)

Definition at line 280 of file MathUtils.cxx.

Referenced by genie::PrimaryLeptonGenerator::AddToEventRecord(), genie::OutgoingDarkGenerator::AddToEventRecord(), and genie::NucBindEnergyAggregator::ProcessEventRecord().

{
// this is used to handle very small numbers in sqrts

  return TMath::Max(0., x);
}

double genie::utils::math::NonNegative

(

float

x

)

Definition at line 287 of file MathUtils.cxx.

{
// this is used to handle very small numbers in sqrts

  return TMath::Max( (float)0., x);
}