yoda is hosted by Hepforge, IPPP Durham
YODA - Yet more Objects for Data Analysis  1.9.0
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 <functional>
24 #include <numeric>
25 #include <cassert>
26 #include <limits>
27 #include <climits>
28 #include <cfloat>
29 
30 namespace YODA {
31 
32 
35  const static double MAXDOUBLE = DBL_MAX; // was std::numeric_limits<double>::max(); -- warns in GCC5
36  const static double MAXINT = INT_MAX; // was std::numeric_limits<int>::max(); -- warns in GCC5
37 
39  static const double PI = M_PI;
40 
42  static const double TWOPI = 2*M_PI;
43 
45  static const double HALFPI = M_PI_2;
46 
48  enum Sign { MINUS = -1, ZERO = 0, PLUS = 1 };
49 
50 
52 
53 
56  inline bool isZero(double val, double tolerance=1E-8) {
57  return (fabs(val) < tolerance);
58  }
59 
65  inline bool isZero(long val, double=1E-8) {
66  return val == 0;
67  }
68 
69 
73  inline bool fuzzyEquals(double a, double b, double tolerance=1E-5) {
74  const double absavg = (fabs(a) + fabs(b))/2.0;
75  const double absdiff = fabs(a - b);
76  const bool rtn = (isZero(a) && isZero(b)) || absdiff < tolerance*absavg;
77  // cout << a << " == " << b << "? " << rtn << endl;
78  return rtn;
79  }
80 
81 
89  inline bool fuzzyEquals(long a, long b, double=1E-5) {
90  return a == b;
91  }
92 
93 
97  inline bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5) {
98  return a > b || fuzzyEquals(a, b, tolerance);
99  }
100 
108  inline bool fuzzyGtrEquals(long a, long b, double=1E-5) {
109  return a >= b;
110  }
111 
112 
116  inline bool fuzzyLessEquals(double a, double b, double tolerance=1E-5) {
117  return a < b || fuzzyEquals(a, b, tolerance);
118  }
119 
127  inline bool fuzzyLessEquals(long a, long b, double=1E-5) {
128  return a <= b;
129  }
130 
132  inline double approx(double a, int n = 5) {
133  double roundTo = pow(10.0,n);
134  a *= roundTo;
135  a = floor(a);
136  return a/roundTo;
137  }
139 
140 
142 
143 
148  enum RangeBoundary { OPEN=0, SOFT=0, CLOSED=1, HARD=1 };
149 
150 
155  template<typename NUM>
156  inline bool inRange(NUM value, NUM low, NUM high,
157  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
158  if (lowbound == OPEN && highbound == OPEN) {
159  return (value > low && value < high);
160  } else if (lowbound == OPEN && highbound == CLOSED) {
161  return (value > low && value <= high);
162  } else if (lowbound == CLOSED && highbound == OPEN) {
163  return (value >= low && value < high);
164  } else { // if (lowbound == CLOSED && highbound == CLOSED) {
165  return (value >= low && value <= high);
166  }
167  }
168 
170  template<typename NUM>
171  inline bool inRange(NUM value, std::pair<NUM, NUM> lowhigh,
172  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
173  return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
174  }
175 
176 
181  inline bool inRange(int value, int low, int high,
182  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=CLOSED) {
183  if (lowbound == OPEN && highbound == OPEN) {
184  return (value > low && value < high);
185  } else if (lowbound == OPEN && highbound == CLOSED) {
186  return (value > low && value <= high);
187  } else if (lowbound == CLOSED && highbound == OPEN) {
188  return (value >= low && value < high);
189  } else { // if (lowbound == CLOSED && highbound == CLOSED) {
190  return (value >= low && value <= high);
191  }
192  }
193 
195  inline bool inRange(int value, std::pair<int, int> lowhigh,
196  RangeBoundary lowbound=CLOSED, RangeBoundary highbound=OPEN) {
197  return inRange(value, lowhigh.first, lowhigh.second, lowbound, highbound);
198  }
199 
201 
202 
204 
205 
207  template <typename NUM>
208  inline NUM sqr(NUM a) {
209  return a*a;
210  }
211 
213  template <typename Num>
214  inline Num add_quad(Num a, Num b) {
215  return sqrt(a*a + b*b);
216  }
217 
219  template <typename Num>
220  inline Num add_quad(Num a, Num b, Num c) {
221  return sqrt(a*a + b*b + c*c);
222  }
223 
225  inline int sign(double val) {
226  if (isZero(val)) return ZERO;
227  const int valsign = (val > 0) ? PLUS : MINUS;
228  return valsign;
229  }
230 
232  inline int sign(int val) {
233  if (val == 0) return ZERO;
234  return (val > 0) ? PLUS : MINUS;
235  }
236 
238  inline int sign(long val) {
239  if (val == 0) return ZERO;
240  return (val > 0) ? PLUS : MINUS;
241  }
242 
244 
245 
247 
248 
253  inline std::vector<double> linspace(size_t nbins, double xmin, double xmax, bool include_end=true) {
254  assert(xmax >= xmin);
255  assert(nbins > 0);
256  std::vector<double> rtn;
257  const double interval = (xmax-xmin)/static_cast<double>(nbins);
258  for (size_t i = 0; i < nbins; ++i) {
259  rtn.push_back(xmin + i*interval);
260  }
261  assert(rtn.size() == nbins);
262  if (include_end) rtn.push_back(xmax); // exact xmax, not result of n * interval
263  return rtn;
264  }
265 
266 
272  inline std::vector<double> logspace(size_t nbins, double xmin, double xmax, bool include_end=true) {
273  assert(xmax >= xmin);
274  assert(xmin > 0);
275  assert(nbins > 0);
276  const double logxmin = std::log(xmin);
277  const double logxmax = std::log(xmax);
278  const std::vector<double> logvals = linspace(nbins, logxmin, logxmax);
279  assert(logvals.size() == nbins+1);
280  std::vector<double> rtn; rtn.reserve(logvals.size());
281  rtn.push_back(xmin);
282  for (size_t i = 1; i < logvals.size()-1; ++i) {
283  rtn.push_back(std::exp(logvals[i]));
284  }
285  assert(rtn.size() == nbins);
286  if (include_end) rtn.push_back(xmax);
287  return rtn;
288  }
289 
290 
292  //inline std::vector<double> fspace(size_t nbins, double xmin, double xmax, std::function<double(double)>& fn) {
293 
294 
314  inline std::vector<double> pdfspace(size_t nbins, double xmin, double xmax, std::function<double(double)>& fn, size_t nsample=10000) {
315  const double dx = (xmax-xmin)/(double)nsample;
316  const std::vector<double> xs = linspace(nsample, xmin, xmax);
317  std::vector<double> ys(0, nsample);
318  auto posfn = [&](double x){return std::max(fn(x), 0.0);};
319  std::transform(xs.begin(), xs.end(), ys.begin(), posfn);
320  std::vector<double> areas; areas.reserve(nsample);
321  double areasum = 0;
322  for (size_t i = 0; i < ys.size()-1; ++i) {
323  const double area = (ys[i] + ys[i+1])*dx/2.0;
324  areas[i] = area;
325  areasum += area;
326  }
327  const double df = areasum/(double)nbins;
328  std::vector<double> xedges{xmin}; xedges.reserve(nbins+1);
329  double fsum = 0;
330  for (size_t i = 0; i < nsample-1; ++i) {
331  fsum += areas[i];
332  if (fsum > df) {
333  fsum = 0;
334  xedges.push_back(xs[i+1]);
335  }
336  }
337  xedges.push_back(xmax);
338  assert(xedges.size() == nbins+1);
339  return xedges;
340  }
341 
342 
346  template <typename NUM>
347  inline int index_between(const NUM& val, const std::vector<NUM>& binedges) {
348  if (!inRange(val, binedges.front(), binedges.back())) return -1; //< Out of histo range
349  int index = -1;
350  for (size_t i = 1; i < binedges.size(); ++i) {
351  if (val < binedges[i]) {
352  index = i-1;
353  break;
354  }
355  }
356  assert(inRange(index, -1, binedges.size()-1));
357  return index;
358  }
359 
361 
362 
364 
365 
367  inline double mean(const std::vector<int>& sample) {
368  double mean = 0.0;
369  for (size_t i=0; i<sample.size(); ++i) {
370  mean += sample[i];
371  }
372  return mean/sample.size();
373  }
374 
375 
377  inline double covariance(const std::vector<int>& sample1, const std::vector<int>& sample2) {
378  const double mean1 = mean(sample1);
379  const double mean2 = mean(sample2);
380  const size_t N = sample1.size();
381  double cov = 0.0;
382  for (size_t i = 0; i < N; i++) {
383  const double cov_i = (sample1[i] - mean1)*(sample2[i] - mean2);
384  cov += cov_i;
385  }
386  if (N > 1) return cov/(N-1);
387  else return 0.0;
388  }
389 
390 
392  inline double correlation(const std::vector<int>& sample1, const std::vector<int>& sample2) {
393  const double cov = covariance(sample1, sample2);
394  const double var1 = covariance(sample1, sample1);
395  const double var2 = covariance(sample2, sample2);
396  const double correlation = cov/sqrt(var1*var2);
397  const double corr_strength = correlation*sqrt(var2/var1);
398  return corr_strength;
399  }
400 
402 
403 
404 }
405 
406 #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:156
bool isZero(double val, double tolerance=1E-8)
Definition: MathUtils.h:56
bool fuzzyGtrEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for >= with a degree of fuzziness.
Definition: MathUtils.h:97
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:347
static const double PI
A pre-defined value of .
Definition: MathUtils.h:39
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:73
static const double MAXINT
Definition: MathUtils.h:36
NUM sqr(NUM a)
Named number-type squaring operation.
Definition: MathUtils.h:208
static const double HALFPI
A pre-defined value of .
Definition: MathUtils.h:45
static const double TWOPI
A pre-defined value of .
Definition: MathUtils.h:42
std::vector< double > logspace(size_t nbins, double xmin, double xmax, bool include_end=true)
Make a list of nbins + 1 values uniformly spaced in log(x) between xmin and xmax inclusive.
Definition: MathUtils.h:272
double covariance(const std::vector< int > &sample1, const std::vector< int > &sample2)
Calculate the covariance (variance) between two samples.
Definition: MathUtils.h:377
Num add_quad(Num a, Num b)
Named number-type addition in quadrature operation.
Definition: MathUtils.h:214
int sign(double val)
Find the sign of a number.
Definition: MathUtils.h:225
std::vector< double > pdfspace(size_t nbins, double xmin, double xmax, std::function< double(double)> &fn, size_t nsample=10000)
Make a list of nbins + 1 values spaced with density ~ f(x) between xmin and end inclusive.
Definition: MathUtils.h:314
double approx(double a, int n=5)
Returns a number floored at the nth decimal place.
Definition: MathUtils.h:132
Sign
Enum for signs of numbers.
Definition: MathUtils.h:48
double mean(const std::vector< int > &sample)
Calculate the mean of a sample.
Definition: MathUtils.h:367
static const double MAXDOUBLE
Definition: MathUtils.h:35
RangeBoundary
Definition: MathUtils.h:148
std::vector< double > linspace(size_t nbins, double xmin, double xmax, bool include_end=true)
Make a list of nbins + 1 values uniformly spaced between xmin and xmax inclusive. ...
Definition: MathUtils.h:253
double correlation(const std::vector< int > &sample1, const std::vector< int > &sample2)
Calculate the correlation strength between two samples.
Definition: MathUtils.h:392
bool fuzzyLessEquals(double a, double b, double tolerance=1E-5)
Compare two floating point numbers for <= with a degree of fuzziness.
Definition: MathUtils.h:116