From 20cc336a7a20ec3f6cec5bfba67242f13342babf Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 21 Jan 2013 13:29:25 +0000 Subject: [PATCH 1/4] added periodic kernels --- GPy/kern/__init__.py | 2 +- GPy/kern/constructors.py | 57 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 1 deletion(-) diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index 4a36d6d0..eae717be 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -2,5 +2,5 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) -from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, linear_ARD, rbf_sympy, sympykern +from constructors import rbf, Matern32, Matern52, exponential, linear, white, bias, finite_dimensional, spline, Brownian, linear_ARD, rbf_sympy, sympykern, periodic_exponential, periodic_Matern32, periodic_Matern52 from kern import kern diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index ab39fcb2..49639711 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -16,6 +16,9 @@ from bias import bias as biaspart from finite_dimensional import finite_dimensional as finite_dimensionalpart from spline import spline as splinepart from Brownian import Brownian as Brownianpart +from periodic_exponential import periodic_exponential as periodic_exponentialpart +from periodic_Matern32 import periodic_Matern32 as periodic_Matern32part +from periodic_Matern52 import periodic_Matern52 as periodic_Matern52part #TODO these s=constructors are not as clean as we'd like. Tidy the code up #using meta-classes to make the objects construct properly wthout them. @@ -196,3 +199,57 @@ def sympykern(D,k): A kernel from a symbolic sympy representation """ return kern(D,[spkern(D,k)]) + +def periodic_exponential(D=1,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): + """ + Construct an periodic exponential kernel + + :param D: dimensionality, only defined for D=1 + :type D: int + :param variance: the variance of the kernel + :type variance: float + :param lengthscale: the lengthscale of the kernel + :type lengthscale: float + :param period: the period + :type period: float + :param n_freq: the number of frequencies considered for the periodic subspace + :type n_freq: int + """ + part = periodic_exponentialpart(D,variance, lengthscale, ARD) + return kern(D, [part]) + +def periodic_Matern32(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): + """ + Construct a periodic Matern 3/2 kernel. + + :param D: dimensionality, only defined for D=1 + :type D: int + :param variance: the variance of the kernel + :type variance: float + :param lengthscale: the lengthscale of the kernel + :type lengthscale: float + :param period: the period + :type period: float + :param n_freq: the number of frequencies considered for the periodic subspace + :type n_freq: int + """ + part = periodic_Matern32part(D,variance, lengthscale, ARD) + return kern(D, [part]) + +def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): + """ + Construct a periodic Matern 5/2 kernel. + + :param D: dimensionality, only defined for D=1 + :type D: int + :param variance: the variance of the kernel + :type variance: float + :param lengthscale: the lengthscale of the kernel + :type lengthscale: float + :param period: the period + :type period: float + :param n_freq: the number of frequencies considered for the periodic subspace + :type n_freq: int + """ + part = periodic_Matern52part(D,variance, lengthscale, ARD) + return kern(D, [part]) From 8d98e09b9eaabbcada161d98d293b3f2f999ba60 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 21 Jan 2013 13:31:27 +0000 Subject: [PATCH 2/4] added some more files for periodic kernels --- GPy/kern/periodic_Matern32.py | 168 +++++++++++++++++++++++++++++ GPy/kern/periodic_Matern52.py | 177 +++++++++++++++++++++++++++++++ GPy/kern/periodic_exponential.py | 164 ++++++++++++++++++++++++++++ 3 files changed, 509 insertions(+) create mode 100644 GPy/kern/periodic_Matern32.py create mode 100644 GPy/kern/periodic_Matern52.py create mode 100644 GPy/kern/periodic_exponential.py diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py new file mode 100644 index 00000000..8f616baf --- /dev/null +++ b/GPy/kern/periodic_Matern32.py @@ -0,0 +1,168 @@ +from kernpart import kernpart +import numpy as np +from GPy.util.linalg import mdot, pdinv + +class periodic_Matern32(kernpart): + """ + Kernel of the periodic subspace (up to a given frequency) of a Matern 3/2 RKHS. Only defined for D=1. + + :param D: the number of input dimensions + :type D: 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 (D,) + :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,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): + assert D==1 + self.name = 'periodic_Mat32' + self.D = D + if lengthscale is not None: + assert lengthscale.shape==(self.D,) + else: + lengthscale = np.ones(self.D) + self.lower,self.upper = lower, upper + self.Nparam = 3 + self.n_freq = n_freq + self.n_basis = 2*n_freq + self.set_param(np.hstack((variance,lengthscale,period))) + + def _cos(self,alpha,omega,phase): + def f(x): + return alpha*np.cos(omega*x+phase) + return f + 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 + 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_param(self): + """return the value of the parameters.""" + return np.hstack((self.variance,self.lengthscale,self.period)) + def set_param(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.""" + 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) + np.add(mdot(FX,self.Gi,FX2.T), target,target) + + def dK_dtheta(self,partial,X,X2,target): + """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" + 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*partial) + #np.add(target[:,:,1],dK_dlen, target[:,:,1]) + target[1] += np.sum(dK_dlen*partial) + #np.add(target[:,:,2],dK_dper, target[:,:,2]) + target[2] += np.sum(dK_dper*partial) + + + diff --git a/GPy/kern/periodic_Matern52.py b/GPy/kern/periodic_Matern52.py new file mode 100644 index 00000000..13540435 --- /dev/null +++ b/GPy/kern/periodic_Matern52.py @@ -0,0 +1,177 @@ +from kernpart import kernpart +import numpy as np +from GPy.util.linalg import mdot, pdinv + +class periodic_Matern52(kernpart): + """ + Kernel of the periodic subspace (up to a given frequency) of a Matern 5/2 RKHS. Only defined for D=1. + + :param D: the number of input dimensions + :type D: 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 (D,) + :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,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): + assert D==1 + self.name = 'periodic_Mat52' + self.D = D + if lengthscale is not None: + assert lengthscale.shape==(self.D,) + else: + lengthscale = np.ones(self.D) + self.lower,self.upper = lower, upper + self.Nparam = 3 + self.n_freq = n_freq + self.n_basis = 2*n_freq + self.set_param(np.hstack((variance,lengthscale,period))) + + def _cos(self,alpha,omega,phase): + def f(x): + return alpha*np.cos(omega*x+phase) + return f + 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 + 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_param(self): + """return the value of the parameters.""" + return np.hstack((self.variance,self.lengthscale,self.period)) + def set_param(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 = [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.)] + + self.basis_alpha = np.ones((2*self.n_freq,)) + 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, 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] + 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.""" + 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) + np.add(mdot(FX,self.Gi,FX2.T), target,target) + + def dK_dtheta(self,partial,X,X2,target): + """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" + 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, 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,FX2.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,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, -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)) + 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)) + 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 = 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*partial) + #np.add(target[:,:,1],dK_dlen, target[:,:,1]) + target[1] += np.sum(dK_dlen*partial) + #np.add(target[:,:,2],dK_dper, target[:,:,2]) + target[2] += np.sum(dK_dper*partial) + + diff --git a/GPy/kern/periodic_exponential.py b/GPy/kern/periodic_exponential.py new file mode 100644 index 00000000..11716e13 --- /dev/null +++ b/GPy/kern/periodic_exponential.py @@ -0,0 +1,164 @@ +from kernpart import kernpart +import numpy as np +from GPy.util.linalg import mdot, pdinv + +class periodic_exponential(kernpart): + """ + Kernel of the periodic subspace (up to a given frequency) of a exponential (Matern 1/2) RKHS. Only defined for D=1. + + :param D: the number of input dimensions + :type D: 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 (D,) + :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,D=1,variance=1.,lengthscale=None,period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): + assert D==1 + self.name = 'periodic_exp' + self.D = D + if lengthscale is not None: + assert lengthscale.shape==(self.D,) + else: + lengthscale = np.ones(self.D) + self.lower,self.upper = lower, upper + self.Nparam = 3 + self.n_freq = n_freq + self.n_basis = 2*n_freq + self.set_param(np.hstack((variance,lengthscale,period))) + + def _cos(self,alpha,omega,phase): + def f(x): + return alpha*np.cos(omega*x+phase) + return f + 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 + 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_param(self): + """return the value of the parameters.""" + return np.hstack((self.variance,self.lengthscale,self.period)) + def set_param(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 = [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)],[])) + 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.""" + 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) + np.add(mdot(FX,self.Gi,FX2.T), target,target) + + def dK_dtheta(self,partial,X,X2,target): + """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" + 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) + + # np.add(target[:,:,0],dK_dvar, target[:,:,0]) + target[0] += np.sum(dK_dvar*partial) + #np.add(target[:,:,1],dK_dlen, target[:,:,1]) + target[1] += np.sum(dK_dlen*partial) + #np.add(target[:,:,2],dK_dper, target[:,:,2]) + target[2] += np.sum(dK_dper*partial) + + From b66eab22fbeac48c2b738845fc5a9125446fc06c Mon Sep 17 00:00:00 2001 From: Nicolas Date: Mon, 21 Jan 2013 18:03:24 +0000 Subject: [PATCH 3/4] few bugs fixed in periodic kernels --- GPy/kern/constructors.py | 6 ++-- GPy/kern/periodic_Matern32.py | 50 +++++++++++++++++-------------- GPy/kern/periodic_Matern52.py | 49 +++++++++++++++++------------- GPy/kern/periodic_exponential.py | 51 ++++++++++++++++++-------------- 4 files changed, 86 insertions(+), 70 deletions(-) diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 49639711..1bd241e7 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -215,7 +215,7 @@ def periodic_exponential(D=1,variance=1., lengthscale=None, period=2*np.pi,n_fre :param n_freq: the number of frequencies considered for the periodic subspace :type n_freq: int """ - part = periodic_exponentialpart(D,variance, lengthscale, ARD) + part = periodic_exponentialpart(D,variance, lengthscale, period, n_freq, lower, upper) return kern(D, [part]) def periodic_Matern32(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): @@ -233,7 +233,7 @@ def periodic_Matern32(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10, :param n_freq: the number of frequencies considered for the periodic subspace :type n_freq: int """ - part = periodic_Matern32part(D,variance, lengthscale, ARD) + part = periodic_Matern32part(D,variance, lengthscale, period, n_freq, lower, upper) return kern(D, [part]) def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10,lower=0.,upper=4*np.pi): @@ -251,5 +251,5 @@ def periodic_Matern52(D,variance=1., lengthscale=None, period=2*np.pi,n_freq=10, :param n_freq: the number of frequencies considered for the periodic subspace :type n_freq: int """ - part = periodic_Matern52part(D,variance, lengthscale, ARD) + part = periodic_Matern52part(D,variance, lengthscale, period, n_freq, lower, upper) return kern(D, [part]) diff --git a/GPy/kern/periodic_Matern32.py b/GPy/kern/periodic_Matern32.py index 8f616baf..3e81044c 100644 --- a/GPy/kern/periodic_Matern32.py +++ b/GPy/kern/periodic_Matern32.py @@ -32,7 +32,7 @@ class periodic_Matern32(kernpart): self.Nparam = 3 self.n_freq = n_freq self.n_basis = 2*n_freq - self.set_param(np.hstack((variance,lengthscale,period))) + self._set_params(np.hstack((variance,lengthscale,period))) def _cos(self,alpha,omega,phase): def f(x): @@ -45,16 +45,16 @@ class periodic_Matern32(kernpart): psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2)) return r,omega[:,0:1], psi 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) ) + 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]) + #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_param(self): + def _get_params(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscale,self.period)) - def set_param(self,x): + def _set_params(self,x): """set the value of the parameters.""" assert x.size==3 self.variance = x[0] @@ -67,15 +67,15 @@ class periodic_Matern32(kernpart): 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): + 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): + 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)) @@ -88,11 +88,18 @@ class periodic_Matern32(kernpart): def K(self,X,X2,target): """Compute the covariance matrix between X and X2.""" - 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) + 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) + def dK_dtheta(self,partial,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X @@ -127,34 +134,34 @@ class periodic_Matern32(kernpart): 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) + 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) + #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]) + #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) + 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] + 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]) @@ -163,6 +170,3 @@ class periodic_Matern32(kernpart): target[1] += np.sum(dK_dlen*partial) #np.add(target[:,:,2],dK_dper, target[:,:,2]) target[2] += np.sum(dK_dper*partial) - - - diff --git a/GPy/kern/periodic_Matern52.py b/GPy/kern/periodic_Matern52.py index 13540435..abdd6a7d 100644 --- a/GPy/kern/periodic_Matern52.py +++ b/GPy/kern/periodic_Matern52.py @@ -32,29 +32,31 @@ class periodic_Matern52(kernpart): self.Nparam = 3 self.n_freq = n_freq self.n_basis = 2*n_freq - self.set_param(np.hstack((variance,lengthscale,period))) + 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 + 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 + 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) ) + 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]) + #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_param(self): + def _get_params(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscale,self.period)) - def set_param(self,x): + def _set_params(self,x): """set the value of the parameters.""" assert x.size==3 self.variance = x[0] @@ -67,15 +69,15 @@ class periodic_Matern52(kernpart): self.basis_alpha = np.ones((2*self.n_freq,)) 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): + 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): + 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)) 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)) @@ -85,16 +87,23 @@ class periodic_Matern52(kernpart): 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] - 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) + 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.""" - 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) + 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) + def dK_dtheta(self,partial,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X @@ -131,25 +140,25 @@ class periodic_Matern52(kernpart): 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) + 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) + #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]) + #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)) 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) + 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 @@ -173,5 +182,3 @@ class periodic_Matern52(kernpart): target[1] += np.sum(dK_dlen*partial) #np.add(target[:,:,2],dK_dper, target[:,:,2]) target[2] += np.sum(dK_dper*partial) - - diff --git a/GPy/kern/periodic_exponential.py b/GPy/kern/periodic_exponential.py index 11716e13..2637afe5 100644 --- a/GPy/kern/periodic_exponential.py +++ b/GPy/kern/periodic_exponential.py @@ -32,29 +32,31 @@ class periodic_exponential(kernpart): self.Nparam = 3 self.n_freq = n_freq self.n_basis = 2*n_freq - self.set_param(np.hstack((variance,lengthscale,period))) + 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 + 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 + 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) ) + 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]) + #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_param(self): + def _get_params(self): """return the value of the parameters.""" return np.hstack((self.variance,self.lengthscale,self.period)) - def set_param(self,x): + def _set_params(self,x): """set the value of the parameters.""" assert x.size==3 self.variance = x[0] @@ -67,32 +69,37 @@ class periodic_exponential(kernpart): 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): + 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): + 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.""" - 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) + 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) + def dK_dtheta(self,partial,X,X2,target): """derivative of the covariance matrix with respect to the parameters (shape is NxMxNparam)""" if X2 is None: X2 = X @@ -125,24 +132,24 @@ class periodic_exponential(kernpart): 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) + 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) + #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]) + #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)) + 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) @@ -151,7 +158,7 @@ class periodic_exponential(kernpart): 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) # np.add(target[:,:,0],dK_dvar, target[:,:,0]) @@ -160,5 +167,3 @@ class periodic_exponential(kernpart): target[1] += np.sum(dK_dlen*partial) #np.add(target[:,:,2],dK_dper, target[:,:,2]) target[2] += np.sum(dK_dper*partial) - - From 4304c1dbe2f627f6adbbd5329a3baa22b9267997 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Thu, 24 Jan 2013 16:55:15 +0000 Subject: [PATCH 4/4] Missing scale and location arguments. --- GPy/inference/likelihoods.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/GPy/inference/likelihoods.py b/GPy/inference/likelihoods.py index 5f0eb7ff..c9b36e10 100644 --- a/GPy/inference/likelihoods.py +++ b/GPy/inference/likelihoods.py @@ -9,7 +9,7 @@ import pylab as pb from ..util.plot import gpplot class likelihood: - def __init__(self,Y): + def __init__(self,Y,location=0,scale=1): """ Likelihood class for doing Expectation propagation @@ -18,6 +18,8 @@ class likelihood: """ self.Y = Y self.N = self.Y.shape[0] + self.location = location + self.scale = scale def plot1Da(self,X_new,Mean_new,Var_new,X_u,Mean_u,Var_u): """