1# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
2# Licensed under the BSD 3-clause license (see LICENSE.txt)
3
4from .kern import Kern
5from ...core.parameterization import Param
6from paramz.transformations import Logexp
7import numpy as np
8from ...util.multioutput import index_to_slices
9
10class ODE_UYC(Kern):
11    def __init__(self, input_dim, variance_U=3., variance_Y=1., lengthscale_U=1., lengthscale_Y=1., ubias =1. ,active_dims=None, name='ode_uyc'):
12        assert input_dim ==2, "only defined for 2 input dims"
13        super(ODE_UYC, self).__init__(input_dim, active_dims, name)
14
15        self.variance_Y = Param('variance_Y', variance_Y, Logexp())
16        self.variance_U = Param('variance_U', variance_U, Logexp())
17        self.lengthscale_Y = Param('lengthscale_Y', lengthscale_Y, Logexp())
18        self.lengthscale_U = Param('lengthscale_U', lengthscale_U, Logexp())
19        self.ubias = Param('ubias', ubias, Logexp())
20
21        self.link_parameters(self.variance_Y, self.variance_U, self.lengthscale_Y, self.lengthscale_U, self.ubias)
22
23    def K(self, X, X2=None):
24        # model :   a * dy/dt + b * y = U
25        #lu=sqrt(3)/theta1  ly=1/theta2  theta2= a/b :thetay   sigma2=1/(2ab) :sigmay
26
27        X,slices = X[:,:-1],index_to_slices(X[:,-1])
28        if X2 is None:
29            X2,slices2 = X,slices
30            K = np.zeros((X.shape[0], X.shape[0]))
31        else:
32            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
33            K = np.zeros((X.shape[0], X2.shape[0]))
34
35        #stop
36        #rdist = X[:,0][:,None] - X2[:,0][:,None].T
37        rdist = X - X2.T
38        ly=1/self.lengthscale_Y
39        lu=np.sqrt(3)/self.lengthscale_U
40        #iu=self.input_lengthU  #dimention of U
41        Vu=self.variance_U
42        Vy=self.variance_Y
43        #Vy=ly/2
44        #stop
45
46
47        # kernel for kuu  matern3/2
48        kuu = lambda dist:Vu * (1 + lu* np.abs(dist)) * np.exp(-lu * np.abs(dist)) +self.ubias
49
50        # kernel for kyy
51        k1 = lambda dist:np.exp(-ly*np.abs(dist))*(2*lu+ly)/(lu+ly)**2
52        k2 = lambda dist:(np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2
53        k3 = lambda dist:np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
54        kyy = lambda dist:Vu*Vy*(k1(dist) + k2(dist) + k3(dist))
55
56
57        # cross covariance function
58        kyu3 = lambda dist:np.exp(-lu*dist)/(lu+ly)*(1+lu*(dist+1/(lu+ly)))
59        #kyu3 = lambda dist: 0
60
61        k1cros = lambda dist:np.exp(ly*dist)/(lu-ly) * ( 1- np.exp( (lu-ly)*dist) + lu* ( dist*np.exp( (lu-ly)*dist ) + (1- np.exp( (lu-ly)*dist ) ) /(lu-ly)   )    )
62        #k1cros = lambda dist:0
63
64        k2cros = lambda dist:np.exp(ly*dist)*( 1/(lu+ly) + lu/(lu+ly)**2 )
65        #k2cros = lambda dist:0
66
67        Vyu=np.sqrt(Vy*ly*2)
68
69        # cross covariance kuy
70        kuyp = lambda dist:Vu*Vyu*(kyu3(dist))       #t>0 kuy
71        kuyn = lambda dist:Vu*Vyu*(k1cros(dist)+k2cros(dist))      #t<0 kuy
72        # cross covariance kyu
73        kyup = lambda dist:Vu*Vyu*(k1cros(-dist)+k2cros(-dist))    #t>0 kyu
74        kyun = lambda dist:Vu*Vyu*(kyu3(-dist))       #t<0 kyu
75
76
77        for i, s1 in enumerate(slices):
78            for j, s2 in enumerate(slices2):
79                for ss1 in s1:
80                    for ss2 in s2:
81                        if i==0 and j==0:
82                            K[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
83                        elif i==0 and j==1:
84                            #K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) )   )
85                            K[ss1,ss2]=  np.where(  rdist[ss1,ss2]>0 , kuyp(rdist[ss1,ss2]), kuyn(rdist[ss1,ss2] )   )
86                        elif i==1 and j==1:
87                            K[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
88                        else:
89                            #K[ss1,ss2]= 0
90                            #K[ss1,ss2]= np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) )   )
91                            K[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(rdist[ss1,ss2]), kyun(rdist[ss1,ss2] )   )
92        return K
93
94
95
96    def Kdiag(self, X):
97        """Compute the diagonal of the covariance matrix associated to X."""
98        Kdiag = np.zeros(X.shape[0])
99        ly=1/self.lengthscale_Y
100        lu=np.sqrt(3)/self.lengthscale_U
101
102        Vu = self.variance_U
103        Vy=self.variance_Y
104
105        k1 = (2*lu+ly)/(lu+ly)**2
106        k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2
107        k3 = 1/(lu+ly) + (lu)/(lu+ly)**2
108
109        slices = index_to_slices(X[:,-1])
110
111        for i, ss1 in enumerate(slices):
112            for s1 in ss1:
113                if i==0:
114                    Kdiag[s1]+= self.variance_U + self.ubias
115                elif i==1:
116                    Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
117                else:
118                    raise ValueError("invalid input/output index")
119        #Kdiag[slices[0][0]]+= self.variance_U   #matern32 diag
120        #Kdiag[slices[1][0]]+= self.variance_U*self.variance_Y*(k1+k2+k3)  #  diag
121        return Kdiag
122
123
124    def update_gradients_full(self, dL_dK, X, X2=None):
125        """derivative of the covariance matrix with respect to the parameters."""
126        X,slices = X[:,:-1],index_to_slices(X[:,-1])
127        if X2 is None:
128            X2,slices2 = X,slices
129        else:
130            X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
131        #rdist = X[:,0][:,None] - X2[:,0][:,None].T
132
133        rdist = X - X2.T
134        ly=1/self.lengthscale_Y
135        lu=np.sqrt(3)/self.lengthscale_U
136
137        Vu=self.variance_U
138        Vy=self.variance_Y
139        Vyu = np.sqrt(Vy*ly*2)
140        dVdly = 0.5/np.sqrt(ly)*np.sqrt(2*Vy)
141        dVdVy = 0.5/np.sqrt(Vy)*np.sqrt(2*ly)
142
143        rd=rdist.shape[0]
144        dktheta1 = np.zeros([rd,rd])
145        dktheta2 = np.zeros([rd,rd])
146        dkUdvar = np.zeros([rd,rd])
147        dkYdvar = np.zeros([rd,rd])
148
149        dkdubias = np.zeros([rd,rd])
150
151        # dk dtheta for UU
152        UUdtheta1 = lambda dist: np.exp(-lu* dist)*dist + (-dist)*np.exp(-lu* dist)*(1+lu*dist)
153        UUdtheta2 = lambda dist: 0
154        #UUdvar = lambda dist: (1 + lu*dist)*np.exp(-lu*dist)
155        UUdvar = lambda dist: (1 + lu* np.abs(dist)) * np.exp(-lu * np.abs(dist))
156
157        # dk dtheta for YY
158
159        dk1theta1 = lambda dist: np.exp(-ly*dist)*2*(-lu)/(lu+ly)**3
160
161        dk2theta1 = lambda dist: (1.0)*(
162            np.exp(-lu*dist)*dist*(-ly+2*lu-lu*ly*dist+dist*lu**2)*(ly-lu)**(-2) + np.exp(-lu*dist)*(-2+ly*dist-2*dist*lu)*(ly-lu)**(-2)
163            +np.exp(-dist*lu)*(ly-2*lu+ly*lu*dist-dist*lu**2)*2*(ly-lu)**(-3)
164            +np.exp(-dist*ly)*2*(ly-lu)**(-2)
165            +np.exp(-dist*ly)*2*(2*lu-ly)*(ly-lu)**(-3)
166            )
167
168        dk3theta1 = lambda dist: np.exp(-dist*lu)*(lu+ly)**(-2)*((2*lu+ly+dist*lu**2+lu*ly*dist)*(-dist-2/(lu+ly))+2+2*lu*dist+ly*dist)
169
170        #dktheta1 = lambda dist: self.variance_U*self.variance_Y*(dk1theta1+dk2theta1+dk3theta1)
171
172
173
174
175        dk1theta2 = lambda dist: np.exp(-ly*dist) * ((lu+ly)**(-2)) * (  (-dist)*(2*lu+ly)  +  1  +  (-2)*(2*lu+ly)/(lu+ly)  )
176
177        dk2theta2 =lambda dist:  1*(
178            np.exp(-dist*lu)*(ly-lu)**(-2) * ( 1+lu*dist+(-2)*(ly-2*lu+lu*ly*dist-dist*lu**2)*(ly-lu)**(-1) )
179            +np.exp(-dist*ly)*(ly-lu)**(-2) * ( (-dist)*(2*lu-ly) -1+(2*lu-ly)*(-2)*(ly-lu)**(-1) )
180            )
181
182        dk3theta2 = lambda dist: np.exp(-dist*lu) * (-3*lu-ly-dist*lu**2-lu*ly*dist)/(lu+ly)**3
183
184        #dktheta2 = lambda dist: self.variance_U*self.variance_Y*(dk1theta2 + dk2theta2 +dk3theta2)
185
186        # kyy kernel
187
188        k1 = lambda dist: np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
189        k2 = lambda dist: (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2
190        k3 = lambda dist: np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
191        #dkdvar = k1+k2+k3
192
193
194
195        # cross covariance function
196        kyu3 = lambda dist:np.exp(-lu*dist)/(lu+ly)*(1+lu*(dist+1/(lu+ly)))
197
198        k1cros = lambda dist:np.exp(ly*dist)/(lu-ly) * ( 1- np.exp( (lu-ly)*dist) + lu* ( dist*np.exp( (lu-ly)*dist ) + (1- np.exp( (lu-ly)*dist ) ) /(lu-ly)   )    )
199
200        k2cros = lambda dist:np.exp(ly*dist)*( 1/(lu+ly) + lu/(lu+ly)**2 )
201        # cross covariance kuy
202        kuyp = lambda dist:(kyu3(dist))       #t>0 kuy
203        kuyn = lambda dist:(k1cros(dist)+k2cros(dist))      #t<0 kuy
204        # cross covariance kyu
205        kyup = lambda dist:(k1cros(-dist)+k2cros(-dist))    #t>0 kyu
206        kyun = lambda dist:(kyu3(-dist))       #t<0 kyu
207
208        # dk dtheta for UY
209
210
211        dkyu3dtheta2 = lambda dist: np.exp(-lu*dist) * ( (-1)*(lu+ly)**(-2)*(1+lu*dist+lu*(lu+ly)**(-1)) + (lu+ly)**(-1)*(-lu)*(lu+ly)**(-2) )
212        dkyu3dtheta1 = lambda dist: np.exp(-lu*dist)*(lu+ly)**(-1)* ( (-dist)*(1+dist*lu+lu*(lu+ly)**(-1)) -\
213         (lu+ly)**(-1)*(1+dist*lu+lu*(lu+ly)**(-1)) +dist+(lu+ly)**(-1)-lu*(lu+ly)**(-2) )
214
215        dkcros2dtheta1 = lambda dist: np.exp(ly*dist)* ( -(ly+lu)**(-2) + (ly+lu)**(-2) + (-2)*lu*(lu+ly)**(-3)  )
216        dkcros2dtheta2 = lambda dist: np.exp(ly*dist)*dist* ( (ly+lu)**(-1) + lu*(lu+ly)**(-2) ) + \
217                                      np.exp(ly*dist)*( -(lu+ly)**(-2) + lu*(-2)*(lu+ly)**(-3)  )
218
219        dkcros1dtheta1 = lambda dist: np.exp(ly*dist)*(     -(lu-ly)**(-2)*(  1-np.exp((lu-ly)*dist) + lu*dist*np.exp((lu-ly)*dist)+ \
220          lu*(1-np.exp((lu-ly)*dist))/(lu-ly)  )  +  (lu-ly)**(-1)*(  -np.exp( (lu-ly)*dist )*dist + dist*np.exp( (lu-ly)*dist)+\
221          lu*dist**2*np.exp((lu-ly)*dist)+(1-np.exp((lu-ly)*dist))/(lu-ly) - lu*np.exp((lu-ly)*dist)*dist/(lu-ly) -\
222          lu*(1-np.exp((lu-ly)*dist))/(lu-ly)**2  )   )
223
224        dkcros1dtheta2 = lambda t: np.exp(ly*t)*t/(lu-ly)*( 1-np.exp((lu-ly)*t) +lu*t*np.exp((lu-ly)*t)+\
225            lu*(1-np.exp((lu-ly)*t))/(lu-ly)  )+\
226            np.exp(ly*t)/(lu-ly)**2* ( 1-np.exp((lu-ly)*t) +lu*t*np.exp((lu-ly)*t) + lu*( 1-np.exp((lu-ly)*t) )/(lu-ly)  )+\
227            np.exp(ly*t)/(lu-ly)*( np.exp((lu-ly)*t)*t -lu*t*t*np.exp((lu-ly)*t) +lu*t*np.exp((lu-ly)*t)/(lu-ly)+\
228            lu*( 1-np.exp((lu-ly)*t) )/(lu-ly)**2 )
229
230        dkuypdtheta1 = lambda dist:(dkyu3dtheta1(dist))       #t>0 kuy
231        dkuyndtheta1 = lambda dist:(dkcros1dtheta1(dist)+dkcros2dtheta1(dist))      #t<0 kuy
232        # cross covariance kyu
233        dkyupdtheta1 = lambda dist:(dkcros1dtheta1(-dist)+dkcros2dtheta1(-dist))    #t>0 kyu
234        dkyundtheta1 = lambda dist:(dkyu3dtheta1(-dist))       #t<0 kyu
235
236        dkuypdtheta2 = lambda dist:(dkyu3dtheta2(dist))       #t>0 kuy
237        dkuyndtheta2 = lambda dist:(dkcros1dtheta2(dist)+dkcros2dtheta2(dist))      #t<0 kuy
238        # cross covariance kyu
239        dkyupdtheta2 = lambda dist:(dkcros1dtheta2(-dist)+dkcros2dtheta2(-dist))    #t>0 kyu
240        dkyundtheta2 = lambda dist:(dkyu3dtheta2(-dist))       #t<0 kyu
241
242
243        for i, s1 in enumerate(slices):
244            for j, s2 in enumerate(slices2):
245                for ss1 in s1:
246                    for ss2 in s2:
247                        if i==0 and j==0:
248                            #target[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
249                            dktheta1[ss1,ss2] = Vu*UUdtheta1(np.abs(rdist[ss1,ss2]))
250                            dktheta2[ss1,ss2] = 0
251                            dkUdvar[ss1,ss2] = UUdvar(np.abs(rdist[ss1,ss2]))
252                            dkYdvar[ss1,ss2] = 0
253                            dkdubias[ss1,ss2] = 1
254                        elif i==0 and j==1:
255                            ########target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
256                            #np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
257                            #dktheta1[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , self.variance_U*self.variance_Y*dkcrtheta1(np.abs(rdist[ss1,ss2])) ,self.variance_U*self.variance_Y*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2])))    )
258                            #dktheta2[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , self.variance_U*self.variance_Y*dkcrtheta2(np.abs(rdist[ss1,ss2])) ,self.variance_U*self.variance_Y*(dk1theta2(np.abs(rdist[ss1,ss2]))+dk2theta2(np.abs(rdist[ss1,ss2])))    )
259                            dktheta1[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vu*Vyu*dkuypdtheta1(rdist[ss1,ss2]),Vu*Vyu*dkuyndtheta1(rdist[ss1,ss2]) )
260                            dkUdvar[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vyu*kuyp(rdist[ss1,ss2]), Vyu* kuyn(rdist[ss1,ss2])  )
261                            dktheta2[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vu*Vyu*dkuypdtheta2(rdist[ss1,ss2])+Vu*dVdly*kuyp(rdist[ss1,ss2]),Vu*Vyu*dkuyndtheta2(rdist[ss1,ss2])+Vu*dVdly*kuyn(rdist[ss1,ss2]) )
262                            dkYdvar[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vu*dVdVy*kuyp(rdist[ss1,ss2]), Vu*dVdVy* kuyn(rdist[ss1,ss2])  )
263                            dkdubias[ss1,ss2] = 0
264                        elif i==1 and j==1:
265                            #target[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
266                            dktheta1[ss1,ss2] = self.variance_U*self.variance_Y*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))+dk3theta1(np.abs(rdist[ss1,ss2])))
267                            dktheta2[ss1,ss2] = self.variance_U*self.variance_Y*(dk1theta2(np.abs(rdist[ss1,ss2])) + dk2theta2(np.abs(rdist[ss1,ss2])) +dk3theta2(np.abs(rdist[ss1,ss2])))
268                            dkUdvar[ss1,ss2] = self.variance_Y*(k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2])) )
269                            dkYdvar[ss1,ss2] = self.variance_U*(k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2])) )
270                            dkdubias[ss1,ss2] = 0
271                        else:
272                            #######target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[s1[0],s2[0]]) )   )
273                            #dktheta1[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 ,self.variance_U*self.variance_Y*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))) , self.variance_U*self.variance_Y*dkcrtheta1(np.abs(rdist[ss1,ss2])) )
274                            #dktheta2[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 ,self.variance_U*self.variance_Y*(dk1theta2(np.abs(rdist[ss1,ss2]))+dk2theta2(np.abs(rdist[ss1,ss2]))) , self.variance_U*self.variance_Y*dkcrtheta2(np.abs(rdist[ss1,ss2])) )
275                            dktheta1[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vu*Vyu*dkyupdtheta1(rdist[ss1,ss2]),Vu*Vyu*dkyundtheta1(rdist[ss1,ss2])  )
276                            dkUdvar[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vyu*kyup(rdist[ss1,ss2]),Vyu*kyun(rdist[ss1,ss2]))
277                            dktheta2[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vu*Vyu*dkyupdtheta2(rdist[ss1,ss2])+Vu*dVdly*kyup(rdist[ss1,ss2]),Vu*Vyu*dkyundtheta2(rdist[ss1,ss2])+Vu*dVdly*kyun(rdist[ss1,ss2])  )
278                            dkYdvar[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , Vu*dVdVy*kyup(rdist[ss1,ss2]), Vu*dVdVy*kyun(rdist[ss1,ss2]))
279                            dkdubias[ss1,ss2] = 0
280        #stop
281        self.variance_U.gradient = np.sum(dkUdvar * dL_dK)     # Vu
282
283        self.variance_Y.gradient = np.sum(dkYdvar * dL_dK)     # Vy
284
285        self.lengthscale_U.gradient = np.sum(dktheta1*(-np.sqrt(3)*self.lengthscale_U**(-2))* dL_dK)     #lu
286
287        self.lengthscale_Y.gradient = np.sum(dktheta2*(-self.lengthscale_Y**(-2)) * dL_dK)              #ly
288
289        self.ubias.gradient = np.sum(dkdubias * dL_dK)
290
291