adding Gaussian Sampling

src/Sample/OpalSample.cpp

@@ -7,6 +7,7 @@
 
 #include "Sample/Uniform.h"
 #include "Sample/SampleSequence.h"
+#include "Sample/SampleGaussianSequence.h"
 #include "Sample/FromFile.h"
 
 
@@ -30,13 +31,13 @@ OpalSample::OpalSample():
 {
     itsAttr[TYPE]       = Attributes::makeString
                           ("TYPE", "UNIFORM_INT, UNIFORM_REAL, SEQUENCE, FROMFILE");
-    
+
     itsAttr[VARIABLE]   = Attributes::makeString
                           ("VARIABLE", "Name of design variable");
-    
+
     itsAttr[SEED]       = Attributes::makeReal
                           ("SEED", "seed for random sampling (default: 42)", 42);
-    
+
     itsAttr[FNAME]      = Attributes::makeString
                           ("FNAME", "File to read from the sampling points");
 
@@ -56,7 +57,7 @@ OpalSample *OpalSample::clone(const std::string &name) {
 
 
 void OpalSample::execute() {
-    
+
 }
 
 
@@ -72,21 +73,23 @@ OpalSample *OpalSample::find(const std::string &name) {
 
 
 void OpalSample::initOpalSample(double lower, double upper, int nSample) {
-    
+
     if ( lower >= upper )
         throw OpalException("OpalSample::initOpalSample()",
                                 "Lower bound >= upper bound.");
-    
+
     std::string type = Util::toUpper(Attributes::getString(itsAttr[TYPE]));
-    
+
     int seed = Attributes::getReal(itsAttr[SEED]);
-        
+
     if (type == "UNIFORM_INT") {
         sampleMethod_m.reset( new Uniform<int>(lower, upper, seed) );
     } else if (type == "UNIFORM_REAL") {
         sampleMethod_m.reset( new Uniform<double>(lower, upper, seed) );
     } else if (type == "SEQUENCE") {
         sampleMethod_m.reset( new SampleSequence(lower, upper, nSample) );
+    } else if (type == "GAUSSIAN_SEQUENCE") {
+        sampleMethod_m.reset( new SampleGaussianSequence(lower, upper, nSample) );
     } else if (type == "FROMFILE") {
         std::string fname = Attributes::getString(itsAttr[FNAME]);
         sampleMethod_m.reset( new FromFile(fname) );
@@ -99,4 +102,4 @@ void OpalSample::initOpalSample(double lower, double upper, int nSample) {
 
 std::string OpalSample::getVariable() const {
     return Attributes::getString(itsAttr[VARIABLE]);
-}
+}
\ No newline at end of file

src/Sample/SampleGaussianSequence.h

+#ifndef OPAL_SAMPLE_GAUSSIAN_SEQUENCE_H
+#define OPAL_SAMPLE_GAUSSIAN_SEQUENCE_H
+
+#include "Sample/SamplingMethod.h"
+
+#ifdef WITH_UNIT_TESTS
+#include <gtest/gtest_prod.h>
+#endif
+
+class SampleGaussianSequence : public SamplingMethod
+{
+
+public:
+
+    SampleGaussianSequence(double lower, double upper, int nSample)
+        : n_m(0)
+    {
+        double mean = 0.5 * (lower + upper);
+        double sigma = (upper - lower) / 10; // +- 5 sigma
+        double dx = 2.0 / nSample;
+        double oldY = -5.0;
+        for (long i = 1; i < nSample; ++ i) {
+            double x = -1.0 + i * dx;
+            double y = erfinv(x);
+            chain_m.push_back(mean + 0.5 * sigma * (y + oldY));
+            oldY = y;
+        }
+
+        chain_m.push_back(mean + 0.5 * sigma * (5.0 + oldY));
+    }
+
+    void create(boost::shared_ptr<SampleIndividual>& ind, size_t i) {
+        ind->genes[i] = getNext();
+    }
+
+    double getNext() {
+        double sample = chain_m[n_m];
+        incrementCounter();
+
+        return sample;
+    }
+
+private:
+#ifdef WITH_UNIT_TESTS
+    FRIEND_TEST(GaussianSampleTest, ChainTest);
+#endif
+    std::vector<double> chain_m;
+    unsigned int n_m;
+
+    void incrementCounter() {
+        ++ 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;
+
+    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
+};
+
+#endif
\ No newline at end of file

tests/opal_src/CMakeLists.txt

@@ -2,7 +2,7 @@ add_subdirectory (Algorithms/StatisticalErrors)
 add_subdirectory (BasicActions)
 add_subdirectory (Distribution)
 add_subdirectory (Elements)
+add_subdirectory (Sample)
 add_subdirectory (Utilities)
-#add_subdirectory (Distribution)
 
 set (TEST_SRCS_LOCAL ${TEST_SRCS_LOCAL} PARENT_SCOPE)
\ No newline at end of file

tests/opal_src/Sample/CMakeLists.txt

+set (_SRCS
+    GaussianSampleTest.cpp
+)
+
+include_directories (
+  ${CMAKE_CURRENT_SOURCE_DIR}
+)
+
+add_sources(${_SRCS})
\ No newline at end of file

tests/opal_src/Sample/GaussianSampleTest.cpp

+#include "gtest/gtest.h"
+
+#include "Sample/SampleGaussianSequence.h"
+#include "opal_test_utilities/SilenceTest.h"
+
+TEST(GaussSampleTest, ChainTest) {
+    OpalTestUtilities::SilenceTest silencer;
+    unsigned int nSample = 101;
+    SampleGaussianSequence seq(-5, 5, nSample);
+
+    for (unsigned int i = 0; i * 2 < nSample - 1; ++ i) {
+        seq.getNext();
+    }
+    double x = seq.getNext();
+
+    EXPECT_NEAR(x, 0.0, 1e-8);
+}
\ No newline at end of file