1# -*- coding: utf-8 -*-
2# Copyright (c) 2015, Alex Grigorevskiy, Arno Solin
3# Licensed under the BSD 3-clause license (see LICENSE.txt)
4"""
5Classes in this module enhance Matern covariance functions with the
6Stochastic Differential Equation (SDE) functionality.
7"""
8from .stationary import Matern32
9from .stationary import Matern52
10import numpy as np
11
12class sde_Matern32(Matern32):
13    """
14
15    Class provide extra functionality to transfer this covariance function into
16    SDE forrm.
17
18    Matern 3/2 kernel:
19
20    .. math::
21
22       k(r) = \sigma^2 (1 + \sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\  \text{ where  } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
23
24    """
25    def sde_update_gradient_full(self, gradients):
26        """
27        Update gradient in the order in which parameters are represented in the
28        kernel
29        """
30
31        self.variance.gradient = gradients[0]
32        self.lengthscale.gradient = gradients[1]
33
34    def sde(self):
35        """
36        Return the state space representation of the covariance.
37        """
38
39        variance = float(self.variance.values)
40        lengthscale = float(self.lengthscale.values)
41
42        foo  = np.sqrt(3.)/lengthscale
43        F    = np.array(((0, 1.0), (-foo**2, -2*foo)))
44        L    = np.array(( (0,), (1.0,) ))
45        Qc   = np.array(((12.*np.sqrt(3) / lengthscale**3 * variance,),))
46        H    = np.array(((1.0, 0),))
47        Pinf = np.array(((variance, 0.0), (0.0, 3.*variance/(lengthscale**2))))
48        P0 = Pinf.copy()
49
50        # Allocate space for the derivatives
51        dF    = np.empty([F.shape[0],F.shape[1],2])
52        dQc   = np.empty([Qc.shape[0],Qc.shape[1],2])
53        dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])
54        # The partial derivatives
55        dFvariance       = np.zeros((2,2))
56        dFlengthscale    = np.array(((0,0), (6./lengthscale**3,2*np.sqrt(3)/lengthscale**2)))
57        dQcvariance      = np.array((12.*np.sqrt(3)/lengthscale**3))
58        dQclengthscale   = np.array((-3*12*np.sqrt(3)/lengthscale**4*variance))
59        dPinfvariance    = np.array(((1,0),(0,3./lengthscale**2)))
60        dPinflengthscale = np.array(((0,0), (0,-6*variance/lengthscale**3)))
61        # Combine the derivatives
62        dF[:,:,0]    = dFvariance
63        dF[:,:,1]    = dFlengthscale
64        dQc[:,:,0]   = dQcvariance
65        dQc[:,:,1]   = dQclengthscale
66        dPinf[:,:,0] = dPinfvariance
67        dPinf[:,:,1] = dPinflengthscale
68        dP0 = dPinf.copy()
69
70        return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)
71
72class sde_Matern52(Matern52):
73    """
74
75    Class provide extra functionality to transfer this covariance function into
76    SDE forrm.
77
78    Matern 5/2 kernel:
79
80    .. math::
81
82       k(r) = \sigma^2 (1 + \sqrt{5} r + \frac{5}{3}r^2) \exp(- \sqrt{5} r) \\ \\ \\ \\  \text{ where  } r = \sqrt{\sum_{i=1}^{input dim} \frac{(x_i-y_i)^2}{\ell_i^2} }
83
84    """
85    def sde_update_gradient_full(self, gradients):
86        """
87        Update gradient in the order in which parameters are represented in the
88        kernel
89        """
90
91        self.variance.gradient = gradients[0]
92        self.lengthscale.gradient = gradients[1]
93
94    def sde(self):
95        """
96        Return the state space representation of the covariance.
97        """
98
99        variance = float(self.variance.values)
100        lengthscale = float(self.lengthscale.values)
101
102        lamda = np.sqrt(5.0)/lengthscale
103        kappa = 5.0/3.0*variance/lengthscale**2
104
105        F = np.array(((0, 1,0), (0, 0, 1), (-lamda**3, -3.0*lamda**2, -3*lamda)))
106        L = np.array(((0,),(0,),(1,)))
107        Qc = np.array((((variance*400.0*np.sqrt(5.0)/3.0/lengthscale**5),),))
108        H = np.array(((1,0,0),))
109
110        Pinf = np.array(((variance,0,-kappa), (0, kappa, 0), (-kappa, 0, 25.0*variance/lengthscale**4)))
111        P0 = Pinf.copy()
112        # Allocate space for the derivatives
113        dF = np.empty((3,3,2))
114        dQc = np.empty((1,1,2))
115        dPinf = np.empty((3,3,2))
116
117         # The partial derivatives
118        dFvariance = np.zeros((3,3))
119        dFlengthscale = np.array(((0,0,0),(0,0,0),(15.0*np.sqrt(5.0)/lengthscale**4,
120                                   30.0/lengthscale**3, 3*np.sqrt(5.0)/lengthscale**2)))
121        dQcvariance = np.array((((400*np.sqrt(5)/3/lengthscale**5,),)))
122        dQclengthscale = np.array((((-variance*2000*np.sqrt(5)/3/lengthscale**6,),)))
123
124        dPinf_variance = Pinf/variance
125        kappa2 = -2.0*kappa/lengthscale
126        dPinf_lengthscale = np.array(((0,0,-kappa2),(0,kappa2,0),(-kappa2,
127                                    0,-100*variance/lengthscale**5)))
128        # Combine the derivatives
129        dF[:,:,0] = dFvariance
130        dF[:,:,1] = dFlengthscale
131        dQc[:,:,0] = dQcvariance
132        dQc[:,:,1] = dQclengthscale
133        dPinf[:,:,0] = dPinf_variance
134        dPinf[:,:,1] = dPinf_lengthscale
135        dP0 = dPinf.copy()
136
137        return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)