expose distance derivatives to enable GPy MaternLLD kernel

davmre · Aug 12, 2015 · 873d607 · 873d607
1 parent e3f824e
commit 873d607
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 1 deletion.
diff --git a/src_c/distances.cc b/src_c/distances.cc
@@ -134,7 +134,7 @@ where
 
   if (deriv_denom == 0) {
     double t2 = pow(sin(d/(2*AVG_EARTH_RADIUS_KM)), 2);
-    printf("WARNING: zero denom deriv in t1 %f t2 %f d %f d/stuff %f dlon %f dlat %f num %f\n", t*1e8, t2*1e8, d*1e8, d/(2*AVG_EARTH_RADIUS_KM) * 1e8, rlon1-rlon2, rlat1-rlat2, deriv_num);
+    //printf("WARNING: zero denom deriv in t1 %f t2 %f d %f d/stuff %f dlon %f dlat %f num %f\n", t*1e8, t2*1e8, d*1e8, d/(2*AVG_EARTH_RADIUS_KM) * 1e8, rlon1-rlon2, rlat1-rlat2, deriv_num);
     return dist_km_deriv_wrt_xi_empirical(p1, p2, i, d);
   }
 

diff --git a/src_c/vector_mult.hpp b/src_c/vector_mult.hpp
@@ -146,6 +146,8 @@ class VectorTree {
   pyublas::numpy_matrix<double> sparse_training_kernel_matrix(const pyublas::numpy_matrix<double> &pts, double max_distance, bool distance_only);
   pyublas::numpy_matrix<double> kernel_deriv_wrt_xi(const pyublas::numpy_matrix<double> &pts1, const pyublas::numpy_matrix<double> &pts2, int i, int k);
   void kernel_deriv_wrt_xi_row(const pyublas::numpy_matrix<double> &pts1, int i, int k, pyublas::numpy_vector<double> K);
+  void dist_deriv_wrt_xi_row(const pyublas::numpy_matrix<double> &pts1, const pyublas::numpy_matrix<double> &pts2, int i, int k, pyublas::numpy_vector<double> D);
+
   pyublas::numpy_matrix<double> kernel_deriv_wrt_i(const pyublas::numpy_matrix<double> &pts1, const pyublas::numpy_matrix<double> &pts2, int param_i,  bool symmetric, const pyublas::numpy_matrix<double> distances);
   pyublas::numpy_vector<double> sparse_kernel_deriv_wrt_i(const pyublas::numpy_matrix<double> &pts1, const pyublas::numpy_matrix<double> &pts2, const pyublas::numpy_vector<int> &nzr, const pyublas::numpy_vector<int> &nzc, int param_i, const pyublas::numpy_vector<double> distance_entries);
 

diff --git a/src_c/vector_mult_py.cc b/src_c/vector_mult_py.cc
@@ -465,6 +465,23 @@ void VectorTree::kernel_deriv_wrt_xi_row(const pyublas::numpy_matrix<double> &pt
 
 }
 
+void VectorTree::dist_deriv_wrt_xi_row(const pyublas::numpy_matrix<double> &pts1, const pyublas::numpy_matrix<double> &pts2, int i, int k, pyublas::numpy_vector<double> D) {
+  // return just the row/col corresponding to the i'th input point (everything else should be zero)
+
+  if (this->ddfn_dx == NULL) {
+    printf("ERROR: gradient not implemented for this distance function.\n");
+    exit(1);
+  }
+  point p1 = {&pts1(i, 0), 0};
+  for (unsigned j = 0; j < pts2.size1 (); ++ j) {
+    point p2 = {&pts2(j, 0), 0};
+    double r = this->dfn(p1, p2, std::numeric_limits< double >::max(), this->dist_params, this->dfn_extra);
+    double dr_dp1 = this->ddfn_dx(p1.p, p2.p, k, r, std::numeric_limits< double >::max(), this->dist_params, this->dfn_extra);
+    D(j) = dr_dp1;
+  }
+
+}
+
 
 pyublas::numpy_matrix<double> VectorTree::kernel_deriv_wrt_i(const pyublas::numpy_matrix<double> &pts1, const pyublas::numpy_matrix<double> &pts2, int param_i, bool symmetric, const pyublas::numpy_matrix<double> distances) {
 
@@ -678,6 +695,7 @@ BOOST_PYTHON_MODULE(cover_tree) {
     .def("sparse_training_kernel_matrix", &VectorTree::sparse_training_kernel_matrix)
     .def("kernel_deriv_wrt_xi", &VectorTree::kernel_deriv_wrt_xi)
     .def("kernel_deriv_wrt_xi_row", &VectorTree::kernel_deriv_wrt_xi_row)
+    .def("dist_deriv_wrt_xi_row", &VectorTree::dist_deriv_wrt_xi_row)
     .def("kernel_deriv_wrt_i", &VectorTree::kernel_deriv_wrt_i)
     .def("sparse_kernel_deriv_wrt_i", &VectorTree::sparse_kernel_deriv_wrt_i)
     .def("quadratic_form_from_dense_hack", &VectorTree::quadratic_form_from_dense_hack)