Mh randutils

src/xmipp/applications/tests/function_tests/test_fringe_processing_main.cpp

@@ -25,6 +25,7 @@
 
 #include <core/multidim_array.h>
 #include <reconstruction/fringe_processing.h>
+#include <utils/prng.h>
 #include <core/xmipp_image.h>
 #include <iostream>
 #include <gtest/gtest.h>
@@ -267,9 +268,10 @@ TEST_F( FringeProcessingTests, unwrapping)
     refPhase.setXmippOrigin();
     im.setXmippOrigin();
 
+    GaussGenerator<double> gauss_gen(0, noiseLevel);
     FOR_ALL_ELEMENTS_IN_ARRAY2D(refPhase)
     {
-        A2D_ELEM(refPhase,i,j) = (50*std::exp(-0.5*(std::pow(i*iMaxDim2,2)+std::pow(j*iMaxDim2,2))))+rnd_gaus(0,noiseLevel);
+        A2D_ELEM(refPhase,i,j) = (50*std::exp(-0.5*(std::pow(i*iMaxDim2,2)+std::pow(j*iMaxDim2,2))))+gauss_gen.rand();
         A2D_ELEM(im,i,j) = std::cos(A2D_ELEM(refPhase,i,j));
     }
     //We wrap the phase inside [0 2pi]

src/xmipp/libraries/classification/kerdensom.cpp

@@ -30,6 +30,7 @@
 //-----------------------------------------------------------------------------
 
 #include <fstream>
+#include <utils/prng.h>
 
 #include "kerdensom.h"
 
@@ -373,16 +374,18 @@ double KerDenSOM::randApproxGVC(const TS* _examples, const FuzzyMap* _som, doubl
             num += VV[j] * VV[j];
     }
 
-    init_random_generator();
+
+    GaussGenerator<double> gauss_gen;
+    gauss_gen.reseed(0);
     for (vv = 0; vv < numVectors; vv++)
     {
         for (j = 0; j < dim; j++)
-            tmpTS.theItems[vv][j] += rnd_gaus() * _dataSD;
+            tmpTS.theItems[vv][j] += gauss_gen.rand() * _dataSD;
     }
     updateV(&tmpSOM, &tmpTS, _reg);
     den = 0.0;
 
-    init_random_generator();
+    gauss_gen.reseed(0);
     for (vv = 0; vv < numVectors; vv++)
     {
         for (j = 0; j < dim; j++)
@@ -394,7 +397,7 @@ double KerDenSOM::randApproxGVC(const TS* _examples, const FuzzyMap* _som, doubl
         }
         r = 0.;
         for (j = 0; j < dim; j++)
-            r += VV[j] * rnd_gaus() * _dataSD;
+            r += VV[j] * gauss_gen.rand() * _dataSD;
         den += r * r;
     }
     if (den != 0)

src/xmipp/libraries/data/normalize.cpp

@@ -31,6 +31,7 @@
 #include "normalize.h"
 #include <core/metadata.h>
 #include <core/xmipp_image_generic.h>
+#include <utils/prng.h>
 
 /* Normalizations ---------------------------------------------------------- */
 void normalize_OldXmipp(MultidimArray<double> &I)
@@ -780,13 +781,14 @@ void ProgNormalize::processImage(const FileName &fnImg, const FileName &fnImgOut
         double avg=0., stddev=0., min=0., max=0., zz;
         img.computeStats(avg, stddev, min, max);
 
+        GaussGenerator<double> gauss_gen(avg, stddev);
         if ((min - avg) / stddev < thresh_black_dust && remove_black_dust)
         {
             FOR_ALL_ELEMENTS_IN_ARRAY3D(img)
             {
                 zz = (A3D_ELEM(img, k, i, j) - avg) / stddev;
                 if (zz < thresh_black_dust)
-                    A3D_ELEM(img, k, i, j) = rnd_gaus(avg,stddev);
+                    A3D_ELEM(img, k, i, j) = gauss_gen.rand();
             }
         }
 
@@ -796,7 +798,7 @@ void ProgNormalize::processImage(const FileName &fnImg, const FileName &fnImgOut
             {
                 zz = (A3D_ELEM(img, k, i, j) - avg) / stddev;
                 if (zz > thresh_white_dust)
-                    A3D_ELEM(img, k, i, j) = rnd_gaus(avg,stddev);
+                    A3D_ELEM(img, k, i, j) = gauss_gen.rand();
             }
         }
     }

src/xmipp/libraries/utils/deviation.h

@@ -0,0 +1,42 @@
+#ifndef DEVIATION_H
+#define DEVIATION_H
+
+#include <random>
+#include <cmath>
+
+/**
+ * Implementation of some new distribution that can be
+ * supplied to RandUtils
+ * Fun fact: it should've been normal distribution,
+ * but something went wrong
+*/
+template< typename T >
+class SomeNewDistribution {
+public:
+
+	SomeNewDistribution( T mean=0, T std=1 )
+	: mean( mean )
+	, std( std ) {}
+
+	template< typename generator >
+	T operator()( generator& gen ) {
+		T u, v, x, y, q;
+		do {
+			u = uniform( gen );
+			v = 1.7156 * ( uniform( gen ) - 0.5 );
+			x = u - 0.449871;
+			y = abs( v ) + 0.386595;
+			q = x * x + y * ( 0.19600 * y - 0.25472 * x );
+		} while( q > 0.27597 && ( q > 0.27846 || v * v > -4. * std::log( u ) * u * u ) );
+		return mean + std * v / u;
+	}
+
+private:
+
+	T mean;
+	T std;
+
+	std::uniform_real_distribution<T> uniform;
+};
+
+#endif //DEVIATION_H

src/xmipp/libraries/utils/main.cpp

@@ -0,0 +1,44 @@
+#include <iostream>
+#include <random>
+
+#include "prng.h"
+#include "deviation.h"
+
+/**
+ * Example program using new RandUtils
+*/
+
+int main() {
+	/**
+	 * You can supply own RNG generator
+	*/
+	RandUtils<double, std::exponential_distribution<>, std::mt19937_64> exp_rand( 2.0 );
+
+	/**
+	 * Some distributions have predefined generators
+	*/
+	ExponentialGenerator< double > eg( 2.0 );
+	GaussGenerator< double > gg( 2.0, 0.5 );
+	/**
+	* You can use RandUtils with own deviation class
+	*/
+	RandUtils< double, SomeNewDistribution< double > > exp_rand2;
+
+	/**
+	* reseed with own seed if determinism is needed
+	*/
+	exp_rand.reseed( 1 );
+
+	double x = 0;
+
+	for ( int i = 0; i < 100000; ++i ) {
+		x += gg.rand();
+	}
+
+	/**
+	 * if the generator is used only once, then just call static function
+	*/
+	x += ExponentialGenerator< double >::rand_once( 2.0 );
+
+	std::cout << "x = " << x << std::endl;
+}

src/xmipp/libraries/utils/prng.h

@@ -0,0 +1,57 @@
+#ifndef PRNG_H
+#define PRNG_H
+
+#include <random>
+
+/**
+ `distribution` must comply to https://en.cppreference.com/w/cpp/named_req/RandomNumberDistribution
+
+  standard library distributions are also on that page
+*/
+template< typename T, typename distribution, typename generator=std::mt19937 >
+class RandUtils {
+public:
+
+	template< typename... Args >
+	RandUtils( Args... args )
+	: distr( args... ) {
+		gen.seed( rd() );
+	}
+
+	/**
+	 * Reseeds generator with `seed`
+	*/
+	void reseed( uint64_t seed ) { gen.seed( seed ); }
+	/**
+	 * Uses random device to reseed generator
+	*/
+	void reseed() { gen.seed( rd() ); }
+
+	T rand() {
+		return distr( gen );
+	}
+
+	/**
+	 * Don't use this for generation of large amount of values !!
+	*/
+	template< typename... Args >
+	static T rand_once( Args... args ) {
+		std::random_device rd;
+		std::mt19937 gen( rd() );
+
+		return distribution( args... )( gen );
+	}
+
+private:
+	std::random_device rd;
+	generator gen;
+	distribution distr;
+};
+
+template< typename T >
+using GaussGenerator = RandUtils< T, std::normal_distribution<> >;
+
+template< typename T >
+using ExponentialGenerator = RandUtils< T, std::exponential_distribution<> >;
+
+#endif // PRNG_H