Fix for gp var + kernel parameters

src/examples/experimental/imgpo.cpp

@@ -193,6 +193,7 @@ struct Params {
 
     struct kernel_squared_exp_ard {
         BO_PARAM(int, k, 0);
+        BO_PARAM(double, sigma_sq, 1);
     };
 
     struct stop_maxiterations {

src/examples/experimental/parego.cpp

@@ -20,7 +20,7 @@ struct Params {
     };
 
     struct kernel_maternfivehalfs {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
         BO_PARAM(double, l, 0.2);
     };
 
@@ -35,6 +35,7 @@ struct Params {
 
     struct kernel_squared_exp_ard {
         BO_PARAM(int, k, 0);
+        BO_PARAM(double, sigma_sq, 1);
     };
 
     struct init_randomsampling {

src/examples/mono_dim.cpp

@@ -28,7 +28,7 @@ BO_PARAMS(std::cout,
               };
 
               struct kernel_maternfivehalfs {
-                  BO_PARAM(double, sigma, 1);
+                  BO_PARAM(double, sigma_sq, 1);
                   BO_PARAM(double, l, 0.2);
               };
 

src/examples/obs_multi.cpp

@@ -22,7 +22,7 @@ struct Params {
     };
 #endif
     struct kernel_maternfivehalfs {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
         BO_PARAM(double, l, 0.2);
     };
 
@@ -30,7 +30,7 @@ struct Params {
         BO_PARAM(bool, stats_enabled, true);
     };
 
-    struct bayes_opt_boptimizer : public defaults::bayes_opt_boptimizer{
+    struct bayes_opt_boptimizer : public defaults::bayes_opt_boptimizer {
         BO_PARAM(double, noise, 0.001);
     };
 

src/examples/obs_multi_auto_mean.cpp

@@ -25,7 +25,7 @@ struct Params {
     };
 
     struct kernel_maternfivehalfs {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
         BO_PARAM(double, l, 0.2);
     };
 

src/limbo/kernel/exp.hpp

@@ -9,20 +9,22 @@ namespace limbo {
     namespace defaults {
         struct kernel_exp {
             /// @ingroup kernel_defaults
-            BO_PARAM(double, sigma, 1);
+            BO_PARAM(double, sigma_sq, 1);
+            BO_PARAM(double, l, 1);
         };
     }
     namespace kernel {
         /**
           @ingroup kernel
           \rst
-          Exponential kernel with a :math:`\sigma` parameter that controls the width (see :cite:`brochu2010tutorial` p. 9).
+          Exponential kernel (see :cite:`brochu2010tutorial` p. 9).
 
           .. math::
-              k(v_1, v_2)  = \exp \Big(-\frac{1}{\sigma^2} ||v_1 - v_2||^2\Big)
+              k(v_1, v_2)  = \sigma^2\exp \Big(-\frac{1}{\l^2} ||v_1 - v_2||^2\Big)
 
           Parameters:
-            - ``double sigma``
+            - ``double sigma_sq`` (signal variance)
+            - ``double l`` (characteristic length scale)
           \endrst
         */
         template <typename Params>
@@ -30,8 +32,8 @@ namespace limbo {
             Exp(size_t dim = 1) {}
             double operator()(const Eigen::VectorXd& v1, const Eigen::VectorXd& v2) const
             {
-                double _sigma = Params::kernel_exp::sigma();
-                return (std::exp(-(1 / (2 * std::pow(_sigma, 2))) * std::pow((v1 - v2).norm(), 2)));
+                double _l = Params::kernel_exp::l();
+                return Params::kernel_exp::sigma_sq() * (std::exp(-(1 / (2 * std::pow(_l, 2))) * std::pow((v1 - v2).norm(), 2)));
             }
         };
     }

src/limbo/kernel/matern_five_halfs.hpp

@@ -9,7 +9,7 @@ namespace limbo {
     namespace defaults {
         struct kernel_maternfivehalfs {
             /// @ingroup kernel_defaults
-            BO_PARAM(double, sigma, 1);
+            BO_PARAM(double, sigma_sq, 1);
             /// @ingroup kernel_defaults
             BO_PARAM(double, l, 1);
         };
@@ -20,7 +20,7 @@ namespace limbo {
 
           \rst
 
-          Matern kernel (TODO: formula)
+          Matern kernel
 
           .. math::
             d = ||v1 - v2||
@@ -31,8 +31,8 @@ namespace limbo {
 
 
           Parameters:
-            - ``double sigma``
-            - ``double l``
+            - ``double sigma_sq`` (signal variance)
+            - ``double l`` (characteristic length scale)
 
           Reference: :cite:`matern1960spatial` & :cite:`brochu2010tutorial` p.10 & https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function
           \endrst
@@ -44,7 +44,7 @@ namespace limbo {
             double operator()(const Eigen::VectorXd& v1, const Eigen::VectorXd& v2) const
             {
                 double d = (v1 - v2).norm();
-                return Params::kernel_maternfivehalfs::sigma() * (1 + std::sqrt(5) * d / Params::kernel_maternfivehalfs::l() + 5 * d * d / (3 * Params::kernel_maternfivehalfs::l() * Params::kernel_maternfivehalfs::l())) * std::exp(-std::sqrt(5) * d / Params::kernel_maternfivehalfs::l());
+                return Params::kernel_maternfivehalfs::sigma_sq() * (1 + std::sqrt(5) * d / Params::kernel_maternfivehalfs::l() + 5 * d * d / (3 * Params::kernel_maternfivehalfs::l() * Params::kernel_maternfivehalfs::l())) * std::exp(-std::sqrt(5) * d / Params::kernel_maternfivehalfs::l());
             }
         };
     }

src/limbo/kernel/matern_three_halfs.hpp

@@ -10,7 +10,7 @@ namespace limbo {
         struct kernel_maternthreehalfs {
             /// @ingroup kernel_defaults
             /// This is is sigma squared!
-            BO_PARAM(double, sigma, 1);
+            BO_PARAM(double, sigma_sq, 1);
             /// @ingroup kernel_defaults
             BO_PARAM(double, l, 1);
         };
@@ -30,8 +30,8 @@ namespace limbo {
 
 
          Parameters:
-          - ``double sigma`` **squared**
-          - ``double l``
+          - ``double sigma_sq`` (signal variance)
+          - ``double l`` (characteristic length scale)
 
         Reference: :cite:`matern1960spatial` & :cite:`brochu2010tutorial` p.10 & https://en.wikipedia.org/wiki/Mat%C3%A9rn_covariance_function
         \endrst
@@ -43,7 +43,7 @@ namespace limbo {
             double operator()(const Eigen::VectorXd& v1, const Eigen::VectorXd& v2) const
             {
                 double d = (v1 - v2).norm();
-                return Params::kernel_maternthreehalfs::sigma() * (1 + std::sqrt(3) * d / Params::kernel_maternthreehalfs::l()) * std::exp(-std::sqrt(3) * d / Params::kernel_maternthreehalfs::l());
+                return Params::kernel_maternthreehalfs::sigma_sq() * (1 + std::sqrt(3) * d / Params::kernel_maternthreehalfs::l()) * std::exp(-std::sqrt(3) * d / Params::kernel_maternthreehalfs::l());
             }
         };
     }

src/limbo/kernel/squared_exp_ard.hpp

@@ -7,6 +7,7 @@ namespace limbo {
     namespace defaults {
         struct kernel_squared_exp_ard {
             BO_PARAM(int, k, 0); //equivalent to the standard exp ARD
+            BO_PARAM(double, sigma_sq, 1);
         };
     }
 
@@ -34,6 +35,7 @@ namespace limbo {
                 Eigen::VectorXd p = Eigen::VectorXd::Zero(_ell.size() + _ell.size() * Params::kernel_squared_exp_ard::k() + 1);
                 p.head(_ell.size()) = Eigen::VectorXd::Ones(_ell.size()) * -1;
                 this->set_h_params(p);
+                _sf2 = Params::kernel_squared_exp_ard::sigma_sq();
             }
 
             size_t h_params_size() const { return _ell.size() + _ell.size() * Params::kernel_squared_exp_ard::k() + 1; }
@@ -48,7 +50,7 @@ namespace limbo {
                 for (size_t j = 0; j < (unsigned int)Params::kernel_squared_exp_ard::k(); ++j)
                     for (size_t i = 0; i < _input_dim; ++i)
                         _A(i, j) = p((j + 1) * _input_dim + i); //can be negative
-                _sf2 = 1; // exp(2 * p(p.size()-1));
+                // _sf2 = 1; // exp(2 * p(p.size()-1));
             }
 
             Eigen::VectorXd grad(const Eigen::VectorXd& x1, const Eigen::VectorXd& x2) const

src/limbo/model/gp.hpp

@@ -136,14 +136,14 @@ 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` (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().
              \\endrst
 	  		*/
             std::tuple<Eigen::VectorXd, double> query(const Eigen::VectorXd& v) const
             {
                 if (_samples.size() == 0 && _bl_samples.size() == 0)
                     return std::make_tuple(_mean_function(v, *this),
-                        sqrt(_kernel_function(v, v)));
+                        _kernel_function(v, v));
 
                 if (_samples.size() == 0)
                     return std::make_tuple(_mean_function(v, *this),
@@ -155,7 +155,7 @@ namespace limbo {
 
             /**
              \\rst
-             return :math:`\\mu` (unormalized). If there is no sample, return the value according to the mean function.
+             return :math:`\mu` (unormalized). If there is no sample, return the value according to the mean function.
              \\endrst
 	  		*/
             Eigen::VectorXd mu(const Eigen::VectorXd& v) const
@@ -167,13 +167,13 @@ namespace limbo {
 
             /**
              \\rst
-             return :math:`\\sigma^2` (unormalized). If there is no sample, return the value according to the mean function.
+             return :math:`\sigma^2` (unormalized). If there is no sample, return the max \sigma^2.
              \\endrst
 	  		*/
             double sigma(const Eigen::VectorXd& v) const
             {
                 if (_samples.size() == 0 && _bl_samples.size() == 0)
-                    return sqrt(_kernel_function(v, v));
+                    return _kernel_function(v, v);
                 return _sigma(v, _compute_k_bl(v, _compute_k(v)));
             }
 

src/tests/all_combinations_template.cpp

@@ -9,16 +9,17 @@ using namespace limbo;
 struct Params {
 
     struct kernel_exp {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
+        BO_PARAM(double, l, 1);
     };
 
     struct kernel_maternthreehalfs {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
         BO_PARAM(double, l, 0.2);
     };
 
     struct kernel_maternfivehalfs {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
         BO_PARAM(double, l, 0.2);
     };
 
@@ -74,7 +75,7 @@ struct Params {
         BO_PARAM(bool, stats_enabled, true);
     };
 
-    struct bayes_opt_boptimizer : public defaults::bayes_opt_boptimizer{
+    struct bayes_opt_boptimizer : public defaults::bayes_opt_boptimizer {
         BO_PARAM(double, noise, 0.001);
     };
 };

src/tests/test_gp.cpp

@@ -5,6 +5,8 @@
 
 #include <limbo/acqui/ucb.hpp>
 #include <limbo/kernel/matern_five_halfs.hpp>
+#include <limbo/kernel/matern_three_halfs.hpp>
+#include <limbo/kernel/exp.hpp>
 #include <limbo/kernel/squared_exp_ard.hpp>
 #include <limbo/mean/constant.hpp>
 #include <limbo/model/gp.hpp>
@@ -31,7 +33,7 @@ struct Params {
     };
 
     struct kernel_maternfivehalfs {
-        BO_PARAM(double, sigma, 1);
+        BO_PARAM(double, sigma_sq, 1);
         BO_PARAM(double, l, 0.25);
     };
 
@@ -272,3 +274,63 @@ BOOST_AUTO_TEST_CASE(test_gp_auto)
     BOOST_CHECK(std::abs((mu(0) - 5)) < 1);
     BOOST_CHECK(sigma < 1e-5);
 }
+
+BOOST_AUTO_TEST_CASE(test_gp_init_variance)
+{
+    using namespace limbo;
+
+    struct Parameters {
+        struct kernel_squared_exp_ard {
+            BO_PARAM(int, k, 0);
+            BO_PARAM(double, sigma_sq, 10);
+        };
+        struct kernel_exp {
+            BO_PARAM(double, sigma_sq, 10);
+            BO_PARAM(double, l, 1);
+        };
+        struct kernel_maternthreehalfs {
+            BO_PARAM(double, sigma_sq, 10);
+            BO_PARAM(double, l, 0.25);
+        };
+        struct kernel_maternfivehalfs {
+            BO_PARAM(double, sigma_sq, 10);
+            BO_PARAM(double, l, 0.25);
+        };
+    };
+
+    // MaternThreeHalfs
+    typedef model::GP<Params, kernel::MaternThreeHalfs<Parameters>, mean::Constant<Params>> GP1_t;
+
+    GP1_t gp1(1, 1);
+
+    double sigma = gp1.sigma(tools::random_vector(1));
+
+    BOOST_CHECK_CLOSE(sigma, 10.0, 1e-5);
+
+    // MaternFiveHalfs
+    typedef model::GP<Params, kernel::MaternFiveHalfs<Parameters>, mean::Constant<Params>> GP2_t;
+
+    GP2_t gp2(1, 1);
+
+    sigma = gp2.sigma(tools::random_vector(1));
+
+    BOOST_CHECK_CLOSE(sigma, 10.0, 1e-5);
+
+    // Exponential
+    typedef model::GP<Params, kernel::Exp<Parameters>, mean::Constant<Params>> GP3_t;
+
+    GP3_t gp3(1, 1);
+
+    sigma = gp3.sigma(tools::random_vector(1));
+
+    BOOST_CHECK_CLOSE(sigma, 10.0, 1e-5);
+
+    // ARD Squared Exponential
+    typedef model::GP<Params, kernel::SquaredExpARD<Parameters>, mean::Constant<Params>> GP4_t;
+
+    GP4_t gp4(1, 1);
+
+    sigma = gp4.sigma(tools::random_vector(1));
+
+    BOOST_CHECK_CLOSE(sigma, 10.0, 1e-5);
+}

src/tests/test_kernel.cpp

@@ -10,6 +10,7 @@ using namespace limbo;
 struct Params {
     struct kernel_squared_exp_ard {
         BO_DYN_PARAM(int, k); //equivalent to the standard exp ARD
+        BO_PARAM(double, sigma_sq, 1);
     };
 };