refactored normal distribution

mlmc/src/generators/algorithms/NormalDistribution.cpp → mlmc/src/generators/NormalDistribution.cpp

@@ -4,27 +4,27 @@
 
 #define PI 3.14159265358979323846
 
+// Todo Ziggurat algorithm
 void NormalDistributionReal::drawSample(const SampleID &id) {
-    // 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);
+  if (internalCounter) {
     updateCounter();
+    return;
+  }
+  double randNum1 = uniformDist.DrawAndEvalSample(id);
+  double randNum2 = uniformDist.DrawAndEvalSample(id);
+  firstSample = sqrt(-2.0 * log(randNum1)) * cos(2.0 * PI * randNum2);
+  secondSample = sqrt(-2.0 * log(randNum1)) * sin(2.0 * PI * randNum2);
+  updateCounter();
 }
 
 void NormalDistributionReal::updateCounter() {
-    internalCounter = internalCounter + 1;
-    internalCounter = internalCounter % 2;
+  internalCounter = internalCounter + 1;
+  internalCounter = internalCounter % 2;
 }
 
 Scalar NormalDistributionReal::EvalSample() {
-    if (internalCounter)
-        return firstSample;
-    else
-        return secondSample;
+  if (internalCounter)
+    return firstSample;
+  else
+    return secondSample;
 }

mlmc/src/generators/algorithms/NormalDistribution.hpp → mlmc/src/generators/NormalDistribution.hpp

 #ifndef NORMALDISTRIBUTION_HPP
 #define NORMALDISTRIBUTION_HPP
 
-#include "../SampleGenerator.hpp"
+#include "SampleGenerator.hpp"
+#include "UniformDistribution.hpp"
 
-//#define USE_MPI
 
-#include "sprng_cpp.h"
-
-
-#define SEED 985456376
-
-/*
- * Todo
- */
-
-//class CollocationUniformDist : public SampleGenerator<Scalar> {
-//  list {
-//
-//
-//
-//  }
-//
-//  drawSample(id)
-//    list(id)
-//
-//};
-
-// Todo use UniformRandomNumber
 class NormalDistributionReal : public SampleGenerator<Scalar> {
 private:
-  Sprng *randNumberGen;
+  UniformDistributionReal uniformDist;
 
   Scalar firstSample;
 
@@ -42,19 +20,8 @@ private:
   void updateCounter();
 
 public:
-//    seed + PPM->Proc(0)
-  NormalDistributionReal(Meshes &meshes,
-                         int streamID = 0, // PPM->Proc(0)
-                         int numStreams = 1, // PPM->Size(0)
-                         int gType = 0,
-                         int seed = SEED) :
-    SampleGenerator(meshes) {
-    this->randNumberGen = SelectType(gType);
-    this->randNumberGen->init_sprng(streamID, numStreams,
-                                    seed, SPRNG_DEFAULT);
-  }
-
-  ~NormalDistributionReal() { randNumberGen->free_sprng(); };
+  NormalDistributionReal(Meshes &meshes) :
+    SampleGenerator(meshes), uniformDist(UniformDistributionReal(meshes, 0.0, 1.0)) {}
 
   Scalar EvalSample() override;
 
@@ -70,10 +37,8 @@ private:
   Complex sample;
 
   void drawSample(const SampleID &id) override {
-    normalDist.DrawSample(id);
-    Scalar real = normalDist.EvalSample();
-    normalDist.DrawSample(id);
-    Scalar imag = normalDist.EvalSample();
+    Scalar real = normalDist.DrawAndEvalSample(id);
+    Scalar imag = normalDist.DrawAndEvalSample(id);
     sample = {real, imag};
   }
 
@@ -98,10 +63,8 @@ private:
   RVector sample;
 
   void drawSample(const SampleID &id) override {
-    for (int i = 0; i < sample.size(); i++) {
-      normalDist.DrawSample(id);
-      sample[i] = normalDist.EvalSample();
-    }
+    for (int i = 0; i < sample.size(); i++)
+      sample[i] = normalDist.DrawAndEvalSample(id);
   }
 
 public:
@@ -134,10 +97,8 @@ private:
   CVector sample;
 
   void drawSample(const SampleID &id) override {
-    for (int i = 0; i < sample.size(); i++) {
-      normalDist.DrawSample(id);
-      sample[i] = normalDist.EvalSample();
-    }
+    for (int i = 0; i < sample.size(); i++)
+      sample[i] = normalDist.DrawAndEvalSample(id);
   }
 
 public:
@@ -174,10 +135,8 @@ private:
   RMatrix sample;
 
   void drawSample(const SampleID &id) override {
-    for (int i = 0; i < sample.rows(); i++) {
-      normalDist.DrawSample(id);
-      sample.InsertRow(normalDist.EvalSample(), i);
-    }
+    for (int i = 0; i < sample.rows(); i++)
+      sample.InsertRow(normalDist.DrawAndEvalSample(id), i);
   }
 
 public:
@@ -210,10 +169,8 @@ private:
   CMatrix sample;
 
   void drawSample(const SampleID &id) override {
-    for (int i = 0; i < sample.rows(); i++) {
-      normalDist.DrawSample(id);
-      sample.InsertRow(normalDist.EvalSample(), i);
-    }
+    for (int i = 0; i < sample.rows(); i++)
+      sample.InsertRow(normalDist.DrawAndEvalSample(id), i);
   }
 
 public: