ALGS/RandomVector.h at master · TNWorkers/ALGS · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#ifndef RANDOMVECTOR
#define RANDOMVECTOR

#include <Eigen/Dense>

//#include <boost/random/mersenne_twister.hpp>
//#include <boost/random/normal_distribution.hpp>
//#include <boost/random/uniform_real_distribution.hpp>
//#include <boost/random/variate_generator.hpp>

//boost::random::mt11213b MtEngine;
//boost::normal_distribution<double> NormDist(0.,1.); // mean=0, sigma=1
//boost::random::uniform_01<double> UniformDist; // lower=0, upper=1
//boost::random::uniform_real_distribution<double> UniformDist0; // lower=-1, upper=1

#include <random>
#include <thread>

template<typename Scalar, typename RealScalar> Scalar threadSafeRandUniform (RealScalar min, RealScalar max, bool SEED=false) {};
template<typename Scalar, typename RealScalar> Scalar threadSafeRandNormal (RealScalar mean, RealScalar sigma) {};

template<>
double threadSafeRandUniform<double,double> (double min, double max, bool SEED)
{
	static thread_local mt19937 generatorUniformReal(random_device{}());
	if (SEED) generatorUniformReal.seed(std::time(0));
	uniform_real_distribution<double> distribution(min, max);
	return distribution(generatorUniformReal);
}

template<>
complex<double> threadSafeRandUniform<complex<double>, double> (double min, double max, bool FIXED_SEED)
{
	static thread_local mt19937 generatorUniformComplex(random_device{}());
	if (FIXED_SEED) generatorUniformComplex.seed(std::time(0));
	uniform_real_distribution<double> distribution(min, max);
	return complex<double>(distribution(generatorUniformComplex), distribution(generatorUniformComplex));
}

template<>
int threadSafeRandUniform<int,int> (int min, int max, bool SEED)
{
	static thread_local mt19937 generatorUniformInt(random_device{}());
	if (SEED) generatorUniformInt.seed(std::time(0));
	uniform_int_distribution<unsigned> distribution(min, max);
	return distribution(generatorUniformInt);
}

template<>
double threadSafeRandNormal<double,double> (double mean, double sigma)
{
	static thread_local mt19937 generatorNormalReal(random_device{}());
	normal_distribution<double> distribution(mean, sigma);
	return distribution(generatorNormalReal);
}

template<>
complex<double> threadSafeRandNormal<complex<double>,double> (double mean, double sigma)
{
	static thread_local mt19937 generatorNormalComplex(random_device{}());
	normal_distribution<double> distribution(mean, sigma);
	return complex<double>(distribution(generatorNormalComplex),distribution(generatorNormalComplex));
}

//template<typename VectorType, typename Scalar>
//struct GaussianRandomVector
//{
//	static void fill (size_t N, VectorType &Vout);
//};

template<typename VectorType, typename Scalar>
struct GaussianRandomVector
{
	static void fill (size_t N, VectorType &Vout)
	{
		Vout.resize(N);
		for (size_t i=0; i<N; ++i) {Vout(i) = threadSafeRandUniform<Scalar,double>(-1.,1.);}
		normalize(Vout);
	}
};

//template<typename VectorType>
//struct GaussianRandomVector<VectorType,complex<double> >
//{
//	static void fill (size_t N, VectorType &Vout)
//	{
//		Vout.resize(N);
//		for (size_t i=0; i<N; ++i) {Vout(i) = complex<double>(threadSafeRandNormal(0.,1.), threadSafeRandNormal(0.,1.));}
//		normalize(Vout);
//	}
//};

Eigen::MatrixXd randOrtho (size_t N)
{
	Eigen::MatrixXd M(N,N);
	for (size_t i=0; i<N; ++i)
	for (size_t j=0; j<N; ++j)
	{
		M(i,j) = threadSafeRandUniform<double,double>(0.,1.);
	}
	Eigen::HouseholderQR<Eigen::MatrixXd> Quirinus(M);
	Eigen::MatrixXd Qmatrix = Eigen::MatrixXd::Identity(N,N);
	Qmatrix = Quirinus.householderQ() * Qmatrix;
	return Qmatrix;
}

#endif