GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MathUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <float.h>
12 
13 #include <TMath.h>
14 
18 
19 //____________________________________________________________________________
20 TMatrixD genie::utils::math::CholeskyDecomposition(const TMatrixD& cov_matrix)
21 {
22 // Perform a Cholesky decomposition of the input covariance matrix and
23 // return the lower triangular matrix
24 //
25  const double epsilon = 1E-12;
26 
27  int ncols = cov_matrix.GetNcols();
28  int nrows = cov_matrix.GetNrows();
29 
30  assert(ncols==nrows);
31 
32  int n = nrows;
33 
34  TMatrixD L(n, n);
35 
36  for (int i = 0; i < n; ++i) {
37 
38  // calculate the diagonal term first
39  L(i,i) = cov_matrix(i,i);
40  for (int k = 0; k < i; ++k) {
41  double tmp = L(k,i);
42  L(i,i) -= tmp*tmp;
43  }//k
44 
45  if(L(i,i) <= 0) {
46  if(fabs(L(i,i)) < epsilon){
47  L(i,i)=epsilon;
48  LOG("Cholesky", pINFO)
49  << "Changed element (" << i << ", " << i << ") to " << L(i,i);
50  }
51  else{
52  LOG("Cholesky", pERROR)
53  << "Decomposed covariance matrix not positive-definite";
54  LOG("Cholesky", pERROR)
55  << "L(" << i << "," << i << ") = " << L(i,i);
56  exit(1);
57  }
58  }
59  L(i,i) = TMath::Sqrt(L(i,i));
60  // then the off-diagonal terms
61  for (int j = i+1; j < n; ++j) {
62  L(i,j) = cov_matrix(i,j);
63  for (int k = 0; k < i; ++k) {
64  L(i,j) -= L(k,i)*L(k,j);
65  }
66  L(i,j) /= L(i,i);
67  }//j
68  }//i
69 
70  // create the transpose of L
71  TMatrixD LT(TMatrixD::kTransposed,L);
72 
73  return LT;
74 }
75 //____________________________________________________________________________
77  const TMatrixD& cholesky_triangular, TVectorD& mean_params)
78 {
79 // Generate a vector of correlated params
80 
81  int ncols = cholesky_triangular.GetNcols();
82  int nrows = cholesky_triangular.GetNrows();
83  int npars = mean_params.GetNrows();
84 
85  if(ncols != nrows) {
86  LOG("Cholesky", pERROR)
87  << "Mismatch between number of columns (" << ncols
88  << ") & rows (" << nrows << ")";
89  exit(1);
90  }
91  if(npars != nrows) {
92  LOG("Cholesky", pERROR)
93  << "Mismatch between number of parameters (" << npars
94  << ") & array size (" << nrows << ")";
95  exit(1);
96  }
97 
98  int n = nrows;
99 
100  // create a vector of unit Gaussian variables
101  // and multiply by Lt to introduce the appropriate correlations
102  TVectorD g(n);
103  for (int k = 0; k < n; ++k) {
104  g(k) = RandomGen::Instance()->RndNum().Gaus();
105  }
106  g *= cholesky_triangular;
107 
108  // add the mean value offsets and store the results
109  TVectorD correlated_params(n);
110  for (int i = 0; i < n; ++i) {
111  double v = mean_params[i];
112  v += g(i);
113  correlated_params[i] = v;
114  }
115 
116  return correlated_params;
117 }
118 //____________________________________________________________________________
120  const TMatrixD& cholesky_triangular, TVectorD& mean_params, TVectorD& g_uncorrelated)
121 {
122 // Generate a vector of correlated params
123 
124  int ncols = cholesky_triangular.GetNcols();
125  int nrows = cholesky_triangular.GetNrows();
126  int npars = mean_params.GetNrows();
127  int nunco = g_uncorrelated.GetNrows();
128 
129  if(ncols != nrows) {
130  LOG("Cholesky", pERROR)
131  << "Mismatch between number of columns (" << ncols
132  << ") & rows (" << nrows << ")";
133  exit(1);
134  }
135  if(npars != nrows) {
136  LOG("Cholesky", pERROR)
137  << "Mismatch between number of parameters (" << npars
138  << ") & array size (" << nrows << ")";
139  exit(1);
140  }
141  if(nunco != nrows) {
142  LOG("Cholesky", pERROR)
143  << "Mismatch between size of uncorrelated parameter vector (" << nunco
144  << ") & array size (" << nrows << ")";
145  exit(1);
146  }
147 
148  int n = nrows;
149 
150  // create a vector of unit Gaussian variables
151  // and multiply by Lt to introduce the appropriate correlations
152  g_uncorrelated *= cholesky_triangular;
153 
154  // add the mean value offsets and store the results
155  TVectorD correlated_params(n);
156  for (int i = 0; i < n; ++i) {
157  double v = mean_params[i];
158  v += g_uncorrelated(i);
159  correlated_params[i] = v;
160  }
161 
162  return correlated_params;
163 }
164 //____________________________________________________________________________
166  const TMatrixD& cholesky_triangular)
167 {
168  int ncols = cholesky_triangular.GetNcols();
169  int nrows = cholesky_triangular.GetNrows();
170 
171  assert(ncols==nrows);
172 
173  int n = nrows;
174 
175  // create a vector of unit Gaussian variables
176  // and multiply by Lt to introduce the appropriate correlations
177  TVectorD g(n);
178  for (int k = 0; k < n; ++k) {
179  g(k) = RandomGen::Instance()->RndNum().Gaus();
180  }
181  g *= cholesky_triangular;
182 
183  return g;
184 }
185 //____________________________________________________________________________
187  const TMatrixD& cholesky_triangular, TVectorD & g_uncorrelated)
188 {
189  int ncols = cholesky_triangular.GetNcols();
190  int nrows = cholesky_triangular.GetNrows();
191  int npars = g_uncorrelated.GetNrows();
192 
193  assert(ncols==nrows);
194  assert(npars==nrows);
195 
196  // create a vector of unit Gaussian variables
197  // and multiply by Lt to introduce the appropriate correlations
198  TVectorD g(g_uncorrelated);
199  g *= cholesky_triangular;
200 
201  return g;
202 }
203 //____________________________________________________________________________
204 double genie::utils::math::KahanSummation(double x[], unsigned int n)
205 {
206 // the Kahan summation algorithm - minimizes the error when adding a sequence
207 // of finite precision floating point numbers (compensated summation)
208 
209  double sum = x[0];
210  double c = 0.0;
211  for(unsigned int i=1; i<n; i++) {
212  double y = x[i]-c;
213  double t = sum+y;
214  c = (t-sum) - y;
215  sum = t;
216  }
217  return sum;
218 }
219 //____________________________________________________________________________
220 double genie::utils::math::KahanSummation(const vector<double> & x)
221 {
222 // the Kahan summation algorithm - minimizes the error when adding a sequence
223 // of finite precision floating point numbers (compensated summation)
224 
225  double sum = x[0];
226  double c = 0.0;
227  for(unsigned int i=1; i<x.size(); i++) {
228  double y = x[i]-c;
229  double t = sum+y;
230  c = (t-sum) - y;
231  sum = t;
232  }
233  return sum;
234 }
235 //____________________________________________________________________________
236 bool genie::utils::math::AreEqual(double x1, double x2)
237 {
238  double err = 0.001*DBL_EPSILON;
239  double dx = TMath::Abs(x1-x2);
240  if(dx<err) {
241  LOG("Math", pINFO) << x1 << " := " << x2;
242  return true;
243  }
244  return false;;
245 }
246 //____________________________________________________________________________
247 bool genie::utils::math::AreEqual(float x1, float x2)
248 {
249  float err = FLT_EPSILON;
250  float dx = TMath::Abs(x1-x2);
251  if(dx<err) {
252  LOG("Math", pINFO) << x1 << " := " << x2;
253  return true;
254  }
255  return false;;
256 }
257 //____________________________________________________________________________
259 {
260  return ( x >= range.min && x <= range.max );
261 }
262 //____________________________________________________________________________
264 {
265  return ( x >= range.min && x <= range.max );
266 }
267 //____________________________________________________________________________
269 {
270  return ( i >= range.min && i <= range.max );
271 }
272 //____________________________________________________________________________
274 {
275 // this is used to handle very small numbers in sqrts
276 
277  return TMath::Max(0., x);
278 }
279 //____________________________________________________________________________
281 {
282 // this is used to handle very small numbers in sqrts
283 
284  return TMath::Max( (float)0., x);
285 }
286 //____________________________________________________________________________
A simple [min,max] interval for integers.
Definition: Range1.h:56
static constexpr double g
Definition: Units.h:144
#define pERROR
Definition: Messenger.h:59
bool AreEqual(double x1, double x2)
Definition: MathUtils.cxx:236
A simple [min,max] interval for doubles.
Definition: Range1.h:42
TVectorD CholeskyGenerateCorrelatedParams(const TMatrixD &Lch, TVectorD &mean)
Definition: MathUtils.cxx:76
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:258
A simple [min,max] interval for floats.
Definition: Range1.h:28
TVectorD CholeskyCalculateCorrelatedParamVariations(const TMatrixD &Lch, TVectorD &g_uncorrelated)
Definition: MathUtils.cxx:186
TMatrixD CholeskyDecomposition(const TMatrixD &cov)
Definition: MathUtils.cxx:20
const double epsilon
double KahanSummation(double x[], unsigned int n)
Definition: MathUtils.cxx:204
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
#define pINFO
Definition: Messenger.h:62
double max
Definition: Range1.h:53
TVectorD CholeskyGenerateCorrelatedParamVariations(const TMatrixD &Lch)
Definition: MathUtils.cxx:165
double NonNegative(double x)
Definition: MathUtils.cxx:273
double min
Definition: Range1.h:52