Added the product of kernels defined on the same space + a few bugs fixed in the prod_orthogonal

GPy/core/parameterised.py

@@ -333,8 +333,7 @@ class parameterised(object):
         header_string = ["{h:^{col}}".format(h = header[i], col = cols[i]) for i in range(len(cols))]
         header_string = map(lambda x: '|'.join(x), [header_string])
         separator = '-'*len(header_string[0])
-        param_string = ["{n:^{c0}}|{v:^{c1}}|{c:^{c2}}|{t:^{c3}}".format(n = names[i], v = values[i], c = constraints[i], t = ties[i],
-                                                                        c0 = cols[0], c1 = cols[1], c2 = cols[2], c3 = cols[3]) for i in range(len(values))]
+        param_string = ["{n:^{c0}}|{v:^{c1}}|{c:^{c2}}|{t:^{c3}}".format(n = names[i], v = values[i], c = constraints[i], t = ties[i], c0 = cols[0], c1 = cols[1], c2 = cols[2], c3 = cols[3]) for i in range(len(values))]
 
 
         return ('\n'.join([header_string[0], separator]+param_string)) + '\n'

GPy/kern/__init__.py

@@ -2,5 +2,5 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product_orthogonal
+from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52, product, product_orthogonal
 from kern import kern

GPy/kern/constructors.py

@@ -18,6 +18,7 @@ from Brownian import Brownian as Brownianpart
 from periodic_exponential import periodic_exponential as periodic_exponentialpart
 from periodic_Matern32 import periodic_Matern32 as periodic_Matern32part
 from periodic_Matern52 import periodic_Matern52 as periodic_Matern52part
+from product import product as productpart
 from product_orthogonal import product_orthogonal as product_orthogonalpart
 #TODO these s=constructors are not as clean as we'd like. Tidy the code up
 #using meta-classes to make the objects construct properly wthout them.
@@ -242,10 +243,21 @@ def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,
     part = periodic_Matern52part(D,variance, lengthscale, period, n_freq, lower, upper)
     return kern(D, [part])
 
+def product(k1,k2):
+    """
+     Construct a product kernel over D from two kernels over D
+
+    :param k1, k2: the kernels to multiply
+    :type k1, k2: kernpart
+    :rtype: kernel object
+    """
+    part = productpart(k1,k2)
+    return kern(k1.D, [part])
+
 def product_orthogonal(k1,k2):
     """
-     Construct a product kernel
-     
+     Construct a product kernel over D1 x D2 from a kernel over D1 and another over D2.
+
     :param k1, k2: the kernels to multiply
     :type k1, k2: kernpart
     :rtype: kernel object

GPy/kern/product.py

@@ -0,0 +1,86 @@
+# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+from kernpart import kernpart
+import numpy as np
+import hashlib
+#from scipy import integrate # This may not be necessary (Nicolas, 20th Feb)
+
+class product(kernpart):
+    """
+    Computes the product of 2 kernels that are defined on the same space
+
+    :param k1, k2: the kernels to multiply
+    :type k1, k2: kernpart
+    :rtype: kernel object
+
+    """
+    def __init__(self,k1,k2):
+        assert k1.D == k2.D, "Error: The input spaces of the kernels to multiply must have the same dimension"
+        self.D = k1.D
+        self.Nparam = k1.Nparam + k2.Nparam
+        self.name = k1.name + '<times>' + k2.name
+        self.k1 = k1
+        self.k2 = k2
+        self._set_params(np.hstack((k1._get_params(),k2._get_params())))
+
+    def _get_params(self):
+        """return the value of the parameters."""
+        return self.params
+
+    def _set_params(self,x):
+        """set the value of the parameters."""
+        self.k1._set_params(x[:self.k1.Nparam])
+        self.k2._set_params(x[self.k1.Nparam:])
+        self.params = x
+
+    def _get_param_names(self):
+        """return parameter names."""
+        return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
+
+    def K(self,X,X2,target):
+        """Compute the covariance matrix between X and X2."""
+        if X2 is None: X2 = X
+        target1 = np.zeros((X.shape[0],X2.shape[0]))
+        target2 = np.zeros((X.shape[0],X2.shape[0]))
+        self.k1.K(X,X2,target1)
+        self.k2.K(X,X2,target2)
+        target += target1 * target2
+
+    def Kdiag(self,X,target):
+        """Compute the diagonal of the covariance matrix associated to X."""
+        target1 = np.zeros((X.shape[0],))
+        target2 = np.zeros((X.shape[0],))
+        self.k1.Kdiag(X,target1)
+        self.k2.Kdiag(X,target2)
+        target += target1 * target2
+
+    def dK_dtheta(self,partial,X,X2,target):
+        """derivative of the covariance matrix with respect to the parameters."""
+        if X2 is None: X2 = X
+        K1 = np.zeros((X.shape[0],X2.shape[0]))
+        K2 = np.zeros((X.shape[0],X2.shape[0]))
+        self.k1.K(X,X2,K1)
+        self.k2.K(X,X2,K2)
+
+        k1_target = np.zeros(self.k1.Nparam)
+        k2_target = np.zeros(self.k2.Nparam)
+        self.k1.dK_dtheta(partial*K2, X, X2, k1_target)
+        self.k2.dK_dtheta(partial*K1, X, X2, k2_target)
+
+        target[:self.k1.Nparam] += k1_target
+        target[self.k1.Nparam:] += k2_target
+
+    def dK_dX(self,partial,X,X2,target):
+        """derivative of the covariance matrix with respect to X."""
+        if X2 is None: X2 = X
+        K1 = np.zeros((X.shape[0],X2.shape[0]))
+        K2 = np.zeros((X.shape[0],X2.shape[0]))
+        self.k1.K(X,X2,K1)
+        self.k2.K(X,X2,K2)
+
+        self.k1.dK_dX(partial*K2, X, X2, target)
+        self.k2.dK_dX(partial*K1, X, X2, target)
+
+    def dKdiag_dX(self,X,target):
+        pass

GPy/kern/product_orthogonal.py

@@ -35,7 +35,7 @@ class product_orthogonal(kernpart):
 
     def _get_param_names(self):
         """return parameter names."""
-        return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k1._get_param_names()]
+        return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
     
     def K(self,X,X2,target):
         """Compute the covariance matrix between X and X2."""
@@ -44,7 +44,7 @@ class product_orthogonal(kernpart):
         target2 = np.zeros((X.shape[0],X2.shape[0]))
         self.k1.K(X[:,0:self.k1.D],X2[:,0:self.k1.D],target1)
         self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],target2)
-        target += self.params[0]*target1 * target2
+        target += target1 * target2
 
     def Kdiag(self,X,target):
         """Compute the diagonal of the covariance matrix associated to X."""
@@ -52,7 +52,7 @@ class product_orthogonal(kernpart):
         target2 = np.zeros((X.shape[0],))
         self.k1.Kdiag(X[:,0:self.k1.D],target1)
         self.k2.Kdiag(X[:,self.k1.D:],target2)
-        target += self.params[0]*target1 * target2
+        target += target1 * target2
 
     def dK_dtheta(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to the parameters."""
@@ -67,13 +67,8 @@ class product_orthogonal(kernpart):
         self.k1.dK_dtheta(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], k1_target)
         self.k2.dK_dtheta(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], k2_target)
 
-        target[0] += np.sum(K1*K2*partial)
-        target[1:self.k1.Nparam] += self.params[0]* k1_target[1:]
-        target[self.k1.Nparam:] += self.params[0]* k2_target[1:]
-
-    def dKdiag_dtheta(self,partial,X,target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters."""
-        target[0] += 1
+        target[:self.k1.Nparam] += k1_target
+        target[self.k1.Nparam:] += k2_target
 
     def dK_dX(self,partial,X,X2,target):
         """derivative of the covariance matrix with respect to X."""