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