[shark] 51/58: rewrote DirectSearch to make more explicitely use of Rng. DirectSearch should now be thread safe when every instance of the algorithms gets its own Rng instead of the globalRng, which is default

Ghislain Vaillant ghisvail-guest at moszumanska.debian.org
Wed Mar 16 10:05:33 UTC 2016


This is an automated email from the git hooks/post-receive script.

ghisvail-guest pushed a commit to branch master
in repository shark.

commit a985b79c034e426f8d7d2338d2c7e0246f62e6a0
Author: Oswin Krause <oswin.krause at di.ku.dk>
Date:   Mon Feb 29 16:44:47 2016 +0100

    rewrote DirectSearch to make more explicitely use of Rng. DirectSearch should now be thread safe when every instance of the algorithms gets its own Rng instead of the globalRng, which is default
---
 .../Operators/Mutation/BitflipMutation.cpp         |  2 +-
 .../DirectSearch/Operators/Selection/Selection.cpp |  4 +--
 Test/Algorithms/Trainers/FisherLDA.cpp             |  2 +-
 Test/Algorithms/Trainers/LDA.cpp                   |  8 ++---
 Test/Algorithms/Trainers/LinearRegression.cpp      |  2 +-
 Test/Algorithms/Trainers/Normalization.cpp         |  6 ++--
 Test/Algorithms/Trainers/PCA.cpp                   |  6 ++--
 Test/ObjectiveFunctions/ErrorFunction.cpp          |  2 +-
 Test/Rng/MultiNomial.cpp                           |  4 +--
 Test/Rng/MultiVariateNormal.cpp                    |  4 +--
 examples/EA/SOO/AckleyES.tpp                       |  4 +--
 examples/EA/SOO/TSP.tpp                            |  4 +--
 examples/Unsupervised/PCA.tpp                      |  2 +-
 include/shark/Algorithms/DirectSearch/CMA.h        |  3 +-
 .../Algorithms/DirectSearch/CMA/CMAIndividual.h    |  5 +--
 .../shark/Algorithms/DirectSearch/CMA/Chromosome.h |  7 ++--
 include/shark/Algorithms/DirectSearch/CMSA.h       | 11 ++++--
 include/shark/Algorithms/DirectSearch/ElitistCMA.h |  4 ++-
 include/shark/Algorithms/DirectSearch/LMCMA.h      |  9 ++---
 include/shark/Algorithms/DirectSearch/MOCMA.h      |  8 ++---
 .../Operators/Mutation/BitflipMutator.h            |  4 +--
 .../Operators/Mutation/PolynomialMutation.h        |  8 ++---
 .../Operators/Recombination/OnePointCrossover.h    |  6 ++--
 .../Recombination/SimulatedBinaryCrossover.h       | 10 +++---
 .../Operators/Recombination/UniformCrossover.h     |  8 ++---
 .../Operators/Selection/EPTournamentSelection.h    | 10 +++---
 .../Operators/Selection/LinearRanking.h            |  3 +-
 .../Operators/Selection/RouletteWheelSelection.h   |  4 +--
 .../Operators/Selection/TournamentSelection.h      |  9 ++---
 .../Operators/Selection/UniformRanking.h           |  3 +-
 .../Algorithms/DirectSearch/RealCodedNSGAII.h      | 16 +++++----
 include/shark/Algorithms/DirectSearch/SMS-EMOA.h   | 21 +++++------
 .../Algorithms/DirectSearch/SteadyStateMOCMA.h     | 10 +++---
 include/shark/Algorithms/DirectSearch/VDCMA.h      |  8 ++---
 include/shark/Data/DataDistribution.h              |  2 +-
 .../Distributions/MultiNomialDistribution.h        | 13 +++----
 .../Distributions/MultiVariateNormalDistribution.h | 41 ++++++++++------------
 src/Algorithms/DirectSearch/CMA.cpp                |  7 ++--
 src/Algorithms/DirectSearch/CMSA.cpp               |  4 +--
 src/Algorithms/DirectSearch/ElitistCMA.cpp         |  4 +--
 40 files changed, 153 insertions(+), 135 deletions(-)

diff --git a/Test/Algorithms/DirectSearch/Operators/Mutation/BitflipMutation.cpp b/Test/Algorithms/DirectSearch/Operators/Mutation/BitflipMutation.cpp
index 27219bd..5c9fd06 100644
--- a/Test/Algorithms/DirectSearch/Operators/Mutation/BitflipMutation.cpp
+++ b/Test/Algorithms/DirectSearch/Operators/Mutation/BitflipMutation.cpp
@@ -21,7 +21,7 @@ BOOST_AUTO_TEST_CASE( BitflipMutation ) {
 	ind1.searchPoint() = std::vector< bool >( n, false );
 	ind2.searchPoint() = ind1.searchPoint();
 	for(std::size_t t = 0; t != trials; ++t){
-		flip(ind2);
+		flip(shark::Rng::globalRng, ind2);
 		for(std::size_t i = 0; i != n; ++i){
 			flipCount[i] += (ind1.searchPoint()[i] == ind2.searchPoint()[i])? 0:1;
 		}
diff --git a/Test/Algorithms/DirectSearch/Operators/Selection/Selection.cpp b/Test/Algorithms/DirectSearch/Operators/Selection/Selection.cpp
index d275a19..9bf82fb 100644
--- a/Test/Algorithms/DirectSearch/Operators/Selection/Selection.cpp
+++ b/Test/Algorithms/DirectSearch/Operators/Selection/Selection.cpp
@@ -36,7 +36,7 @@ BOOST_AUTO_TEST_CASE( Tournament_Selection ) {
 	std::vector<unsigned int> counter(11,0);
 
 	for( unsigned int i = 0; i < 10000; i++ ) {
-		counter[ std::distance( pop, ts( pop, pop + 11 ) ) ]++;
+		counter[ std::distance( pop, ts(shark::Rng::globalRng, pop, pop + 11 ) ) ]++;
 	}
 
 	std::copy( counter.begin(), counter.begin() + 11, std::ostream_iterator<double>( std::cout, "," ) );
@@ -57,7 +57,7 @@ BOOST_AUTO_TEST_CASE( RouletteWheel_Selection ) {
 	
 	shark::RealVector hist(11,0);
 	for( unsigned int i = 0; i < 1000000; i++ )
-		hist[ *ts( pop, pop + 11,prob) ]+=1.0;
+		hist[ *ts(shark::Rng::globalRng, pop, pop + 11,prob) ]+=1.0;
 	
 	hist /=1000000;
 	for(std::size_t i = 0; i != 11; ++i){
diff --git a/Test/Algorithms/Trainers/FisherLDA.cpp b/Test/Algorithms/Trainers/FisherLDA.cpp
index c90115c..8fa9a5d 100644
--- a/Test/Algorithms/Trainers/FisherLDA.cpp
+++ b/Test/Algorithms/Trainers/FisherLDA.cpp
@@ -51,7 +51,7 @@ BOOST_AUTO_TEST_CASE( FISHER_LDA_TEST ){
 	for(size_t i=0;i!=trainExamples;++i) {
 		//create sample
 		target[i] = i % 3;
-		input[i] = dist().first + mean[target[i]];
+		input[i] = dist(Rng::globalRng).first + mean[target[i]];
 	}
 	//statisticalBayesRisk/=trainExamples;
 
diff --git a/Test/Algorithms/Trainers/LDA.cpp b/Test/Algorithms/Trainers/LDA.cpp
index f303482..8eb3799 100644
--- a/Test/Algorithms/Trainers/LDA.cpp
+++ b/Test/Algorithms/Trainers/LDA.cpp
@@ -44,7 +44,7 @@ BOOST_AUTO_TEST_CASE( LDA_TEST_TWOCLASS ){
 	for(size_t i=0;i!=TrainExamples;++i){
 		//create samples. class 1 has double as many elements as class 0
 		target[i]=(i%3 == 0);
-		input[i]=dist().first+mean[target[i]];
+		input[i]=dist(Rng::globalRng).first+mean[target[i]];
 		//calculate bayes Target - the best fit to the distributions
 		RealVector diff=input[i]-mean[0];
 		double dist0=inner_prod(diff,prod(inverse,diff))+prior[0];
@@ -108,7 +108,7 @@ BOOST_AUTO_TEST_CASE( LDA_TEST_TWOCLASS_SINGULAR ){
 	for(size_t i=0;i!=TrainExamples;++i){
 		//create sample
 		target[i]= (i%3 != 0);
-		RealVector vec = dist().first+mean[target[i]];
+		RealVector vec = dist(Rng::globalRng).first+mean[target[i]];
 		//calculate bayes Target - the best fit to the distributions
 		RealVector diff=vec-mean[0];
 		double dist0=inner_prod(diff,prod(inverse,diff)) + prior[0];
@@ -171,7 +171,7 @@ BOOST_AUTO_TEST_CASE( LDA_TEST_MULTICLASS ){
 	for(size_t i=0;i!=TrainExamples;++i){
 		//create sample
 		target[i]=i%classes;
-		input[i]=dist().first+mean[target[i]];
+		input[i]=dist(Rng::globalRng).first+mean[target[i]];
 		//calculate bayes Target - the best fit to the distributions
 		unsigned int bayesTarget = 0;
 		double minDist = 1.e30;
@@ -235,7 +235,7 @@ BOOST_AUTO_TEST_CASE( LDA_TEST_MULTICLASS_WEIGHTING ){
 	for(size_t i=0;i!=TrainExamples;++i){
 		//create sample
 		target[i]=i%classes;
-		input[i]=dist().first+mean[target[i]];
+		input[i]=dist(Rng::globalRng).first+mean[target[i]];
 	}
 	ClassificationDataset dataset = createLabeledDataFromRange(input,target);
 	
diff --git a/Test/Algorithms/Trainers/LinearRegression.cpp b/Test/Algorithms/Trainers/LinearRegression.cpp
index 3fcf147..2f950ed 100644
--- a/Test/Algorithms/Trainers/LinearRegression.cpp
+++ b/Test/Algorithms/Trainers/LinearRegression.cpp
@@ -44,7 +44,7 @@ BOOST_AUTO_TEST_CASE( LinearRegression_TEST ){
 		input[i](0) = uniform();
 		input[i](1) = uniform();
 		testTarget[i] =  model(input[i]);
-		trainTarget[i] = noise().first + testTarget[i];
+		trainTarget[i] = noise(Rng::globalRng).first + testTarget[i];
 	}
 
 	// let the model forget...
diff --git a/Test/Algorithms/Trainers/Normalization.cpp b/Test/Algorithms/Trainers/Normalization.cpp
index 8489a64..e5169b9 100644
--- a/Test/Algorithms/Trainers/Normalization.cpp
+++ b/Test/Algorithms/Trainers/Normalization.cpp
@@ -65,7 +65,7 @@ BOOST_AUTO_TEST_CASE( NORMALIZE_WHITENING)
 	
 	std::vector<RealVector> input(1000,RealVector(3));
 	for(std::size_t i = 0; i != 1000;++i)
-		input[i]=dist().first+mean;
+		input[i]=dist(Rng::globalRng).first+mean;
 
 	UnlabeledData<RealVector> set = createDataFromRange(input);
 	NormalizeComponentsWhitening normalizer(1.5);
@@ -105,7 +105,7 @@ BOOST_AUTO_TEST_CASE( NORMALIZE_WHITENING_RANK_2)
 	
 	std::vector<RealVector> input(1000,RealVector(3));
 	for(std::size_t i = 0; i != 1000;++i)
-		input[i]=dist().first+mean;
+		input[i]=dist(Rng::globalRng).first+mean;
 
 	UnlabeledData<RealVector> set = createDataFromRange(input);
 	NormalizeComponentsWhitening normalizer(1.5);
@@ -149,7 +149,7 @@ BOOST_AUTO_TEST_CASE( NORMALIZE_ZCA)
 	
 	std::vector<RealVector> input(1000,RealVector(3));
 	for(std::size_t i = 0; i != 1000;++i)
-		input[i]=dist().first+mean;
+		input[i]=dist(Rng::globalRng).first+mean;
 
 	UnlabeledData<RealVector> set = createDataFromRange(input);
 	NormalizeComponentsZCA normalizer(1.5);
diff --git a/Test/Algorithms/Trainers/PCA.cpp b/Test/Algorithms/Trainers/PCA.cpp
index c14da41..43edf51 100644
--- a/Test/Algorithms/Trainers/PCA.cpp
+++ b/Test/Algorithms/Trainers/PCA.cpp
@@ -85,7 +85,7 @@ UnlabeledData<RealVector> createData3D()
 	BOOST_FOREACH(RealVector& sample, data)
 	{
 		//first element is the sample, second is the underlying uniform gaussian
-		sample = mean + distribution().first;
+		sample = mean + distribution(Rng::globalRng).first;
 	}
 	return  createDataFromRange(data);
 }
@@ -104,7 +104,7 @@ UnlabeledData<RealVector> createData2D()
 	MultiVariateNormalDistribution distribution(C);
 
 	std::vector<RealVector> v;
-	for(unsigned i=0; i<numberOfExamples; i++) v.push_back(mu + distribution().first);
+	for(unsigned i=0; i<numberOfExamples; i++) v.push_back(mu + distribution(Rng::globalRng).first);
 	return  createDataFromRange(v);
 }
 
@@ -144,7 +144,7 @@ UnlabeledData<RealVector> createDataNotFullRank()
 	BOOST_FOREACH(RealVector& sample, data)
 	{
 		//first element is the sample, second is the underlying uniform gaussian
-		sample = mean + distribution().first;
+		sample = mean + distribution(Rng::globalRng).first;
 	}
 	return  createDataFromRange(data,2);//small batch size to get batching errors
 }
diff --git a/Test/ObjectiveFunctions/ErrorFunction.cpp b/Test/ObjectiveFunctions/ErrorFunction.cpp
index 1791339..d92f9af 100644
--- a/Test/ObjectiveFunctions/ErrorFunction.cpp
+++ b/Test/ObjectiveFunctions/ErrorFunction.cpp
@@ -129,7 +129,7 @@ BOOST_AUTO_TEST_CASE( ObjFunct_ErrorFunction_LinearRegression ){
 		input[i](0) = uniform();
 		input[i](1) = uniform();
 		testTarget[i] =  model(input[i]);
-		RealVector noiseVal = noise().first;
+		RealVector noiseVal = noise(Rng::globalRng).first;
 		trainTarget[i] = noiseVal + testTarget[i];
 		optimalMSE+=norm_sqr(noiseVal);
 	}
diff --git a/Test/Rng/MultiNomial.cpp b/Test/Rng/MultiNomial.cpp
index 38d517e..364a351 100644
--- a/Test/Rng/MultiNomial.cpp
+++ b/Test/Rng/MultiNomial.cpp
@@ -21,7 +21,7 @@ BOOST_AUTO_TEST_CASE( MultiNomial_same_probabilities) {
 	RealVector draws(Dimensions,0.0);
 	RealVector drawsDiscrete(Dimensions,0.0);
 	for(std::size_t s = 0; s != Samples; ++s){
-		MultiNomialDistribution::result_type sample = dist();
+		MultiNomialDistribution::result_type sample = dist(Rng::globalRng);
 		BOOST_REQUIRE(sample < Dimensions);
 		draws[sample]++;
 		drawsDiscrete[Rng::discrete(0,Dimensions-1)]++;
@@ -56,7 +56,7 @@ BOOST_AUTO_TEST_CASE( MultiNomial_different_probabilities) {
 		MultiNomialDistribution dist(probabilities);
 		RealVector draws(Dimensions,0.0);
 		for(std::size_t s = 0; s != Samples; ++s){
-			MultiNomialDistribution::result_type sample = dist();
+			MultiNomialDistribution::result_type sample = dist(Rng::globalRng);
 			BOOST_REQUIRE(sample < Dimensions);
 			draws[sample]++;
 		}
diff --git a/Test/Rng/MultiVariateNormal.cpp b/Test/Rng/MultiVariateNormal.cpp
index 97f4c3c..40d5803 100644
--- a/Test/Rng/MultiVariateNormal.cpp
+++ b/Test/Rng/MultiVariateNormal.cpp
@@ -30,7 +30,7 @@ BOOST_AUTO_TEST_CASE( MULTIVARIATENORMAL_EIGENVALUES ) {
 	Data<RealVector> normalSampleSet(Samples,RealVector(Dimensions));
 	
 	for(std::size_t i = 0; i != Samples; ++i){
-		MultiVariateNormalDistribution::result_type sample = dist();
+		MultiVariateNormalDistribution::result_type sample = dist(Rng::globalRng);
 		sampleSet.element(i) = sample.first;
 		normalSampleSet.element(i) = sample.second;
 	}
@@ -77,7 +77,7 @@ BOOST_AUTO_TEST_CASE( MULTIVARIATENORMAL_Cholesky) {
 	Data<RealVector> normalSampleSet(Samples,RealVector(Dimensions));
 	
 	for(std::size_t i = 0; i != Samples; ++i){
-		MultiVariateNormalDistribution::result_type sample = dist();
+		MultiVariateNormalDistribution::result_type sample = dist(Rng::globalRng);
 		sampleSet.element(i) = sample.first;
 		normalSampleSet.element(i) = sample.second;
 	}
diff --git a/examples/EA/SOO/AckleyES.tpp b/examples/EA/SOO/AckleyES.tpp
index bee36df..b4f9ccd 100644
--- a/examples/EA/SOO/AckleyES.tpp
+++ b/examples/EA/SOO/AckleyES.tpp
@@ -69,9 +69,9 @@ int main( int argc, char ** argv ) {
 			offspring[i].chromosome() *= Rng::logNormal( 0, tau0 + tau1 );
 			
 			// Recombine search points
-			offspring[i].searchPoint() = uniform( mom->searchPoint(), dad->searchPoint() );
+			offspring[i].searchPoint() = uniform(Rng::globalRng, mom->searchPoint(), dad->searchPoint() );
 			// Mutate search point
-			offspring[i].searchPoint() = offspring[i].chromosome() * mutationDistribution().first;
+			offspring[i].searchPoint() = offspring[i].chromosome() * mutationDistribution(Rng::globalRng).first;
 	
 			// Assign fitness
 			offspring[i].unpenalizedFitness() = ackley.eval( offspring[i].searchPoint() );
diff --git a/examples/EA/SOO/TSP.tpp b/examples/EA/SOO/TSP.tpp
index 428bf09..8e00917 100755
--- a/examples/EA/SOO/TSP.tpp
+++ b/examples/EA/SOO/TSP.tpp
@@ -286,8 +286,8 @@ int main( int argc, char ** argv ) {
 			shark::IndividualType, 
 			shark::IndividualType 
 			> result = pmc( 
-				*rws( parents.begin(), parents.end(), selectionProbabilities ),
-				*rws( parents.begin(), parents.end(), selectionProbabilities )
+				*rws( shark::Rng::globalRng, parents.begin(), parents.end(), selectionProbabilities ),
+				*rws( shark::Rng::globalRng, parents.begin(), parents.end(), selectionProbabilities )
 			);
 			offspring[ i ] = result.first;
 			offspring[ i + 1 ] = result.second;
diff --git a/examples/Unsupervised/PCA.tpp b/examples/Unsupervised/PCA.tpp
index 9239c83..62a960a 100644
--- a/examples/Unsupervised/PCA.tpp
+++ b/examples/Unsupervised/PCA.tpp
@@ -53,7 +53,7 @@ UnlabeledData<RealVector> createData()
 	BOOST_FOREACH(RealVector& sample, data)
 	{
 		//first element is the sample, second is the underlying uniform gaussian
-		sample = mean + distribution().first;
+		sample = mean + distribution(Rng::globalRng).first;
 	}
 	return createDataFromRange(data);
 }
diff --git a/include/shark/Algorithms/DirectSearch/CMA.h b/include/shark/Algorithms/DirectSearch/CMA.h
index 4beaed8..b933f28 100644
--- a/include/shark/Algorithms/DirectSearch/CMA.h
+++ b/include/shark/Algorithms/DirectSearch/CMA.h
@@ -69,7 +69,7 @@ public:
 	/**
 	* \brief Default c'tor.
 	*/
-	SHARK_EXPORT_SYMBOL CMA();
+	SHARK_EXPORT_SYMBOL CMA(DefaultRngType& rng = Rng::globalRng);
 
 	/// \brief From INameable: return the class name.
 	std::string name() const
@@ -282,6 +282,7 @@ private:
 	std::size_t m_counter; ///< counter for generations
 
 	MultiVariateNormalDistribution m_mutationDistribution;
+	DefaultRngType* mpe_rng;
 };
 }
 
diff --git a/include/shark/Algorithms/DirectSearch/CMA/CMAIndividual.h b/include/shark/Algorithms/DirectSearch/CMA/CMAIndividual.h
index 448843c..4db0d96 100644
--- a/include/shark/Algorithms/DirectSearch/CMA/CMAIndividual.h
+++ b/include/shark/Algorithms/DirectSearch/CMA/CMAIndividual.h
@@ -64,9 +64,10 @@ public:
 	void updateAsOffspring(){
 		chromosome().updateAsOffspring();
 	}
-	void mutate(){
+	template<class RngType>
+	void mutate(RngType& rng){
 		chromosome().m_mutationDistribution.generate(
-			chromosome().m_lastStep,chromosome().m_lastZ
+			rng, chromosome().m_lastStep,chromosome().m_lastZ
 		);
 		noalias(searchPoint()) += chromosome().m_stepSize * chromosome().m_lastStep;
 	}
diff --git a/include/shark/Algorithms/DirectSearch/CMA/Chromosome.h b/include/shark/Algorithms/DirectSearch/CMA/Chromosome.h
index b6ef8d7..bc2a19c 100644
--- a/include/shark/Algorithms/DirectSearch/CMA/Chromosome.h
+++ b/include/shark/Algorithms/DirectSearch/CMA/Chromosome.h
@@ -46,7 +46,6 @@ struct CMAChromosome{
 		Failure = 3
 	};
 	MultiVariateNormalDistributionCholesky m_mutationDistribution; ///< Models the search distribution using a cholsky matrix
-	//~ RealMatrix m_inverseCholesky;///< inverse cholesky matrix
 
 	RealVector m_evolutionPath; ///< Low-pass filtered accumulation of successful mutative steps.
 	RealVector m_lastStep; ///< The most recent mutative step.
@@ -62,20 +61,18 @@ struct CMAChromosome{
 	double m_covarianceMatrixUnlearningRate; ///< Learning rate (constant) for unlearning unsuccessful directions from the covariance matrix.
 
 	double m_successThreshold; ///< Success threshold \f$p_{\text{thresh}}\f$ for cutting off evolution path updates.
-
+	
 	CMAChromosome(){}
 	CMAChromosome(
 		std::size_t searchSpaceDimension,
 		double successThreshold,
 		double initialStepSize
 	)
-	: m_mutationDistribution(true)//we do a triangular cholesky factorisation
-	, m_stepSize( initialStepSize )
+	: m_stepSize( initialStepSize )
 	, m_covarianceMatrixLearningRate( 0 )
 	, m_successThreshold(successThreshold)
 	{
 		m_mutationDistribution.resize( searchSpaceDimension );
-		//~ m_inverseCholesky = blas::identity_matrix<double>( searchSpaceDimension );
 		m_evolutionPath.resize( searchSpaceDimension );
 		m_lastStep.resize( searchSpaceDimension );
 		m_lastZ.resize( searchSpaceDimension );
diff --git a/include/shark/Algorithms/DirectSearch/CMSA.h b/include/shark/Algorithms/DirectSearch/CMSA.h
index 60ccf8d..5b1cad0 100644
--- a/include/shark/Algorithms/DirectSearch/CMSA.h
+++ b/include/shark/Algorithms/DirectSearch/CMSA.h
@@ -72,7 +72,13 @@ class CMSA : public AbstractSingleObjectiveOptimizer<RealVector > {
 public:
 
 	/// \brief Default c'tor.
-	CMSA() : m_mu( 100 ), m_lambda( 200 ), m_userSetMu(false),m_userSetLambda(false), m_initSigma(-1) {
+	CMSA(DefaultRngType& rng = Rng::globalRng)
+	: m_mu( 100 )
+	, m_lambda( 200 )
+	, m_userSetMu(false)
+	,m_userSetLambda(false)
+	, m_initSigma(-1)
+	, mpe_rng(&rng){
 		m_features |= REQUIRES_VALUE;
 	}
 
@@ -174,7 +180,8 @@ private:
 
 	RealVector m_mean; ///< The current cog of the population.
 
-	shark::MultiVariateNormalDistribution m_mutationDistribution; ///< Multi-variate normal mutation distribution.   
+	MultiVariateNormalDistribution m_mutationDistribution; ///< Multi-variate normal mutation distribution.   
+	DefaultRngType* mpe_rng;
 };
 }
 
diff --git a/include/shark/Algorithms/DirectSearch/ElitistCMA.h b/include/shark/Algorithms/DirectSearch/ElitistCMA.h
index 101692b..4c54269 100644
--- a/include/shark/Algorithms/DirectSearch/ElitistCMA.h
+++ b/include/shark/Algorithms/DirectSearch/ElitistCMA.h
@@ -69,7 +69,7 @@ namespace shark {
 class ElitistCMA : public AbstractSingleObjectiveOptimizer<RealVector >{	    
 public:
 
-	SHARK_EXPORT_SYMBOL ElitistCMA();
+	SHARK_EXPORT_SYMBOL ElitistCMA(DefaultRngType& rng = Rng::globalRng);
 
 	/// \brief From INameable: return the class name.
 	std::string name() const
@@ -125,6 +125,8 @@ private:
 	PenalizingEvaluator m_evaluator;///< evaluates the fitness of the individual and handles constraints
 	std::vector<double> m_ancestralFitness; ///< stores the last k fitness values (by default 5).
 	bool m_activeUpdate;///< Should bad individuals be actively purged from the strategy?
+
+	DefaultRngType* mpe_rng;///< the internal random number generator
 };
 }
 
diff --git a/include/shark/Algorithms/DirectSearch/LMCMA.h b/include/shark/Algorithms/DirectSearch/LMCMA.h
index 7ca25c5..46f0ed8 100644
--- a/include/shark/Algorithms/DirectSearch/LMCMA.h
+++ b/include/shark/Algorithms/DirectSearch/LMCMA.h
@@ -200,7 +200,7 @@ class LMCMA: public AbstractSingleObjectiveOptimizer<RealVector >
 {
 public:
 	/// \brief Default c'tor.
-	LMCMA(){
+	LMCMA(DefaultRngType& rng = Rng::globalRng):mpe_rng(&rng){
 		m_features |= REQUIRES_VALUE;
 	}
 	
@@ -211,7 +211,7 @@ public:
 
 	/// \brief Calculates lambda for the supplied dimensionality n.
 	unsigned suggestLambda( unsigned int dimension ) {
-		unsigned lambda = unsigned( 4. + ::floor( 3. * ::log( static_cast<double>( dimension ) ) ) ); // eq. (44)
+		unsigned lambda = unsigned( 4. + ::floor( 3. * ::log( static_cast<double>( dimension ) ) ) );
 		// heuristic for small search spaces
 		lambda = std::max<unsigned int>( 5, std::min( lambda, dimension ) );
 		return( lambda );
@@ -219,7 +219,7 @@ public:
 
 	/// \brief Calculates mu for the supplied lambda and the recombination strategy.
 	double suggestMu( unsigned int lambda) {
-		return lambda / 2.; // eq. (44)
+		return lambda / 2.;
 	}
 
 	using AbstractSingleObjectiveOptimizer<RealVector >::init;
@@ -369,7 +369,7 @@ private:
 		x.resize(m_numberOfVariables);
 		z.resize(m_numberOfVariables);
 		for(std::size_t i = 0; i != m_numberOfVariables; ++i){
-			z(i) = Rng::gauss(0,1);
+			z(i) = gauss(*mpe_rng,0,1);
 		}
 		m_A.prod(x,z);
 		noalias(x) = sigma()*x +m_mean;
@@ -391,6 +391,7 @@ private:
 
 	RealVector m_evolutionPathC;///< evolution path
 	
+	DefaultRngType* mpe_rng;
 	
 };
 
diff --git a/include/shark/Algorithms/DirectSearch/MOCMA.h b/include/shark/Algorithms/DirectSearch/MOCMA.h
index 579b661..2ba5eef 100644
--- a/include/shark/Algorithms/DirectSearch/MOCMA.h
+++ b/include/shark/Algorithms/DirectSearch/MOCMA.h
@@ -63,7 +63,7 @@ public:
 	};
 public:
 	
-	IndicatorBasedMOCMA() {
+	IndicatorBasedMOCMA(DefaultRngType& rng= Rng::globalRng):mpe_rng(&rng) {
 		m_individualSuccessThreshold = 0.44;
 		initialSigma() = 1.0;
 		mu() = 100;
@@ -191,7 +191,7 @@ protected:
 		}
 		//copy points randomly
 		for(std::size_t i = numPoints; i != mu; ++i){
-			std::size_t index = Rng::discrete(0,initialSearchPoints.size()-1);
+			std::size_t index = discrete(*mpe_rng,0,startingPoints.size()-1);
 			m_parents[i] = IndividualType(noVariables,m_individualSuccessThreshold,m_initialSigma);
 			m_parents[i].searchPoint() = initialSearchPoints[index];
 			m_parents[i].penalizedFitness() = functionValues[index];
@@ -209,7 +209,7 @@ protected:
 		for(std::size_t i = 0; i != mu(); ++i){
 			std::size_t parentId = i;
 			offspring[i] = m_parents[parentId];
-			offspring[i].mutate();
+			offspring[i].mutate(*mpe_rng);
 			offspring[i].parent() = parentId;
 		}
 		return offspring;
@@ -257,7 +257,7 @@ private:
 	NotionOfSuccess m_notionOfSuccess; ///< Flag for deciding whether the improved step-size adaptation shall be used.
 	double m_individualSuccessThreshold;
 	double m_initialSigma;
-
+	DefaultRngType* mpe_rng;
 };
 
 typedef IndicatorBasedMOCMA< HypervolumeIndicator > MOCMA;
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Mutation/BitflipMutator.h b/include/shark/Algorithms/DirectSearch/Operators/Mutation/BitflipMutator.h
index 9bcb1a8..0e9532a 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Mutation/BitflipMutator.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Mutation/BitflipMutator.h
@@ -58,10 +58,10 @@ struct BitflipMutator {
 	/// 
 	/// \param [in,out] ind Individual to be mutated.
 	template<typename IndividualType>
-	void operator()(IndividualType &ind) {
+	void operator()(DefaultRngType& rng, IndividualType &ind) {
 
 		for (unsigned int i = 0; i < ind.searchPoint().size(); i++) {
-			if (Rng::coinToss(m_mutationStrength)) {
+			if (coinToss(rng,m_mutationStrength)) {
 				ind.searchPoint()[ i ] = !ind.searchPoint()[ i ];
 			}
 		}
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Mutation/PolynomialMutation.h b/include/shark/Algorithms/DirectSearch/Operators/Mutation/PolynomialMutation.h
index 051ae20..31fa934 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Mutation/PolynomialMutation.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Mutation/PolynomialMutation.h
@@ -57,14 +57,14 @@ namespace shark {
 		///  for accessing the actual search point.
 		/// \param [in,out] ind Individual to be mutated.
 		template<typename IndividualType>
-		void operator()( IndividualType & ind )const{
+		void operator()(DefaultRngType& rng, IndividualType & ind )const{
 			RealVector& point = ind.searchPoint();
            
 			for( unsigned int i = 0; i < point.size(); i++ ) {
 
-				if( Rng::coinToss( m_prob ) ) {
+				if( coinToss(rng, m_prob ) ) {
 					if( point[i] < m_lower( i ) || point[i] > m_upper( i ) ) { 
-						point[i] = Rng::uni(m_lower(i),m_upper(i));
+						point[i] = uni(rng,m_lower(i),m_upper(i));
 					} else {
 						// Calculate normalized distance from boundaries
 						double delta1 = (m_upper( i ) - point[i]) / (m_upper( i ) - m_lower( i ));
@@ -72,7 +72,7 @@ namespace shark {
 						
 						//compute change in delta
 						double deltaQ=0;
-						double u = Rng::uni(0,1);
+						double u = uni(rng,0,1);
 						if( u <= .5 ) {
 							double delta = std::pow(delta1 , m_nm + 1.);
 							deltaQ =  2.0 * u + (1.0 - 2.0 * u) * delta;
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Recombination/OnePointCrossover.h b/include/shark/Algorithms/DirectSearch/Operators/Recombination/OnePointCrossover.h
index ef949e0..549eef8 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Recombination/OnePointCrossover.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Recombination/OnePointCrossover.h
@@ -42,10 +42,10 @@ namespace shark {
 /// right parent.
 struct OnePointCrossover {
 	/// \brief Performs the one-point crossover
-	template<typename PointType>
-	PointType operator()( const PointType & mom, const PointType & dad ) {
+	template<class RngType, typename PointType>
+	PointType operator()(RngType& rng, const PointType & mom, const PointType & dad ) {
 		SIZE_CHECK(mom.size() == dad.size());
-		std::size_t point = Rng::discrete( 0, mom.size() - 1 );
+		std::size_t point = discrete(rng, 0, mom.size() - 1 );
 	    
 		PointType offspring( mom.size() );
 		std::copy( mom.begin(), mom.begin() + point, offspring.begin() );
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Recombination/SimulatedBinaryCrossover.h b/include/shark/Algorithms/DirectSearch/Operators/Recombination/SimulatedBinaryCrossover.h
index befd3b0..d1a3a60 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Recombination/SimulatedBinaryCrossover.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Recombination/SimulatedBinaryCrossover.h
@@ -57,14 +57,14 @@ namespace shark {
 		/// 
 		/// \param [in,out] i1 Individual to be mated.
 		/// \param [in,out] i2 Individual to be mated.
-		template<typename IndividualType>
-		void operator()( IndividualType & i1, IndividualType & i2 )const{	
+		template<class RngType, typename IndividualType>
+		void operator()(RngType& rng, IndividualType & i1, IndividualType & i2 )const{	
 			RealVector& point1 = i1.searchPoint();
 			RealVector& point2 = i2.searchPoint();
 
 			for( unsigned int i = 0; i < point1.size(); i++ ) {
 
-				if( !Rng::coinToss( m_prob ) )
+				if( !coinToss(rng, m_prob ) )
 					continue;
 				
 				double y1 = 0;
@@ -89,7 +89,7 @@ namespace shark {
 				double alpha1 = 2. - std::pow(beta1 , -expp);
 				double alpha2 = 2. - std::pow(beta2 , -expp);
 
-				double u = Rng::uni( 0., 1. );
+				double u = uni(rng, 0., 1. );
 				alpha1 *=u;
 				alpha2 *=u;
 				if( u > 1. / alpha1 ) {
@@ -105,7 +105,7 @@ namespace shark {
 				point1[i] = 0.5 * ((y1 + y2) - betaQ1 * (y2 - y1));
 				point2[i] = 0.5 * ((y1 + y2) + betaQ2 * (y2 - y1));
 				// randomly swap loci
-				if( Rng::coinToss() ) std::swap(point1[i], point2[i]);
+				if( coinToss(rng,0.5) ) std::swap(point1[i], point2[i]);
 
 
 				//  -> from Deb's implementation, not contained in any paper
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Recombination/UniformCrossover.h b/include/shark/Algorithms/DirectSearch/Operators/Recombination/UniformCrossover.h
index ebc7f9e..59a0f6d 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Recombination/UniformCrossover.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Recombination/UniformCrossover.h
@@ -47,19 +47,19 @@ public:
 	/// \brief Default c'tor, initializes the per element probability.
 	/// 
 	/// \param [in] mixingRatio Mixing ratio between parent individuals.
-	UniformCrossover( double mixingRatio = 0.5 ){
+	UniformCrossover(double mixingRatio = 0.5 ){
 		setMixingRatio(mixingRatio);
 	}
 
 	/// \brief Executes the uniform crossover.
 	///	
 	/// \return The offspring individual.
-	template<typename Point>
-	Point operator()( const Point & mom, const Point & dad ) const {
+	template<class RngType, typename Point>
+	Point operator()(RngType& rng, const Point & mom, const Point & dad ) const {
 		Point result( mom );
 
 		for( std::size_t i = 0; i < std::min( mom.size(), dad.size() ); i++ ) {
-			if( Rng::coinToss( m_mixingRatio ) )
+			if( coinToss(rng, m_mixingRatio ) )
 				result( i ) = dad( i );
 		}
 
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Selection/EPTournamentSelection.h b/include/shark/Algorithms/DirectSearch/Operators/Selection/EPTournamentSelection.h
index d8ff603..d9a938d 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Selection/EPTournamentSelection.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Selection/EPTournamentSelection.h
@@ -52,11 +52,12 @@ struct EPTournamentSelection {
 	/// \param [in] outE Iterator pointing to the first invalid element of the output range.
 	template<typename InIterator,typename OutIterator>
 	void operator()(
+		DefaultRngType& rng,
 		InIterator it, InIterator itE,
 		OutIterator out,  OutIterator outE
 	){
 		std::size_t outputSize = std::distance( out, outE );
-		std::vector<KeyValuePair<int, InIterator> > results = performTournament(it, itE);
+		std::vector<KeyValuePair<int, InIterator> > results = performTournament(rng, it, itE);
 		if(results.size() < outputSize){
 			throw SHARKEXCEPTION("[EPTournamentSelection] Input range must be bigger than output range");
 		}
@@ -74,11 +75,12 @@ struct EPTournamentSelection {
 	/// \param [in] mu number of individuals to select
 	template<typename Population>
 	void operator()(
+		DefaultRngType& rng,
 		Population& population,std::size_t mu
 	){
 		SIZE_CHECK(population.size() >= mu);
 		typedef typename Population::iterator InIterator;
-		std::vector<KeyValuePair<int, InIterator> > results = performTournament(population.begin(),population.end());
+		std::vector<KeyValuePair<int, InIterator> > results = performTournament(rng, population.begin(),population.end());
 		
 		
 		for(std::size_t i = 0; i != mu; ++i){
@@ -95,7 +97,7 @@ private:
 	///Returns a sorted range of pairs indicating, how often every individual won.
 	/// The best individuals are in the front of the range.
 	template<class InIterator>
-	std::vector<KeyValuePair<int, InIterator> > performTournament(InIterator it, InIterator itE){
+	std::vector<KeyValuePair<int, InIterator> > performTournament(DefaultRngType& rng, InIterator it, InIterator itE){
 		std::size_t size = std::distance( it, itE );
 		UIntVector selectionProbability(size,0.0);
 		std::vector<KeyValuePair<int, InIterator> > individualPerformance(size);
@@ -103,7 +105,7 @@ private:
 		for( std::size_t i = 0; i != size(); ++i ) {
 			individualPerformance[i].value = it+i;
 			for( std::size_t round = 0; round < tournamentSize; round++ ) {
-				std::size_t idx = shark::Rng::discrete( 0,size-1 );
+				std::size_t idx = discrete(rng, 0,size-1 );
 				if(e(*it) < e(*(it+idx)){
 					individualPerformance[i].key -= 1;
 				}
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Selection/LinearRanking.h b/include/shark/Algorithms/DirectSearch/Operators/Selection/LinearRanking.h
index ea2b02f..fdf8afc 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Selection/LinearRanking.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Selection/LinearRanking.h
@@ -73,6 +73,7 @@ struct LinearRankingSelection {
 	///
 	template<typename InIterator,typename OutIterator> 
 	void operator()( 
+		DefaultRngType& rng,
 		InIterator individuals,
 		InIterator individualsE,
 		OutIterator out,
@@ -99,7 +100,7 @@ struct LinearRankingSelection {
 
 		RouletteWheelSelection rws;
 		for( ; out != outE; ++out ){
-			InIterator individuals = rws( individualsPerformance.begin(), individualsPerformance.end(), selectionProbability)->value;
+			InIterator individuals = rws(rng, individualsPerformance.begin(), individualsPerformance.end(), selectionProbability)->value;
 			*out = *individuals;
 		}
 	}
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Selection/RouletteWheelSelection.h b/include/shark/Algorithms/DirectSearch/Operators/Selection/RouletteWheelSelection.h
index 058cb51..3101332 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Selection/RouletteWheelSelection.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Selection/RouletteWheelSelection.h
@@ -54,10 +54,10 @@ struct RouletteWheelSelection {
 	* \param [in] probabilities selection probabilities of the individuals
 	*/
 	template<typename Iterator>
-	Iterator operator()(Iterator it, Iterator itE, RealVector const& probabilities) const
+	Iterator operator()(DefaultRngType& rng, Iterator it, Iterator itE, RealVector const& probabilities) const
 	{
 		std::size_t n = probabilities.size();
-		double rnd = Rng::uni(0,1);
+		double rnd = uni(rng, 0,1);
 		double sum = 0;
 		for(std::size_t pos = 0; pos != n; ++pos,++it){
 			sum += probabilities(pos);
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Selection/TournamentSelection.h b/include/shark/Algorithms/DirectSearch/Operators/Selection/TournamentSelection.h
index bfc741e..882a711 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Selection/TournamentSelection.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Selection/TournamentSelection.h
@@ -54,13 +54,14 @@ struct TournamentSelection {
 	
 	template<typename IteratorType1, typename IteratorType2>
 	void operator()(
+		DefaultRngType& rng,
 		IteratorType1 inIt,
 		IteratorType1 inItE,
 		IteratorType2 outIt,
 		IteratorType2 outItE
 	){
 		for(; outIt != outItE; ++outIt ) {
-			*outIt = *(*this)(inIt,inItE);
+			*outIt = *(*this)(rng, inIt,inItE);
 		}
 	}
 	
@@ -69,16 +70,16 @@ struct TournamentSelection {
 	/// \param [in] itE Iterator pointing to the first invalid element.
 	/// \return An iterator pointing to the selected individual.
 	template< typename Iterator>
-	Iterator operator()( Iterator it, Iterator itE) const
+	Iterator operator()(DefaultRngType& rng, Iterator it, Iterator itE) const
 	{
 		std::size_t n = std::distance( it, itE );
 		SHARK_CHECK(tournamentSize > 0, " Tournament size k needs to be larger than 0");
 		SHARK_CHECK(n > tournamentSize, " Size of population needs to be larger than size of tournament");
 		
 		Predicate predicate;
-		Iterator result = it + Rng::discrete( 0, n-1 );
+		Iterator result = it + discrete(rng, 0, n-1 );
 		for( std::size_t i = 1; i < tournamentSize; i++ ) {
-			Iterator itt = it + Rng::discrete(0,n-1);
+			Iterator itt = it + discrete(rng, 0,n-1);
 			if( predicate(*itt, *result) ){
 				result = itt;
 			}
diff --git a/include/shark/Algorithms/DirectSearch/Operators/Selection/UniformRanking.h b/include/shark/Algorithms/DirectSearch/Operators/Selection/UniformRanking.h
index 3ad7f4b..25526fa 100644
--- a/include/shark/Algorithms/DirectSearch/Operators/Selection/UniformRanking.h
+++ b/include/shark/Algorithms/DirectSearch/Operators/Selection/UniformRanking.h
@@ -51,6 +51,7 @@ struct UniformRankingSelection {
 	/// \param [in] outE Iterator pointing to the first invalid element of the output range.
 	template<typename InIterator,typename OutIterator> 
 	void operator()( 
+		DefaultRngType& rng,
 		InIterator individuals,
 		InIterator individualsE,
 		OutIterator out,
@@ -61,7 +62,7 @@ struct UniformRankingSelection {
 		RealVector selectionProbability(size,1.0/size);
 		RouletteWheelSelection rws;
 		for( ; out != outE; ++out ){
-			*out = *rws( individuals, individualsE, selectionProbability);
+			*out = *rws(rng, individuals, individualsE, selectionProbability);
 		}
 	}
 
diff --git a/include/shark/Algorithms/DirectSearch/RealCodedNSGAII.h b/include/shark/Algorithms/DirectSearch/RealCodedNSGAII.h
index 61d675c..cebeaa4 100644
--- a/include/shark/Algorithms/DirectSearch/RealCodedNSGAII.h
+++ b/include/shark/Algorithms/DirectSearch/RealCodedNSGAII.h
@@ -63,7 +63,7 @@ public:
 	/**
 	* \brief Default c'tor.
 	*/
-	IndicatorBasedRealCodedNSGAII(){
+	IndicatorBasedRealCodedNSGAII(DefaultRngType& rng = Rng::globalRng):mpe_rng(&rng){
 		mu() = 100;
 		crossoverProbability() = 0.9;
 		nc() = 20.0;
@@ -214,8 +214,8 @@ protected:
 		}
 		//copy points randomly
 		for(std::size_t i = numPoints; i != mu; ++i){
-			std::size_t index = Rng::discrete(0,initialSearchPoints.size()-1);
-			m_parents[i].searchPoint() = initialSearchPoints[index];
+			std::size_t index = discrete(*mpe_rng, 0,startingPoints.size()-1);
+			m_parents[i].searchPoint() = startingPoints[index];
 			m_parents[i].penalizedFitness() = functionValues[index];
 			m_parents[i].unpenalizedFitness() = functionValues[index];
 		}
@@ -234,6 +234,7 @@ protected:
 		TournamentSelection< IndividualType::RankOrdering > selection;
 		std::vector<IndividualType> offspring(mu());
 		selection(
+			*mpe_rng,
 			m_parents.begin(), 
 			m_parents.end(),
 			offspring.begin(),
@@ -241,12 +242,12 @@ protected:
 		);
 
 		for( std::size_t i = 0; i < mu()-1; i+=2 ) {
-			if( Rng::coinToss( m_crossoverProbability ) ) {
-				m_crossover( offspring[i], offspring[i +1] );
+			if( coinToss(*mpe_rng, m_crossoverProbability ) ) {
+				m_crossover(*mpe_rng,offspring[i], offspring[i +1] );
 			}
 		}
 		for( std::size_t i = 0; i < mu(); i++ ) {
-			m_mutation( offspring[i] );
+			m_mutation(*mpe_rng, offspring[i] );
 		}
 		return offspring;
 	}
@@ -280,7 +281,8 @@ private:
 	SimulatedBinaryCrossover< RealVector > m_crossover; ///< Crossover operator.
 	PolynomialMutator m_mutation; ///< Mutation operator.
 
-	double m_crossoverProbability; ///< Crossover probability.	
+	double m_crossoverProbability; ///< Crossover probability.
+	DefaultRngType* mpe_rng; 
 };
 
 typedef IndicatorBasedRealCodedNSGAII< HypervolumeIndicator > RealCodedNSGAII;
diff --git a/include/shark/Algorithms/DirectSearch/SMS-EMOA.h b/include/shark/Algorithms/DirectSearch/SMS-EMOA.h
index b48014d..9df1865 100644
--- a/include/shark/Algorithms/DirectSearch/SMS-EMOA.h
+++ b/include/shark/Algorithms/DirectSearch/SMS-EMOA.h
@@ -61,7 +61,7 @@ namespace shark {
 */
 class SMSEMOA : public AbstractMultiObjectiveOptimizer<RealVector >{
 public:
-	SMSEMOA() {
+	SMSEMOA(DefaultRngType& rng = Rng::globalRng):mpe_rng(&rng) {
 		m_mu = 100;
 		m_mutator.m_nm = 20.0;
 		m_crossover.m_nc = 20.0;
@@ -193,8 +193,8 @@ protected:
 		}
 		//copy points randomly
 		for(std::size_t i = numPoints; i != mu; ++i){
-			std::size_t index = Rng::discrete(0,initialSearchPoints.size()-1);
-			m_parents[i].searchPoint() = initialSearchPoints[index];
+			std::size_t index = discrete(*mpe_rng,0,startingPoints.size()-1);
+			m_parents[i].searchPoint() = startingPoints[index];
 			m_parents[i].penalizedFitness() = functionValues[index];
 			m_parents[i].unpenalizedFitness() = functionValues[index];
 		}
@@ -242,18 +242,18 @@ private:
 	)const{
 		TournamentSelection< IndividualType::RankOrdering > selection;
 
-		IndividualType mate1( *selection( begin, end ) );
-		IndividualType mate2( *selection( begin, end) );
+		IndividualType mate1( *selection(*mpe_rng, begin, end ) );
+		IndividualType mate2( *selection(*mpe_rng, begin, end) );
 
-		if( Rng::coinToss( m_crossoverProbability ) ) {
-			m_crossover( mate1, mate2 );
+		if( coinToss(*mpe_rng, m_crossoverProbability ) ) {
+			m_crossover(*mpe_rng, mate1, mate2 );
 		}
 
-		if( Rng::coinToss() ) {
-			m_mutator( mate1 );
+		if( coinToss(*mpe_rng,0.5) ) {
+			m_mutator(*mpe_rng, mate1 );
 			return mate1;
 		} else {
-			m_mutator( mate2 );
+			m_mutator(*mpe_rng, mate2 );
 			return mate2;
 		}
 	}
@@ -265,6 +265,7 @@ private:
 	PolynomialMutator m_mutator; ///< Mutation operator.
 
 	double m_crossoverProbability; ///< Crossover probability.
+	DefaultRngType* mpe_rng;
 };
 }
 
diff --git a/include/shark/Algorithms/DirectSearch/SteadyStateMOCMA.h b/include/shark/Algorithms/DirectSearch/SteadyStateMOCMA.h
index 5fc5931..450f121 100644
--- a/include/shark/Algorithms/DirectSearch/SteadyStateMOCMA.h
+++ b/include/shark/Algorithms/DirectSearch/SteadyStateMOCMA.h
@@ -59,7 +59,7 @@ public:
 	};
 public:
 	
-	IndicatorBasedSteadyStateMOCMA() {
+	IndicatorBasedSteadyStateMOCMA(DefaultRngType& rng = Rng::globalRng):mpe_rng(&rng){
 		m_individualSuccessThreshold = 0.44;
 		initialSigma() = 1.0;
 		mu() = 100;
@@ -189,7 +189,7 @@ protected:
 		}
 		//copy points randomly
 		for(std::size_t i = numPoints; i != mu; ++i){
-			std::size_t index = Rng::discrete(0,initialSearchPoints.size()-1);
+			std::size_t index = discrete(*mpe_rng, 0,startingPoints.size()-1);
 			m_parents[i] = IndividualType(noVariables,m_individualSuccessThreshold,m_initialSigma);
 			m_parents[i].searchPoint() = initialSearchPoints[index];
 			m_parents[i].penalizedFitness() = functionValues[index];
@@ -212,10 +212,10 @@ protected:
 				break;
 		}
 		//sample a random parent with rank 1
-		std::size_t parentId = Rng::discrete(0, maxIdx-1);
+		std::size_t parentId = discrete(*mpe_rng, 0, maxIdx-1);
 		std::vector<IndividualType> offspring;
 		offspring.push_back(m_parents[parentId]);
-		offspring[0].mutate();
+		offspring[0].mutate(*mpe_rng);
 		offspring[0].parent() = parentId;
 		return offspring;
 	}
@@ -278,6 +278,8 @@ private:
 			}
 		}
 	}
+	
+	DefaultRngType* mpe_rng;
 };
 
 typedef IndicatorBasedSteadyStateMOCMA< HypervolumeIndicator > SteadyStateMOCMA;
diff --git a/include/shark/Algorithms/DirectSearch/VDCMA.h b/include/shark/Algorithms/DirectSearch/VDCMA.h
index b488c67..7eca1e6 100644
--- a/include/shark/Algorithms/DirectSearch/VDCMA.h
+++ b/include/shark/Algorithms/DirectSearch/VDCMA.h
@@ -58,7 +58,7 @@ private:
 public:
 
 	/// \brief Default c'tor.
-	VDCMA():m_initialSigma(0.0){
+	VDCMA(DefaultRngType& rng = Rng::globalRng):m_initialSigma(0.0), mpe_rng(&rng){
 		m_features |= REQUIRES_VALUE;
 	}
 	
@@ -113,7 +113,7 @@ public:
 		m_mean = blas::repeat(0.0,m_numberOfVariables);
 		m_vn.resize(m_numberOfVariables);
 		for(std::size_t i = 0; i != m_numberOfVariables;++i){
-			m_vn(i) = Rng::uni(0,1.0/m_numberOfVariables);
+			m_vn(i) = uni(*mpe_rng,0,1.0/m_numberOfVariables);
 		}
 		m_normv = norm_2(m_vn);
 		m_vn /= m_normv;
@@ -295,7 +295,7 @@ private:
 		x.resize(m_numberOfVariables);
 		y.resize(m_numberOfVariables);
 		for(std::size_t i = 0; i != m_numberOfVariables; ++i){
-			y(i) = Rng::gauss(0,1);
+			y(i) = gauss(*mpe_rng,0,1);
 		}
 		double a = std::sqrt(1+sqr(m_normv))-1;
 		a *= inner_prod(y,m_vn);
@@ -368,7 +368,7 @@ private:
 
 	unsigned m_counter; ///< counter for generations
 	
-	
+	DefaultRngType* mpe_rng;
 	
 	
 };
diff --git a/include/shark/Data/DataDistribution.h b/include/shark/Data/DataDistribution.h
index 586689c..d2888f8 100644
--- a/include/shark/Data/DataDistribution.h
+++ b/include/shark/Data/DataDistribution.h
@@ -420,7 +420,7 @@ public:
 	void draw(RealVector& input) const{
 		input.resize(m_offset.size());
 		noalias(input) = m_offset;
-		noalias(input) += m_dist().first;
+		noalias(input) += m_dist(Rng::globalRng).first;
 	}
 private:
 	MultiVariateNormalDistributionCholesky m_dist;
diff --git a/include/shark/Statistics/Distributions/MultiNomialDistribution.h b/include/shark/Statistics/Distributions/MultiNomialDistribution.h
index d8f421d..374719d 100644
--- a/include/shark/Statistics/Distributions/MultiNomialDistribution.h
+++ b/include/shark/Statistics/Distributions/MultiNomialDistribution.h
@@ -6,7 +6,7 @@
  * 
  *
  * \author    O.Krause
- * \date        2014
+ * \date        2016
  *
  *
  * \par Copyright 1995-2015 Shark Development Team
@@ -58,7 +58,7 @@ public:
 
 	/// \brief Constructor
 	/// \param [in] probabilities Probability vector
-	MultiNomialDistribution( RealVector const& probabilities ) 
+	MultiNomialDistribution(RealVector const& probabilities ) 
 	: m_probabilities(probabilities){
 		update();
 	}
@@ -87,12 +87,13 @@ public:
 	}
 
 	/// \brief Samples the distribution.
-	result_type operator()() const {
+	template<class RngType>
+	result_type operator()(RngType& rng) const {
 		std::size_t numStates = m_probabilities.size();
  
-		std::size_t index = Rng::discrete(0,numStates-1);
+		std::size_t index = discrete(rng,0,numStates-1);
  
-		if(Rng::coinToss(m_q[index]))
+		if(coinToss(rng, m_q[index]))
 			return index;
 		else
 			return m_J[index];
@@ -136,7 +137,7 @@ public:
 		for(std::size_t i = 0; i != larger.size(); ++i){
 			m_q[larger[i]]=std::min(m_q[larger[i]],1.0);
 		}
-	}			
+	}
 
 private:
 	RealVector m_probabilities; ///< probability of every state.
diff --git a/include/shark/Statistics/Distributions/MultiVariateNormalDistribution.h b/include/shark/Statistics/Distributions/MultiVariateNormalDistribution.h
index 85bca1a..5ce42f3 100644
--- a/include/shark/Statistics/Distributions/MultiVariateNormalDistribution.h
+++ b/include/shark/Statistics/Distributions/MultiVariateNormalDistribution.h
@@ -5,8 +5,8 @@
  * 
  * 
  *
- * \author      T.Voss
- * \date        2010
+ * \author      T.Voss, O.Krause
+ * \date        2016
  *
  *
  * \par Copyright 1995-2015 Shark Development Team
@@ -50,11 +50,12 @@ public:
 
 	/// \brief Constructor
 	/// \param [in] Sigma covariance matrix
-	MultiVariateNormalDistribution( RealMatrix const& Sigma ) {
+	MultiVariateNormalDistribution(RealMatrix const& Sigma ) {
 		m_covarianceMatrix = Sigma;
 		update();
 	}
 	
+	/// \brief Constructor
 	MultiVariateNormalDistribution(){} 
 	
 	/// \brief Stores/Restores the distribution from the supplied archive.
@@ -110,14 +111,13 @@ public:
 	}
 
 	/// \brief Samples the distribution.
-	result_type operator()() const {
+	template<class RngType>
+	result_type operator()(RngType& rng) const {
 		RealVector result( m_eigenValues.size(), 0. );
 		RealVector z( m_eigenValues.size() );
 		
-		SHARK_CRITICAL_REGION{
 		for( std::size_t i = 0; i < result.size(); i++ ) {
-			z( i ) = Rng::gauss( 0., 1. );
-		}
+			z( i ) = gauss(rng, 0., 1. );
 		}
 		for( std::size_t i = 0; i < result.size(); i++ )
 			for( std::size_t j = 0; j < result.size(); j++ )
@@ -129,7 +129,7 @@ public:
 	/// \brief Calculates the evd of the current covariance matrix.
 	void update() {
 		eigensymm( m_covarianceMatrix, m_eigenVectors, m_eigenValues );
-	}			
+	}
 
 private:
 	RealMatrix m_covarianceMatrix; ///< Covariance matrix of the mutation distribution.
@@ -148,14 +148,13 @@ public:
 	typedef std::pair<RealVector,RealVector> result_type;
 
 	/// \brief Constructor
+	/// \param [in] rng the random number generator
 	/// \param [in] covariance covariance matrix
-	/// \param triangular Is the choleksy factor triangular?
-	MultiVariateNormalDistributionCholesky( RealMatrix const& covariance, bool triangular=false ) 
-	:m_triangular(false){
+	MultiVariateNormalDistributionCholesky( RealMatrix const& covariance){
 		setCovarianceMatrix(covariance);
 	}
 	
-	MultiVariateNormalDistributionCholesky(bool triangular=false):m_triangular(triangular){} 
+	MultiVariateNormalDistributionCholesky(){} 
 	
 	/// \brief Stores/Restores the distribution from the supplied archive.
 	///\param [in,out] ar The archive to read from/write to.
@@ -196,17 +195,15 @@ public:
 	}
 	
 	
-	template<class Vector1, class Vector2>
-	void generate(Vector1& y, Vector2& z)const{
+	template<class RngType, class Vector1, class Vector2>
+	void generate(RngType& rng, Vector1& y, Vector2& z)const{
 		z.resize(size());
 		y.resize(size());
-		SHARK_CRITICAL_REGION{
 		for( std::size_t i = 0; i != size(); i++ ) {
-			z( i ) = Rng::gauss( 0, 1 );
-		}
+			z( i ) = gauss(rng, 0, 1 );
 		}
-		if(m_triangular && size() > 400){
-			y=z;
+		if(size() > 400){
+			y = z;
 			blas::triangular_prod<blas::lower>(m_lowerCholesky,y);
 		}else{
 			noalias(y) = prod(m_lowerCholesky,z);
@@ -217,17 +214,17 @@ public:
 	///
 	/// Returns a vector pair (y,z) where  y=Lz and, L is the lower cholesky factor and z is a vector
 	/// of normally distributed numbers. Thus y is the real sampled point.
-	result_type operator()() const {
+	template<class RngType>
+	result_type operator()(RngType& rng) const {
 		result_type result;
 		RealVector& z = result.second;
 		RealVector& y = result.first;
-		generate(y,z);
+		generate(rng,y,z);
 		return result;
 	}
 
 private:
 	blas::matrix<double,blas::column_major> m_lowerCholesky; ///< The lower cholesky factor (actually any is okay as long as it is the left)
-	bool m_triangular;
 };
 
 
diff --git a/src/Algorithms/DirectSearch/CMA.cpp b/src/Algorithms/DirectSearch/CMA.cpp
index c679833..b82cf6b 100644
--- a/src/Algorithms/DirectSearch/CMA.cpp
+++ b/src/Algorithms/DirectSearch/CMA.cpp
@@ -89,7 +89,7 @@ std::size_t CMA::suggestMu( std::size_t lambda, RecombinationType recomb) {
 	return 0;
 }
 
-CMA::CMA()
+CMA::CMA(DefaultRngType& rng)
 : m_userSetMu(false)
 , m_userSetLambda(false)
 , m_initSigma(-1)
@@ -102,7 +102,8 @@ CMA::CMA()
 , m_dSigma( 0 )
 , m_muEff( 0 )
 , m_lowerBound( 1E-20)
-, m_counter( 0 ) {
+, m_counter( 0 )
+, mpe_rng(&rng){
 	m_features |= REQUIRES_VALUE;
 }
 
@@ -270,7 +271,7 @@ void CMA::doInit(
 std::vector<CMA::IndividualType> CMA::generateOffspring( ) const{
 	std::vector< IndividualType > offspring( m_lambda );
 	for( std::size_t i = 0; i < offspring.size(); i++ ) {
-		MultiVariateNormalDistribution::result_type sample = m_mutationDistribution();
+		MultiVariateNormalDistribution::result_type sample = m_mutationDistribution(*mpe_rng);
 		offspring[i].chromosome() = sample.second;
 		offspring[i].searchPoint() = m_mean + m_sigma * sample.first;
 	}
diff --git a/src/Algorithms/DirectSearch/CMSA.cpp b/src/Algorithms/DirectSearch/CMSA.cpp
index 558a6dc..1f2e8ce 100644
--- a/src/Algorithms/DirectSearch/CMSA.cpp
+++ b/src/Algorithms/DirectSearch/CMSA.cpp
@@ -138,8 +138,8 @@ void CMSA::step(ObjectiveFunctionType const& function){
 std::vector<CMSA::IndividualType> CMSA::generateOffspring( ) const{
 	std::vector< IndividualType > offspring( m_lambda );
 	for( std::size_t i = 0; i < offspring.size(); i++ ) {		    
-		MultiVariateNormalDistribution::result_type sample = m_mutationDistribution();
-		offspring[i].chromosome().sigma = m_sigma * ::exp( m_cSigma * Rng::gauss( 0, 1 ) );
+		MultiVariateNormalDistribution::result_type sample = m_mutationDistribution(*mpe_rng);
+		offspring[i].chromosome().sigma = m_sigma * ::exp( m_cSigma * gauss(*mpe_rng, 0, 1 ) );
 		offspring[i].chromosome().step = sample.first;
 		offspring[i].searchPoint() = m_mean + offspring[i].chromosome().sigma * sample.first;
 	}
diff --git a/src/Algorithms/DirectSearch/ElitistCMA.cpp b/src/Algorithms/DirectSearch/ElitistCMA.cpp
index 00cc990..da6ad3b 100644
--- a/src/Algorithms/DirectSearch/ElitistCMA.cpp
+++ b/src/Algorithms/DirectSearch/ElitistCMA.cpp
@@ -34,7 +34,7 @@
 
 using namespace shark;
 
-ElitistCMA::ElitistCMA(): m_activeUpdate(true) {
+ElitistCMA::ElitistCMA(DefaultRngType& rng): m_activeUpdate(true),mpe_rng(&rng){
 	m_features |= REQUIRES_VALUE;
 }
 
@@ -72,7 +72,7 @@ void ElitistCMA::init( ObjectiveFunctionType& function, SearchPointType const& p
 
 void ElitistCMA::step(ObjectiveFunctionType const& function) {
 	//create and evaluate offspring
-	m_individual.mutate();
+	m_individual.mutate(*mpe_rng);
 	m_evaluator( function, m_individual );
 	
 	//evaluate success status of individual

-- 
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/shark.git



More information about the debian-science-commits mailing list