yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.7.2
MathUtils.h
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // This file is part of YODA -- Yet more Objects for Data Analysis
4 // Copyright (C) 2008-2017 The YODA collaboration (see AUTHORS for details)
5 //
6 #ifndef YODA_MathUtils_H
7 #define YODA_MathUtils_H
8 
10 
12 
13 #include <stdexcept>
14 #include <string>
15 #include <ostream>
16 #include <sstream>
17 #include <iostream>
18 #include <cmath>
19 #include <map>
20 #include <vector>
21 #include <utility>
22 #include <algorithm>
23 #include <numeric>
24 #include <cassert>
25 #include <limits>
26 #include <climits>
27 #include <cfloat>
28 
29 namespace YODA {
30 
31 
34  const static double MAXDOUBLE = DBL_MAX; // was std::numeric_limits<double>::max(); -- warns in GCC5
35  const static double MAXINT = INT_MAX; // was std::numeric_limits<int>::max(); -- warns in GCC5
36 
38  static const double PI = M_PI;
39 
41  static const double TWOPI = 2*M_PI;
42 
44  static const double HALFPI = M_PI_2;
45 
47  enum Sign { MINUS = -1, ZERO = 0, PLUS = 1 };
48 
49 
51 
52 
55  inline bool isZero(double val, double tolerance=1E-8) {
56  return (fabs(val) < tolerance);
57  }
58 
64  inline bool isZero(long val, double=1E-8) {
65  return val == 0;
66  }
67 
68 
72  inline bool fuzzyEquals(double a, double b, double tolerance=1E-5) {
73  const double absavg = (fabs(a) + fabs(b))/2.0;
74  const double absdiff = fabs(a - b);
75  const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg;
76  // cout << a << " == " << b << "? " << rtn << endl;
77  return rtn;
78  }
79 
80 
88  inline bool fuzzyEquals(long a, long b, double=1E-5) {
89  return a == b;
90  }
91 
92 
96  inline bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5) {
97  return a > b || fuzzyEquals(a, b, tolerance);
98  }
99 
107  inline bool fuzzyGtrEquals(long a, long b, double=1E-5) {
108  return a >= b;
109  }
110 
111 
115  inline bool fuzzyLessEquals(double a, double b, double tolerance=1E-5) {
116  return a < b || fuzzyEquals(a, b, tolerance);
117  }
118 
126  inline bool fuzzyLessEquals(long a, long b, double=1E-5) {
127  return a <= b;
128  }
129 
131  inline double approx(double a, int n = 5) {
132  double roundTo = pow(10.0,n);
133  a *= roundTo;
134  a = floor(a);
135  return a/roundTo;
136  }
138 
139 
141 
142 
147  enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
148 
149 
154  template<typename NUM>
155  inline bool inRange(NUM value, NUM low, NUM high,
156  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
157  if (lowbound == OPEN && highbound == OPEN) {
158  return (value > low && value < high);
159  } else if (lowbound == OPEN && highbound == CLOSED) {
160  return (value > low && value <= high);
161  } else if (lowbound == CLOSED && highbound == OPEN) {
162  return (value >= low && value < high);
163  } else { // if (lowbound == CLOSED && highbound == CLOSED) {
164  return (value >= low && value <= high);
165  }
166  }
167 
169  template<typename NUM>
170  inline bool inRange(NUM value, std::pair<NUM, NUM> lowhigh,
171  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
172  return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
173  }
174 
175 
180  inline bool inRange(int value, int low, int high,
181  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
182  if (lowbound == OPEN && highbound == OPEN) {
183  return (value > low && value < high);
184  } else if (lowbound == OPEN && highbound == CLOSED) {
185  return (value > low && value <= high);
186  } else if (lowbound == CLOSED && highbound == OPEN) {
187  return (value >= low && value < high);
188  } else { // if (lowbound == CLOSED && highbound == CLOSED) {
189  return (value >= low && value <= high);
190  }
191  }
192 
194  inline bool inRange(int value, std::pair<int, int> lowhigh,
195  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
196  return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
197  }
198 
200 
201 
203 
204 
206  template <typename NUM>
207  inline NUM sqr(NUM a) {
208  return a*a;
209  }
210 
212  template <typename Num>
213  inline Num add_quad(Num a, Num b) {
214  return sqrt(a*a + b*b);
215  }
216 
218  template <typename Num>
219  inline Num add_quad(Num a, Num b, Num c) {
220  return sqrt(a*a + b*b + c*c);
221  }
222 
224  inline int sign(double val) {
225  if (isZero(val)) return ZERO;
226  const int valsign = (val > 0) ? PLUS : MINUS;
227  return valsign;
228  }
229 
231  inline int sign(int val) {
232  if (val == 0) return ZERO;
233  return (val > 0) ? PLUS : MINUS;
234  }
235 
237  inline int sign(long val) {
238  if (val == 0) return ZERO;
239  return (val > 0) ? PLUS : MINUS;
240  }
241 
243 
244 
246 
247 
252  inline std::vector<double> linspace(size_t nbins, double start, double end, bool include_end=true) {
253  assert(end >= start);
254  assert(nbins > 0);
255  std::vector<double> rtn;
256  const double interval = (end-start)/static_cast<double>(nbins);
257  for (size_t i = 0; i < nbins; ++i) {
258  rtn.push_back(start + i*interval);
259  }
260  assert(rtn.size() == nbins);
261  if (include_end) rtn.push_back(end); // exact end, not result of n * interval
262  return rtn;
263  }
264 
265 
271  inline std::vector<double> logspace(size_t nbins, double start, double end, bool include_end=true) {
272  assert(end >= start);
273  assert(start > 0);
274  assert(nbins > 0);
275  const double logstart = std::log(start);
276  const double logend = std::log(end);
277  const std::vector<double> logvals = linspace(nbins, logstart, logend);
278  assert(logvals.size() == nbins+1);
279  std::vector<double> rtn; rtn.reserve(logvals.size());
280  rtn.push_back(start);
281  for (size_t i = 1; i < logvals.size()-1; ++i) {
282  rtn.push_back(std::exp(logvals[i]));
283  }
284  assert(rtn.size() == nbins);
285  if (include_end) rtn.push_back(end);
286  return rtn;
287  }
288 
289 
293  template <typename NUM>
294  inline int index_between(const NUM& val, const std::vector<NUM>& binedges) {
295  if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
296  int index = -1;
297  for (size_t i = 1; i < binedges.size(); ++i) {
298  if (val < binedges[i]) {
299  index = i-1;
300  break;
301  }
302  }
303  assert(inRange(index, -1, binedges.size()-1));
304  return index;
305  }
306 
308 
309 
311 
312 
314  inline double mean(const std::vector<int>& sample) {
315  double mean = 0.0;
316  for (size_t i=0; i<sample.size(); ++i) {
317  mean += sample[i];
318  }
319  return mean/sample.size();
320  }
321 
322 
324  inline double covariance(const std::vector<int>& sample1, const std::vector<int>& sample2) {
325  const double mean1 = mean(sample1);
326  const double mean2 = mean(sample2);
327  const size_t N = sample1.size();
328  double cov = 0.0;
329  for (size_t i = 0; i < N; i++) {
330  const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
331  cov += cov_i;
332  }
333  if (N > 1) return cov/(N-1);
334  else return 0.0;
335  }
336 
337 
339  inline double correlation(const std::vector<int>& sample1, const std::vector<int>& sample2) {
340  const double cov = covariance(sample1, sample2);
341  const double var1 = covariance(sample1, sample1);
342  const double var2 = covariance(sample2, sample2);
343  const double correlation = cov/sqrt(var1*var2);
344  const double corr_strength = correlation*sqrt(var2/var1);
345  return corr_strength;
346  }
347 
349 
350 
351 }
352 
353 #endif
bool inRange(NUM value, NUM low, NUM high, RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN)
Determine if value is in the range low to high, for floating point numbers.
Definition: MathUtils.h:155
bool isZero(double val, double tolerance=1E-8)
Definition: MathUtils.h:55
bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for >= with a degree of fuzziness.
Definition: MathUtils.h:96
int index_between(const NUM &val, const std::vector< NUM > &binedges)
Return the bin index of the given value, val, given a vector of bin edges.
Definition: MathUtils.h:294
std::vector< double > logspace(size_t nbins, double start, double end, bool include_end=true)
Make a list of nbins + 1 values exponentially spaced between start and end inclusive.
Definition: MathUtils.h:271
static const double PI
A pre-defined value of .
Definition: MathUtils.h:38
bool fuzzyEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for equality with a degree of fuzziness.
Definition: MathUtils.h:72
std::vector< double > linspace(size_t nbins, double start, double end, bool include_end=true)
Make a list of nbins + 1 values equally spaced between start and end inclusive.
Definition: MathUtils.h:252
static const double MAXINT
Definition: MathUtils.h:35
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:207
static const double HALFPI
A pre-defined value of .
Definition: MathUtils.h:44
static const double TWOPI
A pre-defined value of .
Definition: MathUtils.h:41
double covariance(const std::vector< int > &sample1, const std::vector< int > &sample2)
Calculate the covariance (variance) between two samples.
Definition: MathUtils.h:324
Num add_quad(Num a, Num b)
Named number-type addition in quadrature operation.
Definition: MathUtils.h:213
int sign(double val)
Find the sign of a number.
Definition: MathUtils.h:224
double approx(double a, int n=5)
Returns a number floored at the nth decimal place.
Definition: MathUtils.h:131
Sign
Enum for signs of numbers.
Definition: MathUtils.h:47
double mean(const std::vector< int > &sample)
Calculate the mean of a sample.
Definition: MathUtils.h:314
static const double MAXDOUBLE
Definition: MathUtils.h:34
RangeBoundary
Definition: MathUtils.h:147
double correlation(const std::vector< int > &sample1, const std::vector< int > &sample2)
Calculate the correlation strength between two samples.
Definition: MathUtils.h:339
bool fuzzyLessEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for <= with a degree of fuzziness.
Definition: MathUtils.h:115