mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-18 13:55:14 +02:00
Merge branch 'devel' of github.com:SheffieldML/GPy into devel
This commit is contained in:
mzwiessele
commit
f26aa7da7d
6 changed files with 730 additions and 3 deletions
|
|
@ -32,7 +32,7 @@ def print_out(len_maxiters, fnow, current_grad, beta, iteration):
|
||||||
sys.stdout.flush()
|
sys.stdout.flush()
|
||||||
|
|
||||||
def exponents(fnow, current_grad):
|
def exponents(fnow, current_grad):
|
||||||
exps = [np.abs(fnow), current_grad]
|
exps = [np.abs(np.float(fnow)), current_grad]
|
||||||
return np.sign(exps) * np.log10(exps).astype(int)
|
return np.sign(exps) * np.log10(exps).astype(int)
|
||||||
|
|
||||||
def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=np.inf, display=True, xtol=None, ftol=None, gtol=None):
|
def SCG(f, gradf, x, optargs=(), maxiters=500, max_f_eval=np.inf, display=True, xtol=None, ftol=None, gtol=None):
|
||||||
|
|
|
||||||
|
|
@ -10,9 +10,11 @@ from _src.independent_outputs import IndependentOutputs, Hierarchical
|
||||||
from _src.coregionalize import Coregionalize
|
from _src.coregionalize import Coregionalize
|
||||||
from _src.ssrbf import SSRBF # TODO: ZD: did you remove this?
|
from _src.ssrbf import SSRBF # TODO: ZD: did you remove this?
|
||||||
from _src.ODE_UY import ODE_UY
|
from _src.ODE_UY import ODE_UY
|
||||||
|
from _src.ODE_UYC import ODE_UYC
|
||||||
|
from _src.ODE_st import ODE_st
|
||||||
|
from _src.ODE_t import ODE_t
|
||||||
from _src.poly import Poly
|
from _src.poly import Poly
|
||||||
#from _src.ODE_UYC import ODE_UYC ADD THIS FILE TO THE REPO!!
|
|
||||||
#from _src.ODE_st import ODE_st
|
|
||||||
# TODO: put this in an init file somewhere
|
# TODO: put this in an init file somewhere
|
||||||
#I'm commenting this out because the files were not added. JH. Remember to add the files before commiting
|
#I'm commenting this out because the files were not added. JH. Remember to add the files before commiting
|
||||||
try:
|
try:
|
||||||
|
|
|
||||||
290
GPy/kern/_src/ODE_UYC.py
Normal file
290
GPy/kern/_src/ODE_UYC.py
Normal file
|
|
@ -0,0 +1,290 @@
|
||||||
|
# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
|
||||||
|
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||||
|
|
||||||
|
from kern import Kern
|
||||||
|
from ...core.parameterization import Param
|
||||||
|
from ...core.parameterization.transformations import Logexp
|
||||||
|
import numpy as np
|
||||||
|
from independent_outputs import index_to_slices
|
||||||
|
|
||||||
|
class ODE_UYC(Kern):
|
||||||
|
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'):
|
||||||
|
assert input_dim ==2, "only defined for 2 input dims"
|
||||||
|
super(ODE_UYC, self).__init__(input_dim, active_dims, name)
|
||||||
|
|
||||||
|
self.variance_Y = Param('variance_Y', variance_Y, Logexp())
|
||||||
|
self.variance_U = Param('variance_U', variance_U, Logexp())
|
||||||
|
self.lengthscale_Y = Param('lengthscale_Y', lengthscale_Y, Logexp())
|
||||||
|
self.lengthscale_U = Param('lengthscale_U', lengthscale_U, Logexp())
|
||||||
|
self.ubias = Param('ubias', ubias, Logexp())
|
||||||
|
|
||||||
|
self.add_parameters(self.variance_Y, self.variance_U, self.lengthscale_Y, self.lengthscale_U, self.ubias)
|
||||||
|
|
||||||
|
def K(self, X, X2=None):
|
||||||
|
# model : a * dy/dt + b * y = U
|
||||||
|
#lu=sqrt(3)/theta1 ly=1/theta2 theta2= a/b :thetay sigma2=1/(2ab) :sigmay
|
||||||
|
|
||||||
|
X,slices = X[:,:-1],index_to_slices(X[:,-1])
|
||||||
|
if X2 is None:
|
||||||
|
X2,slices2 = X,slices
|
||||||
|
K = np.zeros((X.shape[0], X.shape[0]))
|
||||||
|
else:
|
||||||
|
X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
|
||||||
|
K = np.zeros((X.shape[0], X2.shape[0]))
|
||||||
|
|
||||||
|
#stop
|
||||||
|
#rdist = X[:,0][:,None] - X2[:,0][:,None].T
|
||||||
|
rdist = X - X2.T
|
||||||
|
ly=1/self.lengthscale_Y
|
||||||
|
lu=np.sqrt(3)/self.lengthscale_U
|
||||||
|
#iu=self.input_lengthU #dimention of U
|
||||||
|
Vu=self.variance_U
|
||||||
|
Vy=self.variance_Y
|
||||||
|
#Vy=ly/2
|
||||||
|
#stop
|
||||||
|
|
||||||
|
|
||||||
|
# kernel for kuu matern3/2
|
||||||
|
kuu = lambda dist:Vu * (1 + lu* np.abs(dist)) * np.exp(-lu * np.abs(dist)) +self.ubias
|
||||||
|
|
||||||
|
# kernel for kyy
|
||||||
|
k1 = lambda dist:np.exp(-ly*np.abs(dist))*(2*lu+ly)/(lu+ly)**2
|
||||||
|
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
|
||||||
|
k3 = lambda dist:np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
|
||||||
|
kyy = lambda dist:Vu*Vy*(k1(dist) + k2(dist) + k3(dist))
|
||||||
|
|
||||||
|
|
||||||
|
# cross covariance function
|
||||||
|
kyu3 = lambda dist:np.exp(-lu*dist)/(lu+ly)*(1+lu*(dist+1/(lu+ly)))
|
||||||
|
#kyu3 = lambda dist: 0
|
||||||
|
|
||||||
|
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) ) )
|
||||||
|
#k1cros = lambda dist:0
|
||||||
|
|
||||||
|
k2cros = lambda dist:np.exp(ly*dist)*( 1/(lu+ly) + lu/(lu+ly)**2 )
|
||||||
|
#k2cros = lambda dist:0
|
||||||
|
|
||||||
|
Vyu=np.sqrt(Vy*ly*2)
|
||||||
|
|
||||||
|
# cross covariance kuy
|
||||||
|
kuyp = lambda dist:Vu*Vyu*(kyu3(dist)) #t>0 kuy
|
||||||
|
kuyn = lambda dist:Vu*Vyu*(k1cros(dist)+k2cros(dist)) #t<0 kuy
|
||||||
|
# cross covariance kyu
|
||||||
|
kyup = lambda dist:Vu*Vyu*(k1cros(-dist)+k2cros(-dist)) #t>0 kyu
|
||||||
|
kyun = lambda dist:Vu*Vyu*(kyu3(-dist)) #t<0 kyu
|
||||||
|
|
||||||
|
|
||||||
|
for i, s1 in enumerate(slices):
|
||||||
|
for j, s2 in enumerate(slices2):
|
||||||
|
for ss1 in s1:
|
||||||
|
for ss2 in s2:
|
||||||
|
if i==0 and j==0:
|
||||||
|
K[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
|
||||||
|
elif i==0 and j==1:
|
||||||
|
#K[ss1,ss2]= np.where( rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) ) )
|
||||||
|
K[ss1,ss2]= np.where( rdist[ss1,ss2]>0 , kuyp(rdist[ss1,ss2]), kuyn(rdist[ss1,ss2] ) )
|
||||||
|
elif i==1 and j==1:
|
||||||
|
K[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
|
||||||
|
else:
|
||||||
|
#K[ss1,ss2]= 0
|
||||||
|
#K[ss1,ss2]= np.where( rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) ) )
|
||||||
|
K[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , kyup(rdist[ss1,ss2]), kyun(rdist[ss1,ss2] ) )
|
||||||
|
return K
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
def Kdiag(self, X):
|
||||||
|
"""Compute the diagonal of the covariance matrix associated to X."""
|
||||||
|
Kdiag = np.zeros(X.shape[0])
|
||||||
|
ly=1/self.lengthscale_Y
|
||||||
|
lu=np.sqrt(3)/self.lengthscale_U
|
||||||
|
|
||||||
|
Vu = self.variance_U
|
||||||
|
Vy=self.variance_Y
|
||||||
|
|
||||||
|
k1 = (2*lu+ly)/(lu+ly)**2
|
||||||
|
k2 = (ly-2*lu + 2*lu-ly ) / (ly-lu)**2
|
||||||
|
k3 = 1/(lu+ly) + (lu)/(lu+ly)**2
|
||||||
|
|
||||||
|
slices = index_to_slices(X[:,-1])
|
||||||
|
|
||||||
|
for i, ss1 in enumerate(slices):
|
||||||
|
for s1 in ss1:
|
||||||
|
if i==0:
|
||||||
|
Kdiag[s1]+= self.variance_U + self.ubias
|
||||||
|
elif i==1:
|
||||||
|
Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
|
||||||
|
else:
|
||||||
|
raise ValueError, "invalid input/output index"
|
||||||
|
#Kdiag[slices[0][0]]+= self.variance_U #matern32 diag
|
||||||
|
#Kdiag[slices[1][0]]+= self.variance_U*self.variance_Y*(k1+k2+k3) # diag
|
||||||
|
return Kdiag
|
||||||
|
|
||||||
|
|
||||||
|
def update_gradients_full(self, dL_dK, X, X2=None):
|
||||||
|
"""derivative of the covariance matrix with respect to the parameters."""
|
||||||
|
X,slices = X[:,:-1],index_to_slices(X[:,-1])
|
||||||
|
if X2 is None:
|
||||||
|
X2,slices2 = X,slices
|
||||||
|
else:
|
||||||
|
X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
|
||||||
|
#rdist = X[:,0][:,None] - X2[:,0][:,None].T
|
||||||
|
|
||||||
|
rdist = X - X2.T
|
||||||
|
ly=1/self.lengthscale_Y
|
||||||
|
lu=np.sqrt(3)/self.lengthscale_U
|
||||||
|
|
||||||
|
Vu=self.variance_U
|
||||||
|
Vy=self.variance_Y
|
||||||
|
Vyu = np.sqrt(Vy*ly*2)
|
||||||
|
dVdly = 0.5/np.sqrt(ly)*np.sqrt(2*Vy)
|
||||||
|
dVdVy = 0.5/np.sqrt(Vy)*np.sqrt(2*ly)
|
||||||
|
|
||||||
|
rd=rdist.shape[0]
|
||||||
|
dktheta1 = np.zeros([rd,rd])
|
||||||
|
dktheta2 = np.zeros([rd,rd])
|
||||||
|
dkUdvar = np.zeros([rd,rd])
|
||||||
|
dkYdvar = np.zeros([rd,rd])
|
||||||
|
|
||||||
|
dkdubias = np.zeros([rd,rd])
|
||||||
|
|
||||||
|
# dk dtheta for UU
|
||||||
|
UUdtheta1 = lambda dist: np.exp(-lu* dist)*dist + (-dist)*np.exp(-lu* dist)*(1+lu*dist)
|
||||||
|
UUdtheta2 = lambda dist: 0
|
||||||
|
#UUdvar = lambda dist: (1 + lu*dist)*np.exp(-lu*dist)
|
||||||
|
UUdvar = lambda dist: (1 + lu* np.abs(dist)) * np.exp(-lu * np.abs(dist))
|
||||||
|
|
||||||
|
# dk dtheta for YY
|
||||||
|
|
||||||
|
dk1theta1 = lambda dist: np.exp(-ly*dist)*2*(-lu)/(lu+ly)**3
|
||||||
|
|
||||||
|
dk2theta1 = lambda dist: (1.0)*(
|
||||||
|
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)
|
||||||
|
+np.exp(-dist*lu)*(ly-2*lu+ly*lu*dist-dist*lu**2)*2*(ly-lu)**(-3)
|
||||||
|
+np.exp(-dist*ly)*2*(ly-lu)**(-2)
|
||||||
|
+np.exp(-dist*ly)*2*(2*lu-ly)*(ly-lu)**(-3)
|
||||||
|
)
|
||||||
|
|
||||||
|
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)
|
||||||
|
|
||||||
|
#dktheta1 = lambda dist: self.variance_U*self.variance_Y*(dk1theta1+dk2theta1+dk3theta1)
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
dk1theta2 = lambda dist: np.exp(-ly*dist) * ((lu+ly)**(-2)) * ( (-dist)*(2*lu+ly) + 1 + (-2)*(2*lu+ly)/(lu+ly) )
|
||||||
|
|
||||||
|
dk2theta2 =lambda dist: 1*(
|
||||||
|
np.exp(-dist*lu)*(ly-lu)**(-2) * ( 1+lu*dist+(-2)*(ly-2*lu+lu*ly*dist-dist*lu**2)*(ly-lu)**(-1) )
|
||||||
|
+np.exp(-dist*ly)*(ly-lu)**(-2) * ( (-dist)*(2*lu-ly) -1+(2*lu-ly)*(-2)*(ly-lu)**(-1) )
|
||||||
|
)
|
||||||
|
|
||||||
|
dk3theta2 = lambda dist: np.exp(-dist*lu) * (-3*lu-ly-dist*lu**2-lu*ly*dist)/(lu+ly)**3
|
||||||
|
|
||||||
|
#dktheta2 = lambda dist: self.variance_U*self.variance_Y*(dk1theta2 + dk2theta2 +dk3theta2)
|
||||||
|
|
||||||
|
# kyy kernel
|
||||||
|
|
||||||
|
k1 = lambda dist: np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
|
||||||
|
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
|
||||||
|
k3 = lambda dist: np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
|
||||||
|
#dkdvar = k1+k2+k3
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
# cross covariance function
|
||||||
|
kyu3 = lambda dist:np.exp(-lu*dist)/(lu+ly)*(1+lu*(dist+1/(lu+ly)))
|
||||||
|
|
||||||
|
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) ) )
|
||||||
|
|
||||||
|
k2cros = lambda dist:np.exp(ly*dist)*( 1/(lu+ly) + lu/(lu+ly)**2 )
|
||||||
|
# cross covariance kuy
|
||||||
|
kuyp = lambda dist:(kyu3(dist)) #t>0 kuy
|
||||||
|
kuyn = lambda dist:(k1cros(dist)+k2cros(dist)) #t<0 kuy
|
||||||
|
# cross covariance kyu
|
||||||
|
kyup = lambda dist:(k1cros(-dist)+k2cros(-dist)) #t>0 kyu
|
||||||
|
kyun = lambda dist:(kyu3(-dist)) #t<0 kyu
|
||||||
|
|
||||||
|
# dk dtheta for UY
|
||||||
|
|
||||||
|
|
||||||
|
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) )
|
||||||
|
dkyu3dtheta1 = lambda dist: np.exp(-lu*dist)*(lu+ly)**(-1)* ( (-dist)*(1+dist*lu+lu*(lu+ly)**(-1)) -\
|
||||||
|
(lu+ly)**(-1)*(1+dist*lu+lu*(lu+ly)**(-1)) +dist+(lu+ly)**(-1)-lu*(lu+ly)**(-2) )
|
||||||
|
|
||||||
|
dkcros2dtheta1 = lambda dist: np.exp(ly*dist)* ( -(ly+lu)**(-2) + (ly+lu)**(-2) + (-2)*lu*(lu+ly)**(-3) )
|
||||||
|
dkcros2dtheta2 = lambda dist: np.exp(ly*dist)*dist* ( (ly+lu)**(-1) + lu*(lu+ly)**(-2) ) + \
|
||||||
|
np.exp(ly*dist)*( -(lu+ly)**(-2) + lu*(-2)*(lu+ly)**(-3) )
|
||||||
|
|
||||||
|
dkcros1dtheta1 = lambda dist: np.exp(ly*dist)*( -(lu-ly)**(-2)*( 1-np.exp((lu-ly)*dist) + lu*dist*np.exp((lu-ly)*dist)+ \
|
||||||
|
lu*(1-np.exp((lu-ly)*dist))/(lu-ly) ) + (lu-ly)**(-1)*( -np.exp( (lu-ly)*dist )*dist + dist*np.exp( (lu-ly)*dist)+\
|
||||||
|
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) -\
|
||||||
|
lu*(1-np.exp((lu-ly)*dist))/(lu-ly)**2 ) )
|
||||||
|
|
||||||
|
dkcros1dtheta2 = lambda t: np.exp(ly*t)*t/(lu-ly)*( 1-np.exp((lu-ly)*t) +lu*t*np.exp((lu-ly)*t)+\
|
||||||
|
lu*(1-np.exp((lu-ly)*t))/(lu-ly) )+\
|
||||||
|
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) )+\
|
||||||
|
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)+\
|
||||||
|
lu*( 1-np.exp((lu-ly)*t) )/(lu-ly)**2 )
|
||||||
|
|
||||||
|
dkuypdtheta1 = lambda dist:(dkyu3dtheta1(dist)) #t>0 kuy
|
||||||
|
dkuyndtheta1 = lambda dist:(dkcros1dtheta1(dist)+dkcros2dtheta1(dist)) #t<0 kuy
|
||||||
|
# cross covariance kyu
|
||||||
|
dkyupdtheta1 = lambda dist:(dkcros1dtheta1(-dist)+dkcros2dtheta1(-dist)) #t>0 kyu
|
||||||
|
dkyundtheta1 = lambda dist:(dkyu3dtheta1(-dist)) #t<0 kyu
|
||||||
|
|
||||||
|
dkuypdtheta2 = lambda dist:(dkyu3dtheta2(dist)) #t>0 kuy
|
||||||
|
dkuyndtheta2 = lambda dist:(dkcros1dtheta2(dist)+dkcros2dtheta2(dist)) #t<0 kuy
|
||||||
|
# cross covariance kyu
|
||||||
|
dkyupdtheta2 = lambda dist:(dkcros1dtheta2(-dist)+dkcros2dtheta2(-dist)) #t>0 kyu
|
||||||
|
dkyundtheta2 = lambda dist:(dkyu3dtheta2(-dist)) #t<0 kyu
|
||||||
|
|
||||||
|
|
||||||
|
for i, s1 in enumerate(slices):
|
||||||
|
for j, s2 in enumerate(slices2):
|
||||||
|
for ss1 in s1:
|
||||||
|
for ss2 in s2:
|
||||||
|
if i==0 and j==0:
|
||||||
|
#target[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
|
||||||
|
dktheta1[ss1,ss2] = Vu*UUdtheta1(np.abs(rdist[ss1,ss2]))
|
||||||
|
dktheta2[ss1,ss2] = 0
|
||||||
|
dkUdvar[ss1,ss2] = UUdvar(np.abs(rdist[ss1,ss2]))
|
||||||
|
dkYdvar[ss1,ss2] = 0
|
||||||
|
dkdubias[ss1,ss2] = 1
|
||||||
|
elif i==0 and j==1:
|
||||||
|
########target[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) ) )
|
||||||
|
#np.where( rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) ) )
|
||||||
|
#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]))) )
|
||||||
|
#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]))) )
|
||||||
|
dktheta1[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , Vu*Vyu*dkuypdtheta1(rdist[ss1,ss2]),Vu*Vyu*dkuyndtheta1(rdist[ss1,ss2]) )
|
||||||
|
dkUdvar[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , Vyu*kuyp(rdist[ss1,ss2]), Vyu* kuyn(rdist[ss1,ss2]) )
|
||||||
|
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]) )
|
||||||
|
dkYdvar[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , Vu*dVdVy*kuyp(rdist[ss1,ss2]), Vu*dVdVy* kuyn(rdist[ss1,ss2]) )
|
||||||
|
dkdubias[ss1,ss2] = 0
|
||||||
|
elif i==1 and j==1:
|
||||||
|
#target[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
|
||||||
|
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])))
|
||||||
|
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])))
|
||||||
|
dkUdvar[ss1,ss2] = self.variance_Y*(k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2])) )
|
||||||
|
dkYdvar[ss1,ss2] = self.variance_U*(k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2])) )
|
||||||
|
dkdubias[ss1,ss2] = 0
|
||||||
|
else:
|
||||||
|
#######target[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[s1[0],s2[0]]) ) )
|
||||||
|
#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])) )
|
||||||
|
#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])) )
|
||||||
|
dktheta1[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , Vu*Vyu*dkyupdtheta1(rdist[ss1,ss2]),Vu*Vyu*dkyundtheta1(rdist[ss1,ss2]) )
|
||||||
|
dkUdvar[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , Vyu*kyup(rdist[ss1,ss2]),Vyu*kyun(rdist[ss1,ss2]))
|
||||||
|
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]) )
|
||||||
|
dkYdvar[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , Vu*dVdVy*kyup(rdist[ss1,ss2]), Vu*dVdVy*kyun(rdist[ss1,ss2]))
|
||||||
|
dkdubias[ss1,ss2] = 0
|
||||||
|
#stop
|
||||||
|
self.variance_U.gradient = np.sum(dkUdvar * dL_dK) # Vu
|
||||||
|
|
||||||
|
self.variance_Y.gradient = np.sum(dkYdvar * dL_dK) # Vy
|
||||||
|
|
||||||
|
self.lengthscale_U.gradient = np.sum(dktheta1*(-np.sqrt(3)*self.lengthscale_U**(-2))* dL_dK) #lu
|
||||||
|
|
||||||
|
self.lengthscale_Y.gradient = np.sum(dktheta2*(-self.lengthscale_Y**(-2)) * dL_dK) #ly
|
||||||
|
|
||||||
|
self.ubias.gradient = np.sum(dkdubias * dL_dK)
|
||||||
|
|
||||||
267
GPy/kern/_src/ODE_st.py
Normal file
267
GPy/kern/_src/ODE_st.py
Normal file
|
|
@ -0,0 +1,267 @@
|
||||||
|
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
|
||||||
|
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||||
|
from kern import Kern
|
||||||
|
from ...core.parameterization import Param
|
||||||
|
from ...core.parameterization.transformations import Logexp
|
||||||
|
import numpy as np
|
||||||
|
from independent_outputs import index_to_slices
|
||||||
|
|
||||||
|
|
||||||
|
class ODE_st(Kern):
|
||||||
|
"""
|
||||||
|
kernel resultiong from a first order ODE with OU driving GP
|
||||||
|
|
||||||
|
:param input_dim: the number of input dimension, has to be equal to one
|
||||||
|
:type input_dim: int
|
||||||
|
:param varianceU: variance of the driving GP
|
||||||
|
:type varianceU: float
|
||||||
|
:param lengthscaleU: lengthscale of the driving GP (sqrt(3)/lengthscaleU)
|
||||||
|
:type lengthscaleU: float
|
||||||
|
:param varianceY: 'variance' of the transfer function
|
||||||
|
:type varianceY: float
|
||||||
|
:param lengthscaleY: 'lengthscale' of the transfer function (1/lengthscaleY)
|
||||||
|
:type lengthscaleY: float
|
||||||
|
:rtype: kernel object
|
||||||
|
|
||||||
|
"""
|
||||||
|
|
||||||
|
def __init__(self, input_dim, a=1.,b=1., c=1.,variance_Yx=3.,variance_Yt=1.5, lengthscale_Yx=1.5, lengthscale_Yt=1.5, active_dims=None, name='ode_st'):
|
||||||
|
assert input_dim ==3, "only defined for 3 input dims"
|
||||||
|
super(ODE_st, self).__init__(input_dim, active_dims, name)
|
||||||
|
|
||||||
|
self.variance_Yt = Param('variance_Yt', variance_Yt, Logexp())
|
||||||
|
self.variance_Yx = Param('variance_Yx', variance_Yx, Logexp())
|
||||||
|
self.lengthscale_Yt = Param('lengthscale_Yt', lengthscale_Yt, Logexp())
|
||||||
|
self.lengthscale_Yx = Param('lengthscale_Yx', lengthscale_Yx, Logexp())
|
||||||
|
|
||||||
|
self.a= Param('a', a, Logexp())
|
||||||
|
self.b = Param('b', b, Logexp())
|
||||||
|
self.c = Param('c', c, Logexp())
|
||||||
|
|
||||||
|
self.add_parameters(self.a, self.b, self.c, self.variance_Yt, self.variance_Yx, self.lengthscale_Yt,self.lengthscale_Yx)
|
||||||
|
|
||||||
|
|
||||||
|
def K(self, X, X2=None):
|
||||||
|
# model : -a d^2y/dx^2 + b dy/dt + c * y = U
|
||||||
|
# kernel Kyy rbf spatiol temporal
|
||||||
|
# vyt Y temporal variance vyx Y spatiol variance lyt Y temporal lengthscale lyx Y spatiol lengthscale
|
||||||
|
# kernel Kuu doper( doper(Kyy))
|
||||||
|
# a b c lyt lyx vyx*vyt
|
||||||
|
"""Compute the covariance matrix between X and X2."""
|
||||||
|
X,slices = X[:,:-1],index_to_slices(X[:,-1])
|
||||||
|
if X2 is None:
|
||||||
|
X2,slices2 = X,slices
|
||||||
|
K = np.zeros((X.shape[0], X.shape[0]))
|
||||||
|
else:
|
||||||
|
X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
|
||||||
|
K = np.zeros((X.shape[0], X2.shape[0]))
|
||||||
|
|
||||||
|
|
||||||
|
tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
|
||||||
|
xdist = (X[:,1][:,None] - X2[:,1][None,:])**2
|
||||||
|
|
||||||
|
ttdist = (X[:,0][:,None] - X2[:,0][None,:])
|
||||||
|
#rdist = [tdist,xdist]
|
||||||
|
#dist = np.abs(X - X2.T)
|
||||||
|
vyt = self.variance_Yt
|
||||||
|
vyx = self.variance_Yx
|
||||||
|
|
||||||
|
lyt=1/(2*self.lengthscale_Yt)
|
||||||
|
lyx=1/(2*self.lengthscale_Yx)
|
||||||
|
|
||||||
|
a = self.a ## -a is used in the model, negtive diffusion
|
||||||
|
b = self.b
|
||||||
|
c = self.c
|
||||||
|
|
||||||
|
kyy = lambda tdist,xdist: np.exp(-lyt*(tdist) -lyx*(xdist))
|
||||||
|
|
||||||
|
k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )
|
||||||
|
|
||||||
|
k2 = lambda xdist: ( 4*lyx**2 * (xdist) - 2*lyx )
|
||||||
|
|
||||||
|
k3 = lambda xdist: ( 3*4*lyx**2 - 6*8*xdist*lyx**3 + 16*xdist**2*lyx**4 )
|
||||||
|
|
||||||
|
k4 = lambda ttdist: 2*lyt*(ttdist)
|
||||||
|
|
||||||
|
for i, s1 in enumerate(slices):
|
||||||
|
for j, s2 in enumerate(slices2):
|
||||||
|
for ss1 in s1:
|
||||||
|
for ss2 in s2:
|
||||||
|
if i==0 and j==0:
|
||||||
|
K[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
elif i==0 and j==1:
|
||||||
|
K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) + b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
#K[ss1,ss2]= np.where( rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[ss1,ss2]) ) )
|
||||||
|
#K[ss1,ss2]= np.where( rdist[ss1,ss2]>0 , kuyp(rdist[ss1,ss2]), kuyn(rdist[ss1,ss2] ) )
|
||||||
|
elif i==1 and j==1:
|
||||||
|
K[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )* vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
else:
|
||||||
|
K[ss1,ss2] = (-a*k2(xdist[ss1,ss2]) - b*k4(ttdist[ss1,ss2]) + c)*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
#K[ss1,ss2]= np.where( rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[ss1,ss2]) ) )
|
||||||
|
#K[ss1,ss2] = np.where( rdist[ss1,ss2]>0 , kyup(rdist[ss1,ss2]), kyun(rdist[ss1,ss2] ) )
|
||||||
|
|
||||||
|
#stop
|
||||||
|
return K
|
||||||
|
|
||||||
|
def Kdiag(self, X):
|
||||||
|
"""Compute the diagonal of the covariance matrix associated to X."""
|
||||||
|
vyt = self.variance_Yt
|
||||||
|
vyx = self.variance_Yx
|
||||||
|
|
||||||
|
lyt = 1./(2*self.lengthscale_Yt)
|
||||||
|
lyx = 1./(2*self.lengthscale_Yx)
|
||||||
|
|
||||||
|
a = self.a
|
||||||
|
b = self.b
|
||||||
|
c = self.c
|
||||||
|
|
||||||
|
## dk^2/dtdt'
|
||||||
|
k1 = (2*lyt )*vyt*vyx
|
||||||
|
## dk^2/dx^2
|
||||||
|
k2 = ( - 2*lyx )*vyt*vyx
|
||||||
|
## dk^4/dx^2dx'^2
|
||||||
|
k3 = ( 4*3*lyx**2 )*vyt*vyx
|
||||||
|
|
||||||
|
|
||||||
|
Kdiag = np.zeros(X.shape[0])
|
||||||
|
slices = index_to_slices(X[:,-1])
|
||||||
|
|
||||||
|
for i, ss1 in enumerate(slices):
|
||||||
|
for s1 in ss1:
|
||||||
|
if i==0:
|
||||||
|
Kdiag[s1]+= vyt*vyx
|
||||||
|
elif i==1:
|
||||||
|
#i=1
|
||||||
|
Kdiag[s1]+= b**2*k1 - 2*a*c*k2 + a**2*k3 + c**2*vyt*vyx
|
||||||
|
#Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
|
||||||
|
else:
|
||||||
|
raise ValueError, "invalid input/output index"
|
||||||
|
|
||||||
|
return Kdiag
|
||||||
|
|
||||||
|
|
||||||
|
def update_gradients_full(self, dL_dK, X, X2=None):
|
||||||
|
#def dK_dtheta(self, dL_dK, X, X2, target):
|
||||||
|
"""derivative of the covariance matrix with respect to the parameters."""
|
||||||
|
X,slices = X[:,:-1],index_to_slices(X[:,-1])
|
||||||
|
if X2 is None:
|
||||||
|
X2,slices2 = X,slices
|
||||||
|
K = np.zeros((X.shape[0], X.shape[0]))
|
||||||
|
else:
|
||||||
|
X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
|
||||||
|
|
||||||
|
vyt = self.variance_Yt
|
||||||
|
vyx = self.variance_Yx
|
||||||
|
|
||||||
|
lyt = 1./(2*self.lengthscale_Yt)
|
||||||
|
lyx = 1./(2*self.lengthscale_Yx)
|
||||||
|
|
||||||
|
a = self.a
|
||||||
|
b = self.b
|
||||||
|
c = self.c
|
||||||
|
|
||||||
|
tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
|
||||||
|
xdist = (X[:,1][:,None] - X2[:,1][None,:])**2
|
||||||
|
#rdist = [tdist,xdist]
|
||||||
|
ttdist = (X[:,0][:,None] - X2[:,0][None,:])
|
||||||
|
|
||||||
|
rd=tdist.shape[0]
|
||||||
|
|
||||||
|
dka = np.zeros([rd,rd])
|
||||||
|
dkb = np.zeros([rd,rd])
|
||||||
|
dkc = np.zeros([rd,rd])
|
||||||
|
dkYdvart = np.zeros([rd,rd])
|
||||||
|
dkYdvarx = np.zeros([rd,rd])
|
||||||
|
dkYdlent = np.zeros([rd,rd])
|
||||||
|
dkYdlenx = np.zeros([rd,rd])
|
||||||
|
|
||||||
|
|
||||||
|
kyy = lambda tdist,xdist: np.exp(-lyt*(tdist) -lyx*(xdist))
|
||||||
|
#k1 = lambda tdist: (lyt - lyt**2 * (tdist) )
|
||||||
|
#k2 = lambda xdist: ( lyx**2 * (xdist) - lyx )
|
||||||
|
#k3 = lambda xdist: ( 3*lyx**2 - 6*xdist*lyx**3 + xdist**2*lyx**4 )
|
||||||
|
#k4 = lambda tdist: -lyt*np.sqrt(tdist)
|
||||||
|
|
||||||
|
k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )
|
||||||
|
|
||||||
|
k2 = lambda xdist: ( 4*lyx**2 * (xdist) - 2*lyx )
|
||||||
|
|
||||||
|
k3 = lambda xdist: ( 3*4*lyx**2 - 6*8*xdist*lyx**3 + 16*xdist**2*lyx**4 )
|
||||||
|
|
||||||
|
k4 = lambda ttdist: 2*lyt*(ttdist)
|
||||||
|
|
||||||
|
dkyydlyx = lambda tdist,xdist: kyy(tdist,xdist)*(-xdist)
|
||||||
|
dkyydlyt = lambda tdist,xdist: kyy(tdist,xdist)*(-tdist)
|
||||||
|
|
||||||
|
dk1dlyt = lambda tdist: 2. - 4*2.*lyt*tdist
|
||||||
|
dk2dlyx = lambda xdist: (4.*2.*lyx*xdist -2.)
|
||||||
|
dk3dlyx = lambda xdist: (6.*4.*lyx - 18.*8*xdist*lyx**2 + 4*16*xdist**2*lyx**3)
|
||||||
|
|
||||||
|
dk4dlyt = lambda ttdist: 2*(ttdist)
|
||||||
|
|
||||||
|
for i, s1 in enumerate(slices):
|
||||||
|
for j, s2 in enumerate(slices2):
|
||||||
|
for ss1 in s1:
|
||||||
|
for ss2 in s2:
|
||||||
|
if i==0 and j==0:
|
||||||
|
dka[ss1,ss2] = 0
|
||||||
|
dkb[ss1,ss2] = 0
|
||||||
|
dkc[ss1,ss2] = 0
|
||||||
|
dkYdvart[ss1,ss2] = vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdvarx[ss1,ss2] = vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
elif i==0 and j==1:
|
||||||
|
dka[ss1,ss2] = -k2(xdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkb[ss1,ss2] = k4(ttdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkc[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
#dkYdvart[ss1,ss2] = 0
|
||||||
|
#dkYdvarx[ss1,ss2] = 0
|
||||||
|
#dkYdlent[ss1,ss2] = 0
|
||||||
|
#dkYdlenx[ss1,ss2] = 0
|
||||||
|
dkYdvart[ss1,ss2] = (-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdvarx[ss1,ss2] = (-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])* (-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)+\
|
||||||
|
vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*b*dk4dlyt(ttdist[ss1,ss2])
|
||||||
|
dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*k2(xdist[ss1,ss2])+b*k4(ttdist[ss1,ss2])+c)+\
|
||||||
|
vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*dk2dlyx(xdist[ss1,ss2]))
|
||||||
|
elif i==1 and j==1:
|
||||||
|
dka[ss1,ss2] = (2*a*k3(xdist[ss1,ss2]) - 2*c*k2(xdist[ss1,ss2]))*vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkb[ss1,ss2] = 2*b*k1(tdist[ss1,ss2])*vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkc[ss1,ss2] = (-2*a*k2(xdist[ss1,ss2]) + 2*c )*vyt*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdvart[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )*vyx* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdvarx[ss1,ss2] = ( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 )*vyt* kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])*( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 ) +\
|
||||||
|
vyx*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*b**2*dk1dlyt(tdist[ss1,ss2])
|
||||||
|
dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])*( b**2*k1(tdist[ss1,ss2]) - 2*a*c*k2(xdist[ss1,ss2]) + a**2*k3(xdist[ss1,ss2]) + c**2 ) +\
|
||||||
|
vyx*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])* (-2*a*c*dk2dlyx(xdist[ss1,ss2]) + a**2*dk3dlyx(xdist[ss1,ss2]) )
|
||||||
|
else:
|
||||||
|
dka[ss1,ss2] = -k2(xdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkb[ss1,ss2] = -k4(ttdist[ss1,ss2])*vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkc[ss1,ss2] = vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
#dkYdvart[ss1,ss2] = 0
|
||||||
|
#dkYdvarx[ss1,ss2] = 0
|
||||||
|
#dkYdlent[ss1,ss2] = 0
|
||||||
|
#dkYdlenx[ss1,ss2] = 0
|
||||||
|
dkYdvart[ss1,ss2] = (-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdvarx[ss1,ss2] = (-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)*vyt*kyy(tdist[ss1,ss2],xdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*vyx*dkyydlyt(tdist[ss1,ss2],xdist[ss1,ss2])* (-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)+\
|
||||||
|
vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*(-1)*b*dk4dlyt(ttdist[ss1,ss2])
|
||||||
|
dkYdlenx[ss1,ss2] = vyt*vyx*dkyydlyx(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*k2(xdist[ss1,ss2])-b*k4(ttdist[ss1,ss2])+c)+\
|
||||||
|
vyt*vyx*kyy(tdist[ss1,ss2],xdist[ss1,ss2])*(-a*dk2dlyx(xdist[ss1,ss2]))
|
||||||
|
|
||||||
|
self.a.gradient = np.sum(dka * dL_dK)
|
||||||
|
|
||||||
|
self.b.gradient = np.sum(dkb * dL_dK)
|
||||||
|
|
||||||
|
self.c.gradient = np.sum(dkc * dL_dK)
|
||||||
|
|
||||||
|
|
||||||
|
self.variance_Yt.gradient = np.sum(dkYdvart * dL_dK) # Vy
|
||||||
|
|
||||||
|
self.variance_Yx.gradient = np.sum(dkYdvarx * dL_dK)
|
||||||
|
|
||||||
|
self.lengthscale_Yt.gradient = np.sum(dkYdlent*(-0.5*self.lengthscale_Yt**(-2)) * dL_dK) #ly np.sum(dktheta2*(-self.lengthscale_Y**(-2)) * dL_dK)
|
||||||
|
|
||||||
|
self.lengthscale_Yx.gradient = np.sum(dkYdlenx*(-0.5*self.lengthscale_Yx**(-2)) * dL_dK)
|
||||||
|
|
||||||
165
GPy/kern/_src/ODE_t.py
Normal file
165
GPy/kern/_src/ODE_t.py
Normal file
|
|
@ -0,0 +1,165 @@
|
||||||
|
from kern import Kern
|
||||||
|
from ...core.parameterization import Param
|
||||||
|
from ...core.parameterization.transformations import Logexp
|
||||||
|
import numpy as np
|
||||||
|
from independent_outputs import index_to_slices
|
||||||
|
|
||||||
|
|
||||||
|
class ODE_t(Kern):
|
||||||
|
|
||||||
|
def __init__(self, input_dim, a=1., c=1.,variance_Yt=3., lengthscale_Yt=1.5,ubias =1., active_dims=None, name='ode_st'):
|
||||||
|
assert input_dim ==2, "only defined for 2 input dims"
|
||||||
|
super(ODE_t, self).__init__(input_dim, active_dims, name)
|
||||||
|
|
||||||
|
self.variance_Yt = Param('variance_Yt', variance_Yt, Logexp())
|
||||||
|
self.lengthscale_Yt = Param('lengthscale_Yt', lengthscale_Yt, Logexp())
|
||||||
|
|
||||||
|
self.a= Param('a', a, Logexp())
|
||||||
|
self.c = Param('c', c, Logexp())
|
||||||
|
self.ubias = Param('ubias', ubias, Logexp())
|
||||||
|
self.add_parameters(self.a, self.c, self.variance_Yt, self.lengthscale_Yt,self.ubias)
|
||||||
|
|
||||||
|
def K(self, X, X2=None):
|
||||||
|
"""Compute the covariance matrix between X and X2."""
|
||||||
|
X,slices = X[:,:-1],index_to_slices(X[:,-1])
|
||||||
|
if X2 is None:
|
||||||
|
X2,slices2 = X,slices
|
||||||
|
K = np.zeros((X.shape[0], X.shape[0]))
|
||||||
|
else:
|
||||||
|
X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
|
||||||
|
K = np.zeros((X.shape[0], X2.shape[0]))
|
||||||
|
|
||||||
|
tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
|
||||||
|
ttdist = (X[:,0][:,None] - X2[:,0][None,:])
|
||||||
|
|
||||||
|
vyt = self.variance_Yt
|
||||||
|
|
||||||
|
lyt=1/(2*self.lengthscale_Yt)
|
||||||
|
|
||||||
|
a = -self.a
|
||||||
|
c = self.c
|
||||||
|
|
||||||
|
kyy = lambda tdist: np.exp(-lyt*(tdist))
|
||||||
|
|
||||||
|
k1 = lambda tdist: (2*lyt - 4*lyt**2 *(tdist) )
|
||||||
|
|
||||||
|
k4 = lambda tdist: 2*lyt*(tdist)
|
||||||
|
|
||||||
|
for i, s1 in enumerate(slices):
|
||||||
|
for j, s2 in enumerate(slices2):
|
||||||
|
for ss1 in s1:
|
||||||
|
for ss2 in s2:
|
||||||
|
if i==0 and j==0:
|
||||||
|
K[ss1,ss2] = vyt*kyy(tdist[ss1,ss2])
|
||||||
|
elif i==0 and j==1:
|
||||||
|
K[ss1,ss2] = (k4(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
|
||||||
|
#K[ss1,ss2] = (2*lyt*(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
|
||||||
|
elif i==1 and j==1:
|
||||||
|
K[ss1,ss2] = ( k1(tdist[ss1,ss2]) + 1. )*vyt* kyy(tdist[ss1,ss2])+self.ubias
|
||||||
|
else:
|
||||||
|
K[ss1,ss2] = (-k4(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
|
||||||
|
#K[ss1,ss2] = (-2*lyt*(ttdist[ss1,ss2])+1)*vyt*kyy(tdist[ss1,ss2])
|
||||||
|
#stop
|
||||||
|
return K
|
||||||
|
|
||||||
|
|
||||||
|
def Kdiag(self, X):
|
||||||
|
|
||||||
|
vyt = self.variance_Yt
|
||||||
|
lyt = 1./(2*self.lengthscale_Yt)
|
||||||
|
|
||||||
|
a = -self.a
|
||||||
|
c = self.c
|
||||||
|
|
||||||
|
k1 = (2*lyt )*vyt
|
||||||
|
|
||||||
|
Kdiag = np.zeros(X.shape[0])
|
||||||
|
slices = index_to_slices(X[:,-1])
|
||||||
|
|
||||||
|
for i, ss1 in enumerate(slices):
|
||||||
|
for s1 in ss1:
|
||||||
|
if i==0:
|
||||||
|
Kdiag[s1]+= vyt
|
||||||
|
elif i==1:
|
||||||
|
#i=1
|
||||||
|
Kdiag[s1]+= k1 + vyt+self.ubias
|
||||||
|
#Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
|
||||||
|
else:
|
||||||
|
raise ValueError, "invalid input/output index"
|
||||||
|
|
||||||
|
return Kdiag
|
||||||
|
|
||||||
|
def update_gradients_full(self, dL_dK, X, X2=None):
|
||||||
|
"""derivative of the covariance matrix with respect to the parameters."""
|
||||||
|
X,slices = X[:,:-1],index_to_slices(X[:,-1])
|
||||||
|
if X2 is None:
|
||||||
|
X2,slices2 = X,slices
|
||||||
|
K = np.zeros((X.shape[0], X.shape[0]))
|
||||||
|
else:
|
||||||
|
X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
|
||||||
|
|
||||||
|
|
||||||
|
vyt = self.variance_Yt
|
||||||
|
|
||||||
|
lyt = 1./(2*self.lengthscale_Yt)
|
||||||
|
|
||||||
|
tdist = (X[:,0][:,None] - X2[:,0][None,:])**2
|
||||||
|
ttdist = (X[:,0][:,None] - X2[:,0][None,:])
|
||||||
|
#rdist = [tdist,xdist]
|
||||||
|
|
||||||
|
rd=tdist.shape[0]
|
||||||
|
|
||||||
|
dka = np.zeros([rd,rd])
|
||||||
|
dkc = np.zeros([rd,rd])
|
||||||
|
dkYdvart = np.zeros([rd,rd])
|
||||||
|
dkYdlent = np.zeros([rd,rd])
|
||||||
|
|
||||||
|
dkdubias = np.zeros([rd,rd])
|
||||||
|
|
||||||
|
kyy = lambda tdist: np.exp(-lyt*(tdist))
|
||||||
|
dkyydlyt = lambda tdist: kyy(tdist)*(-tdist)
|
||||||
|
|
||||||
|
k1 = lambda tdist: (2*lyt - 4*lyt**2 * (tdist) )
|
||||||
|
|
||||||
|
k4 = lambda ttdist: 2*lyt*(ttdist)
|
||||||
|
|
||||||
|
dk1dlyt = lambda tdist: 2. - 4*2.*lyt*tdist
|
||||||
|
|
||||||
|
dk4dlyt = lambda ttdist: 2*(ttdist)
|
||||||
|
|
||||||
|
for i, s1 in enumerate(slices):
|
||||||
|
for j, s2 in enumerate(slices2):
|
||||||
|
for ss1 in s1:
|
||||||
|
for ss2 in s2:
|
||||||
|
if i==0 and j==0:
|
||||||
|
dkYdvart[ss1,ss2] = kyy(tdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])
|
||||||
|
dkdubias[ss1,ss2] = 0
|
||||||
|
elif i==0 and j==1:
|
||||||
|
dkYdvart[ss1,ss2] = (k4(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
|
||||||
|
#dkYdvart[ss1,ss2] = ((2*lyt*ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (k4(ttdist[ss1,ss2])+1.)+\
|
||||||
|
vyt*kyy(tdist[ss1,ss2])*(dk4dlyt(ttdist[ss1,ss2]))
|
||||||
|
#dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (2*lyt*(ttdist[ss1,ss2])+1.)+\
|
||||||
|
#vyt*kyy(tdist[ss1,ss2])*(2*ttdist[ss1,ss2])
|
||||||
|
dkdubias[ss1,ss2] = 0
|
||||||
|
elif i==1 and j==1:
|
||||||
|
dkYdvart[ss1,ss2] = (k1(tdist[ss1,ss2]) + 1. )* kyy(tdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])*( k1(tdist[ss1,ss2]) + 1. ) +\
|
||||||
|
vyt*kyy(tdist[ss1,ss2])*dk1dlyt(tdist[ss1,ss2])
|
||||||
|
dkdubias[ss1,ss2] = 1
|
||||||
|
else:
|
||||||
|
dkYdvart[ss1,ss2] = (-k4(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
|
||||||
|
#dkYdvart[ss1,ss2] = (-2*lyt*(ttdist[ss1,ss2])+1)*kyy(tdist[ss1,ss2])
|
||||||
|
dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (-k4(ttdist[ss1,ss2])+1.)+\
|
||||||
|
vyt*kyy(tdist[ss1,ss2])*(-dk4dlyt(ttdist[ss1,ss2]) )
|
||||||
|
dkdubias[ss1,ss2] = 0
|
||||||
|
#dkYdlent[ss1,ss2] = vyt*dkyydlyt(tdist[ss1,ss2])* (-2*lyt*(ttdist[ss1,ss2])+1.)+\
|
||||||
|
#vyt*kyy(tdist[ss1,ss2])*(-2)*(ttdist[ss1,ss2])
|
||||||
|
|
||||||
|
|
||||||
|
self.variance_Yt.gradient = np.sum(dkYdvart * dL_dK)
|
||||||
|
|
||||||
|
self.lengthscale_Yt.gradient = np.sum(dkYdlent*(-0.5*self.lengthscale_Yt**(-2)) * dL_dK)
|
||||||
|
|
||||||
|
self.ubias.gradient = np.sum(dkdubias * dL_dK)
|
||||||
|
|
@ -180,6 +180,9 @@ class Hierarchical(CombinationKernel):
|
||||||
def Kdiag(self,X):
|
def Kdiag(self,X):
|
||||||
return np.diag(self.K(X))
|
return np.diag(self.K(X))
|
||||||
|
|
||||||
|
def gradients_X(self, dL_dK, X, X2=None):
|
||||||
|
raise NotImplementedError
|
||||||
|
|
||||||
def update_gradients_full(self,dL_dK,X,X2=None):
|
def update_gradients_full(self,dL_dK,X,X2=None):
|
||||||
slices = [index_to_slices(X[:,i]) for i in self.extra_dims]
|
slices = [index_to_slices(X[:,i]) for i in self.extra_dims]
|
||||||
if X2 is None:
|
if X2 is None:
|
||||||
|
|
|
||||||
Loading…
Add table
Add a link
Reference in a new issue