Functions
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().

244 {
245  double err = 0.001*DBL_EPSILON;
246  double dx = TMath::Abs(x1-x2);
247  if(dx<err) {
248  LOG("Math", pINFO) << x1 << " := " << x2;
249  return true;
250  }
251  return false;;
252 }
Float_t x1[n_points_granero]
Definition: compare.C:5
double dx[NP][NC]
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define pINFO
Definition: Messenger.h:63
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.

255 {
256  float err = FLT_EPSILON;
257  float dx = TMath::Abs(x1-x2);
258  if(dx<err) {
259  LOG("Math", pINFO) << x1 << " := " << x2;
260  return true;
261  }
262  return false;;
263 }
Float_t x1[n_points_granero]
Definition: compare.C:5
double dx[NP][NC]
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define pINFO
Definition: Messenger.h:63
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.

195 {
196  int ncols = cholesky_triangular.GetNcols();
197  int nrows = cholesky_triangular.GetNrows();
198  int npars = g_uncorrelated.GetNrows();
199 
200  assert(ncols==nrows);
201  assert(npars==nrows);
202 
203  // create a vector of unit Gaussian variables
204  // and multiply by Lt to introduce the appropriate correlations
205  TVectorD g(g_uncorrelated);
206  g *= cholesky_triangular;
207 
208  return g;
209 }
assert(nhit_max >=nhit_nbins)
Int_t ncols
Definition: plot.C:53
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.

28 {
29 // Perform a Cholesky decomposition of the input covariance matrix and
30 // return the lower triangular matrix
31 //
32  const double epsilon = 1E-12;
33 
34  int ncols = cov_matrix.GetNcols();
35  int nrows = cov_matrix.GetNrows();
36 
37  assert(ncols==nrows);
38 
39  int n = nrows;
40 
41  TMatrixD L(n, n);
42 
43  for (int i = 0; i < n; ++i) {
44 
45  // calculate the diagonal term first
46  L(i,i) = cov_matrix(i,i);
47  for (int k = 0; k < i; ++k) {
48  double tmp = L(k,i);
49  L(i,i) -= tmp*tmp;
50  }//k
51 
52  if(L(i,i) <= 0) {
53  if(fabs(L(i,i)) < epsilon){
54  L(i,i)=epsilon;
55  LOG("Cholesky", pINFO)
56  << "Changed element (" << i << ", " << i << ") to " << L(i,i);
57  }
58  else{
59  LOG("Cholesky", pERROR)
60  << "Decomposed covariance matrix not positive-definite";
61  LOG("Cholesky", pERROR)
62  << "L(" << i << "," << i << ") = " << L(i,i);
63  exit(1);
64  }
65  }
66  L(i,i) = TMath::Sqrt(L(i,i));
67  // then the off-diagonal terms
68  for (int j = i+1; j < n; ++j) {
69  L(i,j) = cov_matrix(i,j);
70  for (int k = 0; k < i; ++k) {
71  L(i,j) -= L(k,i)*L(k,j);
72  }
73  L(i,j) /= L(i,i);
74  }//j
75  }//i
76 
77  // create the transpose of L
78  TMatrixD LT(TMatrixD::kTransposed,L);
79 
80  return LT;
81 }
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
#define pERROR
Definition: Messenger.h:60
Float_t tmp
Definition: plot.C:36
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
static constexpr double L
Float_t E
Definition: plot.C:20
const double j
Definition: BetheBloch.cxx:29
#define pINFO
Definition: Messenger.h:63
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:324
double epsilon
exit(0)
assert(nhit_max >=nhit_nbins)
Int_t ncols
Definition: plot.C:53
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.

85 {
86 // Generate a vector of correlated params
87 
88  int ncols = cholesky_triangular.GetNcols();
89  int nrows = cholesky_triangular.GetNrows();
90  int npars = mean_params.GetNrows();
91 
92  if(ncols != nrows) {
93  LOG("Cholesky", pERROR)
94  << "Mismatch between number of columns (" << ncols
95  << ") & rows (" << nrows << ")";
96  exit(1);
97  }
98  if(npars != nrows) {
99  LOG("Cholesky", pERROR)
100  << "Mismatch between number of parameters (" << npars
101  << ") & array size (" << nrows << ")";
102  exit(1);
103  }
104 
105  int n = nrows;
106 
107  // create a vector of unit Gaussian variables
108  // and multiply by Lt to introduce the appropriate correlations
109  TVectorD g(n);
110  for (int k = 0; k < n; ++k) {
111  g(k) = RandomGen::Instance()->RndNum().Gaus();
112  }
113  g *= cholesky_triangular;
114 
115  // add the mean value offsets and store the results
116  TVectorD correlated_params(n);
117  for (int i = 0; i < n; ++i) {
118  double v = mean_params[i];
119  v += g(i);
120  correlated_params[i] = v;
121  }
122 
123  return correlated_params;
124 }
#define pERROR
Definition: Messenger.h:60
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
exit(0)
Int_t ncols
Definition: plot.C:53
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.

128 {
129 // Generate a vector of correlated params
130 
131  int ncols = cholesky_triangular.GetNcols();
132  int nrows = cholesky_triangular.GetNrows();
133  int npars = mean_params.GetNrows();
134  int nunco = g_uncorrelated.GetNrows();
135 
136  if(ncols != nrows) {
137  LOG("Cholesky", pERROR)
138  << "Mismatch between number of columns (" << ncols
139  << ") & rows (" << nrows << ")";
140  exit(1);
141  }
142  if(npars != nrows) {
143  LOG("Cholesky", pERROR)
144  << "Mismatch between number of parameters (" << npars
145  << ") & array size (" << nrows << ")";
146  exit(1);
147  }
148  if(nunco != nrows) {
149  LOG("Cholesky", pERROR)
150  << "Mismatch between size of uncorrelated parameter vector (" << nunco
151  << ") & array size (" << nrows << ")";
152  exit(1);
153  }
154 
155  int n = nrows;
156 
157  // create a vector of unit Gaussian variables
158  // and multiply by Lt to introduce the appropriate correlations
159  g_uncorrelated *= cholesky_triangular;
160 
161  // add the mean value offsets and store the results
162  TVectorD correlated_params(n);
163  for (int i = 0; i < n; ++i) {
164  double v = mean_params[i];
165  v += g_uncorrelated(i);
166  correlated_params[i] = v;
167  }
168 
169  return correlated_params;
170 }
#define pERROR
Definition: Messenger.h:60
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
exit(0)
Int_t ncols
Definition: plot.C:53
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.

174 {
175  int ncols = cholesky_triangular.GetNcols();
176  int nrows = cholesky_triangular.GetNrows();
177 
178  assert(ncols==nrows);
179 
180  int n = nrows;
181 
182  // create a vector of unit Gaussian variables
183  // and multiply by Lt to introduce the appropriate correlations
184  TVectorD g(n);
185  for (int k = 0; k < n; ++k) {
186  g(k) = RandomGen::Instance()->RndNum().Gaus();
187  }
188  g *= cholesky_triangular;
189 
190  return g;
191 }
assert(nhit_max >=nhit_nbins)
Int_t ncols
Definition: plot.C:53
bool genie::utils::math::IsWithinLimits ( double  x,
Range1D_t  range 
)
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.

271 {
272  return ( x >= range.min && x <= range.max );
273 }
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.

276 {
277  return ( i >= range.min && i <= range.max );
278 }
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.

212 {
213 // the Kahan summation algorithm - minimizes the error when adding a sequence
214 // of finite precision floating point numbers (compensated summation)
215 
216  double sum = x[0];
217  double c = 0.0;
218  for(unsigned int i=1; i<n; i++) {
219  double y = x[i]-c;
220  double t = sum+y;
221  c = (t-sum) - y;
222  sum = t;
223  }
224  return sum;
225 }
Double_t sum
Definition: plot.C:31
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.

228 {
229 // the Kahan summation algorithm - minimizes the error when adding a sequence
230 // of finite precision floating point numbers (compensated summation)
231 
232  double sum = x[0];
233  double c = 0.0;
234  for(unsigned int i=1; i<x.size(); i++) {
235  double y = x[i]-c;
236  double t = sum+y;
237  c = (t-sum) - y;
238  sum = t;
239  }
240  return sum;
241 }
Double_t sum
Definition: plot.C:31
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().

281 {
282 // this is used to handle very small numbers in sqrts
283 
284  return TMath::Max(0., x);
285 }
double genie::utils::math::NonNegative ( float  x)

Definition at line 287 of file MathUtils.cxx.

288 {
289 // this is used to handle very small numbers in sqrts
290 
291  return TMath::Max( (float)0., x);
292 }