doc style change

GPy/kern/Matern52.py

@@ -31,14 +31,14 @@ class Matern52(kernpart):
         self.ARD = ARD
         if ARD == False:
             self.Nparam = 2
-            self.name = 'Mat32'
+            self.name = 'Mat52'
             if lengthscale is not None:
                 assert lengthscale.shape == (1,)
             else:
                 lengthscale = np.ones(1)
         else:
             self.Nparam = self.D + 1
-            self.name = 'Mat32_ARD'
+            self.name = 'Mat52_ARD'
             if lengthscale is not None:
                 assert lengthscale.shape == (self.D,)
             else:

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
+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 kern import kern

GPy/kern/constructors.py

@@ -18,7 +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_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.
 
@@ -241,3 +241,14 @@ 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_orthogonal(k1,k2):
+    """
+     Construct a product kernel
+     
+    :param k1, k2: the kernels to multiply
+    :type k1, k2: kernpart
+    :rtype: kernel object
+    """
+    part = product_orthogonalpart(k1,k2)
+    return kern(k1.D+k2.D, [part])

GPy/kern/kern.py

@@ -4,9 +4,9 @@
 
 import numpy as np
 from ..core.parameterised import parameterised
-from functools import partial
 from kernpart import kernpart
 import itertools
+from product_orthogonal import product_orthogonal 
 
 class kern(parameterised):
     def __init__(self,D,parts=[], input_slices=None):
@@ -45,7 +45,6 @@ class kern(parameterised):
         for p in self.parts:
             assert isinstance(p,kernpart), "bad kernel part"
 
-
         self.compute_param_slices()
 
         parameterised.__init__(self)
@@ -133,6 +132,67 @@ class kern(parameterised):
         newkern.tied_indices = self.tied_indices + [self.Nparam + x for x in other.tied_indices]
         return newkern
 
+    def __mul__(self,other):
+        """
+        Shortcut for `prod_orthogonal`. Note that `+` assumes that we sum 2 kernels defines on the same space whereas `*` assumes that the kernels are defined on different subspaces.
+        """
+        return self.prod_orthogonal(other)
+
+    def prod_orthogonal(self,other):
+        """
+        multiply two kernels. Both kernels are defined on separate spaces. Note that the constrains on the parameters of the kernels to multiply will be lost.
+        :param other: the other kernel to be added
+        :type other: GPy.kern
+        """
+        K1 = self.copy()
+        K2 = other.copy()
+        K1.unconstrain('')
+        K2.unconstrain('')
+
+        prev_ties = K1.tied_indices + [arr + K1.Nparam for arr in K2.tied_indices]
+        K1.untie_everything()
+        K2.untie_everything()
+
+        D = K1.D + K2.D
+
+        newkernparts = [product_orthogonal(k1,k2).parts[0] for k1, k2 in itertools.product(K1.parts,K2.parts)]
+
+        slices = []
+        for sl1, sl2 in itertools.product(K1.input_slices,K2.input_slices):
+            s1, s2 = [False]*K1.D, [False]*K2.D
+            s1[sl1], s2[sl2] = [True], [True]
+            slices += [s1+s2]
+
+        newkern = kern(D, newkernparts, slices)
+
+        # create the ties
+        K1_param = []
+        n = 0
+        for k1 in K1.parts:
+            K1_param += [range(n,n+k1.Nparam)]
+            n += k1.Nparam
+        n = 0
+        K2_param = []
+        for k2 in K2.parts:
+            K2_param += [range(K1.Nparam+n,K1.Nparam+n+k2.Nparam)]
+            n += k2.Nparam
+        index_param = []
+        for p1 in K1_param:
+            for p2 in K2_param:
+                index_param += [0] + p1[1:] + p2[1:]
+        index_param = np.array(index_param)
+
+        # follow the previous ties
+        for arr in prev_ties:
+            for j in arr:
+                index_param[np.where(index_param==j)[0]] = arr[0]
+
+        # tie
+        for i in np.unique(index_param)[1:]:
+            newkern.tie_param(np.where(index_param==i)[0])
+
+        return newkern
+
     def _get_params(self):
         return np.hstack([p._get_params() for p in self.parts])
 

GPy/kern/product_orthogonal.py

@@ -0,0 +1,91 @@
+# 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
+
+class product_orthogonal(kernpart):
+    """
+    Computes the product of 2 kernels
+
+    :param k1, k2: the kernels to multiply
+    :type k1, k2: kernpart
+    :rtype: kernel object
+
+    """
+    def __init__(self,k1,k2):
+        assert k1._get_param_names()[0] == 'variance' and k2._get_param_names()[0] == 'variance', "Error: The multipication of kernels is only defined when the first parameters of the kernels to multiply is the variance."
+        self.D = k1.D + k2.D
+        self.Nparam = k1.Nparam + k2.Nparam - 1
+        self.name = k1.name + '<times>' + k2.name
+        self.k1 = k1
+        self.k2 = k2
+        self._set_params(np.hstack((k1._get_params()[0]*k2._get_params()[0], k1._get_params()[1:],k2._get_params()[1:])))
+
+    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(np.hstack((1.,x[1:self.k1.Nparam])))
+        self.k2._set_params(np.hstack((1.,x[self.k1.Nparam:])))
+        self.params = x
+
+    def _get_param_names(self):
+        """return parameter names."""
+        return ['variance']+[self.k1.name + '_' + self.k1._get_param_names()[i+1] for i in range(self.k1.Nparam-1)] +  [self.k2.name + '_' + self.k2._get_param_names()[i+1] for i in range(self.k2.Nparam-1)]
+
+    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[:,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
+
+    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[:,0:self.k1.D],target1)
+        self.k2.Kdiag(X[:,self.k1.D:],target2)
+        target += self.params[0]*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[:,0:self.k1.D],X2[:,0:self.k1.D],K1)
+        self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],K2)
+
+        k1_target = np.zeros(self.k1.Nparam)
+        k2_target = np.zeros(self.k2.Nparam)
+        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
+
+    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[:,0:self.k1.D],X2[:,0:self.k1.D],K1)
+        self.k2.K(X[:,self.k1.D:],X2[:,self.k1.D:],K2)
+
+        self.k1.dK_dX(partial*K2, X[:,:self.k1.D], X2[:,:self.k1.D], target)
+        self.k2.dK_dX(partial*K1, X[:,self.k1.D:], X2[:,self.k1.D:], target)
+
+    def dKdiag_dX(self,X,target):
+        pass

GPy/kern/rbf.py

@@ -26,28 +26,29 @@ class rbf(kernpart):
     :type ARD: Boolean
     :rtype: kernel object
 
-    .. Note: for rbf with different lengthscale on each dimension, see rbf_ARD
+    .. Note: this object implements both the ARD and 'spherical' version of the function
     """
 
     def __init__(self,D,variance=1.,lengthscale=None,ARD=False):
         self.D = D
+        self.name = 'rbf'
         self.ARD = ARD
         if not ARD:
             self.Nparam = 2
-            self.name = 'rbf'
             if lengthscale is not None:
-                assert lengthscale.shape == (1,)
+                lengthscale = np.asarray(lengthscale)
+                assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel"
             else:
                 lengthscale = np.ones(1)
         else:
             self.Nparam = self.D + 1
-            self.name = 'rbf_ARD'
             if lengthscale is not None:
-                assert lengthscale.shape == (self.D,)
+                lengthscale = np.asarray(lengthscale)
+                assert lengthscale.size == self.D, "bad number of lengthscales"
             else:
                 lengthscale = np.ones(self.D)
 
-        self._set_params(np.hstack((variance,lengthscale)))
+        self._set_params(np.hstack((variance,lengthscale.flatten())))
 
         #initialize cache
         self._Z, self._mu, self._S = np.empty(shape=(3,1))