1. Benchmarking slow and fast variants of Renyi div. #272

2. Futher steps towards reusing the same thread-local random generator #57

similarity_search/apps/bench_distfunc.cc

@@ -1964,6 +1964,82 @@ void TestSparseJaccardSimilarity(const string& dataFile, size_t N, size_t Rep) {
 
 }
 
+template <class T>
+void TestRenyiDivSlow(size_t N, size_t dim, size_t Rep, T alpha) {
+    T* pArr = new T[N * dim];
+
+    T *p = pArr;
+    for (size_t i = 0; i < N; ++i, p+= dim) {
+        GenRandVect(p, dim, T(RANGE_SMALL), T(1.0), true /* norm. for regular KL */);
+    }
+
+    WallClockTimer  t;
+
+    t.reset();
+
+    T DiffSum = 0;
+
+    T fract = T(1)/N;
+
+    for (size_t i = 0; i < Rep; ++i) {
+        for (size_t j = 1; j < N; ++j) {
+            DiffSum += 0.01f * renyiDivergenceSlow(pArr + j*dim, pArr + (j-1)*dim, dim, alpha) / N;
+        }
+        /* 
+         * Multiplying by 0.01 and dividing the sum by N is to prevent Intel from "cheating":
+         *
+         * http://searchivarius.org/blog/problem-previous-version-intels-library-benchmark
+         */
+        DiffSum *= fract;
+    }
+
+    uint64_t tDiff = t.split();
+
+    LOG(LIB_INFO) << "Ignore: " << DiffSum;
+    LOG(LIB_INFO) << typeid(T).name() << " " << "Elapsed: " << tDiff / 1e3 << " ms " << " # of slow Renyi-div. (alpha=" << alpha << ") per second: " << (1e6/tDiff) * N * Rep ;
+
+    delete [] pArr;
+
+}
+
+template <class T>
+void TestRenyiDivFast(size_t N, size_t dim, size_t Rep, T alpha) {
+    T* pArr = new T[N * dim];
+
+    T *p = pArr;
+    for (size_t i = 0; i < N; ++i, p+= dim) {
+        GenRandVect(p, dim, T(RANGE_SMALL), T(1.0), true /* norm. for regular KL */);
+    }
+
+    WallClockTimer  t;
+
+    t.reset();
+
+    T DiffSum = 0;
+
+    T fract = T(1)/N;
+
+    for (size_t i = 0; i < Rep; ++i) {
+        for (size_t j = 1; j < N; ++j) {
+            DiffSum += 0.01f * renyiDivergenceFast(pArr + j*dim, pArr + (j-1)*dim, dim, alpha) / N;
+        }
+        /* 
+         * Multiplying by 0.01 and dividing the sum by N is to prevent Intel from "cheating":
+         *
+         * http://searchivarius.org/blog/problem-previous-version-intels-library-benchmark
+         */
+        DiffSum *= fract;
+    }
+
+    uint64_t tDiff = t.split();
+
+    LOG(LIB_INFO) << "Ignore: " << DiffSum;
+    LOG(LIB_INFO) << typeid(T).name() << " " << "Elapsed: " << tDiff / 1e3 << " ms " << " # of fast Renyi-div. (alpha=" << alpha << ") per second: " << (1e6/tDiff) * N * Rep ;
+
+    delete [] pArr;
+
+}
+
 }  // namespace similarity
 
 using namespace similarity;
@@ -1977,6 +2053,34 @@ int main(int argc, char* argv[]) {
 
     int dim = 128;
 
+    float alphaStepSlow = 1.0f/4;
+    for (float alpha = alphaStepSlow; alpha <= 2; alpha+= alphaStepSlow) {
+      if (fabs(alpha - 1) < 1e-6) continue;
+      nTest++;
+      TestRenyiDivSlow<float>(1024, dim, 100, alpha);
+    }
+#if TEST_SPEED_DOUBLE
+    for (double alpha = alphaStepSlow; alpha <= 2; alpha+= alphaStepSlow) {
+      if (fabs(alpha - 1) < 1e-6) continue;
+      nTest++;
+      TestRenyiDivSlow<double>(1024, dim, 100, alpha);
+    }
+#endif
+
+    float alphaStepFast = 1.0f/32;
+    for (float alpha = alphaStepFast; alpha <= 2; alpha+= alphaStepFast) {
+      if (fabs(alpha - 1) < 1e-6) continue;
+      nTest++;
+      TestRenyiDivFast<float>(1024, dim, 100, alpha);
+    }
+#if TEST_SPEED_DOUBLE
+    for (double alpha = alphaStepFast; alpha <= 2; alpha+= alphaStepFast) {
+      if (fabs(alpha - 1) < 1e-6) continue;
+      nTest++;
+      TestRenyiDivFast<double>(1024, dim, 100, alpha);
+    }
+#endif
+
 #if defined(WITH_EXTRAS)
     nTest++;
     TestSQFDMinus<float>(2000, 50);

similarity_search/apps/tune_vptree.cc

@@ -113,8 +113,7 @@ void RunExper(unsigned AddRestartQty,
 
   LOG(LIB_INFO) << "We are going to tune parameters for " << MethodName;
 
-  static  std::random_device          rd;
-  static  std::mt19937                engine(rd());
+  static  thread_local auto&          engine(GET_RANDOM_GENERATOR);
   static  std::normal_distribution<>  normGen(0.0f, log(FullFactor));
 
   AnyParamManager pmgr(IndexParams);

similarity_search/include/utils.h

@@ -65,13 +65,13 @@ typedef SSIZE_T ssize_t;
 
 #define FIELD_DELIMITER ':'
 
+#define GET_RANDOM_GENERATOR    getRandomGenerator<std::mt19937>()
+
 namespace similarity {
 
 using std::string;
 using std::vector;
 using std::stringstream;
-using std::random_device;
-using std::mt19937;
 
 using namespace std;
 
@@ -85,6 +85,19 @@ bool DoesFileExist(const char *filename);
 inline bool DoesFileExist(const string &filename) { return DoesFileExist(filename.c_str()); }
 
 extern int randomSeed;
+
+/*
+ * Random number generation is thread safe when respective
+ * objects are not shared among threads. So, we will keep one
+ * random number generator per thread.
+ */ 
+template <class RandGenType>
+inline RandGenType & getRandomGenerator() {
+    static thread_local RandGenType gen(randomSeed);
+
+    return gen;
+}
+
 // random 32-bit integer number
 inline int32_t RandomInt() {
     /*
@@ -93,10 +106,9 @@ inline int32_t RandomInt() {
      * random number generator per thread.
     */ 
     // thread_local is static by default, but let's keep it static for clarity
-    static thread_local mt19937 gen(randomSeed);
     static thread_local std::uniform_int_distribution<int32_t> distr(0, std::numeric_limits<int32_t>::max());
 
-    return distr(gen); 
+    return distr(GET_RANDOM_GENERATOR); 
 }
 
 template <class T>
@@ -108,10 +120,9 @@ inline T RandomReal() {
      * random number generator per thread.
     */ 
     // thread_local is static by default, but let's keep it static for clarity
-    static thread_local mt19937 gen(randomSeed);
     static thread_local std::uniform_real_distribution<T> distr(0, 1);
 
-    return distr(gen); 
+    return distr(GET_RANDOM_GENERATOR); 
 }
 
 void RStrip(char* str);

similarity_search/src/randproj_util.cc

@@ -23,6 +23,7 @@
 #include "randproj_util.h"
 #include "distcomp.h"
 #include "logging.h"
+#include "utils.h"
 
 namespace similarity {
 
@@ -32,16 +33,15 @@ template <class dist_t> void initRandProj(size_t nSrcDim, size_t nDstDim,
                                          bool bDoOrth,
                                          vector<vector<dist_t>>& projMatr) {
   // Static is thread-safe in C++-11
-  static  std::random_device          rd;
-  static  std::mt19937                engine(rd());
+  static thread_local auto&           randGen(GET_RANDOM_GENERATOR);
   static  std::normal_distribution<>  normGen(0.0f, 1.0f);
 
   // 1. Create normally distributed vectors
   projMatr.resize(nDstDim);
   for (size_t i = 0; i < nDstDim; ++i) {
     projMatr[i].resize(nSrcDim);
     for (size_t j = 0; j < nSrcDim; ++j)
-      projMatr[i][j] = normGen(engine);
+      projMatr[i][j] = normGen(randGen);
   }
   /* 
    * 2. If bDoOrth == true, normalize the basis using the numerically stable

similarity_search/src/searchoracle.cc

@@ -205,8 +205,8 @@ void PolynomialPruner<dist_t>::SetIndexTimeParams(AnyParamManager& pmgr) {
     float recall = 0, time_best = 0, impr_best = -1, alpha_left = 0, alpha_right = 0; 
     unsigned exp_left = 0, exp_right = 0;
 
-    static  std::random_device          rd;
-    static  std::mt19937                engine(rd());
+
+    static thread_local auto&           randGen(GET_RANDOM_GENERATOR);
     static  std::normal_distribution<>  normGen(0.0f, log(fullFactor));
 
 
@@ -218,8 +218,8 @@ void PolynomialPruner<dist_t>::SetIndexTimeParams(AnyParamManager& pmgr) {
 
       if (k > 0) {
         // Let's do some random normal fun
-        alpha_left_loc = exp(normGen(engine));
-        alpha_right_loc = exp(normGen(engine));
+        alpha_left_loc = exp(normGen(randGen));
+        alpha_right_loc = exp(normGen(randGen));
         LOG(LIB_INFO) << " RANDOM STARTING POINTS: " << alpha_left_loc << " " << alpha_right_loc;
       }