Added Multioutput derivative kernel for adding derivatives easy and changed it to default kernel in multioutput gp model

GPy/kern/__init__.py

@@ -44,4 +44,5 @@ from .src.sde_static import sde_White, sde_Bias
 from .src.sde_stationary import sde_RBF,sde_Exponential,sde_RatQuad
 from .src.sde_brownian import sde_Brownian
 from .src.multioutput_kern import MultioutputKern
+from .src.multioutput_derivative_kern import MultioutputDerivativeKern
 from .src.diff_kern import DiffKern

GPy/kern/src/multioutput_derivative_kern.py

@@ -0,0 +1,81 @@
+from .kern import Kern, CombinationKernel
+from .multioutput_kern import MultioutputKern, ZeroKern
+import numpy as np
+
+class KernWrapper(Kern):
+    def __init__(self, fk, fug, fg, base_kern):
+        self.fk = fk
+        self.fug = fug
+        self.fg = fg
+        self.base_kern
+        super(KernWrapper, self).__init__(1, None, name='KernWrapper',useGPU=False)
+
+    def K(self, X, X2=None):
+        return self.fk(X,X2=X2)
+    
+    def update_gradients_full(self,dL_dK, X, X2=None):
+        return self.fug(dK_dK, X, X2=X2)
+    
+    def gradients_X(self,dL_dK, X, X2=None):
+        return self.fg(dL_dK, X, X2=X2)
+
+    def get_gradient(self):
+        return self.base_kern.gradient.copy()
+
+    def append_gradient(self, gradient):
+        self.base_kern.gradient += gradient
+
+class MultioutputDerivativeKern(MultioutputKern):
+    """
+    Multioutput derivative kernel is a meta class for combining different kernels for multioutput GPs.
+    Multioutput derivative kernel is only a thin wrapper for Multioutput kernel for user not having to define
+    cross covariances.
+    """
+    def __init__(self, kernels, cross_covariances={}, name='MultioutputDerivativeKern'):
+        #kernels contains a list of kernels as input, 
+        if not isinstance(kernels, list):
+            self.single_kern = True
+            self.kern = kernels
+            kernels = [kernels]
+        else:
+            self.single_kern = False
+            self.kern = kernels
+        # The combination kernel ALLWAYS puts the extra dimension last.
+        # Thus, the index dimension of this kernel is always the last dimension
+        # after slicing. This is why the index_dim is just the last column:
+        self.index_dim = -1
+        super(MultioutputKern, self).__init__(kernels=kernels, extra_dims=[self.index_dim], name=name, link_parameters=False)
+
+        nl = len(kernels)
+         
+        #build covariance structure
+        covariance = [[None for i in range(nl)] for j in range(nl)]
+        linked = []
+        for i in range(0,nl):
+            unique=True
+            for j in range(0,nl):
+                if i==j or (kernels[i] is kernels[j]):
+                    covariance[i][j] = {'kern': kernels[i],'K':kernels[i].K, 'update_gradients_full': kernels[i].update_gradients_full, 'gradients_X': kernels[i].gradients_X}
+                    if i>j:
+                        unique=False
+                elif cross_covariances.get((i,j)) is not None: #cross covariance is given
+                    covariance[i][j] = cross_covariances.get((i,j))
+                elif kernels[i].name == 'diffKern' and kernels[i].base_kern == kernels[j]: # one is derivative of other
+                    kern = KernWrapper(kernels[i].dK_dX_wrap,kernels[i].update_gradients_dK_dX,kernels[i].gradients_X, kernels[j])
+                    covariance[i][j] = {'kern':kern, 'K': kern.K, 'update_gradients_full': kern.update_gradients_full, 'gradients_X': kern.gradients_X}
+                    unique=False
+                elif kernels[j].name == 'diffKern' and kernels[j].base_kern == kernels[i]: # one is derivative of other
+                    kern = KernWrapper(kernels[j].dK_dX2_wrap,kernels[j].update_gradients_dK_dX2,kernels[j].gradients_X2, kernels[i])
+                    covariance[i][j] = {'kern':kern, 'K': kern.K, 'update_gradients_full': kern.update_gradients_full, 'gradients_X': kern.gradients_X}
+                elif kernels[i].name == 'diffKern' and kernels[j].name == 'diffKern' and kernels[i].base_kern == kernels[j].base_kern: #both are partial derivatives
+                    kern = KernWrapper(partial(kernels[i].K, dimX2=kernels[j].dimension), partial(kernels[i].update_gradients_full, dimX2=kernels[j].dimension),None, kernels[i].base_kern)
+                    covariance[i][j] = {'kern':kern, 'K': kern.K, 'update_gradients_full': kern.update_gradients_full, 'gradients_X': kern.gradients_X}
+                    if i>j:
+                        unique=False
+                else:
+                    kern = ZeroKern()
+                    covariance[i][j] = {'kern': kern, 'K': kern.K, 'update_gradients_full': kern.update_gradients_full, 'gradients_X': kern.gradients_X}   
+            if unique is True:
+                linked.append(i)
+        self.covariance = covariance
+        self.link_parameters(*[kernels[i] for i in linked])

GPy/models/multioutput_gp.py

@@ -42,7 +42,7 @@ class MultioutputGP(GP):
         Ny = len(Y_list)
         
         assert isinstance(kernel_list, list)
-        kernel = kern.MultioutputKern(kernels=kernel_list, cross_covariances=kernel_cross_covariances)
+        kernel = kern.MultioutputDerivativeKern(kernels=kernel_list, cross_covariances=kernel_cross_covariances)
 
         assert isinstance(likelihood_list, list)
         likelihood = likelihoods.MultioutputLikelihood(likelihood_list)
@@ -53,8 +53,7 @@ class MultioutputGP(GP):
             else:
                 inference_method = expectation_propagation.EP() 
         
-        super(MultioutputGP, self).__init__(X,Y,kernel,likelihood, Y_metadata={'output_index':self.output_index, 'trials':np.ones(self.output_index.shape)}, inference_method = inference_method)# expectation_propagation.MultioutputEP()) # expectation_propagation.EP())                             
+        super(MultioutputGP, self).__init__(X,Y,kernel,likelihood, Y_metadata={'output_index':self.output_index, 'trials':np.ones(self.output_index.shape)}, inference_method = inference_method)
-                                            #expectation_propagation.MultioutputEP())
 
     def predict_noiseless(self,  Xnew, full_cov=False, Y_metadata=None, kern=None):
         if isinstance(Xnew, list):