Add indices to BaseKernel and fix minor issues with restructure of kernels

src/limbo/kernel/kernel.hpp

@@ -78,18 +78,18 @@ namespace limbo {
                 _noise_p = std::log(std::sqrt(_noise));
             }
 
-            double operator()(const Eigen::VectorXd& v1, const Eigen::VectorXd& v2) const
+            double operator()(const Eigen::VectorXd& v1, const Eigen::VectorXd& v2, int i = -1, int j = -2) const
             {
-                return static_cast<const Kernel*>(this)->kernel(v1, v2) + (((v1 - v2).norm() < 1e-8) ? _noise + 1e-8 : 0.0);
+                return static_cast<const Kernel*>(this)->kernel(v1, v2) + ((i == j) ? _noise + 1e-8 : 0.0);
             }
 
-            Eigen::VectorXd grad(const Eigen::VectorXd& x1, const Eigen::VectorXd& x2) const
+            Eigen::VectorXd grad(const Eigen::VectorXd& x1, const Eigen::VectorXd& x2, int i = -1, int j = -2) const
             {
                 Eigen::VectorXd g = static_cast<const Kernel*>(this)->gradient(x1, x2);
 
                 if (Params::kernel::optimize_noise()) {
                     g.conservativeResize(g.size() + 1);
-                    g(g.size() - 1) = (((x1 - x2).norm() < 1e-8) ? 2.0 * _noise : 0.0);
+                    g(g.size() - 1) = ((i == j) ? 2.0 * _noise : 0.0);
                 }
 
                 return g;

src/limbo/model/gp.hpp

@@ -296,7 +296,7 @@ namespace limbo {
                 // O(n^2) [should be negligible]
                 for (size_t i = 0; i < n; i++)
                     for (size_t j = 0; j <= i; ++j)
-                        _kernel(i, j) = _kernel_function(_samples[i], _samples[j]);
+                        _kernel(i, j) = _kernel_function(_samples[i], _samples[j], i, j);
 
                 for (size_t i = 0; i < n; i++)
                     for (size_t j = 0; j < i; ++j)
@@ -319,7 +319,7 @@ namespace limbo {
                 _kernel.conservativeResize(n, n);
 
                 for (size_t i = 0; i < n; ++i) {
-                    _kernel(i, n - 1) = _kernel_function(_samples[i], _samples[n - 1]);
+                    _kernel(i, n - 1) = _kernel_function(_samples[i], _samples[n - 1], i, n - 1);
                     _kernel(n - 1, i) = _kernel(i, n - 1);
                 }
 

src/limbo/model/gp/kernel_lf_opt.hpp

@@ -111,7 +111,7 @@ namespace limbo {
                         Eigen::VectorXd grad = Eigen::VectorXd::Zero(params.size());
                         for (size_t i = 0; i < n; ++i) {
                             for (size_t j = 0; j <= i; ++j) {
-                                Eigen::VectorXd g = gp.kernel_function().grad(gp.samples()[i], gp.samples()[j]);
+                                Eigen::VectorXd g = gp.kernel_function().grad(gp.samples()[i], gp.samples()[j], i, j);
                                 if (i == j)
                                     grad += w(i, j) * g * 0.5;
                                 else

src/limbo/model/gp/kernel_mean_lf_opt.hpp

@@ -117,7 +117,7 @@ namespace limbo {
                         Eigen::VectorXd grad = Eigen::VectorXd::Zero(params.size());
                         for (size_t i = 0; i < n; ++i) {
                             for (size_t j = 0; j <= i; ++j) {
-                                Eigen::VectorXd g = gp.kernel_function().grad(gp.samples()[i], gp.samples()[j]);
+                                Eigen::VectorXd g = gp.kernel_function().grad(gp.samples()[i], gp.samples()[j], i, j);
                                 if (i == j)
                                     grad.head(gp.kernel_function().h_params_size()) += w(i, j) * g * 0.5;
                                 else

src/tests/test_gp.cpp

@@ -145,7 +145,9 @@ BOOST_AUTO_TEST_CASE(test_gp_check_lf_grad)
 
     for (int i = 0; i < N; i++) {
         samples.push_back(tools::random_vector(4));
-        observations.push_back(tools::random_vector(2));
+        Eigen::VectorXd ob(2);
+        ob << std::cos(samples[i](0)), std::sin(samples[i](1));
+        observations.push_back(ob);
     }
 
     for (int i = 0; i < M; i++) {
@@ -235,7 +237,9 @@ BOOST_AUTO_TEST_CASE(test_gp_check_lf_grad_noise)
 
     for (int i = 0; i < N; i++) {
         samples.push_back(tools::random_vector(4));
-        observations.push_back(tools::random_vector(2));
+        Eigen::VectorXd ob(2);
+        ob << std::cos(samples[i](0)), std::sin(samples[i](1));
+        observations.push_back(ob);
     }
 
     for (int i = 0; i < M; i++) {
@@ -290,7 +294,7 @@ BOOST_AUTO_TEST_CASE(test_gp_dim)
     BOOST_CHECK(std::abs((mu(0) - 5)) < 1);
     BOOST_CHECK(std::abs((mu(1) - 5)) < 1);
 
-    BOOST_CHECK(sigma < 1e-5);
+    BOOST_CHECK(sigma <= Params::kernel::noise() + 1e-8);
 }
 
 BOOST_AUTO_TEST_CASE(test_gp)
@@ -312,15 +316,15 @@ BOOST_AUTO_TEST_CASE(test_gp)
     double sigma;
     std::tie(mu, sigma) = gp.query(make_v1(1));
     BOOST_CHECK(std::abs((mu(0) - 5)) < 1);
-    BOOST_CHECK(sigma < 1e-5);
+    BOOST_CHECK(sigma <= 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 < 1e-5);
+    BOOST_CHECK(sigma <= 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 < 1e-5);
+    BOOST_CHECK(sigma <= Params::kernel::noise() + 1e-8);
 
     for (double x = 0; x < 4; x += 0.05) {
         Eigen::VectorXd mu;
@@ -328,8 +332,8 @@ BOOST_AUTO_TEST_CASE(test_gp)
         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);
-        std::cout << x << " " << mu << " " << mu.array() - sigma << " "
-                  << mu.array() + sigma << std::endl;
+        // std::cout << x << " " << mu << " " << mu.array() - sigma << " "
+        //           << mu.array() + sigma << std::endl;
     }
 }
 
@@ -357,9 +361,9 @@ BOOST_AUTO_TEST_CASE(test_gp_bw_inversion)
         auto t1 = std::chrono::steady_clock::now();
         gp.compute(samples, observations);
         auto time_init = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
-        std::cout.precision(17);
-        std::cout << "Time running first batch: " << time_init << "us" << std::endl
-                  << std::endl;
+        // std::cout.precision(17);
+        // std::cout << "Time running first batch: " << time_init << "us" << std::endl
+        //           << std::endl;
 
         observations.push_back(make_v1(rgen.rand()));
         samples.push_back(make_v1(rgen.rand()));
@@ -367,21 +371,21 @@ BOOST_AUTO_TEST_CASE(test_gp_bw_inversion)
         t1 = std::chrono::steady_clock::now();
         gp.add_sample(samples.back(), observations.back());
         auto time_increment = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
-        std::cout << "Time running increment: " << time_increment << "us" << std::endl
-                  << std::endl;
+        // std::cout << "Time running increment: " << time_increment << "us" << std::endl
+        //           << std::endl;
 
         t1 = std::chrono::steady_clock::now();
         gp.recompute(true);
         auto time_recompute = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
-        std::cout << "Time recomputing: " << time_recompute << "us" << std::endl
-                  << std::endl;
+        // std::cout << "Time recomputing: " << time_recompute << "us" << std::endl
+        //           << std::endl;
 
         GP_t gp2;
         t1 = std::chrono::steady_clock::now();
         gp2.compute(samples, observations);
         auto time_full = std::chrono::duration_cast<std::chrono::microseconds>(std::chrono::steady_clock::now() - t1).count();
-        std::cout << "Time running whole batch: " << time_full << "us" << std::endl
-                  << std::endl;
+        // std::cout << "Time running whole batch: " << time_full << "us" << std::endl
+        //           << std::endl;
 
         Eigen::VectorXd s = make_v1(rgen.rand());
         if ((gp.mu(s) - gp2.mu(s)).norm() >= 1e-5)
@@ -440,23 +444,22 @@ BOOST_AUTO_TEST_CASE(test_gp_auto)
     std::vector<Eigen::VectorXd> observations = {make_v1(5), make_v1(10), make_v1(5)};
     std::vector<Eigen::VectorXd> samples = {make_v1(1), make_v1(2), make_v1(3)};
 
-    gp.compute(samples, observations);
+    gp.compute(samples, observations, false);
     gp.optimize_hyperparams();
-    gp.recompute(false);
 
     Eigen::VectorXd mu;
     double sigma;
     std::tie(mu, sigma) = gp.query(make_v1(1));
     BOOST_CHECK(std::abs((mu(0) - 5)) < 1);
-    BOOST_CHECK(sigma < 1e-5);
+    BOOST_CHECK(sigma <= gp.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 < 1e-5);
+    BOOST_CHECK(sigma <= gp.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 < 1e-5);
+    BOOST_CHECK(sigma <= gp.kernel_function().noise() + 1e-8);
 }
 
 BOOST_AUTO_TEST_CASE(test_gp_init_variance)

src/tutorials/plot_gp.py

@@ -37,7 +37,7 @@ ax.plot(gp_ard[:,0], gp_ard[:,1], linewidth=2, color=colors[1])
 ax.plot(gp[:,0], actual, linewidth=2, linestyle='--', color=colors[3])
 ax.plot(data[:,0], data[:, 1], 'o', color=colors[2])
 
-legend = ax.legend(["GP/exp", "GP/expARD", 'cos(x)', 'Data'], loc=4);
+legend = ax.legend(["GP/exp", "GP/expARD", 'cos(x)', 'Data'], loc=8);
 frame = legend.get_frame()
 frame.set_facecolor('1.0')
 frame.set_edgecolor('1.0')