Commit 8cf1a192 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactoring RandomNumberManager

parent c4364cb2
......@@ -60,7 +60,7 @@ ComplexField1D CirculantEmbedding1D::generateField() {
Xf.resize(numberOfCellsEmbedded);
double z1, z2;
for (int i = 0; i < numberOfCellsEmbedded; i++) {
rnmg->get_normal_rn(z1, z2);
rnmg->DrawNormalRandomNumbers(z1, z2);
normallyDistributed[i] = {z1, z2};
}
......@@ -222,7 +222,7 @@ ComplexField2D CirculantEmbedding2D::generateField() {
Z[i].resize(numberOfCellsEmbedded[0]);
Xf[i].resize(numberOfCellsEmbedded[0]);
for (int j = 0; j < numberOfCellsEmbedded[0]; j++) {
rnmg->get_normal_rn(z1, z2);
rnmg->DrawNormalRandomNumbers(z1, z2);
Z[i][j] = {z1, z2};
}
}
......
......@@ -38,7 +38,7 @@ class CirculantEmbedding {
protected:
int internalCounter = 0;
RNManager *rnmg;
RandomNumberManager *rnmg;
double evtol = 1e-10;
......@@ -59,7 +59,7 @@ public:
config.get("mean", mean);
int streamnum = 0, nstreams = 1, gtype = 0;
rnmg = new RNManager(streamnum, nstreams, gtype);
rnmg = new RandomNumberManager(streamnum, nstreams, gtype);
}
virtual void generateFineSample(SampleID id,
......
#include "RandomNumberManager.hpp"
RNManager::RNManager(int streamnum, int nstreams, int gtype, int seed) {
RandomNumberManager::RandomNumberManager(int streamnum, int nstreams, int gtype, int seed) {
this->streamnum = streamnum;
this->nstreams = nstreams;
this->gtpye = gtype;
......@@ -10,26 +11,30 @@ RNManager::RNManager(int streamnum, int nstreams, int gtype, int seed) {
this->sprng_->init_sprng(streamnum, nstreams, seed, SPRNG_DEFAULT);
}
void RNManager::get_info_sprng() {
void RandomNumberManager::get_info_sprng() {
sprng_->print_sprng();
}
void RNManager::get_uniform_rn(double &u1, double &u2) {
void RandomNumberManager::DrawUniformRandomNumbers(double &u1, double &u2) {
u1 = sprng_->sprng();
u2 = sprng_->sprng();
}
void RNManager::get_uniform_irn(int &i1, int &i2) {
void RandomNumberManager::get_uniform_irn(int &i1, int &i2) {
i1 = sprng_->isprng();
i2 = sprng_->isprng();
}
void RNManager::get_normal_rn(double &z1, double &z2) {
void RandomNumberManager::DrawNormalRandomNumbers(double &z1, double &z2) {
/*
* Alternative: Ziggurat algorithm
*/
double u1, u2;
get_uniform_rn(u1, u2);
DrawUniformRandomNumbers(u1, u2);
z1 = sqrt( - 2.0 * log(u1)) * cos( 2.0 * PI * u2 );
z2 = sqrt( - 2.0 * log(u1)) * sin( 2.0 * PI * u2 );
}
double RandomNumberManager::DrawNormalRandomNumber() {
return 0;
}
......@@ -7,7 +7,7 @@
#define PI 3.14159265358979323846
#define SEED 985456376
class RNManager {
class RandomNumberManager {
public:
Sprng *sprng_;
......@@ -16,17 +16,20 @@ public:
int nstreams;
int seed;
RNManager(int streamnum, int nstreams, int gtype, int seed = SEED);
RandomNumberManager(int streamnum, int nstreams, int gtype, int seed = SEED);
void get_uniform_rn(double &u1, double &u2); // For double in [0, 1) uniformly distributed
// For double in [0, 1) uniformly distributed
void DrawUniformRandomNumbers(double &u1, double &u2);
void get_uniform_irn(int &i, int &u2); // For integer in [0, 2^31) uniformly distributed
void get_normal_rn(double &z1, double &z2);
void DrawNormalRandomNumbers(double &z1, double &z2);
double DrawNormalRandomNumber();
void get_info_sprng();
~RNManager() { sprng_->free_sprng(); };
~RandomNumberManager() { sprng_->free_sprng(); };
};
#endif //MLMC_RNG_UTILITIES_HPP
#include "TestEnvironment.hpp"
#include "utility/M_IOFiles.hpp"
#include "main/MultilevelPlotter.hpp"
#include "matrixgraph/MatrixGraph.h"
#include "Algebra.h"
#include "dof/BasicDoFs.hpp"
class TestMultilevelPlotter : public ::Test {
protected:
std::shared_ptr<Meshes> meshes;
int pLevel;
void SetUp() override {
meshes = std::make_shared<Meshes>("UnitSquare", 1, 4);
plotter = std::make_shared<MultilevelPlotter>(*meshes);
pLevel = meshes->pLevel();
}
void TearDown() override {
// Only master checks if files are there, thus only master deletes
if (PPM->master())
system("rm -rf data/vtk/*");
}
};
TEST_F(TestMultilevelPlotter, TestVertexData1D) {
MatrixGraphs vertexMG(*meshes, dof("vertex", 1));
for (int l = pLevel; l <= meshes->Level(); l++) {
Vector vertexData(vertexMG[l - pLevel]);
vertexData = 1.0;
string name = "VertexData1D_" + to_string(l);
plotter->PlotVector(name, vertexData, 1, l, "VertexData");
if (PPM->master())
ASSERT_TRUE(FileExists("data/vtk/" + name + ".vtk"));
}
}
TEST_F(TestMultilevelPlotter, TestVertexData2D) {
MatrixGraphs vertexMG(*meshes, dof("vertex", 2));
for (int l = meshes->pLevel(); l <= meshes->Level(); l++) {
Vector vertexData(vertexMG[l - pLevel]);
vertexData = 1.0;
string name = "VertexData2D_" + to_string(l);
plotter->PlotVector(name, vertexData, 2, l, "VertexData");
if (PPM->master())
ASSERT_TRUE(FileExists("data/vtk/" + name + ".vtk"));
}
}
TEST_F(TestMultilevelPlotter, TestCellData1D) {
MatrixGraphs cellMG(*meshes, dof("cell", 1));
for (int l = meshes->pLevel(); l <= meshes->Level(); l++) {
Vector cellData(cellMG[l - pLevel]);
cellData = 1.0;
string name = "CellData1D_" + to_string(l);
plotter->PlotVector(name, cellData, 1, l, "CellData");
if (PPM->master())
ASSERT_TRUE(FileExists("data/vtk/" + name + ".vtk"));
}
}
TEST_F(TestMultilevelPlotter, TestCellData2D) {
MatrixGraphs cellMG(*meshes, new CellDoF(2));
for (int l = meshes->pLevel(); l <= meshes->Level(); l++) {
Vector cellData(cellMG[l - pLevel]);
cellData = 1.0;
string name = "CellData2D_" + to_string(l);
plotter->PlotVector(name, cellData, 2, l, "CellData");
if (PPM->master())
ASSERT_TRUE(FileExists("data/vtk/" + name + ".vtk"));
}
}
int main(int argc, char **argv) {
MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().
WithScreenLogging().WithSearchPath("../../mlmc/");
return mppTest.RUN_ALL_MPP_TESTS();
}
\ No newline at end of file
#include "gtest/gtest.h"
#include "stochastic/RNManager.h"
#include "stochastic/CirculantEmbedding.h"
TEST(TestRNMangaer, TestNormalRN) {
int test_samples = 1000000;
auto *rnmg = new RNManager(0, 1, 0);
double z1, z2;
vector<double> rns;
rns.reserve(test_samples);
for (int i = 0; i < test_samples; ++i) {
rnmg->get_normal_rn(z1, z2);
rns.push_back(z1);
rns.push_back(z1);
}
double mean = calc_mean(rns);
double var = calc_var(rns);
ASSERT_NEAR(0.0, mean, 1e-2);
ASSERT_NEAR(1.0, var, 1e-2);
}
TEST(TestRNMangaer, TestUniformRN) {
int test_samples = 1000000;
auto *rnmg = new RNManager(0, 1, 0);
double u1, u2;
vector<double> rns;
rns.reserve(test_samples);
for (int i = 0; i < test_samples; ++i) {
rnmg->get_uniform_rn(u1, u2);
rns.push_back(u1);
rns.push_back(u1);
}
double mean = calc_mean(rns);
double var = calc_var(rns);
ASSERT_NEAR(0.5, mean, 1e-2);
ASSERT_NEAR(0.08333333333333333, var, 1e-2);
}
TEST(TestRNMangaer, TestUniformIRN) {
int test_samples = 1000000;
auto *rnmg = new RNManager(0, 1, 0);
int i1, i2;
vector<int> rns;
rns.reserve(test_samples);
for (int i = 0; i < test_samples; ++i) {
rnmg->get_uniform_irn(i1, i2);
rns.push_back(i1);
rns.push_back(i1);
}
double mean = calc_mean(rns);
double var = calc_var(rns);
ASSERT_NEAR(pow(2, 30), mean, 10000000);
ASSERT_NEAR(pow(2, 61) / 12, var, 2e17);
}
int main(int argc, char **argv) {
testing::InitGoogleTest(&argc, argv);
return RUN_ALL_TESTS();
}
\ No newline at end of file
......@@ -6,16 +6,16 @@
class TestCirculantEmbedding1D : public Test {
protected:
CirculantEmbedding1D *ce;
CirculantEmbedding1D *generator;
Meshes *meshes;
void SetUp() override {
meshes = new Meshes("Line", 0, 2);
ce = new CirculantEmbedding1D(2, *meshes);
meshes = createMeshes("Interval_DirichletBC", 0, 2);
generator = new CirculantEmbedding1D(2, *meshes);
}
void TearDown() override {
delete ce;
delete generator;
delete meshes;
}
};
......@@ -23,7 +23,7 @@ protected:
TEST_F(TestCirculantEmbedding1D, TestCovarianceMatrix) {
ToeplitzRow ar;
ToeplitzColumn ac;
ce->createCovarianceMatrix(ar, ac);
generator->createCovarianceMatrix(ar, ac);
ToeplitzRow expectedAr{1.0, 0.3291930, 0.10836802, 0.0356740};
ToeplitzColumn expectedAc{1.0, 0.3291930, 0.10836802, 0.0356740};
for (int i = 0; i < expectedAr.size(); i++) {
......@@ -43,7 +43,7 @@ TEST_F(TestCirculantEmbedding1D, TestEmbedCovarianceMatrix) {
CirculantRow br;
CirculantRow expectedBr{1.0, 0.3291930, 0.10836802, 0.0356740,
0.0356740, 0.10836802, 0.3291930};
ce->embedCovarianceMatrix(ar, ac, br);
generator->embedCovarianceMatrix(ar, ac, br);
for (int i = 0; i < expectedBr.size(); i++) {
ASSERT_NEAR(expectedBr[i], br[i], 1e-6);
}
......@@ -53,7 +53,7 @@ TEST_F(TestCirculantEmbedding1D, TestEigenValues) {
CirculantRow br{1.0, 0.1888756, 0.03567399, 0.03567399, 0.1888756};
EigenValues1D ev;
EigenValues1D expectedEv{1.44909919, 1.05900981, 0.7164406, 0.7164406, 1.05900981};
ce->computeEV(br, ev);
generator->computeEV(br, ev);
for (int i = 0; i < expectedEv.size(); i++) {
ASSERT_NEAR(expectedEv[i], ev[i], 1e-6);
}
......@@ -62,9 +62,10 @@ TEST_F(TestCirculantEmbedding1D, TestEigenValues) {
TEST_F(TestCirculantEmbedding1D, TestPaddingDummyExample) {
ToeplitzRow ar{1.0, 0.3291930, 0.10836802, 0.0356740};
ToeplitzColumn ac{1.0, 0.3291930, 0.10836802, 0.0356740};
ce->padding(ar, ac);
ToeplitzRow expectedAr {1.0, 0.3291930, 0.10836802, 0.0356740, 0.0117436, 0.0038659};
ToeplitzColumn expectedAc{1.0, 0.3291930, 0.10836802, 0.0356740, 0.0117436, 0.0038659};
generator->padding(ar, ac);
ToeplitzRow expectedAr{1.0, 0.3291930, 0.10836802, 0.0356740, 0.0117436, 0.0038659};
ToeplitzColumn
expectedAc{1.0, 0.3291930, 0.10836802, 0.0356740, 0.0117436, 0.0038659};
for (int i = 0; i < expectedAr.size(); i++) {
ASSERT_NEAR(expectedAr[i], ar[i], 1e-6);
ASSERT_NEAR(expectedAc[i], ac[i], 1e-6);
......@@ -73,16 +74,16 @@ TEST_F(TestCirculantEmbedding1D, TestPaddingDummyExample) {
TEST_F(TestCirculantEmbedding1D, TestGenerateFieldf) {
ComplexField1D Xf;
Xf = ce->generateField();
Xf = generator->generateField();
ASSERT_EQ(Xf.size(), 2 * (*meshes)[2].CellCount() - 1);
}
TEST_F(TestCirculantEmbedding1D, TestDrawSample) {
SampleID id{Level(2, meshes->pLevel()),0,false};
Vector *fineSample;
Vector *coarseSample;
ce->generateFineSample(id, fineSample, coarseSample);
}
//TEST_F(TestCirculantEmbedding1D, TestDrawSample) {
// SampleID id{Level(2, meshes->pLevel()), 0, false};
// Vector *fineSample;
// Vector *coarseSample;
// generator->generateFineSample(id, fineSample, coarseSample);
//}
//TEST_F(TestCirculantEmbedding1D, TestGetFieldf) {
// int N = 128; // level 7
......@@ -135,7 +136,7 @@ protected:
Meshes *meshes;
void SetUp() override {
meshes = new Meshes("Line", 0, 2);
meshes = createMeshes("Square_DirichletBC", 0, 2);
ce = new CirculantEmbedding2D(2, *meshes);
}
......@@ -145,40 +146,40 @@ protected:
}
};
TEST_F(TestCirculantEmbedding2D, TestIsotropicCovarianceMatrix) {
BlockToeplitzRow ar; // configmap1
BlockToeplitzColumn ac;
ce->createCovarianceMatrix(ar, ac);
vector<vector<double>> expectedAr{{1.0, 0.1888756, 0.03567399},
{0.1888756, 0.03567399, 0.00673795},
{0.03567399, 0.00673795, 0.00127263}};
vector<vector<double>> expectedAc{{1.0, 0.1888756, 0.03567399},
{0.1888756, 0.03567399, 0.00673795},
{0.03567399, 0.00673795, 0.00127263}};
for (int i = 0; i < expectedAr.size(); i++) {
for (int j = 0; j < expectedAr[0].size(); j++) {
ASSERT_NEAR(expectedAr[i][j], ar[i][j], 1e-6);
ASSERT_NEAR(expectedAc[i][j], ac[i][j], 1e-6);
}
}
// Has to be since symmetric covariance function
for (int i = 0; i < expectedAr.size(); i++) {
for (int j = 0; j < expectedAr[0].size(); j++) {
ASSERT_NEAR(ar[i][j], ac[i][j], 1e-6);
}
}
// Has to be since isotropic (transpose)
for (int i = 0; i < expectedAr.size(); i++) {
for (int j = 0; j < expectedAr[0].size(); j++) {
ASSERT_NEAR(ar[i][j], ar[j][i], 1e-6);
ASSERT_NEAR(ac[i][j], ac[j][i], 1e-6);
}
}
}
//TEST_F(TestCirculantEmbedding2D, TestIsotropicCovarianceMatrix) {
// BlockToeplitzRow ar; // configmap1
// BlockToeplitzColumn ac;
// ce->createCovarianceMatrix(ar, ac);
//
// vector<vector<double>> expectedAr{{1.0, 0.1888756, 0.03567399},
// {0.1888756, 0.03567399, 0.00673795},
// {0.03567399, 0.00673795, 0.00127263}};
// vector<vector<double>> expectedAc{{1.0, 0.1888756, 0.03567399},
// {0.1888756, 0.03567399, 0.00673795},
// {0.03567399, 0.00673795, 0.00127263}};
//
// for (int i = 0; i < expectedAr.size(); i++) {
// for (int j = 0; j < expectedAr[0].size(); j++) {
// ASSERT_NEAR(expectedAr[i][j], ar[i][j], 1e-6);
// ASSERT_NEAR(expectedAc[i][j], ac[i][j], 1e-6);
// }
// }
//
// // Has to be since symmetric covariance function
// for (int i = 0; i < expectedAr.size(); i++) {
// for (int j = 0; j < expectedAr[0].size(); j++) {
// ASSERT_NEAR(ar[i][j], ac[i][j], 1e-6);
// }
// }
//
// // Has to be since isotropic (transpose)
// for (int i = 0; i < expectedAr.size(); i++) {
// for (int j = 0; j < expectedAr[0].size(); j++) {
// ASSERT_NEAR(ar[i][j], ar[j][i], 1e-6);
// ASSERT_NEAR(ac[i][j], ac[j][i], 1e-6);
// }
// }
//}
TEST_F(TestCirculantEmbedding2D, TestAnisotropicCovarianceMatirx) {
BlockToeplitzRow ar; // configmap2
......@@ -477,7 +478,6 @@ TEST_F(TestCirculantEmbedding2D, TestAnisotropicCovarianceMatirx) {
//}
int main(int argc, char **argv) {
MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().
WithScreenLogging().WithSearchPath("../../mlmc/");
MppTest mppTest = MppTestBuilder(argc, argv).WithPPM().WithScreenLogging();
return mppTest.RUN_ALL_MPP_TESTS();
}
\ No newline at end of file
#include "stochastics/RandomNumberManager.hpp"
#include "TestEnvironment.hpp"
double calc_mean(const std::vector<int> &vec) {
double mean_ = 0.0;
for (auto const &value: vec)
mean_ += value;
return mean_ / vec.size();
}
double calc_var(const std::vector<int> &vec) {
double sqrs_ = 0.0;
for (auto const &value: vec)
sqrs_ += pow(value, 2);
return sqrs_ / vec.size() - pow(calc_mean(vec), 2);
}
double calc_mean(const std::vector<double> &vec) {
double mean_ = 0.0;
for (auto const &value: vec)
mean_ += value;
return mean_ / vec.size();
}
double calc_var(const std::vector<double> &vec) {
double sqrs_ = 0.0;
for (auto const &value: vec)
sqrs_ += pow(value, 2);
return sqrs_ / vec.size() - pow(calc_mean(vec), 2);
}
class TestRandomNumberManager : public Test {
protected:
RandomNumberManager *rnm;
int testSamples = 1e6;
std::vector<double> rns;
void SetUp() override {
rnm = new RandomNumberManager(0, 1, 0);
}
void TearDown() override {
delete rnm;
}
};
TEST_F(TestRandomNumberManager, TestNormalRandomNumber) {
double z1, z2;
for (int i = 0; i < testSamples; i++) {
rnm->DrawNormalRandomNumbers(z1, z2);
rns.push_back(z1);
rns.push_back(z1);
}
double mean = calc_mean(rns);
double var = calc_var(rns);
ASSERT_NEAR(0.0, mean, 1e-2);
ASSERT_NEAR(1.0, var, 1e-2);
}
TEST_F(TestRandomNumberManager, TestUniformlyDistributedRandomNumber) {
double z1, z2;
for (int i = 0; i < testSamples; i++) {
rnm->DrawUniformRandomNumbers(z1, z2);
rns.push_back(z1);
rns.push_back(z1);
}
ASSERT_NEAR(0.5, calc_mean(rns), 1e-2);
ASSERT_NEAR(0.08333333333333333, calc_var(rns), 1e-2);
}
int main(int argc, char **argv) {
MppTest mppTest = MppTestBuilder(argc, argv);
return mppTest.RUN_ALL_MPP_TESTS();
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment