merged static

GPy/kern/__init__.py

@@ -1,34 +1,26 @@
-from _src.rbf import RBF
-from _src.white import White
 from _src.kern import Kern
+from _src.rbf import RBF
 from _src.linear import Linear
-from _src.bias import Bias
+from _src.static import Bias, White
 from _src.brownian import Brownian
-from _src.stationary import Exponential, Matern32, Matern52, ExpQuad
+from _src.stationary import Exponential, Matern32, Matern52, ExpQuad, RatQuad, Cosine
+from _src.mlp import MLP
+from _src.periodic import PeriodicExponential, PeriodicMatern32, PeriodicMatern52
+from _src.independent_outputs import IndependentOutputs
 #import coregionalize
-#import exponential
 #import eq_ode1
 #import finite_dimensional
 #import fixed
 #import gibbs
 #import hetero
 #import hierarchical
-#import independent_outputs
-#import linear
-#import Matern32
-#import Matern52
-#import mlp
 #import ODE_1
 #import periodic_exponential
 #import periodic_Matern32
 #import periodic_Matern52
 #import poly
-#import prod_orthogonal
-#import prod
-#import rational_quadratic
 #import rbfcos
 #import rbf
 #import rbf_inv
 #import spline
 #import symmetric
-#import white

GPy/kern/_src/independent_outputs.py

@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from kernpart import Kernpart
+from kern import Kern
 import numpy as np
 
 def index_to_slices(index):
@@ -31,67 +31,89 @@ def index_to_slices(index):
     [ret[ind_i].append(slice(*indexes_i)) for ind_i,indexes_i in zip(ind[switchpoints[:-1]],zip(switchpoints,switchpoints[1:]))]
     return ret
 
-class IndependentOutputs(Kernpart):
+class IndependentOutputs(Kern):
     """
-    A kernel part shich can reopresent several independent functions.
+    A kernel which can reopresent several independent functions.
     this kernel 'switches off' parts of the matrix where the output indexes are different.
 
     The index of the functions is given by the last column in the input X
-    the rest of the columns of X are passed to the kernel for computation (in blocks).
+    the rest of the columns of X are passed to the underlying kernel for computation (in blocks).
 
     """
-    def __init__(self,k):
+    def __init__(self, kern, name='independ'):
-        self.input_dim = k.input_dim + 1
+        super(IndependentOutputs, self).__init__(kern.input_dim+1, name)
-        self.num_params = k.num_params
+        self.kern = kern
-        self.name = 'iops('+ k.name + ')'
+        self.add_parameters(self.kern)
-        self.k = k
 
-    def _get_params(self):
+    def K(self,X ,X2=None):
-        return self.k._get_params()
+        X, slices = X[:,:-1], index_to_slices(X[:,-1])
+        if X2 is None:
+            target = np.zeros((X.shape[0], X.shape[0]))
+            [[np.copyto(target[s,s], self.kern.K(X[s], None)) for s in slices_i] for slices_i in slices]
+        else:
+            X2, slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
+            target = np.zeros((X.shape[0], X2.shape[0]))
+            [[[np.copyto(target[s, s2], self.kern.K(X[s],X2[s2])) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+        return target
 
-    def _set_params(self,x):
+    def Kdiag(self,X):
-        self.k._set_params(x)
+        X, slices = X[:,:-1], index_to_slices(X[:,-1])
-        self.params = x
+        target = np.zeros(X.shape[0])
+        [[np.copyto(target[s], self.kern.Kdiag(X[s])) for s in slices_i] for slices_i in slices]
+        return target
 
-    def _get_param_names(self):
+    def update_gradients_full(self,dL_dK,X,X2=None):
-        return self.k._get_param_names()
+        target = np.zeros(self.kern.size)
+        def collate_grads(dL, X, X2):
+            self.kern.update_gradients_full(dL,X,X2)
+            self.kern._collect_gradient(target)
 
-    def K(self,X,X2,target):
-        #Sort out the slices from the input data
         X,slices = X[:,:-1],index_to_slices(X[:,-1])
         if X2 is None:
-            X2,slices2 = X,slices
+            [[collate_grads(dL_dK[s,s], X[s], None) for s in slices_i] for slices_i in slices]
         else:
-            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
+            X2, slices2 = X2[:,:-1], index_to_slices(X2[:,-1])
+            [[[collate_grads(dL_dK[s,s2],X[s],X2[s2]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
 
-        [[[self.k.K(X[s],X2[s2],target[s,s2]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+        self.kern._set_gradient(target)
 
-    def Kdiag(self,X,target):
+    def gradients_X(self,dL_dK, X, X2=None):
-        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+        target = np.zeros_like(X)
-        [[self.k.Kdiag(X[s],target[s]) for s in slices_i] for slices_i in slices]
+        X, slices = X[:,:-1],index_to_slices(X[:,-1])
-
-    def _param_grad_helper(self,dL_dK,X,X2,target):
-        X,slices = X[:,:-1],index_to_slices(X[:,-1])
         if X2 is None:
-            X2,slices2 = X,slices
+            [[np.copyto(target[s,:-1], self.kern.gradients_X(dL_dK[s,s],X[s],None)) for s in slices_i] for slices_i in slices]
         else:
             X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
-        [[[self.k._param_grad_helper(dL_dK[s,s2],X[s],X2[s2],target) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+            [[[np.copyto(target[s,:-1], self.kern.gradients_X(dL_dK[s,s2], X[s], X2[s2])) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
+        return target
 
+    def gradients_X_diag(self, dL_dKdiag, X):
+        X, slices = X[:,:-1], index_to_slices(X[:,-1])
+        target = np.zeros(X.shape)
+        [[np.copyto(target[s,:-1], self.kern.gradients_X_diag(dL_dKdiag[s],X[s])) for s in slices_i] for slices_i in slices]
+        return target
 
-    def gradients_X(self,dL_dK,X,X2,target):
+    def update_gradients_diag(self,dL_dKdiag,X,target):
+        target = np.zeros(self.kern.size)
+        def collate_grads(dL, X):
+            self.kern.update_gradients_diag(dL,X)
+            self.kern._collect_gradient(target)
         X,slices = X[:,:-1],index_to_slices(X[:,-1])
-        if X2 is None:
+        [[collate_grads(dL_dKdiag[s], X[s,:]) for s in slices_i] for slices_i in slices]
-            X2,slices2 = X,slices
+        self.kern._set_gradient(target)
-        else:
-            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
-        [[[self.k.gradients_X(dL_dK[s,s2],X[s],X2[s2],target[s,:-1]) for s in slices_i] for s2 in slices_j] for slices_i,slices_j in zip(slices,slices2)]
 
-    def dKdiag_dX(self,dL_dKdiag,X,target):
+def Hierarchical(kern_f, kern_g, name='hierarchy'):
-        X,slices = X[:,:-1],index_to_slices(X[:,-1])
+    """
-        [[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target[s,:-1]) for s in slices_i] for slices_i in slices]
+    A kernel which can reopresent a simple hierarchical model.
 
+    See Hensman et al 2013, "Hierarchical Bayesian modelling of gene expression time
+    series across irregularly sampled replicates and clusters"
+    http://www.biomedcentral.com/1471-2105/14/252
+
+    The index of the functions is given by the last column in the input X
+    the rest of the columns of X are passed to the underlying kernel for computation (in blocks).
+
+    """
+    assert kern_f.input_dim == kern_g.input_dim
+    return kern_f + IndependentOutputs(kern_g)
 
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        X,slices = X[:,:-1],index_to_slices(X[:,-1])
-        [[self.k.dKdiag_dX(dL_dKdiag[s],X[s],target) for s in slices_i] for slices_i in slices]

GPy/kern/_src/kern.py

@@ -105,7 +105,7 @@ class Kern(Parameterized):
         """ Here we overload the '*' operator. See self.prod for more information"""
         return self.prod(other)
 
-    def __pow__(self, other, tensor=False):
+    def __pow__(self, other):
         """
         Shortcut for tensor `prod`.
         """

GPy/kern/_src/mlp.py

@@ -1,11 +1,13 @@
 # Copyright (c) 2013, GPy authors (see AUTHORS.txt).
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
-from kernpart import Kernpart
+from kern import Kern
+from ...core.parameterization import Param
+from ...core.parameterization.transformations import Logexp
 import numpy as np
 four_over_tau = 2./np.pi
 
-class MLP(Kernpart):
+class MLP(Kern):
     """
 
     Multi layer perceptron kernel (also known as arc sine kernel or neural network kernel)
@@ -13,10 +15,10 @@ class MLP(Kernpart):
     .. math::
 
           k(x,y) = \\sigma^{2}\\frac{2}{\\pi }  \\text{asin} \\left ( \\frac{ \\sigma_w^2 x^\\top y+\\sigma_b^2}{\\sqrt{\\sigma_w^2x^\\top x + \\sigma_b^2 + 1}\\sqrt{\\sigma_w^2 y^\\top y \\sigma_b^2 +1}} \\right )
-          
+
 
     :param input_dim: the number of input dimensions
-    :type input_dim: int 
+    :type input_dim: int
     :param variance: the variance :math:`\sigma^2`
     :type variance: float
     :param weight_variance: the vector of the variances of the prior over input weights in the neural network :math:`\sigma^2_w`
@@ -29,85 +31,58 @@ class MLP(Kernpart):
 
     """
 
-    def __init__(self, input_dim, variance=1., weight_variance=None, bias_variance=100., ARD=False):
+    def __init__(self, input_dim, variance=1., weight_variance=1., bias_variance=100., name='mlp'):
-        self.input_dim = input_dim
+        super(MLP, self).__init__(input_dim, name)
-        self.ARD = ARD
+        self.variance = Param('variance', variance, Logexp())
-        if not ARD:
+        self.weight_variance = Param('weight_variance', weight_variance, Logexp())
-            self.num_params=3
+        self.bias_variance = Param('bias_variance', bias_variance, Logexp())
-            if weight_variance is not None:
+        self.add_parameters(self.variance, self.weight_variance, self.bias_variance)
-                weight_variance = np.asarray(weight_variance)
-                assert weight_variance.size == 1, "Only one weight variance needed for non-ARD kernel"
-            else:
-                weight_variance = 100.*np.ones(1)
-        else:
-            self.num_params = self.input_dim + 2
-            if weight_variance is not None:
-                weight_variance = np.asarray(weight_variance)
-                assert weight_variance.size == self.input_dim, "bad number of weight variances"
-            else:
-                weight_variance = np.ones(self.input_dim)
-            raise NotImplementedError
 
-        self.name='mlp'
-        self._set_params(np.hstack((variance, weight_variance.flatten(), bias_variance)))
 
-    def _get_params(self):
+    def K(self, X, X2=None):
-        return np.hstack((self.variance, self.weight_variance.flatten(), self.bias_variance))
-
-    def _set_params(self, x):
-        assert x.size == (self.num_params)
-        self.variance = x[0]
-        self.weight_variance = x[1:-1]
-        self.weight_std = np.sqrt(self.weight_variance)
-        self.bias_variance = x[-1]
-
-    def _get_param_names(self):
-        if self.num_params == 3:
-            return ['variance', 'weight_variance', 'bias_variance']
-        else:
-            return ['variance'] + ['weight_variance_%i' % i for i in range(self.lengthscale.size)] + ['bias_variance']
-
-    def K(self, X, X2, target):
-        """Return covariance between X and X2."""
         self._K_computations(X, X2)
-        target += self.variance*self._K_dvar
+        return self.variance*self._K_dvar
 
-    def Kdiag(self, X, target):
+    def Kdiag(self, X):
         """Compute the diagonal of the covariance matrix for X."""
         self._K_diag_computations(X)
-        target+= self.variance*self._K_diag_dvar
+        return self.variance*self._K_diag_dvar
 
-    def _param_grad_helper(self, dL_dK, X, X2, target):
+    def update_gradients_full(self, dL_dK, X, X2=None):
         """Derivative of the covariance with respect to the parameters."""
         self._K_computations(X, X2)
-        denom3 = self._K_denom*self._K_denom*self._K_denom
+        self.variance.gradient = np.sum(self._K_dvar*dL_dK)
+
+        denom3 = self._K_denom**3
         base = four_over_tau*self.variance/np.sqrt(1-self._K_asin_arg*self._K_asin_arg)
         base_cov_grad = base*dL_dK
 
         if X2 is None:
             vec = np.diag(self._K_inner_prod)
-            target[1] += ((self._K_inner_prod/self._K_denom 
+            self.weight_variance.gradient = ((self._K_inner_prod/self._K_denom
                            -.5*self._K_numer/denom3
-                           *(np.outer((self.weight_variance*vec+self.bias_variance+1.), vec) 
+                           *(np.outer((self.weight_variance*vec+self.bias_variance+1.), vec)
                              +np.outer(vec,(self.weight_variance*vec+self.bias_variance+1.))))*base_cov_grad).sum()
-            target[2] += ((1./self._K_denom 
+            self.bias_variance.gradient = ((1./self._K_denom
-                           -.5*self._K_numer/denom3 
+                           -.5*self._K_numer/denom3
                            *((vec[None, :]+vec[:, None])*self.weight_variance
                            +2.*self.bias_variance + 2.))*base_cov_grad).sum()
         else:
             vec1 = (X*X).sum(1)
             vec2 = (X2*X2).sum(1)
-            target[1] += ((self._K_inner_prod/self._K_denom 
+            self.weight_variance.gradient = ((self._K_inner_prod/self._K_denom
                            -.5*self._K_numer/denom3
                            *(np.outer((self.weight_variance*vec1+self.bias_variance+1.), vec2) + np.outer(vec1, self.weight_variance*vec2 + self.bias_variance+1.)))*base_cov_grad).sum()
-            target[2] += ((1./self._K_denom 
+            self.bias_variance.gradient = ((1./self._K_denom
-                           -.5*self._K_numer/denom3 
+                           -.5*self._K_numer/denom3
                            *((vec1[:, None]+vec2[None, :])*self.weight_variance
                              + 2*self.bias_variance + 2.))*base_cov_grad).sum()
-            
-        target[0] += np.sum(self._K_dvar*dL_dK)
 
-    def gradients_X(self, dL_dK, X, X2, target):
+    def update_gradients_diag(self, X):
+        raise NotImplementedError, "TODO"
+
+
+    def gradients_X(self, dL_dK, X, X2):
         """Derivative of the covariance matrix with respect to X"""
         self._K_computations(X, X2)
         arg = self._K_asin_arg
@@ -116,47 +91,39 @@ class MLP(Kernpart):
         denom3 = denom*denom*denom
         if X2 is not None:
             vec2 = (X2*X2).sum(1)*self.weight_variance+self.bias_variance + 1.
-            target += four_over_tau*self.weight_variance*self.variance*((X2[None, :, :]/denom[:, :, None] - vec2[None, :, None]*X[:, None, :]*(numer/denom3)[:, :, None])*(dL_dK/np.sqrt(1-arg*arg))[:, :, None]).sum(1)
+            return four_over_tau*self.weight_variance*self.variance*((X2[None, :, :]/denom[:, :, None] - vec2[None, :, None]*X[:, None, :]*(numer/denom3)[:, :, None])*(dL_dK/np.sqrt(1-arg*arg))[:, :, None]).sum(1)
         else:
             vec = (X*X).sum(1)*self.weight_variance+self.bias_variance + 1.
-            target += 2*four_over_tau*self.weight_variance*self.variance*((X[None, :, :]/denom[:, :, None] - vec[None, :, None]*X[:, None, :]*(numer/denom3)[:, :, None])*(dL_dK/np.sqrt(1-arg*arg))[:, :, None]).sum(1)
+            return 2*four_over_tau*self.weight_variance*self.variance*((X[None, :, :]/denom[:, :, None] - vec[None, :, None]*X[:, None, :]*(numer/denom3)[:, :, None])*(dL_dK/np.sqrt(1-arg*arg))[:, :, None]).sum(1)
-            
+
     def dKdiag_dX(self, dL_dKdiag, X, target):
         """Gradient of diagonal of covariance with respect to X"""
         self._K_diag_computations(X)
         arg = self._K_diag_asin_arg
         denom = self._K_diag_denom
         numer = self._K_diag_numer
-        target += four_over_tau*2.*self.weight_variance*self.variance*X*(1/denom*(1 - arg)*dL_dKdiag/(np.sqrt(1-arg*arg)))[:, None] 
+        return four_over_tau*2.*self.weight_variance*self.variance*X*(1./denom*(1. - arg)*dL_dKdiag/(np.sqrt(1-arg*arg)))[:, None]
+
 
-    
     def _K_computations(self, X, X2):
         """Pre-computations for the covariance matrix (used for computing the covariance and its gradients."""
-        if self.ARD:
+        if X2 is None:
-            pass
+            self._K_inner_prod = np.dot(X,X.T)
+            self._K_numer = self._K_inner_prod*self.weight_variance + self.bias_variance
+            vec = np.diag(self._K_numer) + 1.
+            self._K_denom = np.sqrt(np.outer(vec,vec))
         else:
-            if X2 is None:
+            self._K_inner_prod = np.dot(X,X2.T)
-                self._K_inner_prod = np.dot(X,X.T)
+            self._K_numer = self._K_inner_prod*self.weight_variance + self.bias_variance
-                self._K_numer = self._K_inner_prod*self.weight_variance+self.bias_variance
+            vec1 = (X*X).sum(1)*self.weight_variance + self.bias_variance + 1.
-                vec = np.diag(self._K_numer) + 1.
+            vec2 = (X2*X2).sum(1)*self.weight_variance + self.bias_variance + 1.
-                self._K_denom = np.sqrt(np.outer(vec,vec))
+            self._K_denom = np.sqrt(np.outer(vec1,vec2))
-                self._K_asin_arg = self._K_numer/self._K_denom
+        self._K_asin_arg = self._K_numer/self._K_denom
-                self._K_dvar = four_over_tau*np.arcsin(self._K_asin_arg)
+        self._K_dvar = four_over_tau*np.arcsin(self._K_asin_arg)
-            else:
-                self._K_inner_prod = np.dot(X,X2.T)
-                self._K_numer = self._K_inner_prod*self.weight_variance + self.bias_variance
-                vec1 = (X*X).sum(1)*self.weight_variance + self.bias_variance + 1.
-                vec2 = (X2*X2).sum(1)*self.weight_variance + self.bias_variance + 1.
-                self._K_denom = np.sqrt(np.outer(vec1,vec2))
-                self._K_asin_arg = self._K_numer/self._K_denom
-                self._K_dvar = four_over_tau*np.arcsin(self._K_asin_arg)
 
     def _K_diag_computations(self, X):
         """Pre-computations concerning the diagonal terms (used for computation of diagonal and its gradients)."""
-        if self.ARD:
+        self._K_diag_numer = (X*X).sum(1)*self.weight_variance + self.bias_variance
-            pass
+        self._K_diag_denom = self._K_diag_numer+1.
-        else:
+        self._K_diag_asin_arg = self._K_diag_numer/self._K_diag_denom
-            self._K_diag_numer = (X*X).sum(1)*self.weight_variance + self.bias_variance
+        self._K_diag_dvar = four_over_tau*np.arcsin(self._K_diag_asin_arg)
-            self._K_diag_denom = self._K_diag_numer+1.
-            self._K_diag_asin_arg = self._K_diag_numer/self._K_diag_denom
-            self._K_diag_dvar = four_over_tau*np.arcsin(self._K_diag_asin_arg)

GPy/kern/_src/periodic.py

@@ -2,12 +2,287 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 
-from kernpart import Kernpart
 import numpy as np
-from GPy.util.linalg import mdot
+from kern import Kern
-from GPy.util.decorators import silence_errors
+from ...util.linalg import mdot
+from ...util.decorators import silence_errors
+from ...core.parameterization.param import Param
+from ...core.parameterization.transformations import Logexp
 
-class PeriodicMatern52(Kernpart):
+class Periodic(Kern):
+    def __init__(self, input_dim, variance, lengthscale, period, n_freq, lower, upper, name):
+        """
+        :type input_dim: int
+        :param variance: the variance of the Matern kernel
+        :type variance: float
+        :param lengthscale: the lengthscale of the Matern kernel
+        :type lengthscale: np.ndarray of size (input_dim,)
+        :param period: the period
+        :type period: float
+        :param n_freq: the number of frequencies considered for the periodic subspace
+        :type n_freq: int
+        :rtype: kernel object
+        """
+
+        assert input_dim==1, "Periodic kernels are only defined for input_dim=1"
+        super(Periodic, self).__init__(input_dim, name)
+        self.input_dim = input_dim
+        self.lower,self.upper = lower, upper
+        self.n_freq = n_freq
+        self.n_basis = 2*n_freq
+        self.variance = Param('variance', np.float64(variance), Logexp())
+        self.lengthscale = Param('lengthscale', np.float64(lengthscale), Logexp())
+        self.period = Param('period', np.float64(period), Logexp())
+        self.add_parameters(self.variance, self.lengthscale, self.period)
+        self.parameters_changed()
+
+    def _cos(self, alpha, omega, phase):
+        def f(x):
+            return alpha*np.cos(omega*x + phase)
+        return f
+
+    @silence_errors
+    def _cos_factorization(self, alpha, omega, phase):
+        r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
+        r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
+        r =  np.sqrt(r1**2 + r2**2)
+        psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
+        return r,omega[:,0:1], psi
+
+    @silence_errors
+    def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
+        Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
+        Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) +  np.cos(phi1-phi2.T)*(self.upper-self.lower)
+        Gint = np.dot(r1,r2.T)/2 * np.where(np.isnan(Gint1),Gint2,Gint1)
+        return Gint
+
+    def K(self, X, X2=None):
+        FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
+        if X2 is None:
+            FX2 = FX
+        else:
+            FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
+        return mdot(FX,self.Gi,FX2.T)
+
+    def Kdiag(self,X):
+        return np.diag(self.K(X))
+
+
+
+
+class PeriodicExponential(Periodic):
+    """
+    Kernel of the periodic subspace (up to a given frequency) of a exponential
+    (Matern 1/2) RKHS.
+
+    Only defined for input_dim=1.
+    """
+
+    def __init__(self, input_dim=1, variance=1., lengthscale=1., period=2.*np.pi, n_freq=10, lower=0., upper=4*np.pi, name='periodic_exponential'):
+        super(PeriodicExponential, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, name)
+
+    def parameters_changed(self):
+        self.a = [1./self.lengthscale, 1.]
+        self.b = [1]
+
+        self.basis_alpha = np.ones((self.n_basis,))
+        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))[:,0]
+        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
+
+        self.G = self.Gram_matrix()
+        self.Gi = np.linalg.inv(self.G)
+
+    def Gram_matrix(self):
+        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega))
+        Lo = np.column_stack((self.basis_omega,self.basis_omega))
+        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2))
+        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
+        Gint = self._int_computation( r,omega,phi, r,omega,phi)
+        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
+        return(self.lengthscale/(2*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T))
+
+    #@silence_errors
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        """derivative of the covariance matrix with respect to the parameters (shape is N x num_inducing x num_params)"""
+        if X2 is None: X2 = X
+        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
+        FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
+
+        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega))
+        Lo = np.column_stack((self.basis_omega,self.basis_omega))
+        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2))
+        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
+        Gint = self._int_computation( r,omega,phi, r,omega,phi)
+
+        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
+
+        #dK_dvar
+        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T)
+
+        #dK_dlen
+        da_dlen = [-1./self.lengthscale**2,0.]
+        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega))
+        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
+        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
+        dGint_dlen = dGint_dlen + dGint_dlen.T
+        dG_dlen = 1./2*Gint + self.lengthscale/2*dGint_dlen
+        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T)
+
+        #dK_dper
+        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
+        dFX2_dper = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X2,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X2)
+
+        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period))
+        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
+        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
+
+        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
+        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
+        # SIMPLIFY!!!       IPPprim1 = (self.upper - self.lower)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
+        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
+        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
+        IPPprim = np.where(np.logical_or(np.isnan(IPPprim1), np.isinf(IPPprim1)), IPPprim2, IPPprim1)
+
+
+        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
+        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
+        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
+        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
+        #IPPint2[0,0] = (self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
+        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
+
+        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period))
+        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2))
+        r2,omega2,phi2 = dLa_dper2.T,Lo[:,0:1],dLp_dper2.T
+
+        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi)
+        dGint_dper = dGint_dper + dGint_dper.T
+
+        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
+
+        dG_dper = 1./self.variance*(self.lengthscale/2*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)))
+
+        dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T)
+
+        self.variance.gradient = np.sum(dK_dvar*dL_dK)
+        self.lengthscale.gradient = np.sum(dK_dlen*dL_dK)
+        self.period.gradient = np.sum(dK_dper*dL_dK)
+
+
+
+class PeriodicMatern32(Periodic):
+    """
+    Kernel of the periodic subspace (up to a given frequency) of a Matern 3/2 RKHS. Only defined for input_dim=1.
+
+    :param input_dim: the number of input dimensions
+    :type input_dim: int
+    :param variance: the variance of the Matern kernel
+    :type variance: float
+    :param lengthscale: the lengthscale of the Matern kernel
+    :type lengthscale: np.ndarray of size (input_dim,)
+    :param period: the period
+    :type period: float
+    :param n_freq: the number of frequencies considered for the periodic subspace
+    :type n_freq: int
+    :rtype: kernel object
+
+    """
+
+    def __init__(self, input_dim=1, variance=1., lengthscale=1., period=2.*np.pi, n_freq=10, lower=0., upper=4*np.pi, name='periodic_Matern32'):
+        super(PeriodicMatern32, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, name)
+    def parameters_changed(self):
+        self.a = [3./self.lengthscale**2, 2*np.sqrt(3)/self.lengthscale, 1.]
+        self.b = [1,self.lengthscale**2/3]
+
+        self.basis_alpha = np.ones((self.n_basis,))
+        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))
+        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
+
+        self.G = self.Gram_matrix()
+        self.Gi = np.linalg.inv(self.G)
+
+    def Gram_matrix(self):
+        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega,self.a[2]*self.basis_omega**2))
+        Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega))
+        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi))
+        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
+        Gint = self._int_computation( r,omega,phi, r,omega,phi)
+
+        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
+        F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
+        return(self.lengthscale**3/(12*np.sqrt(3)*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscale**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
+
+
+    @silence_errors
+    def update_gradients_full(self,dL_dK,X,X2,target):
+        """derivative of the covariance matrix with respect to the parameters (shape is num_data x num_inducing x num_params)"""
+        if X2 is None: X2 = X
+        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
+        FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
+
+        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega,self.a[2]*self.basis_omega**2))
+        Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega))
+        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi))
+        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
+        Gint = self._int_computation( r,omega,phi, r,omega,phi)
+
+        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
+        F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
+
+        #dK_dvar
+        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T)
+
+        #dK_dlen
+        da_dlen = [-6/self.lengthscale**3,-2*np.sqrt(3)/self.lengthscale**2,0.]
+        db_dlen = [0.,2*self.lengthscale/3.]
+        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega,da_dlen[2]*self.basis_omega**2))
+        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
+        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
+        dGint_dlen = dGint_dlen + dGint_dlen.T
+        dG_dlen = self.lengthscale**2/(4*np.sqrt(3))*Gint + self.lengthscale**3/(12*np.sqrt(3))*dGint_dlen + db_dlen[0]*np.dot(Flower,Flower.T) + db_dlen[1]*np.dot(F1lower,F1lower.T)
+        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T)
+
+        #dK_dper
+        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
+        dFX2_dper = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X2,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X2)
+
+        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period, -self.a[2]*self.basis_omega**3/self.period))
+        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi,self.basis_phi+np.pi*3/2))
+        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
+
+        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
+        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
+        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
+        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
+        IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
+
+        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
+        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
+        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
+        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
+        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
+
+        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period, -2*self.a[2]*self.basis_omega**2/self.period))
+        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
+        r2,omega2,phi2 =  self._cos_factorization(dLa_dper2,Lo[:,0:2],dLp_dper2)
+
+        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) +  self._int_computation(r2,omega2,phi2, r,omega,phi)
+        dGint_dper = dGint_dper + dGint_dper.T
+
+        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
+        dF1lower_dper = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega**2/self.period,self.basis_omega,self.basis_phi+np.pi)(self.lower)+self._cos(-self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
+
+        dG_dper = 1./self.variance*(self.lengthscale**3/(12*np.sqrt(3))*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)) + self.b[1]*(np.dot(dF1lower_dper,F1lower.T)+np.dot(F1lower,dF1lower_dper.T)))
+
+        dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T)
+
+        self.variance.gradient = np.sum(dK_dvar*dL_dK)
+        self.lengthscale.gradient = np.sum(dK_dlen*dL_dK)
+        self.period.gradient = np.sum(dK_dper*dL_dK)
+
+
+
+class PeriodicMatern52(Periodic):
     """
     Kernel of the periodic subspace (up to a given frequency) of a Matern 5/2 RKHS. Only defined for input_dim=1.
 
@@ -25,53 +300,10 @@ class PeriodicMatern52(Kernpart):
 
     """
 
-    def __init__(self,input_dim=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi):
+    def __init__(self, input_dim=1, variance=1., lengthscale=1., period=2.*np.pi, n_freq=10, lower=0., upper=4*np.pi, name='periodic_Matern52'):
-        assert input_dim==1, "Periodic kernels are only defined for input_dim=1"
+        super(PeriodicMatern52, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, name)
-        self.name = 'periodic_Mat52'
-        self.input_dim = input_dim
-        if lengthscale is not None:
-            lengthscale = np.asarray(lengthscale)
-            assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
-        else:
-            lengthscale = np.ones(1)
-        self.lower,self.upper = lower, upper
-        self.num_params = 3
-        self.n_freq = n_freq
-        self.n_basis = 2*n_freq
-        self._set_params(np.hstack((variance,lengthscale,period)))
-
-    def _cos(self,alpha,omega,phase):
-        def f(x):
-            return alpha*np.cos(omega*x+phase)
-        return f
-
-    @silence_errors
-    def _cos_factorization(self,alpha,omega,phase):
-        r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
-        r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
-        r =  np.sqrt(r1**2 + r2**2)
-        psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
-        return r,omega[:,0:1], psi
-
-    @silence_errors
-    def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
-        Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
-        Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) +  np.cos(phi1-phi2.T)*(self.upper-self.lower)
-        #Gint2[0,0] = 2.*(self.upper-self.lower)*np.cos(phi1[0,0])*np.cos(phi2[0,0])
-        Gint = np.dot(r1,r2.T)/2 * np.where(np.isnan(Gint1),Gint2,Gint1)
-        return Gint
-
-    def _get_params(self):
-        """return the value of the parameters."""
-        return np.hstack((self.variance,self.lengthscale,self.period))
-
-    def _set_params(self,x):
-        """set the value of the parameters."""
-        assert x.size==3
-        self.variance = x[0]
-        self.lengthscale = x[1]
-        self.period = x[2]
 
+    def parameters_changed(self):
         self.a = [5*np.sqrt(5)/self.lengthscale**3, 15./self.lengthscale**2,3*np.sqrt(5)/self.lengthscale, 1.]
         self.b  = [9./8, 9*self.lengthscale**4/200., 3*self.lengthscale**2/5., 3*self.lengthscale**2/(5*8.), 3*self.lengthscale**2/(5*8.)]
 
@@ -82,10 +314,6 @@ class PeriodicMatern52(Kernpart):
         self.G = self.Gram_matrix()
         self.Gi = np.linalg.inv(self.G)
 
-    def _get_param_names(self):
-        """return parameter names."""
-        return ['variance','lengthscale','period']
-
     def Gram_matrix(self):
         La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)), self.a[1]*self.basis_omega, self.a[2]*self.basis_omega**2, self.a[3]*self.basis_omega**3))
         Lo = np.column_stack((self.basis_omega, self.basis_omega, self.basis_omega, self.basis_omega))
@@ -99,23 +327,8 @@ class PeriodicMatern52(Kernpart):
         lower_terms = self.b[0]*np.dot(Flower,Flower.T) + self.b[1]*np.dot(F2lower,F2lower.T) + self.b[2]*np.dot(F1lower,F1lower.T) + self.b[3]*np.dot(F2lower,Flower.T) + self.b[4]*np.dot(Flower,F2lower.T)
         return(3*self.lengthscale**5/(400*np.sqrt(5)*self.variance) * Gint + 1./self.variance*lower_terms)
 
-    def K(self,X,X2,target):
-        """Compute the covariance matrix between X and X2."""
-        FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        if X2 is None:
-            FX2 = FX
-        else:
-            FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
-        np.add(mdot(FX,self.Gi,FX2.T), target,target)
-
-    def Kdiag(self,X,target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target)
-
     @silence_errors
-    def _param_grad_helper(self,dL_dK,X,X2,target):
+    def update_gradients_full(self, dL_dK, X, X2=None):
-        """derivative of the covariance matrix with respect to the parameters (shape is num_data x num_inducing x num_params)"""
         if X2 is None: X2 = X
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
         FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
@@ -156,14 +369,12 @@ class PeriodicMatern52(Kernpart):
         IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
         IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
         IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
-        #IPPprim2[0,0] = 2*(self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
         IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
 
         IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
         IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
         IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
         IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
-        #IPPint2[0,0] = (self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
         IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
 
         dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period, -2*self.a[2]*self.basis_omega**2/self.period, -3*self.a[3]*self.basis_omega**3/self.period))
@@ -186,81 +397,7 @@ class PeriodicMatern52(Kernpart):
         dG_dper = 1./self.variance*(3*self.lengthscale**5/(400*np.sqrt(5))*dGint_dper + 0.5*dlower_terms_dper)
         dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T)
 
-        # np.add(target[:,:,0],dK_dvar, target[:,:,0])
+        self.variance.gradient = np.sum(dK_dvar*dL_dK)
-        target[0] += np.sum(dK_dvar*dL_dK)
+        self.lengthscale.gradient = np.sum(dK_dlen*dL_dK)
-        #np.add(target[:,:,1],dK_dlen, target[:,:,1])
+        self.period.gradient = np.sum(dK_dper*dL_dK)
-        target[1] += np.sum(dK_dlen*dL_dK)
-        #np.add(target[:,:,2],dK_dper, target[:,:,2])
-        target[2] += np.sum(dK_dper*dL_dK)
 
-    @silence_errors
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters"""
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)), self.a[1]*self.basis_omega, self.a[2]*self.basis_omega**2, self.a[3]*self.basis_omega**3))
-        Lo = np.column_stack((self.basis_omega, self.basis_omega, self.basis_omega, self.basis_omega))
-        Lp = np.column_stack((self.basis_phi, self.basis_phi+np.pi/2, self.basis_phi+np.pi, self.basis_phi+np.pi*3/2))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-        F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-        F2lower = np.array(self._cos(self.basis_alpha*self.basis_omega**2,self.basis_omega,self.basis_phi+np.pi)(self.lower))[:,None]
-
-        #dK_dvar
-        dK_dvar = 1. / self.variance * mdot(FX, self.Gi, FX.T)
-
-        #dK_dlen
-        da_dlen = [-3*self.a[0]/self.lengthscale, -2*self.a[1]/self.lengthscale, -self.a[2]/self.lengthscale, 0.]
-        db_dlen = [0., 4*self.b[1]/self.lengthscale, 2*self.b[2]/self.lengthscale, 2*self.b[3]/self.lengthscale, 2*self.b[4]/self.lengthscale]
-        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)), da_dlen[1]*self.basis_omega, da_dlen[2]*self.basis_omega**2, da_dlen[3]*self.basis_omega**3))
-        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
-        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
-        dGint_dlen = dGint_dlen + dGint_dlen.T
-        dlower_terms_dlen = db_dlen[0]*np.dot(Flower,Flower.T) + db_dlen[1]*np.dot(F2lower,F2lower.T) + db_dlen[2]*np.dot(F1lower,F1lower.T) + db_dlen[3]*np.dot(F2lower,Flower.T) + db_dlen[4]*np.dot(Flower,F2lower.T)
-        dG_dlen = 15*self.lengthscale**4/(400*np.sqrt(5))*Gint + 3*self.lengthscale**5/(400*np.sqrt(5))*dGint_dlen + dlower_terms_dlen
-        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX.T)
-
-        #dK_dper
-        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
-
-        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period, -self.a[2]*self.basis_omega**3/self.period, -self.a[3]*self.basis_omega**4/self.period))
-        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi,self.basis_phi+np.pi*3/2,self.basis_phi))
-        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
-
-        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
-        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
-        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
-        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
-        IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
-
-        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
-        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
-        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + .5*self.upper**2*np.cos(phi-phi1.T)
-        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + .5*self.lower**2*np.cos(phi-phi1.T)
-        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
-
-        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period, -2*self.a[2]*self.basis_omega**2/self.period, -3*self.a[3]*self.basis_omega**3/self.period))
-        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2, self.basis_phi+np.pi, self.basis_phi+np.pi*3/2))
-        r2,omega2,phi2 =  self._cos_factorization(dLa_dper2,Lo[:,0:2],dLp_dper2)
-
-        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) +  self._int_computation(r2,omega2,phi2, r,omega,phi)
-        dGint_dper = dGint_dper + dGint_dper.T
-
-        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-        dF1lower_dper = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega**2/self.period,self.basis_omega,self.basis_phi+np.pi)(self.lower)+self._cos(-self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-        dF2lower_dper = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega**3/self.period,self.basis_omega,self.basis_phi+np.pi*3/2)(self.lower) + self._cos(-2*self.basis_alpha*self.basis_omega**2/self.period,self.basis_omega,self.basis_phi+np.pi)(self.lower))[:,None]
-
-        dlower_terms_dper  = self.b[0] * (np.dot(dFlower_dper,Flower.T) + np.dot(Flower.T,dFlower_dper))
-        dlower_terms_dper += self.b[1] * (np.dot(dF2lower_dper,F2lower.T) + np.dot(F2lower,dF2lower_dper.T)) - 4*self.b[1]/self.period*np.dot(F2lower,F2lower.T)
-        dlower_terms_dper += self.b[2] * (np.dot(dF1lower_dper,F1lower.T) + np.dot(F1lower,dF1lower_dper.T)) - 2*self.b[2]/self.period*np.dot(F1lower,F1lower.T)
-        dlower_terms_dper += self.b[3] * (np.dot(dF2lower_dper,Flower.T) + np.dot(F2lower,dFlower_dper.T)) - 2*self.b[3]/self.period*np.dot(F2lower,Flower.T)
-        dlower_terms_dper += self.b[4] * (np.dot(dFlower_dper,F2lower.T) + np.dot(Flower,dF2lower_dper.T)) - 2*self.b[4]/self.period*np.dot(Flower,F2lower.T)
-
-        dG_dper = 1./self.variance*(3*self.lengthscale**5/(400*np.sqrt(5))*dGint_dper + 0.5*dlower_terms_dper)
-        dK_dper = 2*mdot(dFX_dper,self.Gi,FX.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX.T)
-
-        target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag)
-        target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag)
-        target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag)

GPy/kern/_src/periodic_Matern32.py

@@ -1,248 +0,0 @@
-# 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
-from GPy.util.linalg import mdot
-from GPy.util.decorators import silence_errors
-
-class PeriodicMatern32(Kernpart):
-    """
-    Kernel of the periodic subspace (up to a given frequency) of a Matern 3/2 RKHS. Only defined for input_dim=1.
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int
-    :param variance: the variance of the Matern kernel
-    :type variance: float
-    :param lengthscale: the lengthscale of the Matern kernel
-    :type lengthscale: np.ndarray of size (input_dim,)
-    :param period: the period
-    :type period: float
-    :param n_freq: the number of frequencies considered for the periodic subspace
-    :type n_freq: int
-    :rtype: kernel object
-
-    """
-
-    def __init__(self, input_dim=1, variance=1., lengthscale=None, period=2 * np.pi, n_freq=10, lower=0., upper=4 * np.pi):
-        assert input_dim==1, "Periodic kernels are only defined for input_dim=1"
-        self.name = 'periodic_Mat32'
-        self.input_dim = input_dim
-        if lengthscale is not None:
-            lengthscale = np.asarray(lengthscale)
-            assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
-        else:
-            lengthscale = np.ones(1)
-        self.lower,self.upper = lower, upper
-        self.num_params = 3
-        self.n_freq = n_freq
-        self.n_basis = 2*n_freq
-        self._set_params(np.hstack((variance,lengthscale,period)))
-
-    def _cos(self,alpha,omega,phase):
-        def f(x):
-            return alpha*np.cos(omega*x+phase)
-        return f
-
-    @silence_errors
-    def _cos_factorization(self,alpha,omega,phase):
-        r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
-        r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
-        r =  np.sqrt(r1**2 + r2**2)
-        psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
-        return r,omega[:,0:1], psi
-
-    @silence_errors
-    def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
-        Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
-        Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) +  np.cos(phi1-phi2.T)*(self.upper-self.lower)
-        #Gint2[0,0] = 2.*(self.upper-self.lower)*np.cos(phi1[0,0])*np.cos(phi2[0,0])
-        Gint = np.dot(r1,r2.T)/2 * np.where(np.isnan(Gint1),Gint2,Gint1)
-        return Gint
-
-    def _get_params(self):
-        """return the value of the parameters."""
-        return np.hstack((self.variance,self.lengthscale,self.period))
-
-    def _set_params(self,x):
-        """set the value of the parameters."""
-        assert x.size==3
-        self.variance = x[0]
-        self.lengthscale = x[1]
-        self.period = x[2]
-
-        self.a = [3./self.lengthscale**2, 2*np.sqrt(3)/self.lengthscale, 1.]
-        self.b = [1,self.lengthscale**2/3]
-
-        self.basis_alpha = np.ones((self.n_basis,))
-        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))
-        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
-
-        self.G = self.Gram_matrix()
-        self.Gi = np.linalg.inv(self.G)
-
-    def _get_param_names(self):
-        """return parameter names."""
-        return ['variance','lengthscale','period']
-
-    def Gram_matrix(self):
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega,self.a[2]*self.basis_omega**2))
-        Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega))
-        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-        F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-        return(self.lengthscale**3/(12*np.sqrt(3)*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscale**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
-
-    def K(self,X,X2,target):
-        """Compute the covariance matrix between X and X2."""
-        FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        if X2 is None:
-            FX2 = FX
-        else:
-            FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
-        np.add(mdot(FX,self.Gi,FX2.T), target,target)
-
-    def Kdiag(self,X,target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target)
-
-    @silence_errors
-    def _param_grad_helper(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters (shape is num_data x num_inducing x num_params)"""
-        if X2 is None: X2 = X
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
-
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega,self.a[2]*self.basis_omega**2))
-        Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega))
-        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-        F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-
-        #dK_dvar
-        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T)
-
-        #dK_dlen
-        da_dlen = [-6/self.lengthscale**3,-2*np.sqrt(3)/self.lengthscale**2,0.]
-        db_dlen = [0.,2*self.lengthscale/3.]
-        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega,da_dlen[2]*self.basis_omega**2))
-        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
-        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
-        dGint_dlen = dGint_dlen + dGint_dlen.T
-        dG_dlen = self.lengthscale**2/(4*np.sqrt(3))*Gint + self.lengthscale**3/(12*np.sqrt(3))*dGint_dlen + db_dlen[0]*np.dot(Flower,Flower.T) + db_dlen[1]*np.dot(F1lower,F1lower.T)
-        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T)
-
-        #dK_dper
-        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
-        dFX2_dper = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X2,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X2)
-
-        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period, -self.a[2]*self.basis_omega**3/self.period))
-        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi,self.basis_phi+np.pi*3/2))
-        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
-
-        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
-        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
-        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
-        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
-        #IPPprim2[0,0] = 2*(self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
-        IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
-
-        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
-        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
-        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
-        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
-        #IPPint2[0,0] = (self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
-        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
-
-        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period, -2*self.a[2]*self.basis_omega**2/self.period))
-        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r2,omega2,phi2 =  self._cos_factorization(dLa_dper2,Lo[:,0:2],dLp_dper2)
-
-        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) +  self._int_computation(r2,omega2,phi2, r,omega,phi)
-        dGint_dper = dGint_dper + dGint_dper.T
-
-        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-        dF1lower_dper = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega**2/self.period,self.basis_omega,self.basis_phi+np.pi)(self.lower)+self._cos(-self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-
-        dG_dper = 1./self.variance*(self.lengthscale**3/(12*np.sqrt(3))*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)) + self.b[1]*(np.dot(dF1lower_dper,F1lower.T)+np.dot(F1lower,dF1lower_dper.T)))
-
-        dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T)
-
-        # np.add(target[:,:,0],dK_dvar, target[:,:,0])
-        target[0] += np.sum(dK_dvar*dL_dK)
-        #np.add(target[:,:,1],dK_dlen, target[:,:,1])
-        target[1] += np.sum(dK_dlen*dL_dK)
-        #np.add(target[:,:,2],dK_dper, target[:,:,2])
-        target[2] += np.sum(dK_dper*dL_dK)
-
-    @silence_errors
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        """derivative of the diagonal covariance matrix with respect to the parameters"""
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega, self.a[2]*self.basis_omega**2))
-        Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega))
-        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-        F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-
-        #dK_dvar
-        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX.T)
-
-        #dK_dlen
-        da_dlen = [-6/self.lengthscale**3,-2*np.sqrt(3)/self.lengthscale**2,0.]
-        db_dlen = [0.,2*self.lengthscale/3.]
-        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega,da_dlen[2]*self.basis_omega**2))
-        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
-        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
-        dGint_dlen = dGint_dlen + dGint_dlen.T
-        dG_dlen = self.lengthscale**2/(4*np.sqrt(3))*Gint + self.lengthscale**3/(12*np.sqrt(3))*dGint_dlen + db_dlen[0]*np.dot(Flower,Flower.T) + db_dlen[1]*np.dot(F1lower,F1lower.T)
-        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX.T)
-
-        #dK_dper
-        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
-
-        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period, -self.a[2]*self.basis_omega**3/self.period))
-        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi,self.basis_phi+np.pi*3/2))
-        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
-
-        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
-        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
-        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
-        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
-        IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
-
-        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
-        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
-        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
-        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
-        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
-
-        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period, -2*self.a[2]*self.basis_omega**2/self.period))
-        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r2,omega2,phi2 =  self._cos_factorization(dLa_dper2,Lo[:,0:2],dLp_dper2)
-
-        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) +  self._int_computation(r2,omega2,phi2, r,omega,phi)
-        dGint_dper = dGint_dper + dGint_dper.T
-
-        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-        dF1lower_dper = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega**2/self.period,self.basis_omega,self.basis_phi+np.pi)(self.lower)+self._cos(-self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-
-        dG_dper = 1./self.variance*(self.lengthscale**3/(12*np.sqrt(3))*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)) + self.b[1]*(np.dot(dF1lower_dper,F1lower.T)+np.dot(F1lower,dF1lower_dper.T)))
-
-        dK_dper = 2* mdot(dFX_dper,self.Gi,FX.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX.T)
-
-        target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag)
-        target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag)
-        target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag)

GPy/kern/_src/periodic_exponential.py

@@ -1,237 +0,0 @@
-# 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
-from GPy.util.linalg import mdot
-from GPy.util.decorators import silence_errors
-from GPy.core.parameterization.param import Param
-
-class PeriodicExponential(Kernpart):
-    """
-    Kernel of the periodic subspace (up to a given frequency) of a exponential (Matern 1/2) RKHS. Only defined for input_dim=1.
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int
-    :param variance: the variance of the Matern kernel
-    :type variance: float
-    :param lengthscale: the lengthscale of the Matern kernel
-    :type lengthscale: np.ndarray of size (input_dim,)
-    :param period: the period
-    :type period: float
-    :param n_freq: the number of frequencies considered for the periodic subspace
-    :type n_freq: int
-    :rtype: kernel object
-
-    """
-
-    def __init__(self, input_dim=1, variance=1., lengthscale=None, period=2 * np.pi, n_freq=10, lower=0., upper=4 * np.pi, name='periodic_exp'):
-        super(PeriodicExponential, self).__init__(input_dim, name)
-        assert input_dim==1, "Periodic kernels are only defined for input_dim=1"
-        self.input_dim = input_dim
-        if lengthscale is not None:
-            lengthscale = np.asarray(lengthscale)
-            assert lengthscale.size == 1, "Wrong size: only one lengthscale needed"
-        else:
-            lengthscale = np.ones(1)
-        self.lower,self.upper = lower, upper
-        self.num_params = 3
-        self.n_freq = n_freq
-        self.n_basis = 2*n_freq
-        self.variance = Param('variance', variance)
-        self.lengthscale = Param('lengthscale', lengthscale)
-        self.period = Param('period', period)
-        self.parameters_changed()
-        #self._set_params(np.hstack((variance,lengthscale,period)))
-
-    def _cos(self,alpha,omega,phase):
-        def f(x):
-            return alpha*np.cos(omega*x+phase)
-        return f
-
-    @silence_errors
-    def _cos_factorization(self,alpha,omega,phase):
-        r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None]
-        r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None]
-        r =  np.sqrt(r1**2 + r2**2)
-        psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2))
-        return r,omega[:,0:1], psi
-
-    @silence_errors
-    def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2):
-        Gint1 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) + 1./(omega1-omega2.T)*( np.sin((omega1-omega2.T)*self.upper+phi1-phi2.T) - np.sin((omega1-omega2.T)*self.lower+phi1-phi2.T) )
-        Gint2 = 1./(omega1+omega2.T)*( np.sin((omega1+omega2.T)*self.upper+phi1+phi2.T) - np.sin((omega1+omega2.T)*self.lower+phi1+phi2.T)) +  np.cos(phi1-phi2.T)*(self.upper-self.lower)
-        #Gint2[0,0] = 2.*(self.upper-self.lower)*np.cos(phi1[0,0])*np.cos(phi2[0,0])
-        Gint = np.dot(r1,r2.T)/2 * np.where(np.isnan(Gint1),Gint2,Gint1)
-        return Gint
-
-    #def _get_params(self):
-    #    """return the value of the parameters."""
-    #    return np.hstack((self.variance,self.lengthscale,self.period))
-
-    def parameters_changed(self):
-        """set the value of the parameters."""
-        self.a = [1./self.lengthscale, 1.]
-        self.b = [1]
-
-        self.basis_alpha = np.ones((self.n_basis,))
-        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))[:,0]
-        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
-
-        self.G = self.Gram_matrix()
-        self.Gi = np.linalg.inv(self.G)
-
-    #def _get_param_names(self):
-    #    """return parameter names."""
-    #    return ['variance','lengthscale','period']
-
-    def Gram_matrix(self):
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega))
-        Lo = np.column_stack((self.basis_omega,self.basis_omega))
-        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-        return(self.lengthscale/(2*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T))
-
-    def K(self,X,X2,target):
-        """Compute the covariance matrix between X and X2."""
-        FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        if X2 is None:
-            FX2 = FX
-        else:
-            FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
-        np.add(mdot(FX,self.Gi,FX2.T), target,target)
-
-    def Kdiag(self,X,target):
-        """Compute the diagonal of the covariance matrix associated to X."""
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        np.add(target,np.diag(mdot(FX,self.Gi,FX.T)),target)
-
-    @silence_errors
-    def _param_grad_helper(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters (shape is N x num_inducing x num_params)"""
-        if X2 is None: X2 = X
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-        FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2)
-
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega))
-        Lo = np.column_stack((self.basis_omega,self.basis_omega))
-        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-
-        #dK_dvar
-        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T)
-
-        #dK_dlen
-        da_dlen = [-1./self.lengthscale**2,0.]
-        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega))
-        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
-        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
-        dGint_dlen = dGint_dlen + dGint_dlen.T
-        dG_dlen = 1./2*Gint + self.lengthscale/2*dGint_dlen
-        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T)
-
-        #dK_dper
-        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
-        dFX2_dper = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X2,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X2)
-
-        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period))
-        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
-
-        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
-        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
-        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
-        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
-        #IPPprim2[0,0] = 2*(self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
-        IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
-
-        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
-        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
-        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
-        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
-        #IPPint2[0,0] = (self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0])
-        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
-
-        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period))
-        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2))
-        r2,omega2,phi2 = dLa_dper2.T,Lo[:,0:1],dLp_dper2.T
-
-        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi)
-        dGint_dper = dGint_dper + dGint_dper.T
-
-        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-
-        dG_dper = 1./self.variance*(self.lengthscale/2*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)))
-
-        dK_dper = mdot(dFX_dper,self.Gi,FX2.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX2.T) + mdot(FX,self.Gi,dFX2_dper.T)
-
-        target[0] += np.sum(dK_dvar*dL_dK)
-        target[1] += np.sum(dK_dlen*dL_dK)
-        target[2] += np.sum(dK_dper*dL_dK)
-
-    @silence_errors
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        """derivative of the diagonal of the covariance matrix with respect to the parameters"""
-        FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
-
-        La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega))
-        Lo = np.column_stack((self.basis_omega,self.basis_omega))
-        Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2))
-        r,omega,phi =  self._cos_factorization(La,Lo,Lp)
-        Gint = self._int_computation( r,omega,phi, r,omega,phi)
-
-        Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
-
-        #dK_dvar
-        dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX.T)
-
-        #dK_dlen
-        da_dlen = [-1./self.lengthscale**2,0.]
-        dLa_dlen =  np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega))
-        r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp)
-        dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi)
-        dGint_dlen = dGint_dlen + dGint_dlen.T
-        dG_dlen = 1./2*Gint + self.lengthscale/2*dGint_dlen
-        dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX.T)
-
-        #dK_dper
-        dFX_dper  = self._cos(-self.basis_alpha[None,:]*self.basis_omega[None,:]/self.period*X ,self.basis_omega[None,:],self.basis_phi[None,:]+np.pi/2)(X)
-
-        dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period))
-        dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi))
-        r1,omega1,phi1 =  self._cos_factorization(dLa_dper,Lo,dLp_dper)
-
-        IPPprim1 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi/2))
-        IPPprim1 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  +  1./(omega-omega1.T)*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi/2))
-        IPPprim2 =  self.upper*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi/2)  + self.upper*np.cos(phi-phi1.T))
-        IPPprim2 -= self.lower*(1./(omega+omega1.T)*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi/2)  + self.lower*np.cos(phi-phi1.T))
-        IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1)
-
-        IPPint1 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.upper+phi-phi1.T-np.pi)
-        IPPint1 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  +  1./(omega-omega1.T)**2*np.cos((omega-omega1.T)*self.lower+phi-phi1.T-np.pi)
-        IPPint2 =  1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.upper+phi+phi1.T-np.pi)  + 1./2*self.upper**2*np.cos(phi-phi1.T)
-        IPPint2 -= 1./(omega+omega1.T)**2*np.cos((omega+omega1.T)*self.lower+phi+phi1.T-np.pi)  + 1./2*self.lower**2*np.cos(phi-phi1.T)
-        IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1)
-
-        dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period))
-        dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2))
-        r2,omega2,phi2 = dLa_dper2.T,Lo[:,0:1],dLp_dper2.T
-
-        dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi)
-        dGint_dper = dGint_dper + dGint_dper.T
-
-        dFlower_dper  = np.array(self._cos(-self.lower*self.basis_alpha*self.basis_omega/self.period,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None]
-
-        dG_dper = 1./self.variance*(self.lengthscale/2*dGint_dper + self.b[0]*(np.dot(dFlower_dper,Flower.T)+np.dot(Flower,dFlower_dper.T)))
-
-        dK_dper = 2*mdot(dFX_dper,self.Gi,FX.T) - mdot(FX,self.Gi,dG_dper,self.Gi,FX.T)
-
-        target[0] += np.sum(np.diag(dK_dvar)*dL_dKdiag)
-        target[1] += np.sum(np.diag(dK_dlen)*dL_dKdiag)
-        target[2] += np.sum(np.diag(dK_dper)*dL_dKdiag)

GPy/kern/_src/prod_orthogonal.py

@@ -1,101 +0,0 @@
-# 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 prod_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):
-        self.input_dim = k1.input_dim + k2.input_dim
-        self.num_params = k1.num_params + k2.num_params
-        self.name = k1.name + '<times>' + k2.name
-        self.k1 = k1
-        self.k2 = k2
-        self._X, self._X2, self._params = np.empty(shape=(3,1))
-        self._set_params(np.hstack((k1._get_params(),k2._get_params())))
-
-    def _get_params(self):
-        """return the value of the parameters."""
-        return np.hstack((self.k1._get_params(), self.k2._get_params()))
-
-    def _set_params(self,x):
-        """set the value of the parameters."""
-        self.k1._set_params(x[:self.k1.num_params])
-        self.k2._set_params(x[self.k1.num_params:])
-
-    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):
-        self._K_computations(X,X2)
-        target += self._K1 * self._K2
-
-    def _param_grad_helper(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to the parameters."""
-        self._K_computations(X,X2)
-        if X2 is None:
-            self.k1._param_grad_helper(dL_dK*self._K2, X[:,:self.k1.input_dim], None, target[:self.k1.num_params])
-            self.k2._param_grad_helper(dL_dK*self._K1, X[:,self.k1.input_dim:], None, target[self.k1.num_params:])
-        else:
-            self.k1._param_grad_helper(dL_dK*self._K2, X[:,:self.k1.input_dim], X2[:,:self.k1.input_dim], target[:self.k1.num_params])
-            self.k2._param_grad_helper(dL_dK*self._K1, X[:,self.k1.input_dim:], X2[:,self.k1.input_dim:], target[self.k1.num_params:])
-
-    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[:,:self.k1.input_dim],target1)
-        self.k2.Kdiag(X[:,self.k1.input_dim:],target2)
-        target += target1 * target2
-
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        K1 = np.zeros(X.shape[0])
-        K2 = np.zeros(X.shape[0])
-        self.k1.Kdiag(X[:,:self.k1.input_dim],K1)
-        self.k2.Kdiag(X[:,self.k1.input_dim:],K2)
-        self.k1.dKdiag_dtheta(dL_dKdiag*K2,X[:,:self.k1.input_dim],target[:self.k1.num_params])
-        self.k2.dKdiag_dtheta(dL_dKdiag*K1,X[:,self.k1.input_dim:],target[self.k1.num_params:])
-
-    def gradients_X(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to X."""
-        self._K_computations(X,X2)
-        self.k1.gradients_X(dL_dK*self._K2, X[:,:self.k1.input_dim], X2[:,:self.k1.input_dim], target)
-        self.k2.gradients_X(dL_dK*self._K1, X[:,self.k1.input_dim:], X2[:,self.k1.input_dim:], target)
-
-    def dKdiag_dX(self, dL_dKdiag, X, target):
-        K1 = np.zeros(X.shape[0])
-        K2 = np.zeros(X.shape[0])
-        self.k1.Kdiag(X[:,0:self.k1.input_dim],K1)
-        self.k2.Kdiag(X[:,self.k1.input_dim:],K2)
-
-        self.k1.gradients_X(dL_dKdiag*K2, X[:,:self.k1.input_dim], target)
-        self.k2.gradients_X(dL_dKdiag*K1, X[:,self.k1.input_dim:], target)
-
-    def _K_computations(self,X,X2):
-        if not (np.array_equal(X,self._X) and np.array_equal(X2,self._X2) and np.array_equal(self._params , self._get_params())):
-            self._X = X.copy()
-            self._params == self._get_params().copy()
-            if X2 is None:
-                self._X2 = None
-                self._K1 = np.zeros((X.shape[0],X.shape[0]))
-                self._K2 = np.zeros((X.shape[0],X.shape[0]))
-                self.k1.K(X[:,:self.k1.input_dim],None,self._K1)
-                self.k2.K(X[:,self.k1.input_dim:],None,self._K2)
-            else:
-                self._X2 = X2.copy()
-                self._K1 = np.zeros((X.shape[0],X2.shape[0]))
-                self._K2 = np.zeros((X.shape[0],X2.shape[0]))
-                self.k1.K(X[:,:self.k1.input_dim],X2[:,:self.k1.input_dim],self._K1)
-                self.k2.K(X[:,self.k1.input_dim:],X2[:,self.k1.input_dim:],self._K2)
-

GPy/kern/_src/rational_quadratic.py

@@ -1,82 +0,0 @@
-# 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
-
-class RationalQuadratic(Kernpart):
-    """
-    rational quadratic kernel
-
-    .. math::
-
-       k(r) = \sigma^2 \\bigg( 1 + \\frac{r^2}{2 \ell^2} \\bigg)^{- \\alpha} \ \ \ \ \  \\text{ where  } r^2 = (x-y)^2
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int (input_dim=1 is the only value currently supported)
-    :param variance: the variance :math:`\sigma^2`
-    :type variance: float
-    :param lengthscale: the lengthscale :math:`\ell`
-    :type lengthscale: float
-    :param power: the power :math:`\\alpha`
-    :type power: float
-    :rtype: Kernpart object
-
-    """
-    def __init__(self,input_dim,variance=1.,lengthscale=1.,power=1.):
-        assert input_dim == 1, "For this kernel we assume input_dim=1"
-        self.input_dim = input_dim
-        self.num_params = 3
-        self.name = 'rat_quad'
-        self.variance = variance
-        self.lengthscale = lengthscale
-        self.power = power
-
-    def _get_params(self):
-        return np.hstack((self.variance,self.lengthscale,self.power))
-
-    def _set_params(self,x):
-        self.variance = x[0]
-        self.lengthscale = x[1]
-        self.power = x[2]
-
-    def _get_param_names(self):
-        return ['variance','lengthscale','power']
-
-    def K(self,X,X2,target):
-        if X2 is None: X2 = X
-        dist2 = np.square((X-X2.T)/self.lengthscale)
-        target += self.variance*(1 + dist2/2.)**(-self.power)
-
-    def Kdiag(self,X,target):
-        target += self.variance
-
-    def _param_grad_helper(self,dL_dK,X,X2,target):
-        if X2 is None: X2 = X
-        dist2 = np.square((X-X2.T)/self.lengthscale)
-
-        dvar = (1 + dist2/2.)**(-self.power)
-        dl = self.power * self.variance * dist2 / self.lengthscale * (1 + dist2/2.)**(-self.power-1)
-        dp = - self.variance * np.log(1 + dist2/2.) * (1 + dist2/2.)**(-self.power)
-
-        target[0] += np.sum(dvar*dL_dK)
-        target[1] += np.sum(dl*dL_dK)
-        target[2] += np.sum(dp*dL_dK)
-
-    def dKdiag_dtheta(self,dL_dKdiag,X,target):
-        target[0] += np.sum(dL_dKdiag)
-        # here self.lengthscale and self.power have no influence on Kdiag so target[1:] are unchanged
-
-    def gradients_X(self,dL_dK,X,X2,target):
-        """derivative of the covariance matrix with respect to X."""
-        if X2 is None:
-            dist2 = np.square((X-X.T)/self.lengthscale)
-            dX = -2.*self.variance*self.power * (X-X.T)/self.lengthscale**2 *  (1 + dist2/2./self.lengthscale)**(-self.power-1)
-        else:
-            dist2 = np.square((X-X2.T)/self.lengthscale)
-            dX = -self.variance*self.power * (X-X2.T)/self.lengthscale**2 *  (1 + dist2/2./self.lengthscale)**(-self.power-1)
-        target += np.sum(dL_dK*dX,1)[:,np.newaxis]
-
-    def dKdiag_dX(self,dL_dKdiag,X,target):
-        pass

GPy/kern/_src/rbf.py

@@ -128,7 +128,7 @@ class RBF(Kern):
 
     def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, posterior_variational):
         mu = posterior_variational.mean
-        S = posterior_variational.variance        
+        S = posterior_variational.variance
         self._psi_computations(Z, mu, S)
 
         #contributions from psi0:
@@ -385,4 +385,4 @@ class RBF(Kern):
 
     def input_sensitivity(self):
         if self.ARD: return 1./self.lengthscale
-        else: return (1./self.lengthscale).repeat(self.input_dim)
+        else: return (1./self.lengthscale).repeat(self.input_dim)

GPy/kern/_src/static.py

@@ -7,11 +7,66 @@ from ...core.parameterization import Param
 from ...core.parameterization.transformations import Logexp
 import numpy as np
 
-class Bias(Kern):
+class Static(Kern):
-    def __init__(self,input_dim,variance=1.,name=None):
+    def __init__(self, input_dim, variance, name):
+        super(Static, self).__init__(input_dim, name)
+        self.variance = Param('variance', variance, Logexp())
+        self.add_parameters(self.variance)
+
+    def Kdiag(self, X):
+        ret = np.empty((X.shape[0],), dtype=np.float64)
+        ret[:] = self.variance
+        return ret
+
+    def gradients_X(self, dL_dK, X, X2, target):
+        return np.zeros(X.shape)
+
+    def gradients_X_diag(self, dL_dKdiag, X, target):
+        return np.zeros(X.shape)
+
+    def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
+        return np.zeros(Z.shape)
+
+    def gradients_muS_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
+        return np.zeros(mu.shape), np.zeros(S.shape)
+
+    def psi0(self, Z, mu, S):
+        return self.Kdiag(mu)
+
+    def psi1(self, Z, mu, S, target):
+        return self.K(mu, Z)
+
+    def psi2(Z, mu, S):
+        K = self.K(mu, Z)
+        return K[:,:,None]*K[:,None,:] # NB. more efficient implementations on inherriting classes
+
+
+class White(Static):
+    def __init__(self, input_dim, variance=1., name='white'):
+        super(White, self).__init__(input_dim, name)
+
+    def K(self, X, X2=None):
+        if X2 is None:
+            return np.eye(X.shape[0])*self.variance
+        else:
+            return np.zeros((X.shape[0], X2.shape[0]))
+
+    def psi2(self, Z, mu, S, target):
+        return np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]), dtype=np.float64)
+
+    def update_gradients_full(self, dL_dK, X):
+        self.variance.gradient = np.trace(dL_dK)
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        self.variance.gradient = dL_dKdiag.sum()
+
+    def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
+        self.variance.gradient = np.trace(dL_dKmm) + dL_dpsi0.sum()
+
+
+class Bias(Static):
+    def __init__(self, input_dim, variance=1., name='bias'):
         super(Bias, self).__init__(input_dim, name)
-        self.variance = Param("variance", variance, Logexp())
-        self.add_parameter(self.variance)
 
     def K(self, X, X2=None):
         shape = (X.shape[0], X.shape[0] if X2 is None else X2.shape[0])
@@ -19,33 +74,11 @@ class Bias(Kern):
         ret[:] = self.variance
         return ret
 
-    def Kdiag(self,X):
-        ret = np.empty((X.shape[0],), dtype=np.float64)
-        ret[:] = self.variance
-        return ret
-
     def update_gradients_full(self, dL_dK, X, X2=None):
         self.variance.gradient = dL_dK.sum()
 
     def update_gradients_diag(self, dL_dKdiag, X):
-        self.variance.gradient = dL_dKdiag.sum()
+        self.variance.gradient = dL_dK.sum()
-
-    def gradients_X(self, dL_dK,X, X2):
-        return np.zeros(X.shape)
-
-    def gradients_X_diag(self,dL_dKdiag,X):
-        return np.zeros(X.shape)
-
-
-    #---------------------------------------#
-    #             PSI statistics            #
-    #---------------------------------------#
-
-    def psi0(self, Z, mu, S):
-        return self.Kdiag(mu)
-
-    def psi1(self, Z, mu, S, target):
-        return self.K(mu, S)
 
     def psi2(self, Z, mu, S, target):
         ret = np.empty((mu.shape[0], Z.shape[0], Z.shape[0]), dtype=np.float64)
@@ -55,8 +88,3 @@ class Bias(Kern):
     def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
         self.variance.gradient = dL_dKmm.sum() + dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum()
 
-    def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
-        return np.zeros(Z.shape)
-
-    def gradients_muS_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
-        return np.zeros(mu.shape), np.zeros(S.shape)

GPy/kern/_src/stationary.py

@@ -32,6 +32,13 @@ class Stationary(Kern):
         assert self.variance.size==1
         self.add_parameters(self.variance, self.lengthscale)
 
+    def K_of_r(self, r):
+        raise NotImplementedError, "implement the covaraiance functino and a fn of r to use this class"
+
+    def K(self, X, X2=None):
+        r = self._scaled_dist(X, X2)
+        return self.K_of_r(r)
+
     def _dist(self, X, X2):
         if X2 is None:
             X2 = X
@@ -50,11 +57,11 @@ class Stationary(Kern):
         self.lengthscale.gradient = 0.
 
     def update_gradients_full(self, dL_dK, X, X2=None):
-        K = self.K(X, X2)
+        r = self._scaled_dist(X, X2)
-        self.variance.gradient = np.sum(K * dL_dK)/self.variance
+        K = self.K_of_r(r)
 
         rinv = self._inv_dist(X, X2)
-        dL_dr = self.dK_dr(X, X2) * dL_dK
+        dL_dr = self.dK_dr(r) * dL_dK
         x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3
 
         if self.ARD:
@@ -62,6 +69,8 @@ class Stationary(Kern):
         else:
             self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum()
 
+        self.variance.gradient = np.sum(K * dL_dK)/self.variance
+
     def _inv_dist(self, X, X2=None):
         dist = self._scaled_dist(X, X2)
         if X2 is None:
@@ -72,7 +81,8 @@ class Stationary(Kern):
             return 1./np.where(dist != 0., dist, np.inf)
 
     def gradients_X(self, dL_dK, X, X2=None):
-        dL_dr = self.dK_dr(X, X2) * dL_dK
+        r = self._scaled_dist(X, X2)
+        dL_dr = self.dK_dr(r) * dL_dK
         invdist = self._inv_dist(X, X2)
         ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2
         if X2 is None:
@@ -82,19 +92,65 @@ class Stationary(Kern):
     def gradients_X_diag(self, dL_dKdiag, X):
         return np.zeros(X.shape)
 
+    def add(self, other, tensor=False):
+        if not tensor:
+            return StatAdd(self, other)
+        else:
+            return super(Stationary, self).add(other, tensor)
+
+    def prod(self, other, tensor=False):
+        if not tensor:
+            return StatProd(self, other)
+        else:
+            return super(Stationary, self).prod(other, tensor)
 
 
 
+class StatAdd(Stationary):
+    """
+    Addition of two Stationary kernels on the same space is still stationary.
+
+    If you need to add two (stationary) kernels on separate spaces, use the generic add class.
+    """
+    def __init__(self, k1, k2):
+        assert isinstance(k1, Stationary)
+        assert isinstance(k2, Stationary)
+        self.k1, self.k2 = k1, k2
+        self.add_parameters(k1, k2)
+
+    def K_of_r(self, r):
+        return self.k1.K(r) + self.k2.K(r)
+
+    def dK_dr(self, r):
+        return self.k1.dK_dr + self.k2.dK_dr(r)
+
+class StatProd(Stationary):
+    """
+    Product of two Stationary kernels on the same space is still stationary.
+
+    If you need to multiply two (stationary) kernels on separate spaces, use the generic Prod class.
+    """
+    def __init__(self, k1, k2):
+        assert isinstance(k1, Stationary)
+        assert isinstance(k2, Stationary)
+        self.k1, self.k2 = k1, k2
+        self.add_parameters(k1, k2)
+
+    def K_of_r(self, r):
+        return self.k1.K(r) * self.k2.K(r)
+
+    def dK_dr(self, r):
+        return self.k1.dK_dr(r) * self.k2.K_of_r(r) + self.k2.dK_dr(r) * self.k1.K_of_r(r)
+
 class Exponential(Stationary):
     def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Exponential'):
         super(Exponential, self).__init__(input_dim, variance, lengthscale, ARD, name)
 
     def K(self, X, X2=None):
-        dist = self._scaled_dist(X, X2)
         return self.variance * np.exp(-0.5 * dist)
 
-    def dK_dr(self, X, X2):
+    def dK_dr(self, r):
-        return -0.5*self.K(X, X2)
+        return -0.5*self.K_of_r(r)
 
 class Matern32(Stationary):
     """
@@ -109,13 +165,11 @@ class Matern32(Stationary):
     def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat32'):
         super(Matern32, self).__init__(input_dim, variance, lengthscale, ARD, name)
 
-    def K(self, X, X2=None):
+    def K_of_r(self, r):
-        dist = self._scaled_dist(X, X2)
+        return self.variance * (1. + np.sqrt(3.) * r) * np.exp(-np.sqrt(3.) * r)
-        return self.variance * (1. + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist)
 
-    def dK_dr(self, X, X2):
+    def dK_dr(self,r):
-        dist = self._scaled_dist(X, X2)
+        return -3.*self.variance*r*np.exp(-np.sqrt(3.)*r)
-        return -3.*self.variance*dist*np.exp(-np.sqrt(3.)*dist)
 
     def Gram_matrix(self, F, F1, F2, lower, upper):
         """
@@ -153,12 +207,10 @@ class Matern52(Stationary):
        k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \  \\text{ where  } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
        """
 
-    def K(self, X, X2=None):
+    def K_of_r(self, r):
-        r = self._scaled_dist(X, X2)
         return self.variance*(1+np.sqrt(5.)*r+5./3*r**2)*np.exp(-np.sqrt(5.)*r)
 
-    def dK_dr(self, X, X2):
+    def dK_dr(self, r):
-        r = self._scaled_dist(X, X2)
         return self.variance*(10./3*r -5.*r -5.*np.sqrt(5.)/3*r**2)*np.exp(-np.sqrt(5.)*r)
 
     def Gram_matrix(self,F,F1,F2,F3,lower,upper):
@@ -193,19 +245,55 @@ class Matern52(Stationary):
         return(1./self.variance* (G_coef*G + orig + orig2))
 
 
-
-
 class ExpQuad(Stationary):
     def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='ExpQuad'):
         super(ExpQuad, self).__init__(input_dim, variance, lengthscale, ARD, name)
 
-    def K(self, X, X2=None):
+    def K_of_r(self, r):
-        r = self._scaled_dist(X, X2)
         return self.variance * np.exp(-0.5 * r**2)
 
-    def dK_dr(self, X, X2):
+    def dK_dr(self, r):
-        dist = self._scaled_dist(X, X2)
+        return -r*self.K_of_r(r)
-        return -dist*self.K(X, X2)
+
+class Cosine(Stationary):
+    def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Cosine'):
+        super(Cosine, self).__init__(input_dim, variance, lengthscale, ARD, name)
+
+    def K_of_r(self, r)
+        return self.variance * np.cos(r)
+
+    def dK_dr(self, r):
+        return -self.variance * np.sin(r)
+
+
+class RatQuad(Stationary):
+    """
+    Rational Quadratic Kernel
+
+    .. math::
+
+       k(r) = \sigma^2 \\bigg( 1 + \\frac{r^2}{2} \\bigg)^{- \\alpha}
+
+    """
+
+
+    def __init__(self, input_dim, variance=1., lengthscale=None, power=2., ARD=False, name='ExpQuad'):
+        super(RatQuad, self).__init__(input_dim, variance, lengthscale, ARD, name)
+        self.power = Param('power', power, Logexp())
+        self.add_parameters(self.power)
+
+    def K_of_r(self, r)
+        return self.variance*(1. + r**2/2.)**(-self.power)
+
+    def dK_dr(self, r):
+        return -self.variance*self.power*r*(1. + r**2/2)**(-self.power - 1.)
+
+    def update_gradients_full(self, dL_dK, X, X2=None):
+        super(RatQuad, self).update_gradients_full(dL_dK, X, X2)
+        r = self._scaled_dist(X, X2)
+        r2 = r**2
+        dpow = -2.**self.power*(r2 + 2.)**(-self.power)*np.log(0.5*(r2+2.))
+        self.power.gradient = np.sum(dL_dK*dpow)
 
 
 

GPy/kern/_src/white.py

@@ -1,78 +0,0 @@
-# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-from kern import Kern
-import numpy as np
-from ...core.parameterization import Param
-from ...core.parameterization.transformations import Logexp
-
-class White(Kern):
-    """
-    White noise kernel.
-
-    :param input_dim: the number of input dimensions
-    :type input_dim: int
-    :param variance:
-    :type variance: float
-    """
-    def __init__(self,input_dim,variance=1., name='white'):
-        super(White, self).__init__(input_dim, name)
-        self.input_dim = input_dim
-        self.variance = Param('variance', variance, Logexp())
-        self.add_parameters(self.variance)
-
-    def K(self, X, X2=None):
-        if X2 is None:
-            return np.eye(X.shape[0])*self.variance
-        else:
-            return np.zeros((X.shape[0], X2.shape[0]))
-
-    def Kdiag(self,X):
-        ret = np.ones(X.shape[0])
-        ret[:] = self.variance
-        return ret
-
-    def update_gradients_full(self, dL_dK, X):
-        self.variance.gradient = np.trace(dL_dK)
-
-    def update_gradients_sparse(self, dL_dKmm, dL_dKnm, dL_dKdiag, X, Z):
-        self.variance.gradient = np.trace(dL_dKmm) + np.sum(dL_dKdiag)
-
-    def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z):
-        raise NotImplementedError
-
-    def gradients_X(self,dL_dK,X,X2):
-        return np.zeros_like(X)
-
-    def psi0(self,Z,mu,S,target):
-        pass # target += self.variance
-
-    def dpsi0_dtheta(self,dL_dpsi0,Z,mu,S,target):
-        pass # target += dL_dpsi0.sum()
-
-    def dpsi0_dmuS(self,dL_dpsi0,Z,mu,S,target_mu,target_S):
-        pass
-
-    def psi1(self,Z,mu,S,target):
-        pass
-
-    def dpsi1_dtheta(self,dL_dpsi1,Z,mu,S,target):
-        pass
-
-    def dpsi1_dZ(self,dL_dpsi1,Z,mu,S,target):
-        pass
-
-    def dpsi1_dmuS(self,dL_dpsi1,Z,mu,S,target_mu,target_S):
-        pass
-
-    def psi2(self,Z,mu,S,target):
-        pass
-
-    def dpsi2_dZ(self,dL_dpsi2,Z,mu,S,target):
-        pass
-
-    def dpsi2_dtheta(self,dL_dpsi2,Z,mu,S,target):
-        pass
-
-    def dpsi2_dmuS(self,dL_dpsi2,Z,mu,S,target_mu,target_S):
-        pass

GPy/plotting/matplot_dep/models_plots.py

@@ -20,7 +20,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
       - In higher dimensions, use fixed_inputs to plot the GP  with some of the inputs fixed.
 
     Can plot only part of the data and part of the posterior functions
-    using which_data_rowsm which_data_ycols. 
+    using which_data_rowsm which_data_ycols.
 
     :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
     :type plot_limits: np.array
@@ -56,10 +56,12 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
     if ax is None:
         fig = pb.figure(num=fignum)
         ax = fig.add_subplot(111)
-    
+
-    X, Y, Z = param_to_array(model.X, model.Y, model.Z)
+    X, Y = param_to_array(model.X, model.Y)
-    if model.has_uncertain_inputs(): X_variance = param_to_array(model.q.variance)
+    if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): X_variance = model.X_variance
-    
+
+    if hasattr(model, 'Z'): Z = param_to_array(model.Z)
+
     #work out what the inputs are for plotting (1D or 2D)
     fixed_dims = np.array([i for i,v in fixed_inputs])
     free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims)
@@ -68,8 +70,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
     if len(free_dims) == 1:
 
         #define the frame on which to plot
-        resolution = resolution or 200
+        Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits, resolution=resolution or 200)
-        Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits)
         Xgrid = np.empty((Xnew.shape[0],model.input_dim))
         Xgrid[:,free_dims] = Xnew
         for i,v in fixed_inputs:
@@ -95,7 +96,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
                 ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25)
                 #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs.
 
-        
+
         #add error bars for uncertain (if input uncertainty is being modelled)
         #if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs():
         #    ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(),

GPy/testing/likelihood_tests.py

@@ -10,6 +10,7 @@ from functools import partial
 #np.random.seed(300)
 #np.random.seed(7)
 
+np.seterr(divide='raise')
 def dparam_partial(inst_func, *args):
     """
     If we have a instance method that needs to be called but that doesn't
@@ -149,9 +150,9 @@ class TestNoiseModels(object):
         noise_models = {"Student_t_default": {
                             "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
                                 "vals": [self.var],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 #"constraints": [("t_noise", constrain_positive), ("deg_free", partial(constrain_fixed, value=5))]
                                 },
                             "laplace": True
@@ -159,66 +160,66 @@ class TestNoiseModels(object):
                         "Student_t_1_var": {
                             "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
                                 "vals": [1.0],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 },
                             "laplace": True
                             },
                         "Student_t_small_deg_free": {
                             "model": GPy.likelihoods.StudentT(deg_free=1.5, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
                                 "vals": [self.var],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 },
                             "laplace": True
                             },
                         "Student_t_small_var": {
                             "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
-                                "vals": [0.0001],
+                                "vals": [0.001],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 },
                             "laplace": True
                             },
                         "Student_t_large_var": {
                             "model": GPy.likelihoods.StudentT(deg_free=5, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
                                 "vals": [10.0],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 },
                             "laplace": True
                             },
                         "Student_t_approx_gauss": {
                             "model": GPy.likelihoods.StudentT(deg_free=1000, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
                                 "vals": [self.var],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 },
                             "laplace": True
                             },
                         "Student_t_log": {
                             "model": GPy.likelihoods.StudentT(gp_link=link_functions.Log(), deg_free=5, sigma2=self.var),
                             "grad_params": {
-                                "names": ["t_noise"],
+                                "names": [".*t_noise"],
                                 "vals": [self.var],
-                                "constraints": [("t_noise", constrain_positive), ("deg_free", constrain_fixed)]
+                                "constraints": [(".*t_noise", constrain_positive), (".*deg_free", constrain_fixed)]
                                 },
                             "laplace": True
                             },
                         "Gaussian_default": {
                             "model": GPy.likelihoods.Gaussian(variance=self.var),
                             "grad_params": {
-                                "names": ["variance"],
+                                "names": [".*variance"],
                                 "vals": [self.var],
-                                "constraints": [("variance", constrain_positive)]
+                                "constraints": [(".*variance", constrain_positive)]
                                 },
                             "laplace": True,
-                            "ep": True
+                            "ep": False # FIXME: Should be True when we have it working again
                             },
                         #"Gaussian_log": {
                             #"model": GPy.likelihoods.gaussian(gp_link=link_functions.Log(), variance=self.var, D=self.D, N=self.N),
@@ -252,7 +253,7 @@ class TestNoiseModels(object):
                             "link_f_constraints": [partial(constrain_bounded, lower=0, upper=1)],
                             "laplace": True,
                             "Y": self.binary_Y,
-                            "ep": True
+                            "ep": False # FIXME: Should be True when we have it working again
                             },
                         #"Exponential_default": {
                             #"model": GPy.likelihoods.exponential(),
@@ -515,10 +516,10 @@ class TestNoiseModels(object):
         #Normalize
         Y = Y/Y.max()
         white_var = 1e-6
-        kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
+        kernel = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
         laplace_likelihood = GPy.inference.latent_function_inference.Laplace()
         m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, inference_method=laplace_likelihood)
-        m['white'].constrain_fixed(white_var)
+        m['.*white'].constrain_fixed(white_var)
 
         #Set constraints
         for constrain_param, constraint in constraints:
@@ -541,9 +542,6 @@ class TestNoiseModels(object):
             #NOTE this test appears to be stochastic for some likelihoods (student t?)
             # appears to all be working in test mode right now...
 
-        if not m.checkgrad():
-            import ipdb; ipdb.set_trace()  # XXX BREAKPOINT
-
         assert m.checkgrad(step=step)
 
     ###########
@@ -555,10 +553,10 @@ class TestNoiseModels(object):
         #Normalize
         Y = Y/Y.max()
         white_var = 1e-6
-        kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
+        kernel = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
         ep_inf = GPy.inference.latent_function_inference.EP()
         m = GPy.core.GP(X.copy(), Y.copy(), kernel=kernel, likelihood=model, inference_method=ep_inf)
-        m['white'].constrain_fixed(white_var)
+        m['.*white'].constrain_fixed(white_var)
 
         for param_num in range(len(param_names)):
             name = param_names[param_num]
@@ -634,26 +632,24 @@ class LaplaceTests(unittest.TestCase):
         Y = Y/Y.max()
         #Yc = Y.copy()
         #Yc[75:80] += 1
-        kernel1 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
+        kernel1 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
         #FIXME: Make sure you can copy kernels when params is fixed
         #kernel2 = kernel1.copy()
-        kernel2 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1])
+        kernel2 = GPy.kern.RBF(X.shape[1]) + GPy.kern.White(X.shape[1])
 
         gauss_distr1 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
         exact_inf = GPy.inference.latent_function_inference.ExactGaussianInference()
         m1 = GPy.core.GP(X, Y.copy(), kernel=kernel1, likelihood=gauss_distr1, inference_method=exact_inf)
-        m1['white'].constrain_fixed(1e-6)
+        m1['.*white'].constrain_fixed(1e-6)
-        m1['variance'] = initial_var_guess
+        m1['.*rbf.variance'] = initial_var_guess
-        m1['variance'].constrain_bounded(1e-4, 10)
+        m1['.*rbf.variance'].constrain_bounded(1e-4, 10)
-        m1['rbf'].constrain_bounded(1e-4, 10)
         m1.randomize()
 
         gauss_distr2 = GPy.likelihoods.Gaussian(variance=initial_var_guess)
         laplace_inf = GPy.inference.latent_function_inference.Laplace()
         m2 = GPy.core.GP(X, Y.copy(), kernel=kernel2, likelihood=gauss_distr2, inference_method=laplace_inf)
-        m2['white'].constrain_fixed(1e-6)
+        m2['.*white'].constrain_fixed(1e-6)
-        m2['rbf'].constrain_bounded(1e-4, 10)
+        m2['.*rbf.variance'].constrain_bounded(1e-4, 10)
-        m2['variance'].constrain_bounded(1e-4, 10)
         m2.randomize()
 
         if debug: