diff --git a/docs/tutorials/gp.rst b/docs/tutorials/gp.rst index bab290bc66..62315984f3 100644 --- a/docs/tutorials/gp.rst +++ b/docs/tutorials/gp.rst @@ -13,7 +13,7 @@ We assume that our samples are in a vector called ``samples`` and that our obser .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 79-88 + :lines: 77-86 Basic usage ------------ @@ -23,14 +23,14 @@ We first create a basic GP with an Exponential kernel (``kernel::Exp``) .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 61-74 + :lines: 61-72 The type of the GP is defined by the following lines: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 89-93 + :lines: 87-91 To use the GP, we need : @@ -40,7 +40,7 @@ To use the GP, we need : .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 94-99 + :lines: 92-97 Here we assume that the noise is the same for all samples and that it is equal to 0.01. @@ -57,7 +57,7 @@ To visualize the predictions of the GP, we can query it for many points and reco .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 101-112 + :lines: 99-110 Hyper-parameter optimization @@ -71,7 +71,7 @@ A new GP type is defined as follows: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 114-118 + :lines: 112-116 It uses the default values for the parameters of ``SquaredExpARD``: @@ -85,7 +85,7 @@ After calling the ``compute()`` method, the hyper-parameters can be optimized by .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 121-123 + :lines: 119-121 We can have a look at the difference between the two GPs: @@ -115,7 +115,7 @@ We can also save our optimized GP model: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 140-141 + :lines: 138-139 This will create a directory called ``myGP`` with several files (the GP data, kernel hyperparameters etc.). If we want a binary format (i.e., more compact), we can replace the ``TextArchive`` by ``BinaryArchive``. @@ -124,6 +124,6 @@ To the load a saved model, we can do the following: .. literalinclude:: ../../src/tutorials/gp.cpp :language: c++ :linenos: - :lines: 143-144 + :lines: 141-142 Note that we need to have the same kernel and mean function (i.e., the same GP type) as the one used for saving. \ No newline at end of file diff --git a/src/examples/mono_dim.cpp b/src/examples/mono_dim.cpp index 88dc53d479..95b7846534 100644 --- a/src/examples/mono_dim.cpp +++ b/src/examples/mono_dim.cpp @@ -105,9 +105,6 @@ BO_PARAMS(std::cout, struct opt_rprop : public defaults::opt_rprop { }; - - struct opt_parallelrepeater : public defaults::opt_parallelrepeater { - }; };) struct fit_eval { diff --git a/src/examples/obs_multi_auto_mean.cpp b/src/examples/obs_multi_auto_mean.cpp index 95c8f036eb..a1d0f85093 100644 --- a/src/examples/obs_multi_auto_mean.cpp +++ b/src/examples/obs_multi_auto_mean.cpp @@ -93,9 +93,6 @@ struct Params { struct stop_maxiterations { BO_PARAM(int, iterations, 100); }; - - struct opt_parallelrepeater : defaults::opt_parallelrepeater { - }; }; template diff --git a/src/limbo/model.hpp b/src/limbo/model.hpp index a11dc31454..45c85a2fa3 100644 --- a/src/limbo/model.hpp +++ b/src/limbo/model.hpp @@ -50,6 +50,7 @@ ///@defgroup model_opt_defaults #include +#include #include #include diff --git a/src/limbo/model/gp.hpp b/src/limbo/model/gp.hpp index 4142f10ca9..4a4302b840 100644 --- a/src/limbo/model/gp.hpp +++ b/src/limbo/model/gp.hpp @@ -80,7 +80,7 @@ namespace limbo { /// useful because the model might be created before knowing anything about the process GP() : _dim_in(-1), _dim_out(-1), _inv_kernel_updated(false) {} - /// useful because the model might be created before having samples + /// useful because the model might be created before having samples GP(int dim_in, int dim_out) : _dim_in(dim_in), _dim_out(dim_out), _kernel_function(dim_in), _mean_function(dim_out), _inv_kernel_updated(false) {} @@ -153,7 +153,7 @@ namespace limbo { /** \\rst - return :math:`\mu`, :math:`\sigma^2` (unormalized). If there is no sample, return the value according to the mean function. Using this method instead of separate calls to mu() and sigma() is more efficient because some computations are shared between mu() and sigma(). + return :math:`\mu`, :math:`\sigma^2` (un-normalized). If there is no sample, return the value according to the mean function. Using this method instead of separate calls to mu() and sigma() is more efficient because some computations are shared between mu() and sigma(). \\endrst */ std::tuple query(const Eigen::VectorXd& v) const @@ -193,14 +193,14 @@ namespace limbo { /// return the number of dimensions of the input int dim_in() const { - assert(_dim_in != -1); // need to compute first ! + assert(_dim_in != -1); // need to compute first! return _dim_in; } /// return the number of dimensions of the output int dim_out() const { - assert(_dim_out != -1); // need to compute first ! + assert(_dim_out != -1); // need to compute first! return _dim_out; } diff --git a/src/limbo/model/gp/hp_opt.hpp b/src/limbo/model/gp/hp_opt.hpp index f22c3a6be3..b0e7e16aef 100644 --- a/src/limbo/model/gp/hp_opt.hpp +++ b/src/limbo/model/gp/hp_opt.hpp @@ -48,7 +48,6 @@ #include -#include #include namespace limbo { @@ -56,7 +55,7 @@ namespace limbo { namespace gp { ///@ingroup model_opt ///base class for optimization of the hyper-parameters of a GP - template >> + template > struct HPOpt { public: HPOpt() : _called(false) {} diff --git a/src/limbo/model/gp/kernel_lf_opt.hpp b/src/limbo/model/gp/kernel_lf_opt.hpp index 878c4cd358..75ce44c534 100644 --- a/src/limbo/model/gp/kernel_lf_opt.hpp +++ b/src/limbo/model/gp/kernel_lf_opt.hpp @@ -47,14 +47,13 @@ #define LIMBO_MODEL_GP_KERNEL_LF_OPT_HPP #include -#include namespace limbo { namespace model { namespace gp { ///@ingroup model_opt ///optimize the likelihood of the kernel only - template >> + template > struct KernelLFOpt : public HPOpt { public: template @@ -96,8 +95,8 @@ namespace limbo { const GP& _original_gp; }; }; - } - } -} + } // namespace gp + } // namespace model +} // namespace limbo #endif diff --git a/src/limbo/model/gp/kernel_loo_opt.hpp b/src/limbo/model/gp/kernel_loo_opt.hpp index 59078e5fd6..cb1db84c42 100644 --- a/src/limbo/model/gp/kernel_loo_opt.hpp +++ b/src/limbo/model/gp/kernel_loo_opt.hpp @@ -47,14 +47,13 @@ #define LIMBO_MODEL_GP_KERNEL_LOO_OPT_HPP #include -#include namespace limbo { namespace model { namespace gp { ///@ingroup model_opt ///optimize the likelihood of the kernel only - template >> + template > struct KernelLooOpt : public HPOpt { public: template @@ -96,8 +95,8 @@ namespace limbo { const GP& _original_gp; }; }; - } - } -} + } // namespace gp + } // namespace model +} // namespace limbo #endif diff --git a/src/limbo/model/gp/kernel_mean_lf_opt.hpp b/src/limbo/model/gp/kernel_mean_lf_opt.hpp index 4889a76975..dcf37e4159 100644 --- a/src/limbo/model/gp/kernel_mean_lf_opt.hpp +++ b/src/limbo/model/gp/kernel_mean_lf_opt.hpp @@ -47,14 +47,13 @@ #define LIMBO_MODEL_GP_KERNEL_MEAN_LF_OPT_HPP #include -#include namespace limbo { namespace model { namespace gp { ///@ingroup model_opt ///optimize the likelihood of both the kernel and the mean (try to align the mean function) - template >> + template > struct KernelMeanLFOpt : public HPOpt { public: template @@ -109,8 +108,8 @@ namespace limbo { const GP& _original_gp; }; }; - } - } -} + } // namespace gp + } // namespace model +} // namespace limbo #endif diff --git a/src/limbo/model/gp/mean_lf_opt.hpp b/src/limbo/model/gp/mean_lf_opt.hpp index b6a0032644..9b0e4dfccb 100644 --- a/src/limbo/model/gp/mean_lf_opt.hpp +++ b/src/limbo/model/gp/mean_lf_opt.hpp @@ -47,14 +47,13 @@ #define LIMBO_MODEL_GP_MEAN_LF_OPT_HPP #include -#include namespace limbo { namespace model { namespace gp { ///@ingroup model_opt ///optimize the likelihood of the mean only (try to align the mean function) - template >> + template > struct MeanLFOpt : public HPOpt { public: template @@ -99,8 +98,8 @@ namespace limbo { GP _original_gp; }; }; - } - } -} + } // namespace gp + } // namespace model +} // namespace limbo #endif diff --git a/src/limbo/model/multi_gp.hpp b/src/limbo/model/multi_gp.hpp new file mode 100644 index 0000000000..19029fc198 --- /dev/null +++ b/src/limbo/model/multi_gp.hpp @@ -0,0 +1,285 @@ +//| 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_MODEL_MULTI_GP_HPP +#define LIMBO_MODEL_MULTI_GP_HPP + +#include + +namespace limbo { + namespace model { + /// @ingroup model + /// A wrapper for N-output Gaussian processes. + /// It is parametrized by: + /// - GP class + /// - a kernel function (the same type for all GPs, but can have different parameters) + /// - a mean function (the same type and parameters for all GPs) + /// - [optional] an optimizer for the hyper-parameters + template class GPClass, typename KernelFunction, typename MeanFunction, class HyperParamsOptimizer = limbo::model::gp::NoLFOpt> + class MultiGP { + public: + using GP_t = GPClass, limbo::model::gp::NoLFOpt>; + + /// useful because the model might be created before knowing anything about the process + MultiGP() : _dim_in(-1), _dim_out(-1) {} + + /// useful because the model might be created before having samples + MultiGP(int dim_in, int dim_out) + : _dim_in(dim_in), _dim_out(dim_out), _mean_function(dim_out) + { + // initialize dim_in models with 1 output + _gp_models.resize(_dim_out); + for (int i = 0; i < _dim_out; i++) { + _gp_models[i] = GP_t(_dim_in, 1); + } + } + + /// Compute the GP from samples and observations. This call needs to be explicit! + void compute(const std::vector& samples, + const std::vector& observations, bool compute_kernel = true) + { + assert(samples.size() != 0); + assert(observations.size() != 0); + assert(samples.size() == observations.size()); + + if (_dim_in != samples[0].size()) { + _dim_in = samples[0].size(); + } + + if (_dim_out != observations[0].size()) { + _dim_out = observations[0].size(); + _mean_function = MeanFunction(_dim_out); // the cost of building a functor should be relatively low + } + + if ((int)_gp_models.size() != _dim_out) { + _gp_models.resize(_dim_out); + for (int i = 0; i < _dim_out; i++) + _gp_models[i] = GP_t(_dim_in, 1); + } + + // save observations + // TO-DO: Check how can we improve for not saving observations twice (one here and one for each GP)!? + _observations = observations; + + // compute the new observations for the GPs + std::vector> obs(_dim_out); + + for (size_t j = 0; j < observations.size(); j++) { + Eigen::VectorXd mean_vector = _mean_function(samples[j], *this); + assert(mean_vector.size() == _dim_out); + for (int i = 0; i < _dim_out; i++) { + obs[i].push_back(limbo::tools::make_vector(observations[j][i] - mean_vector[i])); + } + } + + // do the actual computation + limbo::tools::par::loop(0, _dim_out, [&](size_t i) { + _gp_models[i].compute(samples, obs[i], compute_kernel); + }); + } + + /// Do not forget to call this if you use hyper-parameters optimization!! + void optimize_hyperparams() + { + _hp_optimize(*this); + } + + const MeanFunction& mean_function() const { return _mean_function; } + + MeanFunction& mean_function() { return _mean_function; } + + /// add sample and update the GPs. This code uses an incremental implementation of the Cholesky + /// decomposition. It is therefore much faster than a call to compute() + void add_sample(const Eigen::VectorXd& sample, const Eigen::VectorXd& observation) + { + if (_gp_models.size() == 0) { + if (_dim_in != sample.size()) { + _dim_in = sample.size(); + } + if (_dim_out != observation.size()) { + _dim_out = observation.size(); + _gp_models.resize(_dim_out); + for (int i = 0; i < _dim_out; i++) + _gp_models[i] = GP_t(_dim_in, 1); + + _mean_function = MeanFunction(_dim_out); // the cost of building a functor should be relatively low + } + } + else { + assert(sample.size() == _dim_in); + assert(observation.size() == _dim_out); + } + + _observations.push_back(observation); + + Eigen::VectorXd mean_vector = _mean_function(sample, *this); + assert(mean_vector.size() == _dim_out); + + limbo::tools::par::loop(0, _dim_out, [&](size_t i) { + _gp_models[i].add_sample(sample, limbo::tools::make_vector(observation[i] - mean_vector[i])); + }); + } + + /** + \\rst + return :math:`\mu`, :math:`\sigma^2` (un-normalized; this will return a vector --- one for each GP). Using this method instead of separate calls to mu() and sigma() is more efficient because some computations are shared between mu() and sigma(). + \\endrst + */ + std::tuple query(const Eigen::VectorXd& v) const + { + Eigen::VectorXd mu(_dim_out); + Eigen::VectorXd sigma(_dim_out); + + // query the mean function + Eigen::VectorXd mean_vector = _mean_function(v, *this); + + // parallel query of the GPs + limbo::tools::par::loop(0, _dim_out, [&](size_t i) { + Eigen::VectorXd tmp; + std::tie(tmp, sigma(i)) = _gp_models[i].query(v); + mu(i) = tmp(0) + mean_vector(i); + }); + + return std::make_tuple(mu, sigma); + } + + /** + \\rst + return :math:`\mu` (un-normalized). If there is no sample, return the value according to the mean function. + \\endrst + */ + Eigen::VectorXd mu(const Eigen::VectorXd& v) const + { + Eigen::VectorXd mu(_dim_out); + Eigen::VectorXd mean_vector = _mean_function(v, *this); + + limbo::tools::par::loop(0, _dim_out, [&](size_t i) { + mu(i) = _gp_models[i].mu(v)[0] + mean_vector(i); + }); + + return mu; + } + + /** + \\rst + return :math:`\sigma^2` (un-normalized). This returns a vector; one value for each GP. + \\endrst + */ + Eigen::VectorXd sigma(const Eigen::VectorXd& v) const + { + Eigen::VectorXd sigma(_dim_out); + + limbo::tools::par::loop(0, _dim_out, [&](size_t i) { + sigma(i) = _gp_models[i].sigma(v); + }); + + return sigma; + } + + /// return the number of dimensions of the input + int dim_in() const + { + assert(_dim_in != -1); // need to compute first! + return _dim_in; + } + + /// return the number of dimensions of the output + int dim_out() const + { + assert(_dim_out != -1); // need to compute first! + return _dim_out; + } + + /// return the number of samples used to compute the GP + int nb_samples() const + { + return (_gp_models.size() > 0) ? _gp_models[0].nb_samples() : 0; + } + + /// recomputes the GPs + void recompute(bool update_obs_mean = true, bool update_full_kernel = true) + { + // if there are no GPs, there's nothing to recompute + if (_gp_models.size() == 0) + return; + + if (update_obs_mean) // if the mean is updated, we need to fully re-compute + return compute(_gp_models[0].samples(), _observations, update_full_kernel); + else + limbo::tools::par::loop(0, _dim_out, [&](size_t i) { + _gp_models[i].recompute(false, update_full_kernel); + }); + } + + /// return the list of samples that have been tested so far + const std::vector& samples() const + { + assert(_gp_models.size()); + return _gp_models[0].samples(); + } + + /// return the list of GPs + std::vector gp_models() const + { + return _gp_models; + } + + /// return the list of GPs + std::vector& gp_models() + { + return _gp_models; + } + + protected: + std::vector _gp_models; + int _dim_in, _dim_out; + HyperParamsOptimizer _hp_optimize; + MeanFunction _mean_function; + std::vector _observations; + }; + } // namespace model +} // namespace limbo + +#endif \ No newline at end of file diff --git a/src/limbo/model/multi_gp/parallel_lf_opt.hpp b/src/limbo/model/multi_gp/parallel_lf_opt.hpp new file mode 100644 index 0000000000..093f13e1a5 --- /dev/null +++ b/src/limbo/model/multi_gp/parallel_lf_opt.hpp @@ -0,0 +1,75 @@ +//| 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_MODEL_MULTI_GP_PARALLEL_LF_OPT_HPP +#define LIMBO_MODEL_MULTI_GP_PARALLEL_LF_OPT_HPP + +#include + +namespace limbo { + namespace model { + namespace multi_gp { + ///@ingroup model_opt + ///optimize each GP independently in parallel using HyperParamsOptimizer + template > + struct ParallelLFOpt : public limbo::model::gp::HPOpt { + public: + template + void operator()(GP& gp) + { + this->_called = true; + auto& gps = gp.gp_models(); + limbo::tools::par::loop(0, gps.size(), [&](size_t i) { + HyperParamsOptimizer hp_optimize; + hp_optimize(gps[i]); + }); + } + }; + + } // namespace multi_gp + } // namespace model +} // namespace limbo + +#endif \ No newline at end of file diff --git a/src/tests/test_gp.cpp b/src/tests/test_gp.cpp index 4e57bdc94b..e596dbcc74 100644 --- a/src/tests/test_gp.cpp +++ b/src/tests/test_gp.cpp @@ -61,6 +61,8 @@ #include #include #include +#include +#include #include #include #include @@ -119,9 +121,6 @@ struct Params { struct opt_rprop : public defaults::opt_rprop { }; - struct opt_parallelrepeater : public defaults::opt_parallelrepeater { - }; - struct acqui_ucb : public defaults::acqui_ucb { }; @@ -215,9 +214,6 @@ BOOST_AUTO_TEST_CASE(test_gp_check_lf_grad_noise) struct opt_rprop : public defaults::opt_rprop { }; - struct opt_parallelrepeater : public defaults::opt_parallelrepeater { - }; - struct acqui_ucb : public defaults::acqui_ucb { }; @@ -338,9 +334,6 @@ BOOST_AUTO_TEST_CASE(test_gp_check_loo_grad_noise) struct opt_rprop : public defaults::opt_rprop { }; - struct opt_parallelrepeater : public defaults::opt_parallelrepeater { - }; - struct acqui_ucb : public defaults::acqui_ucb { }; @@ -915,3 +908,249 @@ BOOST_AUTO_TEST_CASE(test_sparse_gp_accuracy) BOOST_CHECK(double(failures) / double(N) < 0.1); } + +BOOST_AUTO_TEST_CASE(test_multi_gp_dim) +{ + using namespace limbo; + + using KF_t = kernel::MaternFiveHalves; + using Mean_t = mean::Constant; + using GP_t = model::GP; + + GP_t gp; // no init with dim + + std::vector observations = {make_v2(5, 5), make_v2(10, 10), + make_v2(5, 5)}; + std::vector samples = {make_v2(1, 1), make_v2(2, 2), make_v2(3, 3)}; + + gp.compute(samples, observations); + + using MultiGP_t = model::MultiGP; + MultiGP_t multi_gp; + multi_gp.compute(samples, observations); + + Eigen::VectorXd mu; + double sigma; + std::tie(mu, sigma) = gp.query(make_v2(1, 1)); + BOOST_CHECK(std::abs((mu(0) - 5)) < 1); + BOOST_CHECK(std::abs((mu(1) - 5)) < 1); + + BOOST_CHECK(sigma <= 2. * (Params::kernel::noise() + 1e-8)); + + Eigen::VectorXd mu_multi, sigma_multi; + + std::tie(mu_multi, sigma_multi) = multi_gp.query(make_v2(1, 1)); + BOOST_CHECK(std::abs((mu_multi(0) - 5)) < 1); + BOOST_CHECK(std::abs((mu_multi(1) - 5)) < 1); + + BOOST_CHECK(sigma_multi(0) <= 2. * (Params::kernel::noise() + 1e-8)); + BOOST_CHECK(sigma_multi(1) <= 2. * (Params::kernel::noise() + 1e-8)); + + BOOST_CHECK(std::abs(mu_multi(0) - mu(0)) < 1e-6); + BOOST_CHECK(std::abs(mu_multi(1) - mu(1)) < 1e-6); + BOOST_CHECK(std::abs(sigma_multi(0) - sigma) < 1e-6); + BOOST_CHECK(std::abs(sigma_multi(1) - sigma) < 1e-6); +} + +BOOST_AUTO_TEST_CASE(test_multi_gp) +{ + using namespace limbo; + + using KF_t = kernel::MaternFiveHalves; + using Mean_t = mean::Constant; + using MultiGP_t = model::MultiGP; + + MultiGP_t gp; + std::vector observations = {make_v1(5), make_v1(10), + make_v1(5)}; + std::vector samples = {make_v1(1), make_v1(2), make_v1(3)}; + + gp.compute(samples, observations); + + Eigen::VectorXd mu, sigma; + std::tie(mu, sigma) = gp.query(make_v1(1)); + BOOST_CHECK(std::abs((mu(0) - 5)) < 1); + BOOST_CHECK(sigma(0) <= 2. * (Params::kernel::noise() + 1e-8)); + + std::tie(mu, sigma) = gp.query(make_v1(2)); + BOOST_CHECK(std::abs((mu(0) - 10)) < 1); + BOOST_CHECK(sigma(0) <= 2. * (Params::kernel::noise() + 1e-8)); + + std::tie(mu, sigma) = gp.query(make_v1(3)); + BOOST_CHECK(std::abs((mu(0) - 5)) < 1); + BOOST_CHECK(sigma(0) <= 2. * (Params::kernel::noise() + 1e-8)); + + for (double x = 0; x < 4; x += 0.05) { + Eigen::VectorXd mu, sigma; + std::tie(mu, sigma) = gp.query(make_v1(x)); + BOOST_CHECK(gp.mu(make_v1(x)) == mu); + BOOST_CHECK(gp.sigma(make_v1(x)) == sigma); + } +} + +BOOST_AUTO_TEST_CASE(test_sparse_multi_gp) +{ + using namespace limbo; + constexpr size_t N = 20; + constexpr size_t M = 100; + size_t failures = 0; + + struct SparseParams { + struct mean_constant : public defaults::mean_constant { + }; + struct kernel : public defaults::kernel { + }; + struct kernel_squared_exp_ard : public defaults::kernel_squared_exp_ard { + }; + struct opt_rprop : public defaults::opt_rprop { + }; + struct model_sparse_gp { + BO_PARAM(int, max_points, M / 2); + }; + }; + + using KF_t = kernel::SquaredExpARD; + using MF_t = mean::Constant; + using MultiGP_t = model::MultiGP>>; + using MultiSparseGP_t = model::MultiGP>>; + + for (size_t i = 0; i < N; i++) { + + std::vector observations; + std::vector samples; + tools::rgen_double_t rgen(-2., 2.); + for (size_t i = 0; i < M; i++) { + samples.push_back(make_v1(rgen.rand())); + observations.push_back(make_v1(std::cos(samples.back()[0]))); + } + + std::vector test_observations; + std::vector test_samples; + for (size_t i = 0; i < M; i++) { + test_samples.push_back(make_v1(rgen.rand())); + test_observations.push_back(make_v1(std::cos(test_samples.back()[0]))); + } + + MultiGP_t gp; + gp.compute(samples, observations, false); + gp.optimize_hyperparams(); + + MultiSparseGP_t sgp; + sgp.compute(samples, observations, false); + sgp.optimize_hyperparams(); + + bool failed = false; + + // check if normal GP and sparse GP produce very similar results + // in the learned points + for (size_t i = 0; i < M; i++) { + Eigen::VectorXd gp_val, sgp_val, gp_sigma, sgp_sigma; + std::tie(gp_val, gp_sigma) = gp.query(samples[i]); + std::tie(sgp_val, sgp_sigma) = sgp.query(samples[i]); + + if (std::abs(gp_val[0] - sgp_val[0]) > 1e-2 || std::abs(gp_sigma[0] - sgp_sigma[0]) > 1e-2) + failed = true; + } + + // check if normal GP and sparse GP produce very similar results + // in the test points + for (size_t i = 0; i < M; i++) { + Eigen::VectorXd gp_val, sgp_val, gp_sigma, sgp_sigma; + std::tie(gp_val, gp_sigma) = gp.query(test_samples[i]); + std::tie(sgp_val, sgp_sigma) = sgp.query(test_samples[i]); + + if (std::abs(gp_val[0] - sgp_val[0]) > 1e-2 || std::abs(gp_sigma[0] - sgp_sigma[0]) > 1e-2) + failed = true; + } + + // check if normal GP and sparse GP produce very similar errors + // in the test points + for (size_t i = 0; i < M; i++) { + double gp_val = gp.mu(test_samples[i])[0]; + double sgp_val = sgp.mu(test_samples[i])[0]; + + double gp_error_val = std::abs(gp_val - test_observations[i][0]); + double sgp_error_val = std::abs(sgp_val - test_observations[i][0]); + + if (std::abs(gp_error_val - sgp_error_val) > 1e-2) + failed = true; + } + + if (failed) + failures++; + } + + BOOST_CHECK(double(failures) / double(N) < 0.1); +} + +BOOST_AUTO_TEST_CASE(test_multi_gp_auto) +{ + using KF_t = kernel::SquaredExpARD; + using Mean_t = mean::Constant; + using MultiGP_t = model::MultiGP>>; + + MultiGP_t gp; + std::vector observations = {make_v2(5, 5), make_v2(10, 10), make_v2(5, 5)}; + std::vector samples = {make_v1(1), make_v1(2), make_v1(3)}; + + gp.compute(samples, observations, false); + gp.optimize_hyperparams(); + + Eigen::VectorXd mu, sigma; + std::tie(mu, sigma) = gp.query(make_v1(1)); + BOOST_CHECK(std::abs((mu(0) - 5)) < 1); + BOOST_CHECK(sigma(0) <= 2. * (gp.gp_models()[0].kernel_function().noise() + 1e-8)); + BOOST_CHECK(sigma(1) <= 2. * (gp.gp_models()[1].kernel_function().noise() + 1e-8)); + + std::tie(mu, sigma) = gp.query(make_v1(2)); + BOOST_CHECK(std::abs((mu(0) - 10)) < 1); + BOOST_CHECK(sigma(0) <= 2. * (gp.gp_models()[0].kernel_function().noise() + 1e-8)); + BOOST_CHECK(sigma(1) <= 2. * (gp.gp_models()[1].kernel_function().noise() + 1e-8)); + + std::tie(mu, sigma) = gp.query(make_v1(3)); + BOOST_CHECK(std::abs((mu(0) - 5)) < 1); + BOOST_CHECK(sigma(0) <= 2. * (gp.gp_models()[0].kernel_function().noise() + 1e-8)); + BOOST_CHECK(sigma(1) <= 2. * (gp.gp_models()[1].kernel_function().noise() + 1e-8)); +} + +BOOST_AUTO_TEST_CASE(test_multi_gp_recompute) +{ + using KF_t = kernel::SquaredExpARD; + using Mean_t = mean::Constant; + using MultiGP_t = model::MultiGP; + + MultiGP_t gp; + + gp.add_sample(make_v2(1, 1), make_v1(2)); + gp.add_sample(make_v2(2, 2), make_v1(10)); + + // make sure that the observations are properly passed + BOOST_CHECK(gp._observations[0](0) == 2.); + BOOST_CHECK(gp._observations[1](0) == 10.); + + // make sure the child GPs have the proper observations + BOOST_CHECK(gp.gp_models()[0]._observations.row(0)[0] == (2. - gp.mean_function().h_params()[0])); + BOOST_CHECK(gp.gp_models()[0]._observations.row(1)[0] == (10. - gp.mean_function().h_params()[0])); + + // now change the mean function parameters + gp.mean_function().set_h_params(make_v1(2)); + + // make sure that the observations are properly passed + BOOST_CHECK(gp._observations[0](0) == 2.); + BOOST_CHECK(gp._observations[1](0) == 10.); + + // make sure the child GPs do not have the proper observations + BOOST_CHECK(!(gp.gp_models()[0]._observations.row(0)[0] == (2. - gp.mean_function().h_params()[0]))); + BOOST_CHECK(!(gp.gp_models()[0]._observations.row(1)[0] == (10. - gp.mean_function().h_params()[0]))); + + // recompute the GP + gp.recompute(); + + // make sure that the observations are properly passed + BOOST_CHECK(gp._observations[0](0) == 2.); + BOOST_CHECK(gp._observations[1](0) == 10.); + + // make sure the child GPs have the proper observations + BOOST_CHECK(gp.gp_models()[0]._observations.row(0)[0] == (2. - gp.mean_function().h_params()[0])); + BOOST_CHECK(gp.gp_models()[0]._observations.row(1)[0] == (10. - gp.mean_function().h_params()[0])); +} diff --git a/src/tests/test_serialize.cpp b/src/tests/test_serialize.cpp index d6e8ccbc20..0f94fa4f8f 100644 --- a/src/tests/test_serialize.cpp +++ b/src/tests/test_serialize.cpp @@ -69,8 +69,6 @@ struct Params { }; struct opt_rprop : public limbo::defaults::opt_rprop { }; - struct opt_parallelrepeater : public limbo::defaults::opt_parallelrepeater { - }; struct kernel_maternfivehalves { BO_PARAM(double, sigma_sq, 1); @@ -95,8 +93,6 @@ struct LoadParams { }; struct opt_rprop : public limbo::defaults::opt_rprop { }; - struct opt_parallelrepeater : public limbo::defaults::opt_parallelrepeater { - }; struct kernel_maternfivehalves { BO_PARAM(double, sigma_sq, 2.); diff --git a/src/tutorials/gp.cpp b/src/tutorials/gp.cpp index d9a6b9a49b..605472ab97 100644 --- a/src/tutorials/gp.cpp +++ b/src/tutorials/gp.cpp @@ -69,8 +69,6 @@ struct Params { }; struct opt_rprop : public defaults::opt_rprop { }; - struct opt_parallelrepeater : public defaults::opt_parallelrepeater { - }; }; int main(int argc, char** argv)