Adding LHS initialization

src/limbo/init.hpp

@@ -50,6 +50,7 @@
 ///@defgroup init_defaults
 
 #include <limbo/init/grid_sampling.hpp>
+#include <limbo/init/lhs.hpp>
 #include <limbo/init/no_init.hpp>
 #include <limbo/init/random_sampling.hpp>
 #include <limbo/init/random_sampling_grid.hpp>

src/limbo/init/lhs.hpp

+//| Copyright Inria May 2015
+//| This project has received funding from the European Research Council (ERC) under
+//| the European Union's Horizon 2020 research and innovation programme (grant
+//| agreement No 637972) - see http://www.resibots.eu
+//|
+//| Contributor(s):
+//|   - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr)
+//|   - Antoine Cully (antoinecully@gmail.com)
+//|   - Konstantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
+//|   - Federico Allocati (fede.allocati@gmail.com)
+//|   - Vaios Papaspyros (b.papaspyros@gmail.com)
+//|   - Roberto Rama (bertoski@gmail.com)
+//|
+//| This software is a computer library whose purpose is to optimize continuous,
+//| black-box functions. It mainly implements Gaussian processes and Bayesian
+//| optimization.
+//| Main repository: http://github.com/resibots/limbo
+//| Documentation: http://www.resibots.eu/limbo
+//|
+//| This software is governed by the CeCILL-C license under French law and
+//| abiding by the rules of distribution of free software.  You can  use,
+//| modify and/ or redistribute the software under the terms of the CeCILL-C
+//| license as circulated by CEA, CNRS and INRIA at the following URL
+//| "http://www.cecill.info".
+//|
+//| As a counterpart to the access to the source code and  rights to copy,
+//| modify and redistribute granted by the license, users are provided only
+//| with a limited warranty  and the software's author,  the holder of the
+//| economic rights,  and the successive licensors  have only  limited
+//| liability.
+//|
+//| In this respect, the user's attention is drawn to the risks associated
+//| with loading,  using,  modifying and/or developing or reproducing the
+//| software by the user in light of its specific status of free software,
+//| that may mean  that it is complicated to manipulate,  and  that  also
+//| therefore means  that it is reserved for developers  and  experienced
+//| professionals having in-depth computer knowledge. Users are therefore
+//| encouraged to load and test the software's suitability as regards their
+//| requirements in conditions enabling the security of their systems and/or
+//| data to be ensured and,  more generally, to use and operate it in the
+//| same conditions as regards security.
+//|
+//| The fact that you are presently reading this means that you have had
+//| knowledge of the CeCILL-C license and that you accept its terms.
+//|
+#ifndef LIMBO_INIT_LHS_HPP
+#define LIMBO_INIT_LHS_HPP
+
+#include <Eigen/Core>
+
+#include <limbo/tools/macros.hpp>
+#include <limbo/tools/random_generator.hpp>
+
+namespace limbo {
+    namespace defaults {
+        struct init_lhs {
+            ///@ingroup init_defaults
+            BO_PARAM(int, samples, 10);
+        };
+    } // namespace defaults
+    namespace init {
+        /** @ingroup init
+          \rst
+          LHS sampling in [0, 1]^n
+
+          Parameters:
+            - ``int samples`` (total number of samples)
+          \endrst
+        */
+        template <typename Params>
+        struct LHS {
+            template <typename StateFunction, typename AggregatorFunction, typename Opt>
+            void operator()(const StateFunction& seval, const AggregatorFunction&, Opt& opt) const
+            {
+                assert(Params::bayes_opt_bobase::bounded());
+
+                Eigen::MatrixXd H = tools::random_lhs(StateFunction::dim_in(), Params::init_lhs::samples());
+
+                for (int i = 0; i < Params::init_lhs::samples(); i++) {
+                    opt.eval_and_add(seval, H.row(i));
+                }
+            }
+        };
+    } // namespace init
+} // namespace limbo
+
+#endif

src/limbo/tools/random_generator.hpp

@@ -47,15 +47,15 @@
 #ifndef LIMBO_TOOLS_RANDOM_GENERATOR_HPP
 #define LIMBO_TOOLS_RANDOM_GENERATOR_HPP
 
-#include <cstdlib>
 #include <cmath>
+#include <cstdlib>
 #include <ctime>
+#include <external/rand_utils.hpp>
 #include <list>
-#include <stdlib.h>
+#include <mutex>
 #include <random>
+#include <stdlib.h>
 #include <utility>
-#include <mutex>
-#include <external/rand_utils.hpp>
 
 namespace limbo {
     namespace tools {
@@ -135,7 +135,38 @@ namespace limbo {
                 return random_vector_bounded(size);
             return random_vector_unbounded(size);
         }
-    }
-}
+
+        /// @ingroup tools
+        /// generate random samples with LHS in [0, 1]^n
+        Eigen::MatrixXd random_lhs(int dim, int samples)
+        {
+            Eigen::VectorXd cut = Eigen::VectorXd::LinSpaced(samples + 1, 0., 1.);
+            Eigen::MatrixXd u = Eigen::MatrixXd::Zero(samples, dim);
+
+            for (int i = 0; i < samples; i++) {
+                u.row(i) = tools::random_vector(dim, true);
+            }
+
+            Eigen::VectorXd a = cut.head(samples);
+            Eigen::VectorXd b = cut.tail(samples);
+
+            Eigen::MatrixXd rdpoints = Eigen::MatrixXd::Zero(samples, dim);
+            for (int i = 0; i < dim; i++) {
+                rdpoints.col(i) = u.col(i).array() * (b - a).array() + a.array();
+            }
+
+            Eigen::MatrixXd H = Eigen::MatrixXd::Zero(samples, dim);
+            Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> perm(samples);
+            for (int i = 0; i < dim; i++) {
+                perm.setIdentity();
+                std::random_shuffle(perm.indices().data(), perm.indices().data() + perm.indices().size());
+                Eigen::MatrixXd tmp = perm * rdpoints;
+                H.col(i) = tmp.col(i);
+            }
+
+            return H;
+        }
+    } // namespace tools
+} // namespace limbo
 
 #endif

src/tests/test_init_functions.cpp

@@ -116,6 +116,32 @@ BOOST_AUTO_TEST_CASE(no_init)
     BOOST_CHECK(opt.samples().size() == 0);
 }
 
+BOOST_AUTO_TEST_CASE(random_lhs)
+{
+    std::cout << "LHS" << std::endl;
+    struct MyParams : public Params {
+        struct init_lhs {
+            BO_PARAM(int, samples, 10);
+        };
+    };
+
+    using Init_t = init::LHS<MyParams>;
+    using Opt_t = bayes_opt::BOptimizer<MyParams, initfun<Init_t>>;
+
+    Opt_t opt;
+    opt.optimize(fit_eval());
+    BOOST_CHECK(opt.observations().size() == 10);
+    BOOST_CHECK(opt.samples().size() == 10);
+    for (size_t j = 0; j < opt.samples().size() - 1; ++j) {
+        const Eigen::VectorXd& x = opt.samples()[j];
+        std::cout << x.transpose() << std::endl;
+        for (int i = 0; i < x.size(); ++i) {
+            BOOST_CHECK(x[i] >= 0);
+            BOOST_CHECK(x[i] <= 1);
+        }
+    }
+}
+
 BOOST_AUTO_TEST_CASE(random_sampling)
 {
     std::cout << "RandomSampling" << std::endl;