Fleet  0.0.9
Inference in the LOT
Numerics.h
Go to the documentation of this file.
1 #pragma once
2 
3 // This includes our thread-safe lgamma and must be declared before math.h is imported
4 // This means that we should NOT import math.h anywhere else
5 #define _REENTRANT 1
6 
7 #include <random>
8 #include <iostream>
9 #include <assert.h>
10 #include <math.h>
11 
13 // Constants
15 
16 const double LOG2 = log(2.0); // clang doesn't like constexpr??
17 const double ROOT2 = sqrt(2.0);
19 constexpr double NaN = std::numeric_limits<double>::quiet_NaN();
20 constexpr double pi = M_PI;
21 constexpr double tau = 2*pi;
22 
24 // A fater logarithm
26 
27 //float __int_as_float(int32_t a) { float r; memcpy(&r, &a, sizeof(r)); return r;}
28 //int __float_as_int(float a) { int r; memcpy(&r, &a, sizeof(r)); return r;}
29 //
32 //float fastlogf(float a) {
33 // float m, r, s, t, i, f;
34 // int32_t e;
35 //
36 // e = (__float_as_int(a) - 0x3f2aaaab) & 0xff800000;
37 // m = __int_as_float (__float_as_int(a) - e);
38 // i = (float) e * 1.19209290e-7f; // 0x1.0p-23
39 // /* m in [2/3, 4/3] */
40 // f = m - 1.0f;
41 // s = f * f;
42 // /* Compute log1p(f) for f in [-1/3, 1/3] */
43 // r = fmaf(0.230836749f, f, -0.279208571f); // 0x1.d8c0f0p-3, -0x1.1de8dap-2
44 // t = fmaf(0.331826031f, f, -0.498910338f); // 0x1.53ca34p-2, -0x1.fee25ap-2
45 // r = fmaf(r, s, t);
46 // r = fmaf(r, s, f);
47 // r = fmaf(i, 0.693147182f, r); // 0x1.62e430p-1 // log(2)
48 // return r;
49 //}
50 
51 
53 // Probability
55 
56 
68 template<typename T>
69 T get_log1pexp_breakout_bound(const double precision, T f(T), const T lower=-1e6, const T upper=1e6) {
70 
71  assert(f(lower) > precision or f(upper) > precision);
72 
73  const T m = lower + (upper-lower)/2.0;
74  if (upper-lower < 1e-3) {
75  return m;
76  }
77  else if (f(m) < precision) { // bound on precision here
78  return get_log1pexp_breakout_bound<T>(precision, f, m, upper);
79  }
80  else {
81  return get_log1pexp_breakout_bound<T>(precision, f, lower, m);
82  }
83 }
84 
85 
86 
87 
88 template<typename T>
89 T logplusexp(const T a, const T b) {
90  // An approximate logplusexp -- the check on z<-25 makes it much faster since it saves many log,exp operations
91  // It is easy to derive a good polynomial approximation that is a bit faster (using sollya) on [-25,0] but that appears
92  // not to be worth it at this point.
93 
94  if (a == -infinity) return b;
95  else if(b == -infinity) return a;
96 
97  T mx = std::max(a,b);
98 
99  // compute a bound here at such that log1p(exp(z)) doesn't affect the result much (up to 1e-6)
100  // for floats, this ends up being about -13.8
101  // TODO: I Wish we could make this constexpr, but the math functions are not constexpr. Static is fine though
102  // it should only be run once at the start.
103  static const float breakout = get_log1pexp_breakout_bound<T>(1e-6, +[](T x) -> T { return log1p(exp(x)); });
104 
105  T z = std::min(a,b)-mx;
106  if(z < breakout) {
107  return mx; // save us from having to do anything
108  }
109  else {
110  return mx + log1p(exp(z)); // mx + log(exp(a-mx)+exp(b-mx));
111  }
112 }
113 
114 
121 template<typename t>
122 double logsumexp(const t& v) {
123  double lse = -infinity;
124  for(auto& x : v) {
125  lse = logplusexp(lse,x);
126  }
127  return lse;
128 }
129 
130 template<typename t>
131 double logsumexp(const std::vector<t>& v, double f(const t&) ) {
132  double lse = -infinity;
133  for(auto& x : v) {
134  lse = logplusexp(lse,f(x));
135  }
136  return lse;
137 }
138 
140 // Numerical functions
142 
143 double mylgamma(double v) {
144  // thread safe version gives our own place to store sgn
145  int sgn = 1;
146  auto ret = lgamma_r(v, &sgn);
147  return ret;
148 }
149 
150 double mygamma(double v) {
151  // thread safe usage
152  int sgn = 1;
153  auto ret = lgamma_r(v, &sgn);
154  return sgn*exp(ret);
155 }
156 
157 double lfactorial(double x) {
158  return mylgamma(x+1);
159 }
constexpr double pi
Definition: Numerics.h:20
const double ROOT2
Definition: Numerics.h:17
double mylgamma(double v)
Definition: Numerics.h:143
constexpr double infinity
Definition: Numerics.h:18
T logplusexp(const T a, const T b)
Definition: Numerics.h:89
double logsumexp(const t &v)
Compute log(sum(exp(v)). For now, this just unrolls logplusexp, but we might consider the faster (sta...
Definition: Numerics.h:122
const double LOG2
Definition: Numerics.h:16
constexpr double tau
Definition: Numerics.h:21
double lfactorial(double x)
Definition: Numerics.h:157
T get_log1pexp_breakout_bound(const double precision, T f(T), const T lower=-1e6, const T upper=1e6)
Tis takes a function f(x) and finds a bound B so that f(x) < precision for all x < B...
Definition: Numerics.h:69
constexpr double NaN
Definition: Numerics.h:19
double mygamma(double v)
Definition: Numerics.h:150