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