All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
kaldi-math.h
Go to the documentation of this file.
1 // base/kaldi-math.h
2 
3 // Copyright 2009-2011 Ondrej Glembek; Microsoft Corporation; Yanmin Qian;
4 // Jan Silovsky; Saarland University
5 //
6 // See ../../COPYING for clarification regarding multiple authors
7 //
8 // Licensed under the Apache License, Version 2.0 (the "License");
9 // you may not use this file except in compliance with the License.
10 // You may obtain a copy of the License at
11 //
12 // http://www.apache.org/licenses/LICENSE-2.0
13 //
14 // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
15 // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
16 // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
17 // MERCHANTABLITY OR NON-INFRINGEMENT.
18 // See the Apache 2 License for the specific language governing permissions and
19 // limitations under the License.
20 
21 #ifndef KALDI_BASE_KALDI_MATH_H_
22 #define KALDI_BASE_KALDI_MATH_H_ 1
23 
24 #ifdef _MSC_VER
25 #include <float.h>
26 #endif
27 
28 #include <cmath>
29 #include <limits>
30 #include <vector>
31 
32 #include "base/kaldi-types.h"
33 #include "base/kaldi-common.h"
34 
35 
36 #ifndef DBL_EPSILON
37 #define DBL_EPSILON 2.2204460492503131e-16
38 #endif
39 #ifndef FLT_EPSILON
40 #define FLT_EPSILON 1.19209290e-7f
41 #endif
42 
43 #ifndef M_PI
44 #define M_PI 3.1415926535897932384626433832795
45 #endif
46 
47 #ifndef M_SQRT2
48 #define M_SQRT2 1.4142135623730950488016887
49 #endif
50 
51 #ifndef M_2PI
52 #define M_2PI 6.283185307179586476925286766559005
53 #endif
54 
55 #ifndef M_SQRT1_2
56 #define M_SQRT1_2 0.7071067811865475244008443621048490
57 #endif
58 
59 #ifndef M_LOG_2PI
60 #define M_LOG_2PI 1.8378770664093454835606594728112
61 #endif
62 
63 #ifndef M_LN2
64 #define M_LN2 0.693147180559945309417232121458
65 #endif
66 
67 #ifndef M_LN10
68 #define M_LN10 2.302585092994045684017991454684
69 #endif
70 
71 
72 #define KALDI_ISNAN std::isnan
73 #define KALDI_ISINF std::isinf
74 #define KALDI_ISFINITE(x) std::isfinite(x)
75 
76 #if !defined(KALDI_SQR)
77 # define KALDI_SQR(x) ((x) * (x))
78 #endif
79 
80 namespace kaldi {
81 
82 #if !defined(_MSC_VER) || (_MSC_VER >= 1900)
83 inline double Exp(double x) { return exp(x); }
84 #ifndef KALDI_NO_EXPF
85 inline float Exp(float x) { return expf(x); }
86 #else
87 inline float Exp(float x) { return exp(static_cast<double>(x)); }
88 #endif // KALDI_NO_EXPF
89 #else
90 inline double Exp(double x) { return exp(x); }
91 #if !defined(__INTEL_COMPILER) && _MSC_VER == 1800 && defined(_M_X64)
92 // Microsoft CL v18.0 buggy 64-bit implementation of
93 // expf() incorrectly returns -inf for exp(-inf).
94 inline float Exp(float x) { return exp(static_cast<double>(x)); }
95 #else
96 inline float Exp(float x) { return expf(x); }
97 #endif // !defined(__INTEL_COMPILER) && _MSC_VER == 1800 && defined(_M_X64)
98 #endif // !defined(_MSC_VER) || (_MSC_VER >= 1900)
99 
100 inline double Log(double x) { return log(x); }
101 inline float Log(float x) { return logf(x); }
102 
103 #if !defined(_MSC_VER) || (_MSC_VER >= 1700)
104 inline double Log1p(double x) { return log1p(x); }
105 inline float Log1p(float x) { return log1pf(x); }
106 #else
107 inline double Log1p(double x) {
108  const double cutoff = 1.0e-08;
109  if (x < cutoff)
110  return x - 2 * x * x;
111  else
112  return Log(1.0 + x);
113 }
114 
115 inline float Log1p(float x) {
116  const float cutoff = 1.0e-07;
117  if (x < cutoff)
118  return x - 2 * x * x;
119  else
120  return Log(1.0 + x);
121 }
122 #endif
123 
124 static const double kMinLogDiffDouble = Log(DBL_EPSILON); // negative!
125 static const float kMinLogDiffFloat = Log(FLT_EPSILON); // negative!
126 
127 // -infinity
128 const float kLogZeroFloat = -std::numeric_limits<float>::infinity();
129 const double kLogZeroDouble = -std::numeric_limits<double>::infinity();
130 const BaseFloat kLogZeroBaseFloat = -std::numeric_limits<BaseFloat>::infinity();
131 
132 // Returns a random integer between 0 and RAND_MAX, inclusive
133 int Rand(struct RandomState* state = NULL);
134 
135 // State for thread-safe random number generator
136 struct RandomState {
137  RandomState();
138  unsigned seed;
139 };
140 
141 // Returns a random integer between min and max inclusive.
142 int32 RandInt(int32 min, int32 max, struct RandomState* state = NULL);
143 
144 // Returns true with probability "prob",
145 bool WithProb(BaseFloat prob, struct RandomState* state = NULL);
146 // with 0 <= prob <= 1 [we check this].
147 // Internally calls Rand(). This function is carefully implemented so
148 // that it should work even if prob is very small.
149 
151 inline float RandUniform(struct RandomState* state = NULL) {
152  return static_cast<float>((Rand(state) + 1.0) / (RAND_MAX+2.0));
153 }
154 
155 inline float RandGauss(struct RandomState* state = NULL) {
156  return static_cast<float>(sqrtf (-2 * Log(RandUniform(state)))
157  * cosf(2*M_PI*RandUniform(state)));
158 }
159 
160 // Returns poisson-distributed random number. Uses Knuth's algorithm.
161 // Take care: this takes time proportinal
162 // to lambda. Faster algorithms exist but are more complex.
163 int32 RandPoisson(float lambda, struct RandomState* state = NULL);
164 
165 // Returns a pair of gaussian random numbers. Uses Box-Muller transform
166 void RandGauss2(float *a, float *b, RandomState *state = NULL);
167 void RandGauss2(double *a, double *b, RandomState *state = NULL);
168 
169 // Also see Vector<float,double>::RandCategorical().
170 
171 // This is a randomized pruning mechanism that preserves expectations,
172 // that we typically use to prune posteriors.
173 template<class Float>
174 inline Float RandPrune(Float post, BaseFloat prune_thresh,
175  struct RandomState* state = NULL) {
176  KALDI_ASSERT(prune_thresh >= 0.0);
177  if (post == 0.0 || std::abs(post) >= prune_thresh)
178  return post;
179  return (post >= 0 ? 1.0 : -1.0) *
180  (RandUniform(state) <= fabs(post)/prune_thresh ? prune_thresh : 0.0);
181 }
182 
183 
184 inline double LogAdd(double x, double y) {
185  double diff;
186  if (x < y) {
187  diff = x - y;
188  x = y;
189  } else {
190  diff = y - x;
191  }
192  // diff is negative. x is now the larger one.
193 
194  if (diff >= kMinLogDiffDouble) {
195  double res;
196  res = x + Log1p(Exp(diff));
197  return res;
198  } else {
199  return x; // return the larger one.
200  }
201 }
202 
203 
204 inline float LogAdd(float x, float y) {
205  float diff;
206  if (x < y) {
207  diff = x - y;
208  x = y;
209  } else {
210  diff = y - x;
211  }
212  // diff is negative. x is now the larger one.
213 
214  if (diff >= kMinLogDiffFloat) {
215  float res;
216  res = x + Log1p(Exp(diff));
217  return res;
218  } else {
219  return x; // return the larger one.
220  }
221 }
222 
223 
224 // returns exp(x) - exp(y).
225 inline double LogSub(double x, double y) {
226  if (y >= x) { // Throws exception if y>=x.
227  if (y == x)
228  return kLogZeroDouble;
229  else
230  KALDI_ERR << "Cannot subtract a larger from a smaller number.";
231  }
232 
233  double diff = y - x; // Will be negative.
234  double res = x + Log(1.0 - Exp(diff));
235 
236  // res might be NAN if diff ~0.0, and 1.0-exp(diff) == 0 to machine precision
237  if (KALDI_ISNAN(res))
238  return kLogZeroDouble;
239  return res;
240 }
241 
242 
243 // returns exp(x) - exp(y).
244 inline float LogSub(float x, float y) {
245  if (y >= x) { // Throws exception if y>=x.
246  if (y == x)
247  return kLogZeroDouble;
248  else
249  KALDI_ERR << "Cannot subtract a larger from a smaller number.";
250  }
251 
252  float diff = y - x; // Will be negative.
253  float res = x + Log(1.0f - Exp(diff));
254 
255  // res might be NAN if diff ~0.0, and 1.0-exp(diff) == 0 to machine precision
256  if (KALDI_ISNAN(res))
257  return kLogZeroFloat;
258  return res;
259 }
260 
262 static inline bool ApproxEqual(float a, float b,
263  float relative_tolerance = 0.001) {
264  // a==b handles infinities.
265  if (a == b) return true;
266  float diff = std::abs(a-b);
267  if (diff == std::numeric_limits<float>::infinity()
268  || diff != diff) return false; // diff is +inf or nan.
269  return (diff <= relative_tolerance*(std::abs(a)+std::abs(b)));
270 }
271 
273 static inline void AssertEqual(float a, float b,
274  float relative_tolerance = 0.001) {
275  // a==b handles infinities.
276  KALDI_ASSERT(ApproxEqual(a, b, relative_tolerance));
277 }
278 
279 
280 // RoundUpToNearestPowerOfTwo does the obvious thing. It crashes if n <= 0.
281 int32 RoundUpToNearestPowerOfTwo(int32 n);
282 
284 static inline int32 DivideRoundingDown(int32 a, int32 b) {
285  KALDI_ASSERT(b != 0);
286  if (a * b >= 0)
287  return a / b;
288  else if (a < 0)
289  return (a - b + 1) / b;
290  else
291  return (a - b - 1) / b;
292 }
293 
294 template<class I> I Gcd(I m, I n) {
295  if (m == 0 || n == 0) {
296  if (m == 0 && n == 0) { // gcd not defined, as all integers are divisors.
297  KALDI_ERR << "Undefined GCD since m = 0, n = 0.";
298  }
299  return (m == 0 ? (n > 0 ? n : -n) : ( m > 0 ? m : -m));
300  // return absolute value of whichever is nonzero
301  }
302  // could use compile-time assertion
303  // but involves messing with complex template stuff.
304  KALDI_ASSERT(std::numeric_limits<I>::is_integer);
305  while (1) {
306  m %= n;
307  if (m == 0) return (n > 0 ? n : -n);
308  n %= m;
309  if (n == 0) return (m > 0 ? m : -m);
310  }
311 }
312 
315 template<class I> I Lcm(I m, I n) {
316  KALDI_ASSERT(m > 0 && n > 0);
317  I gcd = Gcd(m, n);
318  return gcd * (m/gcd) * (n/gcd);
319 }
320 
321 
322 template<class I> void Factorize(I m, std::vector<I> *factors) {
323  // Splits a number into its prime factors, in sorted order from
324  // least to greatest, with duplication. A very inefficient
325  // algorithm, which is mainly intended for use in the
326  // mixed-radix FFT computation (where we assume most factors
327  // are small).
328  KALDI_ASSERT(factors != NULL);
329  KALDI_ASSERT(m >= 1); // Doesn't work for zero or negative numbers.
330  factors->clear();
331  I small_factors[10] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };
332 
333  // First try small factors.
334  for (I i = 0; i < 10; i++) {
335  if (m == 1) return; // We're done.
336  while (m % small_factors[i] == 0) {
337  m /= small_factors[i];
338  factors->push_back(small_factors[i]);
339  }
340  }
341  // Next try all odd numbers starting from 31.
342  for (I j = 31;; j += 2) {
343  if (m == 1) return;
344  while (m % j == 0) {
345  m /= j;
346  factors->push_back(j);
347  }
348  }
349 }
350 
351 inline double Hypot(double x, double y) { return hypot(x, y); }
352 inline float Hypot(float x, float y) { return hypotf(x, y); }
353 
354 
355 
356 
357 } // namespace kaldi
358 
359 
360 #endif // KALDI_BASE_KALDI_MATH_H_
Relabels neural network egs with the read pdf-id alignments.
Definition: chain.dox:20
double Exp(double x)
Definition: kaldi-math.h:83
static int32 DivideRoundingDown(int32 a, int32 b)
Returns a / b, rounding towards negative infinity in all cases.
Definition: kaldi-math.h:284
#define DBL_EPSILON
Definition: kaldi-math.h:37
float RandUniform(struct RandomState *state=NULL)
Returns a random number strictly between 0 and 1.
Definition: kaldi-math.h:151
#define M_PI
Definition: kaldi-math.h:44
bool WithProb(BaseFloat prob, struct RandomState *state)
Definition: kaldi-math.cc:77
Float RandPrune(Float post, BaseFloat prune_thresh, struct RandomState *state=NULL)
Definition: kaldi-math.h:174
I Gcd(I m, I n)
Definition: kaldi-math.h:294
float RandGauss(struct RandomState *state=NULL)
Definition: kaldi-math.h:155
int32 RoundUpToNearestPowerOfTwo(int32 n)
Definition: kaldi-math.cc:31
I Lcm(I m, I n)
Returns the least common multiple of two integers.
Definition: kaldi-math.h:315
int32 RandPoisson(float lambda, struct RandomState *state)
Definition: kaldi-math.cc:133
float BaseFloat
Definition: kaldi-types.h:29
double Log(double x)
Definition: kaldi-math.h:100
double LogSub(double x, double y)
Definition: kaldi-math.h:225
static const double kMinLogDiffDouble
Definition: kaldi-math.h:124
void RandGauss2(float *a, float *b, RandomState *state)
Definition: kaldi-math.cc:146
struct rnnlm::@11::@12 n
#define KALDI_ERR
Definition: kaldi-error.h:127
int Rand(struct RandomState *state)
Definition: kaldi-math.cc:46
static const float kMinLogDiffFloat
Definition: kaldi-math.h:125
double LogAdd(double x, double y)
Definition: kaldi-math.h:184
#define FLT_EPSILON
Definition: kaldi-math.h:40
#define KALDI_ISNAN
Definition: kaldi-math.h:72
#define KALDI_ASSERT(cond)
Definition: kaldi-error.h:169
static void AssertEqual(float a, float b, float relative_tolerance=0.001)
assert abs(a - b) <= relative_tolerance * (abs(a)+abs(b))
Definition: kaldi-math.h:273
const float kLogZeroFloat
Definition: kaldi-math.h:128
double Log1p(double x)
Definition: kaldi-math.h:104
const BaseFloat kLogZeroBaseFloat
Definition: kaldi-math.h:130
const double kLogZeroDouble
Definition: kaldi-math.h:129
static bool ApproxEqual(float a, float b, float relative_tolerance=0.001)
return abs(a - b) <= relative_tolerance * (abs(a)+abs(b)).
Definition: kaldi-math.h:262
int32 RandInt(int32 min_val, int32 max_val, struct RandomState *state)
Definition: kaldi-math.cc:100
double Hypot(double x, double y)
Definition: kaldi-math.h:351
void Factorize(I m, std::vector< I > *factors)
Definition: kaldi-math.h:322