fixing type of member; renaming members; adding comments

src/Classic/Utilities/Util.cpp

@@ -16,6 +16,87 @@ namespace Util {
         return std::string(GIT_VERSION);
     }
 
+#define erfinv_a3 -0.140543331
+#define erfinv_a2 0.914624893
+#define erfinv_a1 -1.645349621
+#define erfinv_a0 0.886226899
+
+#define erfinv_b4 0.012229801
+#define erfinv_b3 -0.329097515
+#define erfinv_b2 1.442710462
+#define erfinv_b1 -2.118377725
+#define erfinv_b0 1
+
+#define erfinv_c3 1.641345311
+#define erfinv_c2 3.429567803
+#define erfinv_c1 -1.62490649
+#define erfinv_c0 -1.970840454
+
+#define erfinv_d2 1.637067800
+#define erfinv_d1 3.543889200
+#define erfinv_d0 1
+
+    double erfinv (double x) // inverse error function
+    {
+        double x2, r, y;
+        int  sign_x;
+
+        if (x < -1 || x > 1)
+            return NAN;
+
+        if (x == 0)
+            return 0;
+
+        if (x > 0)
+            sign_x = 1;
+        else {
+            sign_x = -1;
+            x = -x;
+        }
+
+        if (x <= 0.7) {
+
+            x2 = x * x;
+            r =
+                x * (((erfinv_a3 * x2 + erfinv_a2) * x2 + erfinv_a1) * x2 + erfinv_a0);
+            r /= (((erfinv_b4 * x2 + erfinv_b3) * x2 + erfinv_b2) * x2 +
+                  erfinv_b1) * x2 + erfinv_b0;
+        }
+        else {
+            y = sqrt (-log ((1 - x) / 2));
+            r = (((erfinv_c3 * y + erfinv_c2) * y + erfinv_c1) * y + erfinv_c0);
+            r /= ((erfinv_d2 * y + erfinv_d1) * y + erfinv_d0);
+        }
+
+        r = r * sign_x;
+        x = x * sign_x;
+
+        r -= (erf (r) - x) / (2 / sqrt (M_PI) * exp (-r * r));
+        r -= (erf (r) - x) / (2 / sqrt (M_PI) * exp (-r * r));
+
+        return r;
+    }
+
+#undef erfinv_a3
+#undef erfinv_a2
+#undef erfinv_a1
+#undef erfinv_a0
+
+#undef erfinv_b4
+#undef erfinv_b3
+#undef erfinv_b2
+#undef erfinv_b1
+#undef erfinv_b0
+
+#undef erfinv_c3
+#undef erfinv_c2
+#undef erfinv_c1
+#undef erfinv_c0
+
+#undef erfinv_d2
+#undef erfinv_d1
+#undef erfinv_d0
+
     Vector_t getTaitBryantAngles(Quaternion rotation, const std::string &elementName) {
         Quaternion rotationBAK = rotation;
 

src/Classic/Utilities/Util.h

@@ -16,6 +16,8 @@
 namespace Util {
     std::string getGitRevision();
 
+    double erfinv(double x);
+
     inline
     double getGamma(Vector_t p) {
         return sqrt(dot(p, p) + 1.0);

src/Sample/SampleGaussianSequence.h

@@ -2,6 +2,7 @@
 #define OPAL_SAMPLE_GAUSSIAN_SEQUENCE_H
 
 #include "Sample/SamplingMethod.h"
+#include "Utilities/Util.h"
 
 #ifdef WITH_UNIT_TESTS
 #include <gtest/gtest_prod.h>
@@ -9,28 +10,32 @@
 
 class SampleGaussianSequence : public SamplingMethod
 {
+    // provides a sequence of sampling points that have a Gaussian distribution
+    // with
+    //      mean = 0.5 * (upper + lower)
+    //      sigma = (upper - lower) / 10
+    // This can be achieved if the integral of the Gaussian between the sampling
+    // points are all equal. The sampling points are therefore computed using
+    // the inverse error function at equally distributed arguments between
+    // -1 and 1.
 
 public:
 
     SampleGaussianSequence(double lower, double upper, size_t modulo, int nSample)
-        : n_m(0)
-        , counter_m(0)
-        , mod_m(modulo)
+        : sampleNr_m(0)
+        , numSamples_m(nSample)
+        , volumeLowerDimensions_m(modulo)
+        , individualCounter_m(0)
     {
         double mean = 0.5 * (lower + upper);
         double sigma = (upper - lower) / 10; // +- 5 sigma
         double factor = sigma / sqrt(2);
         double dx = 2.0 / nSample;
-        double oldY = -2.5;
-        for (long i = 1; i < nSample; ++ i) {
-            double x = -1.0 + i * dx;
-            double y = erfinv(x);
-            chain_m.push_back(mean + factor * (y + oldY));
-            oldY = y;
+        for (long i = 0; i < nSample; ++ i) {
+            double x = -1.0 + (i + 0.5) * dx;
+            double y = Util::erfinv(x);
+            sampleChain_m.push_back(mean + factor * y);
         }
-
-        chain_m.push_back(mean + factor * (2.5 + oldY));
-
     }
 
     void create(boost::shared_ptr<SampleIndividual>& ind, size_t i) {
@@ -38,7 +43,7 @@ public:
     }
 
     double getNext() {
-        double sample = chain_m[n_m];
+        double sample = sampleChain_m[sampleNr_m];
         incrementCounter();
 
         return sample;
@@ -48,99 +53,19 @@ private:
 #ifdef WITH_UNIT_TESTS
     FRIEND_TEST(GaussianSampleTest, ChainTest);
 #endif
-    std::vector<double> chain_m;
-    unsigned int n_m;
-    size_t counter_m;
-    size_t mod_m;
+    std::vector<double> sampleChain_m;
+    unsigned int sampleNr_m;
+    unsigned int numSamples_m; // size of this "dimension"
+    size_t volumeLowerDimensions_m; // the "volume" of the sampling space of the lower "dimensions"
+    size_t individualCounter_m; // counts how many "individuals" have been created
 
     void incrementCounter() {
-        ++ counter_m;
-        if (counter_m % mod_m == 0)
-            ++ n_m;
-        if (n_m >= chain_m.size())
-            n_m = 0;
-    }
-
-#define erfinv_a3 -0.140543331
-#define erfinv_a2 0.914624893
-#define erfinv_a1 -1.645349621
-#define erfinv_a0 0.886226899
-
-#define erfinv_b4 0.012229801
-#define erfinv_b3 -0.329097515
-#define erfinv_b2 1.442710462
-#define erfinv_b1 -2.118377725
-#define erfinv_b0 1
-
-#define erfinv_c3 1.641345311
-#define erfinv_c2 3.429567803
-#define erfinv_c1 -1.62490649
-#define erfinv_c0 -1.970840454
-
-#define erfinv_d2 1.637067800
-#define erfinv_d1 3.543889200
-#define erfinv_d0 1
-
-double erfinv (double x)
-{
-    double x2, r, y;
-    int  sign_x;
-
-    if (x < -1 || x > 1)
-        return NAN;
-
-    if (x == 0)
-        return 0;
+        ++ individualCounter_m;
+        if (individualCounter_m % volumeLowerDimensions_m == 0)
+            ++ sampleNr_m;
 
-    if (x > 0)
-        sign_x = 1;
-    else {
-        sign_x = -1;
-        x = -x;
+        sampleNr_m = sampleNr_m % numSamples_m;
     }
-
-    if (x <= 0.7) {
-
-        x2 = x * x;
-        r =
-            x * (((erfinv_a3 * x2 + erfinv_a2) * x2 + erfinv_a1) * x2 + erfinv_a0);
-        r /= (((erfinv_b4 * x2 + erfinv_b3) * x2 + erfinv_b2) * x2 +
-              erfinv_b1) * x2 + erfinv_b0;
-    }
-    else {
-        y = sqrt (-log ((1 - x) / 2));
-        r = (((erfinv_c3 * y + erfinv_c2) * y + erfinv_c1) * y + erfinv_c0);
-        r /= ((erfinv_d2 * y + erfinv_d1) * y + erfinv_d0);
-    }
-
-    r = r * sign_x;
-    x = x * sign_x;
-
-    r -= (erf (r) - x) / (2 / sqrt (M_PI) * exp (-r * r));
-    r -= (erf (r) - x) / (2 / sqrt (M_PI) * exp (-r * r));
-
-    return r;
-}
-
-#undef erfinv_a3
-#undef erfinv_a2
-#undef erfinv_a1
-#undef erfinv_a0
-
-#undef erfinv_b4
-#undef erfinv_b3
-#undef erfinv_b2
-#undef erfinv_b1
-#undef erfinv_b0
-
-#undef erfinv_c3
-#undef erfinv_c2
-#undef erfinv_c1
-#undef erfinv_c0
-
-#undef erfinv_d2
-#undef erfinv_d1
-#undef erfinv_d0
 };
 
 #endif
\ No newline at end of file

src/Sample/SampleSequence.h

@@ -6,39 +6,43 @@
 template <typename T>
 class SampleSequence : public SamplingMethod
 {
+    // provides a sequence of equidistant sampling points. It
+    // can't be garanteed that the sampling is equidistant if
+    // an integer type is chosen and the difference between
+    // the upper and lower limit isn't divisible by the number
+    // of sampling points.
 
 public:
 
     SampleSequence(T lower, T upper, size_t modulo, int nSample)
-        : lower_m(lower)
-        , step_m( (upper - lower) / double(nSample - 1) )
-        , n_m(0)
-        , size_m(nSample)
-        , counter_m(0)
-        , mod_m(modulo)
-
+        : lowerLimit_m(lower)
+        , stepSize_m( (upper - lower) / double(nSample - 1) )
+        , sampleNr_m(0)
+        , numSamples_m(nSample)
+        , volumeLowerDimensions_m(modulo)
+        , individualCounter_m(0)
     { }
 
     void create(boost::shared_ptr<SampleIndividual>& ind, size_t i) {
-        ind->genes[i] = static_cast<T>(lower_m + step_m * n_m);
+        ind->genes[i] = static_cast<T>(lowerLimit_m + stepSize_m * sampleNr_m);
         incrementCounter();
     }
 
 private:
     void incrementCounter() {
-        ++ counter_m;
-        if (counter_m % mod_m == 0)
-            ++ n_m;
-        if (n_m >= size_m)
-            n_m = 0;
+        ++ individualCounter_m;
+        if (individualCounter_m % volumeLowerDimensions_m == 0)
+            ++ sampleNr_m;
+
+        sampleNr_m = sampleNr_m % numSamples_m;
     }
 
-    T lower_m;
-    double step_m;
-    unsigned int n_m;
-    int size_m;
-    size_t counter_m;
-    size_t mod_m;
+    T lowerLimit_m;
+    double stepSize_m;
+    unsigned int sampleNr_m;
+    unsigned int numSamples_m; // size of this "dimension"
+    size_t volumeLowerDimensions_m; // the "volume" of the sampling space of the lower "dimensions"
+    size_t individualCounter_m; // counts how many "individuals" have been created
 };
 
 #endif
\ No newline at end of file