WIP: Add a DefaultInducingPointStrategy

include/albatross/CovarianceFunctions

@@ -18,6 +18,7 @@
 #include <albatross/src/core/declarations.hpp>
 #include <albatross/src/core/traits.hpp>
 #include <albatross/src/core/priors.hpp>
+#include <albatross/src/core/concatenate.hpp>
 #include <albatross/src/core/parameter_handling_mixin.hpp>
 #include <albatross/src/core/parameter_macros.hpp>
 

include/albatross/Dataset

@@ -15,7 +15,6 @@
 
 #include "Distribution"
 
-#include <albatross/src/core/concatenate.hpp>
 #include <albatross/src/core/dataset.hpp>
 
 #endif

include/albatross/src/core/traits.hpp

@@ -15,18 +15,34 @@
 
 namespace albatross {
 
-/*
- * This determines whether or not a class, T, has a method,
- *   `std::string T.name() const`
- */
-template <typename T> class has_name {
-  template <typename C,
-            typename ReturnType = decltype(std::declval<const C>().name())>
-  static typename std::enable_if<std::is_same<std::string, ReturnType>::value,
-                                 std::true_type>::type
-  test(int);
+HAS_METHOD_WITH_RETURN_TYPE(name);
+
+template <typename T>
+class has_name : public has_name_with_return_type<T, std::string> {};
+
+HAS_METHOD_WITH_RETURN_TYPE(gridded_features);
+
+template <typename T, typename FeatureType> class ssr_feature_type {
+  template <
+      typename C,
+      typename ReturnType = decltype(std::declval<const C>()._ssr_features(
+          std::declval<const std::vector<FeatureType>>()))>
+  static typename std::enable_if<is_vector<ReturnType>::value,
+                                 typename ReturnType::value_type>::type
+  test(C *);
+  template <typename> static void test(...);
+
+public:
+  typedef decltype(test<T>(0)) type;
+};
 
-  template <typename C> static std::false_type test(...);
+template <typename T, typename FeatureType> class has_valid_ssr_features {
+  template <typename C, typename ReturnType =
+                            typename ssr_feature_type<C, FeatureType>::type>
+  static typename std::enable_if<!std::is_same<ReturnType, void>::value,
+                                 std::true_type>::type
+  test(C *);
+  template <typename> static std::false_type test(...);
 
 public:
   static constexpr bool value = decltype(test<T>(0))::value;

include/albatross/src/covariance_functions/covariance_function.hpp

@@ -195,6 +195,17 @@ class CovarianceFunction : public ParameterHandlingMixin {
     return diag;
   }
 
+  /*
+   * Cross covariance between two vectors of (possibly) different types.
+   */
+  template <
+      typename FeatureType,
+      typename std::enable_if<
+          has_valid_ssr_features<Derived, FeatureType>::value, int>::type = 0>
+  auto get_ssr_features(const std::vector<FeatureType> &features) const {
+    return derived()._ssr_features(features);
+  }
+
   /*
    * Stubs to catch the case where a covariance function was called
    * with arguments that aren't supported.
@@ -299,6 +310,35 @@ class SumOfCovarianceFunctions
     return this->rhs_(x, y);
   }
 
+  /*
+   * Manage the concatenation of ssr features.
+   */
+  template <typename X,
+            typename std::enable_if<has_valid_ssr_features<LHS, X>::value &&
+                                        !has_valid_ssr_features<RHS, X>::value,
+                                    int>::type = 0>
+  auto _ssr_features(const std::vector<X> &features) const {
+    return this->lhs_._ssr_features(features);
+  }
+
+  template <typename X,
+            typename std::enable_if<!has_valid_ssr_features<LHS, X>::value &&
+                                        has_valid_ssr_features<RHS, X>::value,
+                                    int>::type = 0>
+  auto _ssr_features(const std::vector<X> &features) const {
+    return this->rhs_._ssr_features(features);
+  }
+
+  template <typename X,
+            typename std::enable_if<has_valid_ssr_features<LHS, X>::value &&
+                                        has_valid_ssr_features<RHS, X>::value,
+                                    int>::type = 0>
+  auto _ssr_features(const std::vector<X> &features) const {
+    const auto lhs = this->lhs_._ssr_features(features);
+    const auto rhs = this->rhs_._ssr_features(features);
+    return concatenate(lhs, rhs);
+  }
+
 protected:
   LHS lhs_;
   RHS rhs_;
@@ -372,6 +412,35 @@ class ProductOfCovarianceFunctions
     return this->rhs_(x, y);
   }
 
+  /*
+   * Manage the concatenation of ssr features.
+   */
+  template <typename X,
+            typename std::enable_if<has_valid_ssr_features<LHS, X>::value &&
+                                        !has_valid_ssr_features<RHS, X>::value,
+                                    int>::type = 0>
+  auto _ssr_features(const std::vector<X> &features) const {
+    return this->lhs_._ssr_features(features);
+  }
+
+  template <typename X,
+            typename std::enable_if<!has_valid_ssr_features<LHS, X>::value &&
+                                        has_valid_ssr_features<RHS, X>::value,
+                                    int>::type = 0>
+  auto _ssr_features(const std::vector<X> &features) const {
+    return this->rhs_._ssr_features(features);
+  }
+
+  template <typename X,
+            typename std::enable_if<has_valid_ssr_features<LHS, X>::value &&
+                                        has_valid_ssr_features<RHS, X>::value,
+                                    int>::type = 0>
+  auto _ssr_features(const std::vector<X> &features) const {
+    const auto lhs = this->lhs_._ssr_features(features);
+    const auto rhs = this->rhs_._ssr_features(features);
+    return concatenate(lhs, rhs);
+  }
+
 protected:
   LHS lhs_;
   RHS rhs_;

include/albatross/src/covariance_functions/distance_metrics.hpp

@@ -17,6 +17,16 @@ namespace albatross {
 
 constexpr double EPSILON = 1e-16;
 
+std::vector<double> inline linspace(double a, double b, std::size_t n) {
+  double h = (b - a) / static_cast<double>(n - 1);
+  std::vector<double> xs(n);
+  typename std::vector<double>::iterator x;
+  double val;
+  for (x = xs.begin(), val = a; x != xs.end(); ++x, val += h)
+    *x = val;
+  return xs;
+}
+
 class DistanceMetric : public ParameterHandlingMixin {
 public:
   DistanceMetric(){};
@@ -42,6 +52,24 @@ class EuclideanDistance : public DistanceMetric {
                     const Eigen::Matrix<_Scalar, _Rows, 1> &y) const {
     return (x - y).norm();
   }
+
+  std::vector<double> gridded_features(const std::vector<double> &features,
+                                       const double &spacing) const {
+    assert(features.size() > 0);
+    double min = *std::min_element(features.begin(), features.end());
+    double max = *std::max_element(features.begin(), features.end());
+
+    double count = (max - min) / spacing;
+    count = ceil(count);
+    assert(count >= 0.);
+    assert(count < std::numeric_limits<std::size_t>::max());
+    std::size_t n = static_cast<std::size_t>(count);
+    if (n == 1) {
+      return {min, max};
+    } else {
+      return linspace(min, max, n);
+    }
+  }
 };
 
 template <typename _Scalar, int _Rows>

include/albatross/src/covariance_functions/radial.hpp

@@ -71,6 +71,22 @@ class SquaredExponential
         sigma_squared_exponential.value);
   }
 
+  template <
+      typename FeatureType,
+      typename std::enable_if<has_gridded_features_with_return_type<
+                                  DistanceMetricType, std::vector<FeatureType>,
+                                  std::vector<FeatureType>, double>::value,
+                              int>::type = 0>
+  std::vector<FeatureType>
+  _ssr_features(const std::vector<FeatureType> &features) const {
+    //  std::vector<double> _ssr_features(const std::vector<double> &features)
+    //  const {
+    // This is the spacing for which the correlation between points drops
+    // by one percent
+    double spacing = 0.1 * squared_exponential_length_scale.value;
+    return distance_metric_.gridded_features(features, spacing);
+  }
+
   DistanceMetricType distance_metric_;
 };
 

include/albatross/src/models/sparse_gp.hpp

@@ -21,33 +21,32 @@ const std::string measurement_nugget_name = "measurement_nugget";
 const std::string inducing_nugget_name = "inducing_nugget";
 } // namespace details
 
-template <typename CovFunc, typename InducingPointStrategy,
-          typename IndexingFunction>
-class SparseGaussianProcessRegression;
-
-std::vector<double> inline linspace(double a, double b, std::size_t n) {
-  double h = (b - a) / static_cast<double>(n - 1);
-  std::vector<double> xs(n);
-  typename std::vector<double>::iterator x;
-  double val;
-  for (x = xs.begin(), val = a; x != xs.end(); ++x, val += h)
-    *x = val;
-  return xs;
-}
+template <typename CovFunc> struct DefaultInducingPointStrategy;
 
-struct UniformlySpacedInducingPoints {
+template <typename CovFunc, typename IndexingFunction,
+          typename InducingPointStrategy =
+              DefaultInducingPointStrategy<CovFunc>>
+class SparseGaussianProcessRegression;
 
-  UniformlySpacedInducingPoints(std::size_t num_points_ = 10)
-      : num_points(num_points_) {}
+/*
+ * This uses the state space representation provided by the
+ * covariance function as the inducing points.  The nice part
+ * about this is that you can have the inducing points change
+ * as a function of the parameters of the covariance function,
+ * the downside is that you have to globally modify the behavior
+ * of a covariance function to get it to behave for a sparse GP.
+ */
+template <typename CovFunc> struct DefaultInducingPointStrategy {
 
-  std::vector<double> operator()(const std::vector<double> &features) const {
-    double min = *std::min_element(features.begin(), features.end());
-    double max = *std::max_element(features.begin(), features.end());
+  DefaultInducingPointStrategy(const CovFunc &cov_func_)
+      : cov_func(cov_func_){};
 
-    return linspace(min, max, num_points);
+  template <typename FeatureType>
+  auto operator()(const std::vector<FeatureType> &features) const {
+    return cov_func.get_ssr_features(features);
   }
 
-  std::size_t num_points;
+  CovFunc cov_func;
 };
 
 /*
@@ -101,35 +100,35 @@ struct UniformlySpacedInducingPoints {
  *     - https://bwengals.github.io/pymc3-fitcvfe-implementation-notes.html
  *     - https://github.com/SheffieldML/GPy see fitc.py
  */
-template <typename CovFunc, typename InducingPointStrategy,
-          typename IndexingFunction>
+template <typename CovFunc, typename IndexingFunction,
+          typename InducingPointStrategy>
 class SparseGaussianProcessRegression
     : public GaussianProcessBase<
-          CovFunc, SparseGaussianProcessRegression<
-                       CovFunc, InducingPointStrategy, IndexingFunction>> {
+          CovFunc, SparseGaussianProcessRegression<CovFunc, IndexingFunction,
+                                                   InducingPointStrategy>> {
 
 public:
   using Base = GaussianProcessBase<
-      CovFunc, SparseGaussianProcessRegression<CovFunc, InducingPointStrategy,
-                                               IndexingFunction>>;
+      CovFunc, SparseGaussianProcessRegression<CovFunc, IndexingFunction,
+                                               InducingPointStrategy>>;
 
   SparseGaussianProcessRegression() : Base() { initialize_params(); };
-  SparseGaussianProcessRegression(CovFunc &covariance_function)
+  SparseGaussianProcessRegression(const CovFunc &covariance_function)
       : Base(covariance_function) {
     initialize_params();
   };
   SparseGaussianProcessRegression(
-      CovFunc &covariance_function,
-      InducingPointStrategy &inducing_point_strategy_,
-      IndexingFunction &independent_group_indexing_function_,
+      const CovFunc &covariance_function,
+      const IndexingFunction &independent_group_indexing_function_,
+      const InducingPointStrategy &inducing_point_strategy_,
       const std::string &model_name)
       : Base(covariance_function, model_name),
-        inducing_point_strategy(inducing_point_strategy_),
         independent_group_indexing_function(
-            independent_group_indexing_function_) {
+            independent_group_indexing_function_),
+        inducing_point_strategy(inducing_point_strategy_) {
     initialize_params();
   };
-  SparseGaussianProcessRegression(CovFunc &covariance_function,
+  SparseGaussianProcessRegression(const CovFunc &covariance_function,
                                   const std::string &model_name)
       : Base(covariance_function, model_name) {
     initialize_params();
@@ -294,20 +293,31 @@ class SparseGaussianProcessRegression
 
   Parameter measurement_nugget_;
   Parameter inducing_nugget_;
-  InducingPointStrategy inducing_point_strategy;
   IndexingFunction independent_group_indexing_function;
+  InducingPointStrategy inducing_point_strategy;
 };
 
-template <typename CovFunc, typename InducingPointStrategy,
-          typename IndexingFunction>
+template <typename CovFunc, typename IndexingFunction,
+          typename InducingPointStrategy>
 auto sparse_gp_from_covariance(CovFunc covariance_function,
+                               IndexingFunction &index_function,
                                InducingPointStrategy &strategy,
+                               const std::string &model_name) {
+  return SparseGaussianProcessRegression<CovFunc, IndexingFunction,
+                                         InducingPointStrategy>(
+      covariance_function, index_function, strategy, model_name);
+};
+
+template <typename CovFunc, typename IndexingFunction>
+auto sparse_gp_from_covariance(CovFunc covariance_function,
                                IndexingFunction &index_function,
                                const std::string &model_name) {
-  return SparseGaussianProcessRegression<CovFunc, InducingPointStrategy,
-                                         IndexingFunction>(
-      covariance_function, strategy, index_function, model_name);
+  const DefaultInducingPointStrategy<CovFunc> strategy(covariance_function);
+  return SparseGaussianProcessRegression<CovFunc, IndexingFunction,
+                                         DefaultInducingPointStrategy<CovFunc>>(
+      covariance_function, index_function, strategy, model_name);
 };
+
 } // namespace albatross
 
 #endif /* INCLUDE_ALBATROSS_MODELS_SPARSE_GP_H_ */