random_generator.hpp 7.08 KB
Newer Older
1
2
3
4
//| 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
5
//|
6
7
8
//| Contributor(s):
//|   - Jean-Baptiste Mouret (jean-baptiste.mouret@inria.fr)
//|   - Antoine Cully (antoinecully@gmail.com)
9
//|   - Konstantinos Chatzilygeroudis (konstantinos.chatzilygeroudis@inria.fr)
10
11
//|   - Federico Allocati (fede.allocati@gmail.com)
//|   - Vaios Papaspyros (b.papaspyros@gmail.com)
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
12
//|   - Roberto Rama (bertoski@gmail.com)
13
//|
14
15
16
17
18
//| 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
19
//|
20
21
22
23
24
//| 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".
25
//|
26
27
28
29
30
//| 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.
31
//|
32
33
34
35
36
37
38
39
40
41
//| 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.
42
//|
43
44
//| 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.
45
//|
46
47
48
49
50

#ifndef LIMBO_TOOLS_RANDOM_GENERATOR_HPP
#define LIMBO_TOOLS_RANDOM_GENERATOR_HPP

#include <cmath>
51
#include <cstdlib>
52
#include <ctime>
53
#include <external/rand_utils.hpp>
54
#include <list>
55
#include <mutex>
56
#include <random>
57
#include <stdlib.h>
58
59
60
61
#include <utility>

namespace limbo {
    namespace tools {
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
62
63
64
65
        /// @ingroup tools
        /// a mt19937-based random generator (mutex-protected)
        ///
        /// usage :
66
        /// - RandomGenerator<dist<double>>(0.0, 1.0);
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
67
        /// - double r = rgen.rand();
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
68
        template <typename D>
69
70
        class RandomGenerator {
        public:
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
71
            using result_type = typename D::result_type;
72
73
74
75
76
            RandomGenerator(result_type a, result_type b, int seed = -1) : _dist(a, b) { this->seed(seed); }

            result_type rand() { return _dist(_rgen); }

            void seed(int seed = -1)
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
77
            {
78
79
80
81
                if (seed >= 0)
                    _rgen.seed(seed);
                else
                    _rgen.seed(randutils::auto_seed_128{}.base());
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
82
            }
83

84
85
86
87
            void reset() { _dist.reset(); }

            void param(const typename D::param_type& param) { _dist.param(param); }

88
        private:
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
89
90
            D _dist;
            std::mt19937 _rgen;
91
        };
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
92
93

        /// @ingroup tools
94
        using rdist_double_t = std::uniform_real_distribution<double>;
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
95
        /// @ingroup tools
96
        using rdist_int_t = std::uniform_int_distribution<int>;
97
98
        /// @ingroup tools
        using rdist_gauss_t = std::normal_distribution<>;
99

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
100
101
        /// @ingroup tools
        /// Double random number generator
102
        using rgen_double_t = RandomGenerator<rdist_double_t>;
103

104
105
106
107
        /// @ingroup tools
        /// Double random number generator (gaussian)
        using rgen_gauss_t = RandomGenerator<rdist_gauss_t>;

108
109
        /// @ingroup tools
        /// integer random number generator
110
        using rgen_int_t = RandomGenerator<rdist_int_t>;
111
112

        /// @ingroup tools
113
114
115
116
117
118
119
120
121
122
123
124
        /// random vector by providing custom RandomGenerator
        template <typename Rng>
        inline Eigen::VectorXd random_vec(int size, Rng& rng)
        {
            Eigen::VectorXd res(size);
            for (int i = 0; i < size; ++i)
                res[i] = rng.rand();
            return res;
        }

        /// @ingroup tools
        /// random vector in [0, 1]
125
        ///
126
        /// - this function is thread safe because we use a random generator for each thread
127
        /// - we use a C++11 random number generator
128
        inline Eigen::VectorXd random_vector_bounded(int size)
129
        {
130
            static thread_local rgen_double_t rgen(0.0, 1.0);
131
            return random_vec(size, rgen);
132
        }
133
134

        /// @ingroup tools
135
        /// random vector generated with a normal distribution centered on 0, with standard deviation of 10
136
        ///
137
        /// - this function is thread safe because we use a random generator for each thread
138
        /// - we use a C++11 random number generator
139
        inline Eigen::VectorXd random_vector_unbounded(int size)
140
        {
141
            static thread_local rgen_gauss_t rgen(0.0, 10.0);
142
            return random_vec(size, rgen);
143
144
145
146
        }

        /// @ingroup tools
        /// random vector wrapper for both bounded and unbounded versions
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
147
        inline Eigen::VectorXd random_vector(int size, bool bounded = true)
148
149
150
151
152
        {
            if (bounded)
                return random_vector_bounded(size);
            return random_vector_unbounded(size);
        }
153
154

        /// @ingroup tools
155
        /// generate n random samples with Latin Hypercube Sampling (LHS) in [0, 1]^dim
156
        inline Eigen::MatrixXd random_lhs(int dim, int n)
157
        {
158
159
            Eigen::VectorXd cut = Eigen::VectorXd::LinSpaced(n + 1, 0., 1.);
            Eigen::MatrixXd u = Eigen::MatrixXd::Zero(n, dim);
160

161
            for (int i = 0; i < n; i++) {
162
163
164
                u.row(i) = tools::random_vector(dim, true);
            }

165
166
            Eigen::VectorXd a = cut.head(n);
            Eigen::VectorXd b = cut.tail(n);
167

168
            Eigen::MatrixXd rdpoints = Eigen::MatrixXd::Zero(n, dim);
169
170
171
172
            for (int i = 0; i < dim; i++) {
                rdpoints.col(i) = u.col(i).array() * (b - a).array() + a.array();
            }

173
174
175
            Eigen::MatrixXd H = Eigen::MatrixXd::Zero(n, dim);
            Eigen::PermutationMatrix<Eigen::Dynamic, Eigen::Dynamic> perm(n);
            static thread_local std::mt19937 rgen(randutils::auto_seed_128{}.base());
176
177
            for (int i = 0; i < dim; i++) {
                perm.setIdentity();
178
                std::shuffle(perm.indices().data(), perm.indices().data() + perm.indices().size(), rgen);
179
180
181
182
183
184
185
186
                Eigen::MatrixXd tmp = perm * rdpoints;
                H.col(i) = tmp.col(i);
            }

            return H;
        }
    } // namespace tools
} // namespace limbo
187
188

#endif