1# Copyright (c) 2012, GPy authors (see AUTHORS.txt). 2# Licensed under the BSD 3-clause license (see LICENSE.txt) 3 4 5import numpy as np 6from .kern import Kern 7from ...util.linalg import mdot 8from ...util.decorators import silence_errors 9from ...core.parameterization.param import Param 10from paramz.transformations import Logexp 11 12class Periodic(Kern): 13 def __init__(self, input_dim, variance, lengthscale, period, n_freq, lower, upper, active_dims, name): 14 """ 15 :type input_dim: int 16 :param variance: the variance of the Matern kernel 17 :type variance: float 18 :param lengthscale: the lengthscale of the Matern kernel 19 :type lengthscale: np.ndarray of size (input_dim,) 20 :param period: the period 21 :type period: float 22 :param n_freq: the number of frequencies considered for the periodic subspace 23 :type n_freq: int 24 :rtype: kernel object 25 """ 26 27 assert input_dim==1, "Periodic kernels are only defined for input_dim=1" 28 super(Periodic, self).__init__(input_dim, active_dims, name) 29 self.input_dim = input_dim 30 self.lower,self.upper = lower, upper 31 self.n_freq = n_freq 32 self.n_basis = 2*n_freq 33 self.variance = Param('variance', np.float64(variance), Logexp()) 34 self.lengthscale = Param('lengthscale', np.float64(lengthscale), Logexp()) 35 self.period = Param('period', np.float64(period), Logexp()) 36 self.link_parameters(self.variance, self.lengthscale, self.period) 37 38 def _cos(self, alpha, omega, phase): 39 def f(x): 40 return alpha*np.cos(omega*x + phase) 41 return f 42 43 @silence_errors 44 def _cos_factorization(self, alpha, omega, phase): 45 r1 = np.sum(alpha*np.cos(phase),axis=1)[:,None] 46 r2 = np.sum(alpha*np.sin(phase),axis=1)[:,None] 47 r = np.sqrt(r1**2 + r2**2) 48 psi = np.where(r1 != 0, (np.arctan(r2/r1) + (r1<0.)*np.pi),np.arcsin(r2)) 49 return r,omega[:,0:1], psi 50 51 @silence_errors 52 def _int_computation(self,r1,omega1,phi1,r2,omega2,phi2): 53 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) ) 54 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) 55 Gint = np.dot(r1,r2.T)/2 * np.where(np.isnan(Gint1),Gint2,Gint1) 56 return Gint 57 58 def K(self, X, X2=None): 59 FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) 60 if X2 is None: 61 FX2 = FX 62 else: 63 FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2) 64 return mdot(FX,self.Gi,FX2.T) 65 66 def Kdiag(self,X): 67 return np.diag(self.K(X)) 68 69 70class PeriodicExponential(Periodic): 71 """ 72 Kernel of the periodic subspace (up to a given frequency) of a exponential 73 (Matern 1/2) RKHS. 74 75 Only defined for input_dim=1. 76 """ 77 78 def __init__(self, input_dim=1, variance=1., lengthscale=1., period=2.*np.pi, n_freq=10, lower=0., upper=4*np.pi, active_dims=None, name='periodic_exponential'): 79 super(PeriodicExponential, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, active_dims, name) 80 81 def parameters_changed(self): 82 self.a = [1./self.lengthscale, 1.] 83 self.b = [1] 84 85 self.basis_alpha = np.ones((self.n_basis,)) 86 self.basis_omega = (2*np.pi*np.arange(1,self.n_freq+1)/self.period).repeat(2) 87 self.basis_phi = np.zeros(self.n_freq * 2) 88 self.basis_phi[::2] = -np.pi/2 89 90 self.G = self.Gram_matrix() 91 self.Gi = np.linalg.inv(self.G) 92 93 def Gram_matrix(self): 94 La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega)) 95 Lo = np.column_stack((self.basis_omega,self.basis_omega)) 96 Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2)) 97 r,omega,phi = self._cos_factorization(La,Lo,Lp) 98 Gint = self._int_computation( r,omega,phi, r,omega,phi) 99 Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None] 100 return(self.lengthscale/(2*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T)) 101 102 @silence_errors 103 def update_gradients_full(self, dL_dK, X, X2=None): 104 """derivative of the covariance matrix with respect to the parameters (shape is N x num_inducing x num_params)""" 105 if X2 is None: X2 = X 106 FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) 107 FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2) 108 109 La = np.column_stack((self.a[0]*np.ones((self.n_basis,1)),self.a[1]*self.basis_omega)) 110 Lo = np.column_stack((self.basis_omega,self.basis_omega)) 111 Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2)) 112 r,omega,phi = self._cos_factorization(La,Lo,Lp) 113 Gint = self._int_computation( r,omega,phi, r,omega,phi) 114 115 Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None] 116 117 #dK_dvar 118 dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T) 119 120 #dK_dlen 121 da_dlen = [-1./self.lengthscale**2,0.] 122 dLa_dlen = np.column_stack((da_dlen[0]*np.ones((self.n_basis,1)),da_dlen[1]*self.basis_omega)) 123 r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp) 124 dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi) 125 dGint_dlen = dGint_dlen + dGint_dlen.T 126 dG_dlen = 1./2*Gint + self.lengthscale/2*dGint_dlen 127 dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T) 128 129 #dK_dper 130 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) 131 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) 132 133 dLa_dper = np.column_stack((-self.a[0]*self.basis_omega/self.period, -self.a[1]*self.basis_omega**2/self.period)) 134 dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi)) 135 r1,omega1,phi1 = self._cos_factorization(dLa_dper,Lo,dLp_dper) 136 137 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)) 138 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)) 139 # 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)) 140 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)) 141 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)) 142 IPPprim = np.where(np.logical_or(np.isnan(IPPprim1), np.isinf(IPPprim1)), IPPprim2, IPPprim1) 143 144 145 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) 146 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) 147 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) 148 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) 149 #IPPint2[0,0] = (self.upper**2 - self.lower**2)*np.cos(phi[0,0])*np.cos(phi1[0,0]) 150 IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1) 151 152 dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period)) 153 dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2)) 154 r2,omega2,phi2 = dLa_dper2.T,Lo[:,0:1],dLp_dper2.T 155 156 dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi) 157 dGint_dper = dGint_dper + dGint_dper.T 158 159 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] 160 161 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))) 162 163 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) 164 165 self.variance.gradient = np.sum(dK_dvar*dL_dK) 166 self.lengthscale.gradient = np.sum(dK_dlen*dL_dK) 167 self.period.gradient = np.sum(dK_dper*dL_dK) 168 169 170 171class PeriodicMatern32(Periodic): 172 """ 173 Kernel of the periodic subspace (up to a given frequency) of a Matern 3/2 RKHS. Only defined for input_dim=1. 174 175 :param input_dim: the number of input dimensions 176 :type input_dim: int 177 :param variance: the variance of the Matern kernel 178 :type variance: float 179 :param lengthscale: the lengthscale of the Matern kernel 180 :type lengthscale: np.ndarray of size (input_dim,) 181 :param period: the period 182 :type period: float 183 :param n_freq: the number of frequencies considered for the periodic subspace 184 :type n_freq: int 185 :rtype: kernel object 186 187 """ 188 189 def __init__(self, input_dim=1, variance=1., lengthscale=1., period=2.*np.pi, n_freq=10, lower=0., upper=4*np.pi, active_dims=None, name='periodic_Matern32'): 190 super(PeriodicMatern32, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, active_dims, name) 191 def parameters_changed(self): 192 self.a = [3./self.lengthscale**2, 2*np.sqrt(3)/self.lengthscale, 1.] 193 self.b = [1,self.lengthscale**2/3] 194 195 self.basis_alpha = np.ones((self.n_basis,)) 196 self.basis_omega = (2*np.pi*np.arange(1,self.n_freq+1)/self.period).repeat(2) 197 self.basis_phi = np.zeros(self.n_freq * 2) 198 self.basis_phi[::2] = -np.pi/2 199 200 self.G = self.Gram_matrix() 201 self.Gi = np.linalg.inv(self.G) 202 203 def Gram_matrix(self): 204 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)) 205 Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega)) 206 Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi)) 207 r,omega,phi = self._cos_factorization(La,Lo,Lp) 208 Gint = self._int_computation( r,omega,phi, r,omega,phi) 209 210 Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None] 211 F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None] 212 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)) 213 214 215 @silence_errors 216 def update_gradients_full(self,dL_dK,X,X2): 217 """derivative of the covariance matrix with respect to the parameters (shape is num_data x num_inducing x num_params)""" 218 if X2 is None: X2 = X 219 FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) 220 FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2) 221 222 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)) 223 Lo = np.column_stack((self.basis_omega,self.basis_omega,self.basis_omega)) 224 Lp = np.column_stack((self.basis_phi,self.basis_phi+np.pi/2,self.basis_phi+np.pi)) 225 r,omega,phi = self._cos_factorization(La,Lo,Lp) 226 Gint = self._int_computation( r,omega,phi, r,omega,phi) 227 228 Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None] 229 F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None] 230 231 #dK_dvar 232 dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T) 233 234 #dK_dlen 235 da_dlen = [-6/self.lengthscale**3,-2*np.sqrt(3)/self.lengthscale**2,0.] 236 db_dlen = [0.,2*self.lengthscale/3.] 237 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)) 238 r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp) 239 dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi) 240 dGint_dlen = dGint_dlen + dGint_dlen.T 241 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) 242 dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T) 243 244 #dK_dper 245 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) 246 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) 247 248 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)) 249 dLp_dper = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi,self.basis_phi+np.pi*3/2)) 250 r1,omega1,phi1 = self._cos_factorization(dLa_dper,Lo,dLp_dper) 251 252 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)) 253 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)) 254 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)) 255 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)) 256 IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1) 257 258 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) 259 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) 260 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) 261 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) 262 IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1) 263 264 dLa_dper2 = np.column_stack((-self.a[1]*self.basis_omega/self.period, -2*self.a[2]*self.basis_omega**2/self.period)) 265 dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2,self.basis_phi+np.pi)) 266 r2,omega2,phi2 = self._cos_factorization(dLa_dper2,Lo[:,0:2],dLp_dper2) 267 268 dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi) 269 dGint_dper = dGint_dper + dGint_dper.T 270 271 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] 272 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] 273 274 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))) 275 276 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) 277 278 self.variance.gradient = np.sum(dK_dvar*dL_dK) 279 self.lengthscale.gradient = np.sum(dK_dlen*dL_dK) 280 self.period.gradient = np.sum(dK_dper*dL_dK) 281 282 283 284class PeriodicMatern52(Periodic): 285 """ 286 Kernel of the periodic subspace (up to a given frequency) of a Matern 5/2 RKHS. Only defined for input_dim=1. 287 288 :param input_dim: the number of input dimensions 289 :type input_dim: int 290 :param variance: the variance of the Matern kernel 291 :type variance: float 292 :param lengthscale: the lengthscale of the Matern kernel 293 :type lengthscale: np.ndarray of size (input_dim,) 294 :param period: the period 295 :type period: float 296 :param n_freq: the number of frequencies considered for the periodic subspace 297 :type n_freq: int 298 :rtype: kernel object 299 300 """ 301 302 def __init__(self, input_dim=1, variance=1., lengthscale=1., period=2.*np.pi, n_freq=10, lower=0., upper=4*np.pi, active_dims=None, name='periodic_Matern52'): 303 super(PeriodicMatern52, self).__init__(input_dim, variance, lengthscale, period, n_freq, lower, upper, active_dims, name) 304 305 def parameters_changed(self): 306 self.a = [5*np.sqrt(5)/self.lengthscale**3, 15./self.lengthscale**2,3*np.sqrt(5)/self.lengthscale, 1.] 307 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.)] 308 309 self.basis_alpha = np.ones((2*self.n_freq,)) 310 self.basis_omega = (2*np.pi*np.arange(1,self.n_freq+1)/self.period).repeat(2) 311 self.basis_phi = np.zeros(self.n_freq * 2) 312 self.basis_phi[::2] = -np.pi/2 313 314 self.G = self.Gram_matrix() 315 self.Gi = np.linalg.inv(self.G) 316 317 def Gram_matrix(self): 318 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)) 319 Lo = np.column_stack((self.basis_omega, self.basis_omega, self.basis_omega, self.basis_omega)) 320 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)) 321 r,omega,phi = self._cos_factorization(La,Lo,Lp) 322 Gint = self._int_computation( r,omega,phi, r,omega,phi) 323 324 Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None] 325 F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None] 326 F2lower = np.array(self._cos(self.basis_alpha*self.basis_omega**2,self.basis_omega,self.basis_phi+np.pi)(self.lower))[:,None] 327 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) 328 return(3*self.lengthscale**5/(400*np.sqrt(5)*self.variance) * Gint + 1./self.variance*lower_terms) 329 330 @silence_errors 331 def update_gradients_full(self, dL_dK, X, X2=None): 332 if X2 is None: X2 = X 333 FX = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X) 334 FX2 = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X2) 335 336 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)) 337 Lo = np.column_stack((self.basis_omega, self.basis_omega, self.basis_omega, self.basis_omega)) 338 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)) 339 r,omega,phi = self._cos_factorization(La,Lo,Lp) 340 Gint = self._int_computation( r,omega,phi, r,omega,phi) 341 342 Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None] 343 F1lower = np.array(self._cos(self.basis_alpha*self.basis_omega,self.basis_omega,self.basis_phi+np.pi/2)(self.lower))[:,None] 344 F2lower = np.array(self._cos(self.basis_alpha*self.basis_omega**2,self.basis_omega,self.basis_phi+np.pi)(self.lower))[:,None] 345 346 #dK_dvar 347 dK_dvar = 1./self.variance*mdot(FX,self.Gi,FX2.T) 348 349 #dK_dlen 350 da_dlen = [-3*self.a[0]/self.lengthscale, -2*self.a[1]/self.lengthscale, -self.a[2]/self.lengthscale, 0.] 351 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] 352 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)) 353 r1,omega1,phi1 = self._cos_factorization(dLa_dlen,Lo,Lp) 354 dGint_dlen = self._int_computation(r1,omega1,phi1, r,omega,phi) 355 dGint_dlen = dGint_dlen + dGint_dlen.T 356 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) 357 dG_dlen = 15*self.lengthscale**4/(400*np.sqrt(5))*Gint + 3*self.lengthscale**5/(400*np.sqrt(5))*dGint_dlen + dlower_terms_dlen 358 dK_dlen = -mdot(FX,self.Gi,dG_dlen/self.variance,self.Gi,FX2.T) 359 360 #dK_dper 361 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) 362 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) 363 364 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)) 365 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)) 366 r1,omega1,phi1 = self._cos_factorization(dLa_dper,Lo,dLp_dper) 367 368 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)) 369 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)) 370 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)) 371 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)) 372 IPPprim = np.where(np.isnan(IPPprim1),IPPprim2,IPPprim1) 373 374 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) 375 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) 376 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) 377 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) 378 IPPint = np.where(np.isnan(IPPint1),IPPint2,IPPint1) 379 380 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)) 381 dLp_dper2 = np.column_stack((self.basis_phi+np.pi/2, self.basis_phi+np.pi, self.basis_phi+np.pi*3/2)) 382 r2,omega2,phi2 = self._cos_factorization(dLa_dper2,Lo[:,0:2],dLp_dper2) 383 384 dGint_dper = np.dot(r,r1.T)/2 * (IPPprim - IPPint) + self._int_computation(r2,omega2,phi2, r,omega,phi) 385 dGint_dper = dGint_dper + dGint_dper.T 386 387 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] 388 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] 389 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] 390 391 dlower_terms_dper = self.b[0] * (np.dot(dFlower_dper,Flower.T) + np.dot(Flower.T,dFlower_dper)) 392 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) 393 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) 394 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) 395 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) 396 397 dG_dper = 1./self.variance*(3*self.lengthscale**5/(400*np.sqrt(5))*dGint_dper + 0.5*dlower_terms_dper) 398 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) 399 400 self.variance.gradient = np.sum(dK_dvar*dL_dK) 401 self.lengthscale.gradient = np.sum(dK_dlen*dL_dK) 402 self.period.gradient = np.sum(dK_dper*dL_dK) 403