kernel.py 14.7 KB
Newer Older
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
1
2
3
4
5
import numpy as np
from abc import ABC, abstractmethod

from scipy.spatial.distance import pdist, cdist, squareform

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
6
7
from time import time

8
class Kernel(ABC):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
9
10
11
12
    def __init__(self):
        self._theta = None

    @abstractmethod
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
13
    def __call__(self):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
14
15
        pass

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
16
17
18
    @abstractmethod
    def kernel(self):
        pass
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
19

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
20
21
22
    @abstractmethod
    def kernel_vector(self):
        pass
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

    @abstractmethod
    def kernel_jacobian(self):
        pass

    @abstractmethod
    def kernel_hyperparameter_gradient(self):
        pass

    @property
    def theta(self):
        return self._theta

    @theta.setter
    def theta(self, theta):
        self._theta = theta
    
    def numerical_jacobian(self,x,y, dx=1.e-5):
        if np.ndim(y) == 1:
            y = y.reshape((1,-1))
        nx = len(x)
        ny = y.shape[0]
        f_jac = np.zeros((ny,nx))
        for i in range(nx):
            x_up = np.copy(x)
            x_down = np.copy(x)
            x_up[i] += 0.5*dx
            x_down[i] -= 0.5*dx
            
            f_up = self.kernel(x_up,y)
            f_down = self.kernel(x_down,y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
54

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
55
56
57
58
            f_jac[:,i] = (f_up - f_down)/dx
        return f_jac

    def numerical_hyperparameter_gradient(self,X, dx=1.e-5):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
59
60
61
        """Calculates the numerical derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
62
63
64
65
66
67
68
69
70
71
72
        N_data = X.shape[0]
        theta = np.copy(self.theta)
        N_hyper = len(theta)
        dK_dTheta = np.zeros((N_hyper, N_data, N_data))
        for i in range(N_hyper):
            theta_up = np.copy(theta)
            theta_down = np.copy(theta)
            theta_up[i] += 0.5*dx
            theta_down[i] -= 0.5*dx
            
            self.theta = theta_up
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
73
            K_up = self(X, eval_gradient=False)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
74
            self.theta = theta_down
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
75
            K_down = self(X, eval_gradient=False)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
76

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
77
            dK_dTheta[i,:,:] = (K_up - K_down)/dx
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
78
        self.theta = theta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
79
80
81
        return dK_dTheta


82
class GaussKernel(Kernel):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
83
84
    def __init__(self, amplitude=100.0, amplitude_bounds=(1e0, 1e5),
                 length_scale=10.0, length_scale_bounds=(1e0, 1e3),
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
85
                 noise=1e-5, noise_bounds=None,
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
86
87
                 eta=1, eta_bounds=(0.1,10),
                 Nsplit_eta=None):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
88
89
        self.amplitude = amplitude
        self.length_scale = length_scale
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
90

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
91
92
93
        self.amplitude_bounds = amplitude_bounds
        self.length_scale_bounds = length_scale_bounds

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
94
95
        self.theta_bounds = np.array([amplitude_bounds, length_scale_bounds])

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
96
97
98
99
100
101
102
103
        self.amplitude = amplitude
        self.length_scale = length_scale
        self.noise = noise
        self.eta = eta
        self.Nsplit_eta = Nsplit_eta

        self.amplitude_bounds = amplitude_bounds
        self.length_scale_bounds = length_scale_bounds
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
104
        self.noise_bounds = (noise,noise)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
105
106
107
108
109
        if self.Nsplit_eta is None:
            self.eta_bounds = (eta,eta)
        else:
            self.eta_bounds = eta_bounds

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
110
        self.theta_bounds = np.log(np.array([amplitude_bounds, length_scale_bounds, self.noise_bounds, self.eta_bounds]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
111

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
112
113
114
115
116
117
118
119
120
121
122
123
    def __call__(self, X, eval_gradient=False):
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        d = cdist(X / self.length_scale,
                  X / self.length_scale, metric='sqeuclidean')
        K = self.amplitude * np.exp(-0.5 * d)
        
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
124
125
126
127
128
129
130
131
132
133
134
135

    def __call__(self, X, eval_gradient=False):
        K = self.kernel(X, with_noise=True)
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K

    def kernel(self, X, Y=None, with_noise=False):
        if with_noise:
            assert Y is None
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
136
137
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
138
139
140
        if Y is None:
            Y = X
        elif np.ndim(Y) == 1:
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
141
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
142
143
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
144
145
        d = cdist(X / self.length_scale,
                  Y / self.length_scale, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
146
147
148
149
        if with_noise:
            K = self.amplitude * (np.exp(-0.5 * d) + self.noise*np.eye(X.shape[0]))
        else:
            K = self.amplitude * np.exp(-0.5 * d)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
150
151
        return K

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
152
153
154
155
    def kernel_value(self, x,y):
        K = self.kernel(x,y)
        return np.asscalar(K)
    
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
156
    def kernel_vector(self, x,Y):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
157
        K = self.kernel(x,Y).reshape(-1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
158
159
160
        return K

    def kernel_jacobian(self, X,Y):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
161
162
163
164
        """ Jacobian of the kernel with respect to X
        """
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
165
166
167
168
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        if np.ndim(Y) == 1:
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
169
170
171
172
173
        d = cdist(X / self.length_scale,
                  Y / self.length_scale, metric='sqeuclidean')
        dK_dd = -self.amplitude * 1/(2*self.length_scale**2) * np.exp(-0.5 * d)
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
174
175
        dd_df = 2*(X - Y)

176
177
        #dk_df = np.multiply(dK_dd, dd_df)
        dk_df = np.multiply(dK_dd.T, dd_df)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
178
179
        return dk_df

180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
    def kernel_hessian(self, X,Y):
        """ Jacobian of the kernel with respect to X
        """
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        if np.ndim(Y) == 1:
            Y = Y.reshape((1,-1))
        d = cdist(X / self.length_scale,
                  Y / self.length_scale, metric='sqeuclidean')
        u = 1/(2*self.length_scale**2)
        K = self.amplitude * np.exp(-0.5 * d)
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)

        N_X, Nf = X.shape
        N_Y = Y.shape[0]

        X_rep = np.tile(X, (N_Y,1,1)).swapaxes(0,1)
        Y_rep = np.tile(Y, (N_X,1,1))

        dd_dX = -2*(Y_rep-X_rep)  # shape: (N_X, N_Y, Nf)
        dd_dY = -dd_dX
        dd_dX_dd_dY = np.einsum('nmi,nmj->nmij',dd_dX, dd_dY)  # shape: (N_X, N_Y, Nf, Nf)
        if self.Nsplit_eta is not None:
            diag = np.diag(self.Nsplit_eta*[1]+(Nf-self.Nsplit_eta)*[self.eta**2])
        else:
            diag = np.identity(Nf)
        d2d_dXdY = np.tile(-2*diag, (N_X, N_Y,1,1))  # shape: (N_X, N_Y, Nf, Nf)
        hess = -u*K.reshape(N_X,N_Y,1,1) * (u*dd_dX_dd_dY - d2d_dXdY)

        return hess

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
214
215
    @property
    def theta(self):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
216
217
        """Returns the log-transformed hyperparameters of the kernel.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
218
        self._theta = np.array([self.amplitude, self.length_scale, self.noise, self.eta]) 
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
219
        return np.log(self._theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
220
221
222

    @theta.setter
    def theta(self, theta):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
223
224
225
226
227
        """Sets the hyperparameters of the kernel.

        theta: log-transformed hyperparameters
        """
        self._theta = np.exp(theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
228
229
        self.amplitude = self._theta[0]
        self.length_scale = self._theta[1]
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
230
231
232
233
234
235
236
237
238
239
240
241
        self.noise = self._theta[2]
        self.eta = self._theta[3]

    def apply_eta(self, X):
        Xeta = np.copy(X)
        if self.Nsplit_eta is not None:
            if np.ndim(X) == 1:
                Xeta[self.Nsplit_eta:] *= self.eta
            else:
                Xeta[:,self.Nsplit_eta:] *= self.eta
        return Xeta

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
242
243
244
    def dK_da(self, X):
        d = cdist(X / self.length_scale,
                  X / self.length_scale, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
245
        dK_da = self.amplitude * (np.exp(-0.5 * d) + self.noise*np.eye(X.shape[0]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
246
247
248
249
250
        return dK_da
        
    def dK_dl(self, X):
        d = cdist(X / self.length_scale,
                  X / self.length_scale, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
251
        dK_dl = self.amplitude * d * np.exp(-0.5 * d)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
252
253
        return dK_dl

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
    def dK_dn(self, X):
        dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
        return dK_dn

    def dK_deta(self, X, dx=1e-5):
        N_data = X.shape[0]
        theta = np.copy(self.theta)
        dK_deta = np.zeros((N_data, N_data))

        theta_up = np.copy(theta)
        theta_down = np.copy(theta)
        theta_up[-1] += 0.5*dx
        theta_down[-1] -= 0.5*dx
        
        self.theta = theta_up
        K_up = self(X, eval_gradient=False)
        self.theta = theta_down
        K_down = self(X, eval_gradient=False)
        dK_dTheta = (K_up - K_down)/dx

        self.theta = theta
        return dK_dTheta

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
277
    def kernel_hyperparameter_gradient(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
278
279
280
        """Calculates the derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
281
282
283
284
        dK_deta = self.dK_deta(X)
        X = self.apply_eta(X)
        return np.array([self.dK_da(X), self.dK_dl(X), self.dK_dn(X), dK_deta])

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
285
286
    

287
class DoubleGaussKernel(Kernel):
288
    def __init__(self, amplitude=100., amplitude_bounds=(1e1,1e5),
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
289
290
                 length_scale1=10.0, length_scale1_bounds=(1e0, 1e3),
                 length_scale2=10.0, length_scale2_bounds=(1e0, 1e3),
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
291
292
                 weight=0.01, weight_bounds=None,
                 noise=1e-5, noise_bounds=None,
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
293
294
                 eta=1, eta_bounds=(0.1,10),
                 Nsplit_eta=None):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
295
296
297
298
        self.amplitude = amplitude
        self.length_scale1 = length_scale1
        self.length_scale2 = length_scale2
        self.weight = weight
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
299
        self.noise = noise
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
300
301
        self.eta = eta
        self.Nsplit_eta = Nsplit_eta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
302
303
304
305

        self.amplitude_bounds = amplitude_bounds
        self.length_scale1_bounds = length_scale1_bounds
        self.length_scale2_bounds = length_scale2_bounds
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
306
307
        self.weight_bounds = (weight, weight)
        self.noise_bounds = (noise, noise)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
308
309
310
311
        if self.Nsplit_eta is None:
            self.eta_bounds = (eta,eta)
        else:
            self.eta_bounds = eta_bounds
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
312

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
313
        self.theta_bounds = np.log(np.array([amplitude_bounds, length_scale1_bounds, length_scale2_bounds, self.weight_bounds, self.noise_bounds, self.eta_bounds]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
314

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
315
    def __call__(self, X, eval_gradient=False):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
316
        K = self.kernel(X, with_noise=True)
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
317
318
319
320
321
        if eval_gradient:
            K_gradient = self.kernel_hyperparameter_gradient(X)
            return K, K_gradient
        else:
            return K
322

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
323
324
325
    def kernel(self, X, Y=None, with_noise=False):
        if with_noise:
            assert Y is None
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
326
327
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
328
329
330
        if Y is None:
            Y = X
        elif np.ndim(Y) == 1:
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
331
            Y = Y.reshape((1,-1))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
332
333
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
334
335
336
337
        d1 = cdist(X / self.length_scale1,
                  Y / self.length_scale1, metric='sqeuclidean')
        d2 = cdist(X / self.length_scale2,
                  Y / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
338
339
340
341
        if with_noise:
            K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2) + self.noise*np.eye(X.shape[0]))
        else:
            K = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
342
343
        return K

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
344
345
346
347
    def kernel_value(self, x,y):
        K = self.kernel(x,y)
        return np.asscalar(K)
    
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
348
    def kernel_vector(self, x,Y):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
349
        K = self.kernel(x,Y).reshape(-1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
350
351
352
353
354
        return K

    def kernel_jacobian(self, X,Y):
        """ Jacobian of the kernel with respect to X
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
355
356
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
357
358
359
360
361
362
363
364
365
366
367
        if np.ndim(X) == 1:
            X = X.reshape((1,-1))
        if np.ndim(Y) == 1:
            Y = Y.reshape((1,-1))
        d1 = cdist(X / self.length_scale1,
                   Y / self.length_scale1, metric='sqeuclidean')
        d2 = cdist(X / self.length_scale2,
                   Y / self.length_scale2, metric='sqeuclidean')
        dK1_dd1 = -1/(2*self.length_scale1**2) * np.exp(-0.5 * d1)
        dK2_dd2 = -1/(2*self.length_scale2**2) * np.exp(-0.5 * d2)
        dK_dd = self.amplitude * (dK1_dd1 + self.weight*dK2_dd2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
368
369
        X = self.apply_eta(X)
        Y = self.apply_eta(Y)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
370
371
372
373
374
375
376
        dd_df = 2*(X - Y)

        dk_df = np.multiply(dK_dd.T, dd_df)
        return dk_df

    @property
    def theta(self):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
377
378
        """Returns the log-transformed hyperparameters of the kernel.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
379
        self._theta = np.array([self.amplitude, self.length_scale1, self.length_scale2, self.weight, self.noise, self.eta])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
380
        return np.log(self._theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
381
382
383

    @theta.setter
    def theta(self, theta):
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
384
385
386
387
388
        """Sets the hyperparameters of the kernel.

        theta: log-transformed hyperparameters
        """
        self._theta = np.exp(theta)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
389
390
391
392
        self.amplitude = self._theta[0]
        self.length_scale1 = self._theta[1]
        self.length_scale2 = self._theta[2]
        self.weight = self._theta[3]
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
393
        self.noise = self._theta[4]
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
394
395
396
397
398
399
400
401
402
403
        self.eta = self._theta[5]

    def apply_eta(self, X):
        Xeta = np.copy(X)
        if self.Nsplit_eta is not None:
            if np.ndim(X) == 1:
                Xeta[self.Nsplit_eta:] *= self.eta
            else:
                Xeta[:,self.Nsplit_eta:] *= self.eta
        return Xeta
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
404
405
406
407
408
409

    def dK_da(self, X):
        d1 = cdist(X / self.length_scale1,
                   X / self.length_scale1, metric='sqeuclidean')
        d2 = cdist(X / self.length_scale2,
                   X / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
410
        dK_da = self.amplitude * (np.exp(-0.5 * d1) + self.weight*np.exp(-0.5 * d2) + self.noise*np.eye(X.shape[0]))
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
411
412
413
414
415
        return dK_da
        
    def dK_dl1(self, X):
        d1 = cdist(X / self.length_scale1,
                   X / self.length_scale1, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
416
        dK_dl1 = self.amplitude*d1 * np.exp(-0.5 * d1)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
417
418
419
420
421
        return dK_dl1

    def dK_dl2(self, X):
        d2 = cdist(X / self.length_scale2,
                   X / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
422
        dK_dl2 = self.amplitude*self.weight*d2 * np.exp(-0.5 * d2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
423
424
425
426
427
        return dK_dl2

    def dK_dw(self, X):
        d2 = cdist(X / self.length_scale2,
                   X / self.length_scale2, metric='sqeuclidean')
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
428
        dK_dl2 = self.amplitude*self.weight*np.exp(-0.5 * d2)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
429
430
        return dK_dl2

Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
431
    def dK_dn(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
432
        dK_dn = self.amplitude * self.noise * np.eye(X.shape[0])
Malthe Kjær Bisbo's avatar
update    
Malthe Kjær Bisbo committed
433
434
        return dK_dn

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
    def dK_deta(self, X, dx=1e-5):
        N_data = X.shape[0]
        theta = np.copy(self.theta)
        dK_deta = np.zeros((N_data, N_data))

        theta_up = np.copy(theta)
        theta_down = np.copy(theta)
        theta_up[-1] += 0.5*dx
        theta_down[-1] -= 0.5*dx
        
        self.theta = theta_up
        K_up = self(X, eval_gradient=False)
        self.theta = theta_down
        K_down = self(X, eval_gradient=False)
        dK_dTheta = (K_up - K_down)/dx

        self.theta = theta
        return dK_dTheta

Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
454
    def kernel_hyperparameter_gradient(self, X):
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
455
456
457
        """Calculates the derivative of the kernel with respect to the
        log transformed hyperparameters.
        """
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
458
        dK_deta = self.dK_deta(X)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
459
        X = self.apply_eta(X)
Malthe Kjær Bisbo's avatar
Malthe Kjær Bisbo committed
460
        return np.array([self.dK_da(X), self.dK_dl1(X), self.dK_dl2(X), self.dK_dw(X), self.dK_dn(X), dK_deta])