Latin Hypercube Sampling

Here is Omar's implementation (to check!).
#ifndef INIT_FUNCTIONS_HPP_
#define INIT_FUNCTIONS_HPP_

#include <vector>
#include <iostream>
#include <Eigen/Core>
#include <typeinfo>

#include <limits>

#include <stdlib.h>     /* srand, rand */
#include <fstream>
#include <string>       // std::string
#include <iostream>     // std::cout
#include <sstream>      // std::stringstream

#include <algorithm>    // std::shuffle
#include <random>       // std::default_random_engine
#include <chrono>       // std::chrono::system_clock

namespace limbo {
  namespace init_functions {

    // params is here only to make it easy to switch
    // from/to the other init functions
    template<typename Params>
    struct NoInit {
      template<typename F, typename Opt>
      void operator()(const F& f, Opt& opt) const {}
    };

    // initialize in [0,1] !
    // params: init::nb_samples
    template<typename Params>
    struct RandomSampling {
      template<typename F, typename Opt>
      void operator()(const F& feval, Opt& opt) const {
        for (int i = 0; i < Params::init::nb_samples(); i++) {
          Eigen::VectorXd new_sample(F::dim);
          for (size_t i = 0; i < F::dim; i++)
            ///////////////////////////////////////////////////////////////////////////////
            new_sample[i] = misc::rand<double>(0, 1); //Check the type of the random function here!
            // I checked the rand.hpp library. This function generates a "uniform" random points based
            // on the built-in function rand from random library I guess.
            ///////////////////////////////////////////////////////////////////////////////
          std::cout << "random sample:" << new_sample.transpose() << std::endl;
          opt.add_new_sample(new_sample, feval(new_sample));
          std::cout << "Current sample number =" << i << std::endl;
        }
      }
    };

    template<typename Params>
    struct NonFillingLHS {
      template<typename F, typename Opt>
      void operator()(const F& feval, Opt& opt) const {
        double step = (0.0 + 1.0) / float(Params::init::nb_samples());
        for (double point_iter = 0; point_iter < (1.0 - step) ; point_iter = point_iter + step)
        {
          Eigen::VectorXd new_sample(F::dim);
          for (int dim_iter = 0; dim_iter < F::dim; dim_iter++){
            // std::cout << "Min = " << point_iter << ", max = " << point_iter+step << std::endl;
            new_sample[dim_iter] = misc::rand<double>(point_iter, point_iter+step);
          }
          std::cout << "NonFillingLHS -- random sample:" << new_sample.transpose() << std::endl;
          opt.add_new_sample(new_sample, feval(new_sample));
        }
      }
    };

    template<typename Params>
    struct FillingLHS {
      template<typename F, typename Opt>
      void operator()(const F& feval, Opt& opt) const {
        /*
        This is a corrected implementation for LHS sampling based on the book 'Design and analysis of computer experiments',
        page 128. I assume here that the all dimensions have similar continuous range.
        */
        int problem_dim = F::dim;
        Eigen::MatrixXd random_points(problem_dim, Params::init::nb_samples()); // This is an NxD matrix, where N is the number of samples, D is the number of dimensions
        Eigen::VectorXd new_sample(problem_dim);

        // Generate an array of length Params::init::nb_samples()
        std::vector<int> samples_index_list(Params::init::nb_samples());
        for (int i = 0; i < Params::init::nb_samples(); i++){
          samples_index_list[i] = i + 1;
        }
        // Generate elements of the nXd matrix
        auto engine = std::default_random_engine{};
        for (int row_iter = 0; row_iter < problem_dim; row_iter++){
          std::shuffle (std::begin(samples_index_list), std::end(samples_index_list), engine);
          for (int col_iter = 0; col_iter < Params::init::nb_samples(); col_iter++){
            random_points(row_iter,col_iter) = samples_index_list[col_iter];
          }
        }
        // Get the matrix transpose
        random_points = random_points.transpose().eval(); // Now, its dimensions are nb_samples * F::dim
        // Make the random pairings according to the number of points you want
        for (int sample_index = 0; sample_index < Params::init::nb_samples(); sample_index++){
          for (int dim_index = 0; dim_index < problem_dim; dim_index++){
            double random_U = misc::rand<double>(0, 1); //Check the type of the random function here!
            double random_point = (random_points(sample_index, dim_index) - 1.0 + random_U) / Params::init::nb_samples();
            new_sample(dim_index) = random_point;
          }
          std::cout << "FillingLHS -- random sample:" << new_sample.transpose() << std::endl;
          opt.add_new_sample(new_sample, feval(new_sample));
        }
        exit (EXIT_FAILURE);
      }
    };
    // params:
    //  -init::nb_bins
    //  - init::nb_samples
    template<typename Params>
    struct RandomSamplingGrid {
      template<typename F, typename Opt>
      void operator()(const F& feval, Opt& opt) const {
        for (int i = 0; i < Params::init::nb_samples(); i++) {
          Eigen::VectorXd new_sample(F::dim);
          for (size_t i = 0; i < F::dim; i++)
            new_sample[i] =
              int(((double) (Params::init::nb_bins() + 1) * rand())
                  / (RAND_MAX + 1.0)) / double(Params::init::nb_bins());
          opt.add_new_sample(new_sample, feval(new_sample));
        }
      }
    };

    // params:
    //  -init::nb_bins
    template<typename Params>
    struct GridSampling {
      template<typename F, typename Opt>
      void operator()(const F& feval, Opt& opt) const {
        _explore(0, feval, Eigen::VectorXd::Constant(F::dim, 0), opt);
      }
     private:
      //recursively explore all the dimensions
      template<typename F, typename Opt>
      void _explore(int dim, const F& feval, const Eigen::VectorXd& current, Opt& opt) const {
        for (double x = 0; x <= 1.0f; x += 1.0f / (double)Params::init::nb_bins()) {
          Eigen::VectorXd point = current;
          point[dim] = x;
          if (dim == current.size() - 1) {
            opt.add_new_sample(point, feval(point));
          } else {
            _explore(dim + 1, feval, point, opt);
          }
        }
      }
    };
  }
}
#endif
resibots / limbo

Latin Hypercube Sampling #34
