kaldi-math.cc 5.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162
  1. // base/kaldi-math.cc
  2. // Copyright 2009-2011 Microsoft Corporation; Yanmin Qian;
  3. // Saarland University; Jan Silovsky
  4. // See ../../COPYING for clarification regarding multiple authors
  5. //
  6. // Licensed under the Apache License, Version 2.0 (the "License");
  7. // you may not use this file except in compliance with the License.
  8. // You may obtain a copy of the License at
  9. //
  10. // http://www.apache.org/licenses/LICENSE-2.0
  11. //
  12. // THIS CODE IS PROVIDED *AS IS* BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
  13. // KIND, EITHER EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION ANY IMPLIED
  14. // WARRANTIES OR CONDITIONS OF TITLE, FITNESS FOR A PARTICULAR PURPOSE,
  15. // MERCHANTABLITY OR NON-INFRINGEMENT.
  16. // See the Apache 2 License for the specific language governing permissions and
  17. // limitations under the License.
  18. #include "base/kaldi-math.h"
  19. #ifndef _MSC_VER
  20. #include <stdlib.h>
  21. #include <unistd.h>
  22. #endif
  23. #include <string>
  24. #include <mutex>
  25. namespace kaldi {
  26. // These routines are tested in matrix/matrix-test.cc
  27. int32 RoundUpToNearestPowerOfTwo(int32 n) {
  28. KALDI_ASSERT(n > 0);
  29. n--;
  30. n |= n >> 1;
  31. n |= n >> 2;
  32. n |= n >> 4;
  33. n |= n >> 8;
  34. n |= n >> 16;
  35. return n+1;
  36. }
  37. static std::mutex _RandMutex;
  38. int Rand(struct RandomState* state) {
  39. #if !defined(_POSIX_THREAD_SAFE_FUNCTIONS)
  40. // On Windows and Cygwin, just call Rand()
  41. return rand();
  42. #else
  43. if (state) {
  44. return rand_r(&(state->seed));
  45. } else {
  46. std::lock_guard<std::mutex> lock(_RandMutex);
  47. return rand();
  48. }
  49. #endif
  50. }
  51. RandomState::RandomState() {
  52. // we initialize it as Rand() + 27437 instead of just Rand(), because on some
  53. // systems, e.g. at the very least Mac OSX Yosemite and later, it seems to be
  54. // the case that rand_r when initialized with rand() will give you the exact
  55. // same sequence of numbers that rand() will give if you keep calling rand()
  56. // after that initial call. This can cause problems with repeated sequences.
  57. // For example if you initialize two RandomState structs one after the other
  58. // without calling rand() in between, they would give you the same sequence
  59. // offset by one (if we didn't have the "+ 27437" in the code). 27437 is just
  60. // a randomly chosen prime number.
  61. seed = unsigned(Rand()) + 27437;
  62. }
  63. bool WithProb(BaseFloat prob, struct RandomState* state) {
  64. KALDI_ASSERT(prob >= 0 && prob <= 1.1); // prob should be <= 1.0,
  65. // but we allow slightly larger values that could arise from roundoff in
  66. // previous calculations.
  67. KALDI_COMPILE_TIME_ASSERT(RAND_MAX > 128 * 128);
  68. if (prob == 0) return false;
  69. else if (prob == 1.0) return true;
  70. else if (prob * RAND_MAX < 128.0) {
  71. // prob is very small but nonzero, and the "main algorithm"
  72. // wouldn't work that well. So: with probability 1/128, we
  73. // return WithProb (prob * 128), else return false.
  74. if (Rand(state) < RAND_MAX / 128) { // with probability 128...
  75. // Note: we know that prob * 128.0 < 1.0, because
  76. // we asserted RAND_MAX > 128 * 128.
  77. return WithProb(prob * 128.0);
  78. } else {
  79. return false;
  80. }
  81. } else {
  82. return (Rand(state) < ((RAND_MAX + static_cast<BaseFloat>(1.0)) * prob));
  83. }
  84. }
  85. int32 RandInt(int32 min_val, int32 max_val, struct RandomState* state) {
  86. // This is not exact.
  87. KALDI_ASSERT(max_val >= min_val);
  88. if (max_val == min_val) return min_val;
  89. #ifdef _MSC_VER
  90. // RAND_MAX is quite small on Windows -> may need to handle larger numbers.
  91. if (RAND_MAX > (max_val-min_val)*8) {
  92. // *8 to avoid large inaccuracies in probability, from the modulus...
  93. return min_val +
  94. ((unsigned int)Rand(state) % (unsigned int)(max_val+1-min_val));
  95. } else {
  96. if ((unsigned int)(RAND_MAX*RAND_MAX) >
  97. (unsigned int)((max_val+1-min_val)*8)) {
  98. // *8 to avoid inaccuracies in probability, from the modulus...
  99. return min_val + ( (unsigned int)( (Rand(state)+RAND_MAX*Rand(state)))
  100. % (unsigned int)(max_val+1-min_val));
  101. } else {
  102. KALDI_ERR << "rand_int failed because we do not support such large "
  103. "random numbers. (Extend this function).";
  104. }
  105. }
  106. #else
  107. return min_val +
  108. (static_cast<int32>(Rand(state)) % static_cast<int32>(max_val+1-min_val));
  109. #endif
  110. }
  111. // Returns poisson-distributed random number.
  112. // Take care: this takes time proportional
  113. // to lambda. Faster algorithms exist but are more complex.
  114. int32 RandPoisson(float lambda, struct RandomState* state) {
  115. // Knuth's algorithm.
  116. KALDI_ASSERT(lambda >= 0);
  117. float L = expf(-lambda), p = 1.0;
  118. int32 k = 0;
  119. do {
  120. k++;
  121. float u = RandUniform(state);
  122. p *= u;
  123. } while (p > L);
  124. return k-1;
  125. }
  126. void RandGauss2(float *a, float *b, RandomState *state) {
  127. KALDI_ASSERT(a);
  128. KALDI_ASSERT(b);
  129. float u1 = RandUniform(state);
  130. float u2 = RandUniform(state);
  131. u1 = sqrtf(-2.0f * logf(u1));
  132. u2 = 2.0f * M_PI * u2;
  133. *a = u1 * cosf(u2);
  134. *b = u1 * sinf(u2);
  135. }
  136. void RandGauss2(double *a, double *b, RandomState *state) {
  137. KALDI_ASSERT(a);
  138. KALDI_ASSERT(b);
  139. float a_float, b_float;
  140. // Just because we're using doubles doesn't mean we need super-high-quality
  141. // random numbers, so we just use the floating-point version internally.
  142. RandGauss2(&a_float, &b_float, state);
  143. *a = a_float;
  144. *b = b_float;
  145. }
  146. } // end namespace kaldi