random_generator.hpp 6.62 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
            RandomGenerator(result_type a, result_type b) : _dist(a, b), _rgen(randutils::auto_seed_128{}.base()) {}
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
73
74
75
76
            result_type rand()
            {
                return _dist(_rgen);
            }
77
        private:
Konstantinos Chatzilygeroudis's avatar
Konstantinos Chatzilygeroudis committed
78
79
            D _dist;
            std::mt19937 _rgen;
80
        };
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
81
82

        /// @ingroup tools
83
        using rdist_double_t = std::uniform_real_distribution<double>;
Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
84
        /// @ingroup tools
85
        using rdist_int_t = std::uniform_int_distribution<int>;
86
87
        /// @ingroup tools
        using rdist_gauss_t = std::normal_distribution<>;
88

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
89
90
        /// @ingroup tools
        /// Double random number generator
91
        using rgen_double_t = RandomGenerator<rdist_double_t>;
92

93
94
95
96
        /// @ingroup tools
        /// Double random number generator (gaussian)
        using rgen_gauss_t = RandomGenerator<rdist_gauss_t>;

Jean-Baptiste Mouret's avatar
Jean-Baptiste Mouret committed
97
98
        ///@ingroup tools
        ///integer random number generator
99
        using rgen_int_t = RandomGenerator<rdist_int_t>;
100
101

        /// @ingroup tools
102
        /// random vector in [0, 1.0]
103
        ///
104
        /// - this function is thread safe because we use a random generator for each thread
105
        /// - we use a C++11 random number generator
106
        Eigen::VectorXd random_vector_bounded(int size)
107
        {
108
            static thread_local rgen_double_t rgen(0.0, 1.0);
109
110
            Eigen::VectorXd res(size);
            for (int i = 0; i < size; ++i)
111
                res[i] = rgen.rand();
112
113
            return res;
        }
114
115

        /// @ingroup tools
116
        /// random vector generated with a normal distribution centered on 0, with standard deviation of 10.0
117
        ///
118
        /// - this function is thread safe because we use a random generator for each thread
119
120
121
        /// - we use a C++11 random number generator
        Eigen::VectorXd random_vector_unbounded(int size)
        {
122
            static thread_local rgen_gauss_t rgen(0.0, 10.0);
123
124
125
126
127
128
129
130
131
132
133
134
135
136
            Eigen::VectorXd res(size);
            for (int i = 0; i < size; ++i)
                res[i] = rgen.rand();
            return res;
        }

        /// @ingroup tools
        /// random vector wrapper for both bounded and unbounded versions
        Eigen::VectorXd random_vector(int size, bool bounded = true)
        {
            if (bounded)
                return random_vector_bounded(size);
            return random_vector_unbounded(size);
        }
137
138

        /// @ingroup tools
139
140
        /// generate n random samples with Latin Hypercube Sampling (LHS) in [0, 1]^dim
        Eigen::MatrixXd random_lhs(int dim, int n)
141
        {
142
143
            Eigen::VectorXd cut = Eigen::VectorXd::LinSpaced(n + 1, 0., 1.);
            Eigen::MatrixXd u = Eigen::MatrixXd::Zero(n, dim);
144

145
            for (int i = 0; i < n; i++) {
146
147
148
                u.row(i) = tools::random_vector(dim, true);
            }

149
150
            Eigen::VectorXd a = cut.head(n);
            Eigen::VectorXd b = cut.tail(n);
151

152
            Eigen::MatrixXd rdpoints = Eigen::MatrixXd::Zero(n, dim);
153
154
155
156
            for (int i = 0; i < dim; i++) {
                rdpoints.col(i) = u.col(i).array() * (b - a).array() + a.array();
            }

157
158
159
            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());
160
161
            for (int i = 0; i < dim; i++) {
                perm.setIdentity();
162
                std::shuffle(perm.indices().data(), perm.indices().data() + perm.indices().size(), rgen);
163
164
165
166
167
168
169
170
                Eigen::MatrixXd tmp = perm * rdpoints;
                H.col(i) = tmp.col(i);
            }

            return H;
        }
    } // namespace tools
} // namespace limbo
171
172

#endif