fix halton randomness

IOHprofiler · Feb 13, 2024 · c539faf · c539faf
1 parent 1271da0
commit c539faf
Show file tree

Hide file tree

Showing 7 changed files with 77 additions and 71 deletions.
diff --git a/include/common.hpp b/include/common.hpp
@@ -21,6 +21,7 @@
 using Matrix = Eigen::MatrixXd;
 using Vector = Eigen::VectorXd;
 using Array = Eigen::ArrayXd;
+using size_to = std::optional<size_t>;
 
 template <typename T>
 std::ostream &operator<<(std::ostream &os, const std::vector<T> &x);

diff --git a/include/parameters.hpp b/include/parameters.hpp
@@ -10,8 +10,6 @@
 #include "selection.hpp"
 #include "weights.hpp"
 
-using size_to = std::optional<size_t>;
-
 namespace parameters
 {
     struct Parameters

diff --git a/include/sampling.hpp b/include/sampling.hpp
@@ -133,12 +133,12 @@ namespace sampling
      */
     struct Halton : Sampler
     {
-        Halton(const size_t d, const size_t i = 1);
+        Halton(const size_t d, const size_to i = std::nullopt);
 
         [[nodiscard]] Vector operator()() override;
 
     private:
-        size_t i;
+        size_t seed;
         std::vector<int> primes;
 
         static double next(int index, int base);

diff --git a/src/common.cpp b/src/common.cpp
@@ -1,93 +1,96 @@
 #include "common.hpp"
 
 template <typename T>
-std::ostream& operator<<(std::ostream& os, const std::vector<T>& x)
+std::ostream &operator<<(std::ostream &os, const std::vector<T> &x)
 {
-    for (auto& xi : x)
+    for (auto &xi : x)
         os << xi << ' ';
     return os;
 }
 
-namespace utils {
-std::vector<size_t> sort_indexes(const Vector& v)
-{ 
-    std::vector<size_t> idx(v.size());
-    std::iota(idx.begin(), idx.end(), 0);
+namespace utils
+{
+    std::vector<size_t> sort_indexes(const Vector &v)
+    {
+        std::vector<size_t> idx(v.size());
+        std::iota(idx.begin(), idx.end(), 0);
 
-    std::stable_sort(idx.begin(), idx.end(),
-        [&v](size_t i1, size_t i2) { return v[i1] < v[i2]; });
+        std::stable_sort(idx.begin(), idx.end(),
+                         [&v](size_t i1, size_t i2)
+                         { return v[i1] < v[i2]; });
 
-    return idx;
-}
+        return idx;
+    }
 
-std::vector<size_t> sort_indexes(const std::vector<size_t>& v)
-{
-    std::vector<size_t> idx(v.size());
-    std::iota(idx.begin(), idx.end(), 0);
+    std::vector<size_t> sort_indexes(const std::vector<size_t> &v)
+    {
+        std::vector<size_t> idx(v.size());
+        std::iota(idx.begin(), idx.end(), 0);
 
-    std::stable_sort(idx.begin(), idx.end(),
-        [&v](size_t i1, size_t i2) { return v[i1] < v[i2]; });
+        std::stable_sort(idx.begin(), idx.end(),
+                         [&v](size_t i1, size_t i2)
+                         { return v[i1] < v[i2]; });
 
-    return idx;
-}
+        return idx;
+    }
 
-void hstack(Matrix& X, const Matrix& Y)
-{
-    X.conservativeResize(Eigen::NoChange, X.cols() + Y.cols());
-    X.rightCols(Y.cols()) = Y;
-}
+    void hstack(Matrix &X, const Matrix &Y)
+    {
+        X.conservativeResize(Eigen::NoChange, X.cols() + Y.cols());
+        X.rightCols(Y.cols()) = Y;
+    }
 
-void vstack(Matrix& X, const Matrix& Y)
-{
-    X.conservativeResize(X.rows() + Y.rows(), Eigen::NoChange);
-    X.bottomRows(Y.rows()) = Y;
-}
+    void vstack(Matrix &X, const Matrix &Y)
+    {
+        X.conservativeResize(X.rows() + Y.rows(), Eigen::NoChange);
+        X.bottomRows(Y.rows()) = Y;
+    }
 
-void concat(Vector& x, const Vector& y)
-{
-    x.conservativeResize(x.rows() + y.rows(), Eigen::NoChange);
-    x.bottomRows(y.rows()) = y;
+    void concat(Vector &x, const Vector &y)
+    {
+        x.conservativeResize(x.rows() + y.rows(), Eigen::NoChange);
+        x.bottomRows(y.rows()) = y;
+    }
+
+    std::pair<double, size_t> compute_ert(const std::vector<size_t> &running_times, const size_t budget)
+    {
+        size_t successfull_runs = 0, total_rt = 0;
+
+        for (const auto &rt : running_times)
+        {
+            if (rt < budget)
+                successfull_runs++;
+            total_rt += rt;
+        }
+        return {static_cast<double>(total_rt) / successfull_runs, successfull_runs};
+    }
 }
 
-std::pair<double, size_t> compute_ert(const std::vector<size_t>& running_times, const size_t budget)
+namespace rng
 {
-    size_t successfull_runs = 0, total_rt = 0;
+    int SEED = std::random_device()();
+    std::mt19937 GENERATOR(SEED);
 
-    for (const auto& rt : running_times) {
-        if (rt < budget)
-            successfull_runs++;
-        total_rt += rt;
+    void set_seed(const int seed)
+    {
+        SEED = seed;
+        GENERATOR.seed(seed);
     }
-    return { static_cast<double>(total_rt) / successfull_runs, successfull_runs };
-}
-}
 
-namespace rng {
-int SEED = 0;
-std::mt19937 GENERATOR(SEED);
-
-void set_seed(const int seed)
-{
-    SEED = seed;
-    GENERATOR.seed(seed);
+    int random_integer(int l, int h)
+    {
+        std::uniform_int_distribution<> distrib(l, h);
+        return distrib(GENERATOR);
+    }
 }
 
-int random_integer(int l, int h)
+namespace functions
 {
-    std::uniform_int_distribution<> distrib(l, h);
-    return distrib(GENERATOR);
-}
-}
-
-
-namespace functions {
 
-
-
-    double sphere(const Vector& x)
+    double sphere(const Vector &x)
     {
         double res = 0;
-        for (auto& xi : x)
+        for (auto &xi : x)
             res += xi * xi;
         return res;
     }

diff --git a/src/interface.cpp b/src/interface.cpp
@@ -117,7 +117,7 @@ void define_samplers(py::module &main)
         .def("__call__", &Sobol::operator());
 
     py::class_<Halton, Sampler, std::shared_ptr<Halton>>(m, "Halton")
-        .def(py::init<size_t, size_t>(), py::arg("d"), py::arg("i") = 1)
+        .def(py::init<size_t, size_to>(), py::arg("d"), py::arg("i") = std::nullopt)
         .def("__call__", &Halton::operator());
 
     py::class_<Mirrored, Sampler, std::shared_ptr<Mirrored>>(m, "Mirrored")

diff --git a/src/sampling.cpp b/src/sampling.cpp
@@ -42,7 +42,7 @@ namespace sampling
         return samples.col(current++);
     }
 
-    Halton::Halton(const size_t d, const size_t i) : Sampler(d), i(i)
+    Halton::Halton(const size_t d, const size_to i) : Sampler(d), seed(i.value_or( rng::random_integer(1, INT_MAX)))
     {
         primes = sieve(std::max(6, static_cast<int>(d)));
         while (primes.size() < d)
@@ -54,8 +54,8 @@ namespace sampling
     {
         Vector res(d);
         for (size_t j = 0; j < d; ++j)
-            res(j) = ppf(next(static_cast<int>(i), primes[j]));
-        i++;
+            res(j) = ppf(next(static_cast<int>(seed), primes[j]));
+        seed++;
         return res;
     }
 

diff --git a/tests/test_c_sampling.py b/tests/test_c_sampling.py
@@ -22,13 +22,17 @@ def test_base_sampler_gauss(self):
         self.sampler_test(sampler,  1.107994899)
 
     def test_base_sampler_halton(self):
-        sampler = sampling.Halton(5)
+        sampler = sampling.Halton(5, 1)
         self.sampler_test(sampler, -3.675096792)
 
     def test_base_sampler_sobol(self):
         sampler = sampling.Sobol(5)
         self.sampler_test(sampler, -1.651717819)
 
+    def test_samplers_are_random(self):
+        for sampler in (sampling.Halton, sampling.Sobol):
+            samples = np.array([sampler(5)() for _ in range(5)])
+            self.assertGreater(abs(samples[0] - samples).sum(), 1e-10)
 
     def test_mirrored(self):
         sampler = sampling.Mirrored(sampling.Gaussian(5))