merged and updated slicing operations

GPy/kern/__init__.py

@@ -1,6 +1,6 @@
 from _src.kern import Kern
 from _src.rbf import RBF
-from _src.linear import Linear
+from _src.linear import Linear, LinearFull
 from _src.static import Bias, White
 from _src.brownian import Brownian
 from _src.sympykern import Sympykern

GPy/kern/_src/kern.py

@@ -11,7 +11,7 @@ from ...util.caching import Cache_this
 
 class Kern(Parameterized):
     #===========================================================================
-    # This adds input slice support. The rather ugly code for slicing can be 
+    # This adds input slice support. The rather ugly code for slicing can be
     # found in kernel_slice_operations
     __metaclass__ = KernCallsViaSlicerMeta
     #===========================================================================

GPy/kern/_src/kernel_slice_operations.py

@@ -5,30 +5,27 @@ Created on 11 Mar 2014
 '''
 from ...core.parameterization.parameterized import ParametersChangedMeta
 import numpy as np
-import functools
 
+def put_clean(dct, name, *args, **kw):
+    if name in dct:
+        dct['_clean_{}'.format(name)] = dct[name]
+        dct[name] = _slice_wrapper(None, dct[name], *args, **kw)
+    
 class KernCallsViaSlicerMeta(ParametersChangedMeta):
-    def __call__(self, *args, **kw):
-        instance = super(ParametersChangedMeta, self).__call__(*args, **kw)
-        instance.K = _Slice_wrapper(instance, instance.K)
-        instance.Kdiag = _Slice_wrapper_diag(instance, instance.Kdiag)
-
-        instance.update_gradients_full = _Slice_wrapper_derivative(instance, instance.update_gradients_full)
-        instance.update_gradients_diag = _Slice_wrapper_diag_derivative(instance, instance.update_gradients_diag)
-
-        instance.gradients_X = _Slice_wrapper_grad_X(instance, instance.gradients_X)
-        instance.gradients_X_diag = _Slice_wrapper_grad_X_diag(instance, instance.gradients_X_diag)
-
-        instance.psi0 = _Slice_wrapper(instance, instance.psi0)
-        instance.psi1 = _Slice_wrapper(instance, instance.psi1)
-        instance.psi2 = _Slice_wrapper(instance, instance.psi2)
-
-        instance.update_gradients_expectations = _Slice_wrapper_psi_stat_derivative_no_ret(instance, instance.update_gradients_expectations)
-        instance.gradients_Z_expectations = _Slice_wrapper_psi_stat_derivative_Z(instance, instance.gradients_Z_expectations)
-        instance.gradients_qX_expectations = _Slice_wrapper_psi_stat_derivative(instance, instance.gradients_qX_expectations)
-        instance.parameters_changed()
-        return instance
-
+    def __new__(cls, name, bases, dct):
+        put_clean(dct, 'K')
+        put_clean(dct, 'Kdiag', diag=True)
+        put_clean(dct, 'update_gradients_full', diag=False, derivative=True)
+        put_clean(dct, 'gradients_X', diag=False, derivative=True, ret_X=True)
+        put_clean(dct, 'gradients_X_diag', diag=True, derivative=True, ret_X=True)
+        put_clean(dct, 'psi0', diag=False, derivative=False)
+        put_clean(dct, 'psi1', diag=False, derivative=False)
+        put_clean(dct, 'psi2', diag=False, derivative=False)
+        put_clean(dct, 'update_gradients_expectations', derivative=True, psi_stat=True)
+        put_clean(dct, 'gradients_Z_expectations', derivative=True, psi_stat_Z=True, ret_X=True)
+        put_clean(dct, 'gradients_qX_expectations', derivative=True, psi_stat=True, ret_X=True)
+        return super(KernCallsViaSlicerMeta, cls).__new__(cls, name, bases, dct)
+    
 class _Slice_wrap(object):
     def __init__(self, instance, f):
         self.k = instance

GPy/kern/_src/linear.py

@@ -313,3 +313,47 @@ class Linear(Kern):
 
     def input_sensitivity(self):
         return np.ones(self.input_dim) * self.variances
+
+class LinearFull(Kern):
+    def __init__(self, input_dim, rank, W=None, kappa=None, active_dims=None, name='linear_full'):
+        super(LinearFull, self).__init__(input_dim, active_dims, name)
+        if W is None:
+            W = np.ones((input_dim, rank))
+        if kappa is None:
+            kappa = np.ones(input_dim)
+        assert W.shape == (input_dim, rank)
+        assert kappa.shape == (input_dim,)
+
+        self.W = Param('W', W)
+        self.kappa = Param('kappa', kappa, Logexp())
+        self.add_parameters(self.W, self.kappa)
+
+    def K(self, X, X2=None):
+        P = np.dot(self.W, self.W.T) + np.diag(self.kappa)
+        return np.einsum('ij,jk,lk->il', X, P, X if X2 is None else X2)
+
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        self.kappa.gradient = np.einsum('ij,ik,kj->j', X, dL_dK, X if X2 is None else X2)
+        self.W.gradient = np.einsum('ij,kl,ik,lm->jm', X, X if X2 is None else X2, dL_dK, self.W)
+        self.W.gradient += np.einsum('ij,kl,ik,jm->lm', X, X if X2 is None else X2, dL_dK, self.W)
+
+    def Kdiag(self, X):
+        P = np.dot(self.W, self.W.T) + np.diag(self.kappa)
+        return np.einsum('ij,jk,ik->i', X, P, X)
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        self.kappa.gradient = np.einsum('ij,i->j', np.square(X), dL_dKdiag)
+        self.W.gradient = 2.*np.einsum('ij,ik,jl,i->kl', X, X, self.W, dL_dKdiag)
+
+    def gradients_X(self, dL_dK, X, X2=None):
+        P = np.dot(self.W, self.W.T) + np.diag(self.kappa)
+        if X2 is None:
+            return 2.*np.einsum('ij,jk,kl->il', dL_dK, X, P)
+        else:
+            return np.einsum('ij,jk,kl->il', dL_dK, X2, P)
+
+    def gradients_X_diag(self, dL_dKdiag, X):
+        P = np.dot(self.W, self.W.T) + np.diag(self.kappa)
+        return 2.*np.einsum('jk,i,ij->ik', P, dL_dKdiag, X)
+
+

GPy/kern/_src/rbf.py

@@ -64,7 +64,7 @@ class RBF(Stationary):
             if self.ARD:
                 self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).reshape(-1,self.input_dim).sum(axis=0)
             else:
-                self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).sum()  
+                self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).sum()
 
             #from psi2
             self.variance.gradient += (dL_dpsi2 * _dpsi2_dvariance).sum()