Better MeanFunctionARD gradients and fixes in tests

src/limbo/mean/function_ard.hpp

@@ -53,40 +53,59 @@ namespace limbo {
 
         /// Functor used to optimize the mean function using the maximum likelihood principle
         ///
-        /// For the moment, it ignores the hyperparameters of the underlying mean function, if any
+        /// It incorporates the hyperparameters of the underlying mean function, if any
         /// @see limbo::model::gp::KernelMeanLFOpt, limbo::model::gp::MeanLFOpt
         template <typename Params, typename MeanFunction>
         struct FunctionARD : public BaseMean<Params> {
             FunctionARD(size_t dim_out = 1)
                 : _mean_function(dim_out), _tr(dim_out, dim_out + 1)
             {
-                Eigen::VectorXd h = Eigen::VectorXd::Zero(dim_out * (dim_out + 1));
+                Eigen::VectorXd h = Eigen::VectorXd::Zero(dim_out * (dim_out + 1) + _mean_function.h_params_size());
                 for (size_t i = 0; i < dim_out; i++)
                     h[i * (dim_out + 2)] = 1;
+                if (_mean_function.h_params_size() > 0)
+                    h.tail(_mean_function.h_params_size()) = _mean_function.h_params();
                 this->set_h_params(h);
             }
 
-            size_t h_params_size() const { return _tr.rows() * _tr.cols(); }
+            size_t h_params_size() const { return _tr.rows() * _tr.cols() + _mean_function.h_params_size(); }
 
-            Eigen::VectorXd h_params() const { return _h_params; }
+            Eigen::VectorXd h_params() const
+            {
+                Eigen::VectorXd params(h_params_size());
+                params.head(_tr.rows() * _tr.cols()) = _h_params;
+                if (_mean_function.h_params_size() > 0)
+                    params.tail(_mean_function.h_params_size()) = _mean_function.h_params();
+
+                return params;
+            }
 
             void set_h_params(const Eigen::VectorXd& p)
             {
-                _h_params = p;
+                _h_params = p.head(_tr.rows() * _tr.cols());
                 for (int c = 0; c < _tr.cols(); c++)
                     for (int r = 0; r < _tr.rows(); r++)
                         _tr(r, c) = p[r * _tr.cols() + c];
+
+                if (_mean_function.h_params_size() > 0)
+                    _mean_function.set_h_params(p.tail(_mean_function.h_params_size()));
             }
 
             template <typename GP>
             Eigen::MatrixXd grad(const Eigen::VectorXd& x, const GP& gp) const
             {
-                Eigen::MatrixXd grad = Eigen::MatrixXd::Zero(_tr.rows(), _h_params.size());
+                Eigen::MatrixXd grad = Eigen::MatrixXd::Zero(_tr.rows(), h_params_size());
                 Eigen::VectorXd m = _mean_function(x, gp);
                 for (int i = 0; i < _tr.rows(); i++) {
                     grad.block(i, i * _tr.cols(), 1, _tr.cols() - 1) = m.transpose();
                     grad(i, (i + 1) * _tr.cols() - 1) = 1;
                 }
+                if (_mean_function.h_params_size() > 0) {
+                    Eigen::MatrixXd m_grad = Eigen::MatrixXd::Zero(_tr.rows() + 1, _mean_function.h_params_size());
+                    m_grad.block(0, 0, _tr.rows(), _mean_function.h_params_size()) = _mean_function.grad(x, gp);
+                    Eigen::MatrixXd gg = _tr * m_grad;
+                    grad.block(0, h_params_size() - _mean_function.h_params_size(), _tr.rows(), _mean_function.h_params_size()) = gg;
+                }
                 return grad;
             }
 

src/tests/test_gp.cpp

@@ -50,15 +50,15 @@
 #include <boost/test/unit_test.hpp>
 
 #include <limbo/acqui/ucb.hpp>
+#include <limbo/kernel/exp.hpp>
 #include <limbo/kernel/matern_five_halves.hpp>
 #include <limbo/kernel/matern_three_halves.hpp>
-#include <limbo/kernel/exp.hpp>
 #include <limbo/kernel/squared_exp_ard.hpp>
 #include <limbo/mean/constant.hpp>
 #include <limbo/mean/function_ard.hpp>
 #include <limbo/model/gp.hpp>
-#include <limbo/model/gp/kernel_loo_opt.hpp>
 #include <limbo/model/gp/kernel_lf_opt.hpp>
+#include <limbo/model/gp/kernel_loo_opt.hpp>
 #include <limbo/model/gp/kernel_mean_lf_opt.hpp>
 #include <limbo/model/gp/mean_lf_opt.hpp>
 #include <limbo/opt/grid_search.hpp>
@@ -152,9 +152,9 @@ BOOST_AUTO_TEST_CASE(test_gp_check_lf_grad)
     }
 
     for (int i = 0; i < M; i++) {
-        test_samples.push_back(tools::random_vector(4 + 1));
-        test_samples_mean.push_back(tools::random_vector(6));
-        test_samples_kernel_mean.push_back(tools::random_vector(6 + 4 + 1));
+        test_samples.push_back(tools::random_vector(gp.kernel_function().h_params_size()));
+        test_samples_mean.push_back(tools::random_vector(gp.mean_function().h_params_size()));
+        test_samples_kernel_mean.push_back(tools::random_vector(gp.kernel_function().h_params_size() + gp.mean_function().h_params_size()));
     }
 
     gp.compute(samples, observations);
@@ -244,8 +244,8 @@ BOOST_AUTO_TEST_CASE(test_gp_check_lf_grad_noise)
     }
 
     for (int i = 0; i < M; i++) {
-        test_samples.push_back(tools::random_vector(4 + 2));
-        test_samples_kernel_mean.push_back(tools::random_vector(6 + 4 + 2));
+        test_samples.push_back(tools::random_vector(gp.kernel_function().h_params_size()));
+        test_samples_kernel_mean.push_back(tools::random_vector(gp.kernel_function().h_params_size() + gp.mean_function().h_params_size()));
     }
 
     gp.compute(samples, observations);
@@ -297,7 +297,7 @@ BOOST_AUTO_TEST_CASE(test_gp_check_loo_grad)
     }
 
     for (int i = 0; i < M; i++) {
-        test_samples.push_back(tools::random_vector(4 + 1));
+        test_samples.push_back(tools::random_vector(gp.kernel_function().h_params_size()));
     }
 
     gp.compute(samples, observations);
@@ -367,7 +367,7 @@ BOOST_AUTO_TEST_CASE(test_gp_check_loo_grad_noise)
     }
 
     for (int i = 0; i < M; i++) {
-        test_samples.push_back(tools::random_vector(4 + 2));
+        test_samples.push_back(tools::random_vector(gp.kernel_function().h_params_size()));
     }
 
     gp.compute(samples, observations);

src/tests/test_kernel.cpp

@@ -157,7 +157,7 @@ void check_kernel(size_t N, size_t K)
 
 BOOST_AUTO_TEST_CASE(test_grad_exp)
 {
-    for (int i = 1; i < 10; i++) {
+    for (int i = 1; i <= 10; i++) {
         check_kernel<kernel::Exp<Params>>(i, 100);
         check_kernel<kernel::Exp<ParamsNoise>>(i, 100);
     }
@@ -165,7 +165,7 @@ BOOST_AUTO_TEST_CASE(test_grad_exp)
 
 BOOST_AUTO_TEST_CASE(test_grad_matern_three)
 {
-    for (int i = 1; i < 10; i++) {
+    for (int i = 1; i <= 10; i++) {
         check_kernel<kernel::MaternThreeHalves<Params>>(i, 100);
         check_kernel<kernel::MaternThreeHalves<ParamsNoise>>(i, 100);
     }
@@ -173,7 +173,7 @@ BOOST_AUTO_TEST_CASE(test_grad_matern_three)
 
 BOOST_AUTO_TEST_CASE(test_grad_matern_five)
 {
-    for (int i = 1; i < 10; i++) {
+    for (int i = 1; i <= 10; i++) {
         check_kernel<kernel::MaternFiveHalves<Params>>(i, 100);
         check_kernel<kernel::MaternFiveHalves<ParamsNoise>>(i, 100);
     }
@@ -182,15 +182,15 @@ BOOST_AUTO_TEST_CASE(test_grad_matern_five)
 BOOST_AUTO_TEST_CASE(test_grad_SE_ARD)
 {
     Params::kernel_squared_exp_ard::set_k(0);
-    for (int i = 1; i < 10; i++) {
+    for (int i = 1; i <= 10; i++) {
         check_kernel<kernel::SquaredExpARD<Params>>(i, 100);
         check_kernel<kernel::SquaredExpARD<ParamsNoise>>(i, 100);
     }
 
     // THIS TEST FAILS!
     // Params::kernel_squared_exp_ard::set_k(1);
-    // for (int i = 1; i < 2; i++) {
-    //     check_kernel<kernel::SquaredExpARD<Params>>(i, 10);
+    // for (int i = 1; i <= 10; i++) {
+    //     check_kernel<kernel::SquaredExpARD<Params>>(i, 100);
     // }
 }
 

src/tests/test_mean.cpp

@@ -52,6 +52,7 @@
 
 #include <limbo/mean/constant.hpp>
 #include <limbo/mean/function_ard.hpp>
+#include <limbo/mean/null_function.hpp>
 #include <limbo/tools/macros.hpp>
 #include <limbo/tools/random_generator.hpp>
 
@@ -103,21 +104,32 @@ void check_mean(size_t N, size_t K)
         Eigen::VectorXd v = tools::random_vector(N).array() * 10. - 5.;
 
         std::tie(error, analytic, finite_diff) = check_grad(mean, hp, v);
-        // std::cout << error << ": " << analytic.transpose() << " vs " << finite_diff.transpose() << std::endl;
+        // std::cout << error << ": " << analytic << " vs " << finite_diff << std::endl;
         BOOST_CHECK(error < 1e-6);
     }
 }
 
 BOOST_AUTO_TEST_CASE(test_mean_constant)
 {
-    for (int i = 1; i < 10; i++) {
+    for (int i = 1; i <= 10; i++) {
         check_mean<mean::Constant<Params>>(i, 100);
     }
 }
 
 BOOST_AUTO_TEST_CASE(test_mean_function_ard)
 {
-    for (int i = 1; i < 10; i++) {
+    // This test checks the gradients computation of FunctionARD when the base mean function
+    // also has tunable parameters
+    for (int i = 1; i <= 10; i++) {
         check_mean<mean::FunctionARD<Params, mean::Constant<Params>>>(i, 100);
     }
 }
+
+BOOST_AUTO_TEST_CASE(test_mean_function_ard_dummy)
+{
+    // This test checks the gradients computation of FunctionARD when the base mean function
+    // has no tunable parameters
+    for (int i = 1; i <= 10; i++) {
+        check_mean<mean::FunctionARD<Params, mean::NullFunction<Params>>>(i, 100);
+    }
+}