diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index d858ad5b..f91f5ac6 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -3,8 +3,9 @@ from _src.kern import Kern from _src.linear import Linear from _src.static import Bias, White from _src.brownian import Brownian -from _src.stationary import Exponential, Matern32, Matern52, ExpQuad, RatQuad +from _src.stationary import Exponential, Matern32, Matern52, ExpQuad, RatQuad, Cosine from _src.mlp import MLP +from _src.periodic import PeriodicExponential, PeriodicMatern32, PeriodicMatern52 #import coregionalize #import eq_ode1 #import finite_dimensional @@ -18,9 +19,6 @@ from _src.mlp import MLP #import periodic_Matern32 #import periodic_Matern52 #import poly -#import prod_orthogonal -#import prod -#import rational_quadratic #import rbfcos #import rbf #import rbf_inv diff --git a/GPy/kern/_src/mlp.py b/GPy/kern/_src/mlp.py index f2f40e62..85792acd 100644 --- a/GPy/kern/_src/mlp.py +++ b/GPy/kern/_src/mlp.py @@ -32,10 +32,10 @@ class MLP(Kern): """ def __init__(self, input_dim, variance=1., weight_variance=1., bias_variance=100., name='mlp'): - super(Linear, self).__init__(input_dim, name) - self.variance = Param('variance', variance, Logexp) - self.weight_variance = Param('weight_variance', weight_variance, Logexp) - self.bias_variance = Param('bias_variance', bias_variance, Logexp) + super(MLP, self).__init__(input_dim, name) + self.variance = Param('variance', variance, Logexp()) + self.weight_variance = Param('weight_variance', weight_variance, Logexp()) + self.bias_variance = Param('bias_variance', bias_variance, Logexp()) self.add_parameters(self.variance, self.weight_variance, self.bias_variance) @@ -109,14 +109,15 @@ class MLP(Kern): """Pre-computations for the covariance matrix (used for computing the covariance and its gradients.""" if X2 is None: 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: 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_numer = self._K_inner_prod*self.weight_variance + self.bias_variance self._K_asin_arg = self._K_numer/self._K_denom self._K_dvar = four_over_tau*np.arcsin(self._K_asin_arg) diff --git a/GPy/kern/_src/periodic_Matern52.py b/GPy/kern/_src/periodic.py similarity index 53% rename from GPy/kern/_src/periodic_Matern52.py rename to GPy/kern/_src/periodic.py index 1f9d90b3..e4e659a2 100644 --- a/GPy/kern/_src/periodic_Matern52.py +++ b/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 GPy.util.decorators import silence_errors +from kern import Kern +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): - assert input_dim==1, "Periodic kernels are only defined for input_dim=1" - 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 __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'): + super(PeriodicMatern52, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, name) + 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): - """derivative of the covariance matrix with respect to the parameters (shape is num_data x num_inducing x num_params)""" + def update_gradients_full(self, dL_dK, X, X2=None): 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]) - 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) + 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) - @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) diff --git a/GPy/kern/_src/periodic_Matern32.py b/GPy/kern/_src/periodic_Matern32.py deleted file mode 100644 index 24ec45f9..00000000 --- a/GPy/kern/_src/periodic_Matern32.py +++ /dev/null @@ -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) diff --git a/GPy/kern/_src/periodic_exponential.py b/GPy/kern/_src/periodic_exponential.py deleted file mode 100644 index 4562cd56..00000000 --- a/GPy/kern/_src/periodic_exponential.py +++ /dev/null @@ -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) diff --git a/GPy/kern/_src/prod_orthogonal.py b/GPy/kern/_src/prod_orthogonal.py deleted file mode 100644 index e7dd1fdc..00000000 --- a/GPy/kern/_src/prod_orthogonal.py +++ /dev/null @@ -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 + '' + 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) - diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index e47b2b63..86db393a 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -205,6 +205,19 @@ class ExpQuad(Stationary): dist = self._scaled_dist(X, X2) 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(self, X, X2=None): + r = self._scaled_dist(X, X2) + return self.variance * np.cos(r) + + def dK_dr(self, X, X2): + r = self._scaled_dist(X, X2) + return -self.variance * np.sin(r) + + class RatQuad(Stationary): """ Rational Quadratic Kernel diff --git a/GPy/kern/_src/white.py b/GPy/kern/_src/white.py deleted file mode 100644 index 1fc022f5..00000000 --- a/GPy/kern/_src/white.py +++ /dev/null @@ -1,8 +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 - diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 4d49dd12..974b5740 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -68,8 +68,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) + Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits, resolution=resolution or 200) Xgrid = np.empty((Xnew.shape[0],model.input_dim)) Xgrid[:,free_dims] = Xnew for i,v in fixed_inputs: