MathUtils.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 8, 2016 - CA
14  Added Cholesky's method functions
15 */
16 //____________________________________________________________________________
17 
18 #include <float.h>
19 
20 #include <TMath.h>
21 
25 
26 //____________________________________________________________________________
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 }
82 //____________________________________________________________________________
84  const TMatrixD& cholesky_triangular, TVectorD& mean_params)
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 }
125 //____________________________________________________________________________
127  const TMatrixD& cholesky_triangular, TVectorD& mean_params, TVectorD& g_uncorrelated)
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 }
171 //____________________________________________________________________________
173  const TMatrixD& cholesky_triangular)
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 }
192 //____________________________________________________________________________
194  const TMatrixD& cholesky_triangular, TVectorD & g_uncorrelated)
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 }
210 //____________________________________________________________________________
211 double genie::utils::math::KahanSummation(double x[], unsigned int n)
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 }
226 //____________________________________________________________________________
227 double genie::utils::math::KahanSummation(const vector<double> & x)
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 }
242 //____________________________________________________________________________
243 bool genie::utils::math::AreEqual(double x1, double x2)
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 }
253 //____________________________________________________________________________
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 }
264 //____________________________________________________________________________
266 {
267  return ( x >= range.min && x <= range.max );
268 }
269 //____________________________________________________________________________
271 {
272  return ( x >= range.min && x <= range.max );
273 }
274 //____________________________________________________________________________
276 {
277  return ( i >= range.min && i <= range.max );
278 }
279 //____________________________________________________________________________
281 {
282 // this is used to handle very small numbers in sqrts
283 
284  return TMath::Max(0., x);
285 }
286 //____________________________________________________________________________
288 {
289 // this is used to handle very small numbers in sqrts
290 
291  return TMath::Max( (float)0., x);
292 }
293 //____________________________________________________________________________
294 
A simple [min,max] interval for integers.
Definition: Range1.h:57
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
#define pERROR
Definition: Messenger.h:60
bool AreEqual(double x1, double x2)
Definition: MathUtils.cxx:243
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Float_t x1[n_points_granero]
Definition: compare.C:5
TVectorD CholeskyGenerateCorrelatedParams(const TMatrixD &Lch, TVectorD &mean)
Definition: MathUtils.cxx:83
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:265
A simple [min,max] interval for floats.
Definition: Range1.h:29
Float_t tmp
Definition: plot.C:36
TVectorD CholeskyCalculateCorrelatedParamVariations(const TMatrixD &Lch, TVectorD &g_uncorrelated)
Definition: MathUtils.cxx:193
TMatrixD CholeskyDecomposition(const TMatrixD &cov)
Definition: MathUtils.cxx:27
double KahanSummation(double x[], unsigned int n)
Definition: MathUtils.cxx:211
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
static constexpr double L
Float_t E
Definition: plot.C:20
const double j
Definition: BetheBloch.cxx:29
#define pINFO
Definition: Messenger.h:63
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double epsilon
exit(0)
assert(nhit_max >=nhit_nbins)
TVectorD CholeskyGenerateCorrelatedParamVariations(const TMatrixD &Lch)
Definition: MathUtils.cxx:172
double NonNegative(double x)
Definition: MathUtils.cxx:280
double min
Definition: Range1.h:53
Double_t sum
Definition: plot.C:31
Int_t ncols
Definition: plot.C:53