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 - 0.5 * 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 - 0.5 * 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 first and last inclusive.
142 int32 RandInt(int32 first, int32 last, 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 proportional
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 // returns log(exp(x) + exp(y)).
184 inline double LogAdd(double x, double y) {
185  double diff;
186 
187  if (x < y) {
188  diff = x - y;
189  x = y;
190  } else {
191  diff = y - x;
192  }
193  // diff is negative. x is now the larger one.
194 
195  if (diff >= kMinLogDiffDouble) {
196  double res;
197  res = x + Log1p(Exp(diff));
198  return res;
199  } else {
200  return x; // return the larger one.
201  }
202 }
203 
204 
205 // returns log(exp(x) + exp(y)).
206 inline float LogAdd(float x, float y) {
207  float diff;
208 
209  if (x < y) {
210  diff = x - y;
211  x = y;
212  } else {
213  diff = y - x;
214  }
215  // diff is negative. x is now the larger one.
216 
217  if (diff >= kMinLogDiffFloat) {
218  float res;
219  res = x + Log1p(Exp(diff));
220  return res;
221  } else {
222  return x; // return the larger one.
223  }
224 }
225 
226 
227 // returns log(exp(x) - exp(y)).
228 inline double LogSub(double x, double y) {
229  if (y >= x) { // Throws exception if y>=x.
230  if (y == x)
231  return kLogZeroDouble;
232  else
233  KALDI_ERR << "Cannot subtract a larger from a smaller number.";
234  }
235 
236  double diff = y - x; // Will be negative.
237  double res = x + Log(1.0 - Exp(diff));
238 
239  // res might be NAN if diff ~0.0, and 1.0-exp(diff) == 0 to machine precision
240  if (KALDI_ISNAN(res))
241  return kLogZeroDouble;
242  return res;
243 }
244 
245 
246 // returns log(exp(x) - exp(y)).
247 inline float LogSub(float x, float y) {
248  if (y >= x) { // Throws exception if y>=x.
249  if (y == x)
250  return kLogZeroDouble;
251  else
252  KALDI_ERR << "Cannot subtract a larger from a smaller number.";
253  }
254 
255  float diff = y - x; // Will be negative.
256  float res = x + Log(1.0f - Exp(diff));
257 
258  // res might be NAN if diff ~0.0, and 1.0-exp(diff) == 0 to machine precision
259  if (KALDI_ISNAN(res))
260  return kLogZeroFloat;
261  return res;
262 }
263 
265 static inline bool ApproxEqual(float a, float b,
266  float relative_tolerance = 0.001) {
267  // a==b handles infinities.
268  if (a == b) return true;
269  float diff = std::abs(a-b);
270  if (diff == std::numeric_limits<float>::infinity()
271  || diff != diff) return false; // diff is +inf or nan.
272  return (diff <= relative_tolerance*(std::abs(a)+std::abs(b)));
273 }
274 
276 static inline void AssertEqual(float a, float b,
277  float relative_tolerance = 0.001) {
278  // a==b handles infinities.
279  KALDI_ASSERT(ApproxEqual(a, b, relative_tolerance));
280 }
281 
282 
283 // RoundUpToNearestPowerOfTwo does the obvious thing. It crashes if n <= 0.
285 
287 static inline int32 DivideRoundingDown(int32 a, int32 b) {
288  KALDI_ASSERT(b != 0);
289  if (a * b >= 0)
290  return a / b;
291  else if (a < 0)
292  return (a - b + 1) / b;
293  else
294  return (a - b - 1) / b;
295 }
296 
297 template<class I> I Gcd(I m, I n) {
298  if (m == 0 || n == 0) {
299  if (m == 0 && n == 0) { // gcd not defined, as all integers are divisors.
300  KALDI_ERR << "Undefined GCD since m = 0, n = 0.";
301  }
302  return (m == 0 ? (n > 0 ? n : -n) : ( m > 0 ? m : -m));
303  // return absolute value of whichever is nonzero
304  }
305  // could use compile-time assertion
306  // but involves messing with complex template stuff.
307  KALDI_ASSERT(std::numeric_limits<I>::is_integer);
308  while (1) {
309  m %= n;
310  if (m == 0) return (n > 0 ? n : -n);
311  n %= m;
312  if (n == 0) return (m > 0 ? m : -m);
313  }
314 }
315 
318 template<class I> I Lcm(I m, I n) {
319  KALDI_ASSERT(m > 0 && n > 0);
320  I gcd = Gcd(m, n);
321  return gcd * (m/gcd) * (n/gcd);
322 }
323 
324 
325 template<class I> void Factorize(I m, std::vector<I> *factors) {
326  // Splits a number into its prime factors, in sorted order from
327  // least to greatest, with duplication. A very inefficient
328  // algorithm, which is mainly intended for use in the
329  // mixed-radix FFT computation (where we assume most factors
330  // are small).
331  KALDI_ASSERT(factors != NULL);
332  KALDI_ASSERT(m >= 1); // Doesn't work for zero or negative numbers.
333  factors->clear();
334  I small_factors[10] = { 2, 3, 5, 7, 11, 13, 17, 19, 23, 29 };
335 
336  // First try small factors.
337  for (I i = 0; i < 10; i++) {
338  if (m == 1) return; // We're done.
339  while (m % small_factors[i] == 0) {
340  m /= small_factors[i];
341  factors->push_back(small_factors[i]);
342  }
343  }
344  // Next try all odd numbers starting from 31.
345  for (I j = 31;; j += 2) {
346  if (m == 1) return;
347  while (m % j == 0) {
348  m /= j;
349  factors->push_back(j);
350  }
351  }
352 }
353 
354 inline double Hypot(double x, double y) { return hypot(x, y); }
355 inline float Hypot(float x, float y) { return hypotf(x, y); }
356 
357 
358 
359 
360 } // namespace kaldi
361 
362 
363 #endif // KALDI_BASE_KALDI_MATH_H_
This code computes Goodness of Pronunciation (GOP) and extracts phone-level pronunciation feature for...
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:287
#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:72
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:297
float RandGauss(struct RandomState *state=NULL)
Definition: kaldi-math.h:155
kaldi::int32 int32
int32 RoundUpToNearestPowerOfTwo(int32 n)
Definition: kaldi-math.cc:32
I Lcm(I m, I n)
Returns the least common multiple of two integers.
Definition: kaldi-math.h:318
int32 RandPoisson(float lambda, struct RandomState *state)
Definition: kaldi-math.cc:126
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:228
static const double kMinLogDiffDouble
Definition: kaldi-math.h:124
void RandGauss2(float *a, float *b, RandomState *state)
Definition: kaldi-math.cc:139
struct rnnlm::@11::@12 n
#define KALDI_ERR
Definition: kaldi-error.h:147
int Rand(struct RandomState *state)
Definition: kaldi-math.cc:45
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:185
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:276
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:265
int32 RandInt(int32 min_val, int32 max_val, struct RandomState *state)
Definition: kaldi-math.cc:95
double Hypot(double x, double y)
Definition: kaldi-math.h:354
void Factorize(I m, std::vector< I > *factors)
Definition: kaldi-math.h:325