removed random number manager,

only generators building on each other will be used

mlmc/src/basics/Utilities.cpp

@@ -28,4 +28,32 @@ std::string vec2str<std::string>(std::vector<std::string> vec) {
         else str += vec[i];
     }
     return (str + "]");
+}
+
+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);
 }
\ No newline at end of file

mlmc/src/basics/Utilities.hpp

@@ -3,9 +3,26 @@
 
 #include <vector>
 #include <string>
+#include <math.h>
 
 
 template<typename T>
 std::string vec2str(std::vector<T> vec);
 
+double calc_mean(const std::vector<double> &vec);
+
+double calc_mean(const std::vector<std::vector<double>> &vec);
+
+double calc_mean(const std::vector<int> &vec);
+
+double calc_mean(const std::vector<std::vector<int>> &vec);
+
+double calc_var(const std::vector<double> &vec);
+
+double calc_var(const std::vector<std::vector<double>> &vec);
+
+double calc_var(const std::vector<int> &vec);
+
+double calc_var(const std::vector<std::vector<int>> &vec);
+
 #endif //TOOLS_HPP

mlmc/src/problems/CMakeLists.txt

@@ -4,6 +4,7 @@ add_library(PROBLEMS STATIC
         StochasticTransportProblem.cpp
         generators/CirculantEmbedding.cpp
         generators/HybridFluxGenerator.cpp
-        generators/RandomNumberManager.cpp
+        generators/NormalDistribution.cpp
+        generators/UniformDistribution.cpp
         )
 target_link_libraries(PROBLEMS BASICS fftw3 sprng ${MPP_LIBRARIES})

mlmc/src/problems/generators/CirculantEmbedding.cpp

 #include "CirculantEmbedding.hpp"
 
 
-void CirculantEmbedding1D::generateFineSample(SampleID id,
+void CirculantEmbedding1D::generateFineSample(const SampleID &id,
                                               Vector *&fineSample,
                                               Vector *&coarseSample) {
     fineSample = new Vector(cellMGraphs[id.level.mGfine]);
 
     if (internalCounter == 0)
-        fineComplexField = generateField();
+        fineComplexField = generateField(id);
 
     for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
         int i = floor(c()[0] * fineSample->size());
@@ -30,7 +30,7 @@ void CirculantEmbedding1D::generateFineSample(SampleID id,
     coarseSample = nullptr;
 }
 
-void CirculantEmbedding1D::generateCoarseSample(SampleID id,
+void CirculantEmbedding1D::generateCoarseSample(const SampleID &id,
                                                 Vector *&fineSample,
                                                 Vector *&coarseSample) {
     coarseSample = new Vector(cellMGraphs[id.level.mGcoarse]);
@@ -52,16 +52,15 @@ void CirculantEmbedding1D::generateCoarseSample(SampleID id,
     fineSample = nullptr;
 }
 
-ComplexField1D CirculantEmbedding1D::generateField() {
+ComplexField1D CirculantEmbedding1D::generateField(const SampleID &id) {
     ComplexField1D Xf;
     ComplexField1D normallyDistributed;
     normallyDistributed.resize(numberOfCellsEmbedded);
     fineComplexField.clear();
     Xf.resize(numberOfCellsEmbedded);
-    double z1, z2;
     for (int i = 0; i < numberOfCellsEmbedded; i++) {
-        rnmg->DrawNormalRandomNumbers(z1, z2);
-        normallyDistributed[i] = {z1, z2};
+        normalDist.DrawSample(id);
+        normallyDistributed[i] = normalDist.EvalComplexSample();
     }
 
     ComplexField1D product;
@@ -157,13 +156,13 @@ SqrtEigenValues1D CirculantEmbedding1D::computeSqrtEV() {
     return sqrt_ev_return;
 }
 
-void CirculantEmbedding2D::generateFineSample(SampleID id,
+void CirculantEmbedding2D::generateFineSample(const SampleID &id,
                                               Vector *&fineSample,
                                               Vector *&coarseSample) {
     fineSample = new Vector(cellMGraphs[id.level.mGfine]);
 
     if (internalCounter == 0)
-        fineComplexField = generateField();
+        fineComplexField = generateField(id);
 
     for (cell c = fineSample->cells(); c != fineSample->cells_end(); c++) {
         int i = floor(c()[0] * sqrt(fineSample->size()));
@@ -189,7 +188,7 @@ void CirculantEmbedding2D::generateFineSample(SampleID id,
     coarseSample = nullptr;
 }
 
-void CirculantEmbedding2D::generateCoarseSample(SampleID id,
+void CirculantEmbedding2D::generateCoarseSample(const SampleID &id,
                                                 Vector *&fineSample,
                                                 Vector *&coarseSample) {
     coarseSample = new Vector(cellMGraphs[id.level.mGcoarse]);
@@ -211,19 +210,18 @@ void CirculantEmbedding2D::generateCoarseSample(SampleID id,
     fineSample = nullptr;
 }
 
-ComplexField2D CirculantEmbedding2D::generateField() {
+ComplexField2D CirculantEmbedding2D::generateField(const SampleID &id) {
     ComplexField2D Xf;
     ComplexField2D Z;
     Z.resize(numberOfCellsEmbedded[1]);
     fineComplexField.clear();
     Xf.resize(numberOfCellsEmbedded[1]);
-    double z1, z2;
     for (int i = 0; i < numberOfCellsEmbedded[1]; i++) {
         Z[i].resize(numberOfCellsEmbedded[0]);
         Xf[i].resize(numberOfCellsEmbedded[0]);
         for (int j = 0; j < numberOfCellsEmbedded[0]; j++) {
-            rnmg->DrawNormalRandomNumbers(z1, z2);
-            Z[i][j] = {z1, z2};
+            normalDist.DrawSample(id);
+            Z[i][j] = normalDist.EvalComplexSample();
         }
     }
 

mlmc/src/problems/generators/CirculantEmbedding.hpp

 #ifndef M_CIRCULANTEMBEDDING_H
 #define M_CIRCULANTEMBEDDING_H
 
-#include "RandomNumberManager.hpp"
 #include "CovarianceFunction.hpp"
 #include "basics/Utilities.hpp"
 #include "Algebra.hpp"
 #include "SampleGenerator.hpp"
 #include "dof/BasicDoFs.hpp"
-//#include "main/MultilevelPlotter.hpp"
-
+#include "NormalDistribution.hpp"
 
 typedef std::vector<std::complex<double>> ComplexField1D;
 
@@ -38,7 +36,7 @@ class CirculantEmbedding {
 protected:
     int internalCounter = 0;
 
-    RandomNumberManager *rnmg;
+    ComplexNormalDistribution normalDist;
 
     double evtol = 1e-10;
 
@@ -52,27 +50,24 @@ protected:
 
 public:
     explicit CirculantEmbedding(Meshes &meshes) :
-        meshes(meshes), cellMGraphs(meshes, dof(new CellDoF(1))) {
+        meshes(meshes),
+        normalDist(meshes),
+        cellMGraphs(meshes, dof(new CellDoF(1))) {
 
         config.get("evtol", evtol);
         config.get("StochasticField", fieldType);
         config.get("mean", mean);
-
-        int streamnum = 0, nstreams = 1, gtype = 0;
-        rnmg = new RandomNumberManager(streamnum, nstreams, gtype);
     }
 
-    virtual void generateFineSample(SampleID id,
+    virtual void generateFineSample(const SampleID &id,
                                     Vector *&fineSample,
                                     Vector *&coarseSample) = 0;
 
-    virtual void generateCoarseSample(SampleID id,
+    virtual void generateCoarseSample(const SampleID &id,
                                       Vector *&fineSample,
                                       Vector *&coarseSample) = 0;
 
-    virtual ~CirculantEmbedding() {
-        delete rnmg;
-    };
+    virtual ~CirculantEmbedding() {};
 };
 
 class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
@@ -94,16 +89,16 @@ public:
         numberOfCellsEmbedded = sqrtEigenvalues.size();
     }
 
-    void generateFineSample(SampleID id,
+    void generateFineSample(const SampleID &id,
                             Vector *&fineSample,
                             Vector *&coarseSample) override;
 
-    void generateCoarseSample(SampleID id,
+    void generateCoarseSample(const SampleID &id,
                               Vector *&fineSample,
                               Vector *&coarseSample) override;
 
 //private:
-    ComplexField1D generateField();
+    ComplexField1D generateField(const SampleID &id);
 
     void createCovarianceMatrix(ToeplitzRow &ar,
                                 ToeplitzColumn &ac);
@@ -144,16 +139,16 @@ public:
         numberOfCellsEmbedded[1] = sqrtEigenvalues[0].size();
     }
 
-    void generateFineSample(SampleID id,
+    void generateFineSample(const SampleID &id,
                             Vector *&fineSample,
                             Vector *&coarseSample) override;
 
-    void generateCoarseSample(SampleID id,
+    void generateCoarseSample(const SampleID &id,
                               Vector *&fineSample,
                               Vector *&coarseSample) override;
 
 //private:
-    ComplexField2D generateField();
+    ComplexField2D generateField(const SampleID &id);
 
     void createCovarianceMatrix(BlockToeplitzRow &ar,
                                 BlockToeplitzColumn &ac);

mlmc/src/problems/generators/CovarianceFunction.hpp

@@ -42,11 +42,11 @@ class CovarianceFunction1D : CovarianceFunction {
 public:
     double covariance_fct(double *tau) override {
         tau[0] /= lambda[0];
-        return pow(sigma, 2) * exp(-pow(norm(tau), smoothing));
+        return std::pow(sigma, 2) * std::exp(-std::pow(norm(tau), smoothing));
     }
 
     double norm(const double *x) override {
-        return abs(x[0]);
+        return std::abs(x[0]);
     }
 };
 
@@ -55,14 +55,14 @@ public:
     double covariance_fct(double *tau) override {
         tau[0] /= lambda[0];
         tau[1] /= lambda[1];
-        return pow(sigma, 2) * exp(-pow(norm(tau), smoothing));
+        return std::pow(sigma, 2) * std::exp(-std::pow(norm(tau), smoothing));
     }
 
     double norm(const double *x) override {
         if (norm_p == 1)
-            return abs(x[0]) + abs(x[1]);
+            return std::abs(x[0]) + std::abs(x[1]);
         if (norm_p == 2)
-            return sqrt(pow(x[0], 2) + pow(x[1], 2));
+            return sqrt(std::pow(x[0], 2) + std::pow(x[1], 2));
         else Exit("Choose p=1 or p=2")
     }
 };

mlmc/src/problems/generators/NormalDistribution.cpp

+#include "NormalDistribution.hpp"

mlmc/src/problems/generators/NormalDistribution.hpp

+#ifndef NORMALDISTRIBUTION_HPP
+#define NORMALDISTRIBUTION_HPP
+
+#include "SampleGenerator.hpp"
+#include "sprng_cpp.h"
+#include <math.h>
+
+
+#define PI 3.14159265358979323846
+#define SEED 985456376
+
+class NormalDistribution : public SampleGenerator {
+private:
+    Sprng *randNumberGen;
+
+    Scalar firstSample;
+
+    Scalar secondSample;
+
+    int internalCounter = 0;
+
+    void drawSample(const SampleID &id) override {
+        // Todo Alternative: Ziggurat algorithm
+        if (internalCounter) {
+            updateCounter();
+            return;
+        }
+        double randNum1 = randNumberGen->sprng();
+        double randNum2 = randNumberGen->sprng();
+        firstSample = sqrt(-2.0 * log(randNum1)) * cos(2.0 * PI * randNum2);
+        secondSample = sqrt(-2.0 * log(randNum1)) * sin(2.0 * PI * randNum2);
+        updateCounter();
+    }
+
+    void updateCounter() {
+        internalCounter = internalCounter + 1;
+        internalCounter = internalCounter % 2;
+    };
+
+public:
+    NormalDistribution(Meshes &meshes,
+                       int streamID = 0,
+                       int numStreams = 1,
+                       int gType = 0,
+                       int seed = SEED) :
+        SampleGenerator(meshes) {
+        this->randNumberGen = SelectType(gType);
+        this->randNumberGen->init_sprng(streamID, numStreams,
+                                        seed, SPRNG_DEFAULT);
+    }
+
+    ~NormalDistribution() { randNumberGen->free_sprng(); };
+
+    Scalar EvalScalarSample() override {
+        if (internalCounter)
+            return firstSample;
+        else
+            return secondSample;
+    }
+
+    string Name() const override {
+        return "NormalDistribution";
+    };
+};
+
+class ComplexNormalDistribution : public SampleGenerator {
+private:
+    NormalDistribution normalDist;
+
+    Complex sample;
+
+    void drawSample(const SampleID &id) override {
+        normalDist.DrawSample(id);
+        Scalar real = normalDist.EvalScalarSample();
+        normalDist.DrawSample(id);
+        Scalar imag = normalDist.EvalScalarSample();
+        sample = {real, imag};
+    }
+
+public:
+    ComplexNormalDistribution(Meshes &meshes) :
+        SampleGenerator(meshes),
+        normalDist(meshes) {}
+
+    Complex EvalComplexSample() override {
+        return sample;
+    }
+
+    string Name() const override {
+        return "ComplexNormalDistribution";
+    };
+};
+
+#endif //NORMALDISTRIBUTION_HPP

mlmc/src/problems/generators/RandomNumberManager.cpp

-#include "RandomNumberManager.hpp"
-
-
-RandomNumberManager::RandomNumberManager(int streamnum, int nstreams, int gtype, int seed) {
-    this->streamnum = streamnum;
-    this->nstreams = nstreams;
-    this->gtpye = gtype;
-    this->seed = seed;
-
-    this->sprng_ = SelectType(gtype);
-    this->sprng_->init_sprng(streamnum, nstreams, seed, SPRNG_DEFAULT);
-}
-
-void RandomNumberManager::get_info_sprng() {
-    sprng_->print_sprng();
-}
-
-void RandomNumberManager::DrawUniformRandomNumbers(double &u1, double &u2) {
-    u1 = sprng_->sprng();
-    u2 = sprng_->sprng();
-}
-
-void RandomNumberManager::get_uniform_irn(int &i1, int &i2) {
-    i1 = sprng_->isprng();
-    i2 = sprng_->isprng();
-}
-
-void RandomNumberManager::DrawNormalRandomNumbers(double &z1, double &z2) {
-    /*
-     * Alternative: Ziggurat algorithm
-     */
-    double 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;
-}

mlmc/src/problems/generators/RandomNumberManager.hpp

-#ifndef MLMC_RNG_UTILITIES_HPP
-#define MLMC_RNG_UTILITIES_HPP
-
-#include "SampleGenerator.hpp"
-
-#include "sprng_cpp.h"
-#include <math.h>
-
-#define PI 3.14159265358979323846
-#define SEED 985456376
-
-class RandomNumberManager {
-public:
-    Sprng *sprng_;
-
-    int gtpye;
-    int streamnum;
-    int nstreams;
-    int seed;
-
-    RandomNumberManager(int streamnum, int nstreams, int gtype, int seed = SEED);
-
-    // 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 DrawNormalRandomNumbers(double &z1, double &z2);
-
-    double DrawNormalRandomNumber();
-
-    void get_info_sprng();
-
-    ~RandomNumberManager() { sprng_->free_sprng(); };
-};
-
-class NormalDistribution : public SampleGenerator {
-
-    virtual string Name() const {
-        return "NormalDistribution";
-    };
-};
-
-class UniformDistribution : public SampleGenerator {
-
-    virtual string Name() const {
-        return "UniformDistribution";
-    };
-};
-
-#endif //MLMC_RNG_UTILITIES_HPP

mlmc/src/problems/generators/SampleGenerator.hpp

@@ -6,6 +6,8 @@
 #include "mesh/Meshes.hpp"
 
 
+typedef std::complex<double> Complex;
+
 class SampleGenerator {
 protected:
     int plotting = 0;
@@ -36,6 +38,10 @@ public:
     /*
      * Scalar return values
      */
+    virtual Scalar EvalScalarSample() {
+        Exit("Not implemented")
+    }
+
     virtual Scalar EvalScalarSample(const Point &) {
         Exit("Not implemented")
     }
@@ -60,10 +66,44 @@ public:
         Exit("Not implemented")
     }
 
+    /*
+     * Complex return values
+     */
+    virtual Complex EvalComplexSample() {
+        Exit("Not implemented")
+    }
+
+    virtual Complex EvalComplexSample(const Point &) {
+        Exit("Not implemented")
+    }
+
+    virtual Complex EvalComplexSample(int, const Point &) {
+        Exit("Not implemented")
+    }
+
+    virtual Complex EvalComplexSample(double, const Point &) {
+        Exit("Not implemented")
+    }
+
+    virtual Complex EvalComplexSample(const cell &) {
+        Exit("Not implemented")
+    }
+
+    virtual Complex EvalComplexSample(int, const cell &) {
+        Exit("Not implemented")
+    }
+
+    virtual Complex EvalComplexSample(double, const cell &) {
+        Exit("Not implemented")
+    }
+
     /*
      * VectorField return values
-     * Todo abstract class as interface Generator spezifisch
      */
+    virtual VectorField EvalVectorFieldSample() {
+        Exit("Not implemented")
+    }
+
     virtual VectorField EvalVectorFieldSample(const Point &) {
         Exit("Not implemented")
     }
@@ -91,6 +131,10 @@ public:
     /*
      * Tensor return values
      */
+    virtual Tensor EvalTensorSample() {
+        Exit("Not implemented")
+    }
+
     virtual Tensor EvalTensorSample(const Point &) {
         Exit("Not implemented")
     }

mlmc/src/problems/generators/UniformDistribution.cpp

+#include "UniformDistribution.hpp"

mlmc/src/problems/generators/UniformDistribution.hpp

+#ifndef UNIFORMDISTRIBUTION_HPP
+#define UNIFORMDISTRIBUTION_HPP
+
+#include "SampleGenerator.hpp"
+
+
+class UniformDistribution : public SampleGenerator {
+
+public:
+    std::string Name() const override {
+        return "UniformDistribution";
+    };
+};
+
+#endif //UNIFORMDISTRIBUTION_HPP

mlmc/tests/CMakeLists.txt

@@ -23,8 +23,9 @@ endmacro()
 add_mlmc_test(basics/TestLevel)
 add_mlmc_test(basics/TestPlotMap)
 
-add_mlmc_test(problems/generators/TestRandomNumberManager)
+add_mlmc_test(problems/generators/TestNormalDistribution)
 add_mlmc_test(problems/generators/TestCirculantEmbedding)
+
 add_mlmc_test(montecarlo/TestMonteCarlo)
 add_mlmc_test(pdesolver/TestPDESolver)
 add_mlmc_test(montecarlo/datastructure/TestEmpiricMeasureLevelMaps)

mlmc/tests/problems/generators/TestCirculantEmbedding.cpp

@@ -74,7 +74,7 @@ TEST_F(TestCirculantEmbedding1D, TestPaddingDummyExample) {
 
 TEST_F(TestCirculantEmbedding1D, TestGenerateFieldf) {