mirror of
https://github.com/SheffieldML/GPy.git
synced 2026-05-21 14:05:14 +02:00
a little merge
This commit is contained in:
Zhenwen Dai
commit
3ee2dc920c
50 changed files with 2255 additions and 10642 deletions
|
|
@ -3,15 +3,18 @@ from _src.rbf import RBF
|
|||
from _src.linear import Linear, LinearFull
|
||||
from _src.static import Bias, White
|
||||
from _src.brownian import Brownian
|
||||
from _src.stationary import Exponential, Matern32, Matern52, ExpQuad, RatQuad, Cosine
|
||||
from _src.stationary import Exponential, OU, Matern32, Matern52, ExpQuad, RatQuad, Cosine
|
||||
from _src.mlp import MLP
|
||||
from _src.periodic import PeriodicExponential, PeriodicMatern32, PeriodicMatern52
|
||||
from _src.independent_outputs import IndependentOutputs, Hierarchical
|
||||
from _src.coregionalize import Coregionalize
|
||||
from _src.ssrbf import SSRBF # TODO: ZD: did you remove this?
|
||||
from _src.ODE_UY import ODE_UY
|
||||
#from _src.ODE_UYC import ODE_UYC ADD THIS FILE TO THE REPO!!
|
||||
#from _src.ODE_st import ODE_st
|
||||
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
|
||||
|
||||
# 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
|
||||
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)
|
||||
|
|
@ -141,10 +141,10 @@ class Add(CombinationKernel):
|
|||
from static import White, Bias
|
||||
target_mu = np.zeros(variational_posterior.shape)
|
||||
target_S = np.zeros(variational_posterior.shape)
|
||||
for p1 in self._parameters_:
|
||||
for p1 in self.parameters:
|
||||
#compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2!
|
||||
eff_dL_dpsi1 = dL_dpsi1.copy()
|
||||
for p2 in self._parameters_:
|
||||
for p2 in self.parameters:
|
||||
if p2 is p1:
|
||||
continue
|
||||
if isinstance(p2, White):
|
||||
|
|
@ -160,7 +160,7 @@ class Add(CombinationKernel):
|
|||
|
||||
def add(self, other, name='sum'):
|
||||
if isinstance(other, Add):
|
||||
other_params = other._parameters_[:]
|
||||
other_params = other.parameters[:]
|
||||
for p in other_params:
|
||||
other.remove_parameter(p)
|
||||
self.add_parameters(*other_params)
|
||||
|
|
@ -170,7 +170,4 @@ class Add(CombinationKernel):
|
|||
return self
|
||||
|
||||
def input_sensitivity(self):
|
||||
in_sen = np.zeros(self.input_dim)
|
||||
for i, p in enumerate(self.parts):
|
||||
in_sen[p.active_dims] += p.input_sensitivity()
|
||||
return in_sen
|
||||
return reduce(np.add, [k.input_sensitivity() for k in self.parts])
|
||||
|
|
|
|||
|
|
@ -32,7 +32,7 @@ def index_to_slices(index):
|
|||
[ret[ind_i].append(slice(*indexes_i)) for ind_i,indexes_i in zip(ind[switchpoints[:-1]],zip(switchpoints,switchpoints[1:]))]
|
||||
return ret
|
||||
|
||||
class IndependentOutputs(Kern):
|
||||
class IndependentOutputs(CombinationKernel):
|
||||
"""
|
||||
A kernel which can represent several independent functions. this kernel
|
||||
'switches off' parts of the matrix where the output indexes are different.
|
||||
|
|
@ -180,6 +180,9 @@ class Hierarchical(CombinationKernel):
|
|||
def Kdiag(self,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):
|
||||
slices = [index_to_slices(X[:,i]) for i in self.extra_dims]
|
||||
if X2 is None:
|
||||
|
|
|
|||
|
|
@ -34,40 +34,28 @@ class Kern(Parameterized):
|
|||
|
||||
is the active_dimensions of inputs X we will work on.
|
||||
All kernels will get sliced Xes as inputs, if active_dims is not None
|
||||
Only positive integers are allowed in active_dims!
|
||||
if active_dims is None, slicing is switched off and all X will be passed through as given.
|
||||
|
||||
:param int input_dim: the number of input dimensions to the function
|
||||
:param array-like|slice|None active_dims: list of indices on which dimensions this kernel works on, or none if no slicing
|
||||
:param array-like|None active_dims: list of indices on which dimensions this kernel works on, or none if no slicing
|
||||
|
||||
Do not instantiate.
|
||||
"""
|
||||
super(Kern, self).__init__(name=name, *a, **kw)
|
||||
try:
|
||||
self.input_dim = int(input_dim)
|
||||
self.active_dims = active_dims# if active_dims is not None else slice(0, input_dim, 1)
|
||||
except TypeError:
|
||||
# input_dim is something else then an integer
|
||||
self.input_dim = input_dim
|
||||
if active_dims is not None:
|
||||
print "WARNING: given input_dim={} is not an integer and active_dims={} is given, switching off slicing"
|
||||
self.active_dims = None
|
||||
self.input_dim = int(input_dim)
|
||||
|
||||
if active_dims is None:
|
||||
active_dims = np.arange(input_dim)
|
||||
|
||||
self.active_dims = np.array(active_dims, dtype=int)
|
||||
|
||||
assert self.active_dims.size == self.input_dim, "input_dim={} does not match len(active_dim)={}, active_dims={}".format(self.input_dim, self.active_dims.size, self.active_dims)
|
||||
|
||||
if self.active_dims is not None and self.input_dim is not None:
|
||||
assert isinstance(self.active_dims, (slice, list, tuple, np.ndarray)), 'active_dims needs to be an array-like or slice object over dimensions, {} given'.format(self.active_dims.__class__)
|
||||
if isinstance(self.active_dims, slice):
|
||||
self.active_dims = slice(self.active_dims.start or 0, self.active_dims.stop or self.input_dim, self.active_dims.step or 1)
|
||||
active_dim_size = int(np.round((self.active_dims.stop-self.active_dims.start)/self.active_dims.step))
|
||||
elif isinstance(self.active_dims, np.ndarray):
|
||||
#assert np.all(self.active_dims >= 0), 'active dimensions need to be positive. negative indexing is not allowed'
|
||||
assert self.active_dims.ndim == 1, 'only flat indices allowed, given active_dims.shape={}, provide only indexes to the dimensions (columns) of the input'.format(self.active_dims.shape)
|
||||
active_dim_size = self.active_dims.size
|
||||
else:
|
||||
active_dim_size = len(self.active_dims)
|
||||
assert active_dim_size == self.input_dim, "input_dim={} does not match len(active_dim)={}, active_dims={}".format(self.input_dim, active_dim_size, self.active_dims)
|
||||
self._sliced_X = 0
|
||||
self.useGPU = self._support_GPU and useGPU
|
||||
|
||||
@Cache_this(limit=10)
|
||||
@Cache_this(limit=20)
|
||||
def _slice_X(self, X):
|
||||
return X[:, self.active_dims]
|
||||
|
||||
|
|
@ -176,8 +164,8 @@ class Kern(Parameterized):
|
|||
"""
|
||||
Shortcut for tensor `prod`.
|
||||
"""
|
||||
assert self.active_dims == range(self.input_dim), "Can only use kernels, which have their input_dims defined from 0"
|
||||
assert other.active_dims == range(other.input_dim), "Can only use kernels, which have their input_dims defined from 0"
|
||||
assert np.all(self.active_dims == range(self.input_dim)), "Can only use kernels, which have their input_dims defined from 0"
|
||||
assert np.all(other.active_dims == range(other.input_dim)), "Can only use kernels, which have their input_dims defined from 0"
|
||||
other.active_dims += self.input_dim
|
||||
return self.prod(other)
|
||||
|
||||
|
|
@ -195,19 +183,19 @@ class Kern(Parameterized):
|
|||
assert isinstance(other, Kern), "only kernels can be added to kernels..."
|
||||
from prod import Prod
|
||||
#kernels = []
|
||||
#if isinstance(self, Prod): kernels.extend(self._parameters_)
|
||||
#if isinstance(self, Prod): kernels.extend(self.parameters)
|
||||
#else: kernels.append(self)
|
||||
#if isinstance(other, Prod): kernels.extend(other._parameters_)
|
||||
#if isinstance(other, Prod): kernels.extend(other.parameters)
|
||||
#else: kernels.append(other)
|
||||
return Prod([self, other], name)
|
||||
|
||||
def _check_input_dim(self, X):
|
||||
assert X.shape[1] == self.input_dim, "You did not specify active_dims and X has wrong shape: X_dim={}, whereas input_dim={}".format(X.shape[1], self.input_dim)
|
||||
|
||||
def _check_active_dims(self, X):
|
||||
assert X.shape[1] >= len(np.r_[self.active_dims]), "At least {} dimensional X needed, X.shape={!s}".format(len(np.r_[self.active_dims]), X.shape)
|
||||
assert X.shape[1] == self.input_dim, "{} did not specify active_dims and X has wrong shape: X_dim={}, whereas input_dim={}".format(self.name, X.shape[1], self.input_dim)
|
||||
|
||||
def _check_active_dims(self, X):
|
||||
assert X.shape[1] >= len(self.active_dims), "At least {} dimensional X needed, X.shape={!s}".format(len(self.active_dims), X.shape)
|
||||
|
||||
|
||||
|
||||
class CombinationKernel(Kern):
|
||||
"""
|
||||
Abstract super class for combination kernels.
|
||||
|
|
@ -222,9 +210,10 @@ class CombinationKernel(Kern):
|
|||
|
||||
:param list kernels: List of kernels to combine (can be only one element)
|
||||
:param str name: name of the combination kernel
|
||||
:param array-like|slice extra_dims: if needed extra dimensions for the combination kernel to work on
|
||||
:param array-like extra_dims: if needed extra dimensions for the combination kernel to work on
|
||||
"""
|
||||
assert all([isinstance(k, Kern) for k in kernels])
|
||||
extra_dims = np.array(extra_dims, dtype=int)
|
||||
input_dim, active_dims = self.get_input_dim_active_dims(kernels, extra_dims)
|
||||
# initialize the kernel with the full input_dim
|
||||
super(CombinationKernel, self).__init__(input_dim, active_dims, name)
|
||||
|
|
@ -233,21 +222,23 @@ class CombinationKernel(Kern):
|
|||
|
||||
@property
|
||||
def parts(self):
|
||||
return self._parameters_
|
||||
return self.parameters
|
||||
|
||||
def get_input_dim_active_dims(self, kernels, extra_dims = None):
|
||||
#active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int))
|
||||
#active_dims = np.array(np.concatenate((active_dims, extra_dims if extra_dims is not None else [])), dtype=int)
|
||||
input_dim = np.array([k.input_dim for k in kernels])
|
||||
if np.all(input_dim[0]==input_dim):
|
||||
input_dim = input_dim[0]
|
||||
active_dims = None
|
||||
input_dim = reduce(max, (k.active_dims.max() for k in kernels)) + 1
|
||||
|
||||
if extra_dims is not None:
|
||||
input_dim += extra_dims.size
|
||||
|
||||
active_dims = np.arange(input_dim)
|
||||
return input_dim, active_dims
|
||||
|
||||
def input_sensitivity(self):
|
||||
raise NotImplementedError("Choose the kernel you want to get the sensitivity for. You need to override the default behaviour for getting the input sensitivity to be able to get the input sensitivity. For sum kernel it is the sum of all sensitivities, TODO: product kernel? Other kernels?, also TODO: shall we return all the sensitivities here in the combination kernel? So we can combine them however we want? This could lead to just plot all the sensitivities here...")
|
||||
|
||||
def _check_input_dim(self, X):
|
||||
def _check_active_dims(self, X):
|
||||
return
|
||||
|
||||
def _check_input_dim(self, X):
|
||||
|
|
|
|||
|
|
@ -12,6 +12,7 @@ from ...core.parameterization.transformations import Logexp
|
|||
from ...util.caching import Cache_this
|
||||
from ...core.parameterization import variational
|
||||
from psi_comp import linear_psi_comp
|
||||
from ...util.config import *
|
||||
|
||||
class Linear(Kern):
|
||||
"""
|
||||
|
|
@ -224,12 +225,23 @@ class Linear(Kern):
|
|||
AZZA = ZA.T[:, None, :, None] * ZA[None, :, None, :]
|
||||
AZZA = AZZA + AZZA.swapaxes(1, 2)
|
||||
AZZA_2 = AZZA/2.
|
||||
if config.getboolean('parallel', 'openmp'):
|
||||
pragma_string = '#pragma omp parallel for private(m,mm,q,qq,factor,tmp)'
|
||||
header_string = '#include <omp.h>'
|
||||
weave_options = {'headers' : ['<omp.h>'],
|
||||
'extra_compile_args': ['-fopenmp -O3'],
|
||||
'extra_link_args' : ['-lgomp'],
|
||||
'libraries': ['gomp']}
|
||||
else:
|
||||
pragma_string = ''
|
||||
header_string = ''
|
||||
weave_options = {'extra_compile_args': ['-O3']}
|
||||
|
||||
#Using weave, we can exploit the symmetry of this problem:
|
||||
code = """
|
||||
int n, m, mm,q,qq;
|
||||
double factor,tmp;
|
||||
#pragma omp parallel for private(m,mm,q,qq,factor,tmp)
|
||||
%s
|
||||
for(n=0;n<N;n++){
|
||||
for(m=0;m<num_inducing;m++){
|
||||
for(mm=0;mm<=m;mm++){
|
||||
|
|
@ -253,26 +265,36 @@ class Linear(Kern):
|
|||
}
|
||||
}
|
||||
}
|
||||
"""
|
||||
""" % pragma_string
|
||||
support_code = """
|
||||
#include <omp.h>
|
||||
%s
|
||||
#include <math.h>
|
||||
"""
|
||||
weave_options = {'headers' : ['<omp.h>'],
|
||||
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
|
||||
'extra_link_args' : ['-lgomp']}
|
||||
""" % header_string
|
||||
mu = vp.mean
|
||||
N,num_inducing,input_dim,mu = mu.shape[0],Z.shape[0],mu.shape[1],param_to_array(mu)
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
weave.inline(code, support_code=support_code,
|
||||
arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'],
|
||||
type_converters=weave.converters.blitz,**weave_options)
|
||||
|
||||
|
||||
def _weave_dpsi2_dZ(self, dL_dpsi2, Z, vp, target):
|
||||
AZA = self.variances*self._ZAinner(vp, Z)
|
||||
|
||||
if config.getboolean('parallel', 'openmp'):
|
||||
pragma_string = '#pragma omp parallel for private(n,mm,q)'
|
||||
header_string = '#include <omp.h>'
|
||||
weave_options = {'headers' : ['<omp.h>'],
|
||||
'extra_compile_args': ['-fopenmp -O3'],
|
||||
'extra_link_args' : ['-lgomp'],
|
||||
'libraries': ['gomp']}
|
||||
else:
|
||||
pragma_string = ''
|
||||
header_string = ''
|
||||
weave_options = {'extra_compile_args': ['-O3']}
|
||||
|
||||
code="""
|
||||
int n,m,mm,q;
|
||||
#pragma omp parallel for private(n,mm,q)
|
||||
%s
|
||||
for(m=0;m<num_inducing;m++){
|
||||
for(q=0;q<input_dim;q++){
|
||||
for(mm=0;mm<num_inducing;mm++){
|
||||
|
|
@ -282,18 +304,15 @@ class Linear(Kern):
|
|||
}
|
||||
}
|
||||
}
|
||||
"""
|
||||
""" % pragma_string
|
||||
support_code = """
|
||||
#include <omp.h>
|
||||
%s
|
||||
#include <math.h>
|
||||
"""
|
||||
weave_options = {'headers' : ['<omp.h>'],
|
||||
'extra_compile_args': ['-fopenmp -O3'], #-march=native'],
|
||||
'extra_link_args' : ['-lgomp']}
|
||||
""" % header_string
|
||||
|
||||
N,num_inducing,input_dim = vp.mean.shape[0],Z.shape[0],vp.mean.shape[1]
|
||||
mu = param_to_array(vp.mean)
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
weave.inline(code, support_code=support_code,
|
||||
arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'],
|
||||
type_converters=weave.converters.blitz,**weave_options)
|
||||
|
||||
|
|
|
|||
42
GPy/kern/_src/poly.py
Normal file
42
GPy/kern/_src/poly.py
Normal file
|
|
@ -0,0 +1,42 @@
|
|||
# Copyright (c) 2014, James Hensman
|
||||
# Licensed under the BSD 3-clause license (see LICENSE.txt)
|
||||
|
||||
import numpy as np
|
||||
from kern import Kern
|
||||
from ...util.misc import param_to_array
|
||||
from ...core.parameterization import Param
|
||||
from ...core.parameterization.transformations import Logexp
|
||||
class Poly(Kern):
|
||||
"""
|
||||
Polynomial kernel
|
||||
"""
|
||||
|
||||
def __init__(self, input_dim, variance=1., order=3., active_dims=None, name='poly'):
|
||||
super(Poly, self).__init__(input_dim, active_dims, name)
|
||||
self.variance = Param('variance', variance, Logexp())
|
||||
self.add_parameter(self.variance)
|
||||
self.order=order
|
||||
|
||||
def K(self, X, X2=None):
|
||||
return (self._dot_product(X, X2) + 1.)**self.order * self.variance
|
||||
|
||||
def _dot_product(self, X, X2=None):
|
||||
if X2 is None:
|
||||
return np.dot(X, X.T)
|
||||
else:
|
||||
return np.dot(X, X2.T)
|
||||
|
||||
def Kdiag(self, X):
|
||||
return self.variance*(np.square(X).sum(1) + 1.)**self.order
|
||||
|
||||
def update_gradients_full(self, dL_dK, X, X2=None):
|
||||
self.variance.gradient = np.sum(dL_dK * (self._dot_product(X, X2) + 1.)**self.order)
|
||||
|
||||
def update_gradients_diag(self, dL_dKdiag, X):
|
||||
raise NotImplementedError
|
||||
|
||||
def gradients_X(self, dL_dK, X, X2=None):
|
||||
raise NotImplementedError
|
||||
|
||||
def gradients_X_diag(self, dL_dKdiag, X):
|
||||
raise NotImplementedError
|
||||
|
|
@ -10,6 +10,7 @@ from GPy.util.caching import Cache_this
|
|||
from ...core.parameterization import variational
|
||||
from psi_comp import ssrbf_psi_comp
|
||||
from psi_comp.ssrbf_psi_gpucomp import PSICOMP_SSRBF
|
||||
from ...util.config import *
|
||||
|
||||
class RBF(Stationary):
|
||||
"""
|
||||
|
|
@ -243,6 +244,16 @@ class RBF(Stationary):
|
|||
|
||||
@Cache_this(limit=1)
|
||||
def _psi2computations(self, Z, vp):
|
||||
|
||||
if config.getboolean('parallel', 'openmp'):
|
||||
pragma_string = '#pragma omp parallel for private(tmp, exponent_tmp)'
|
||||
header_string = '#include <omp.h>'
|
||||
libraries = ['gomp']
|
||||
else:
|
||||
pragma_string = ''
|
||||
header_string = ''
|
||||
libraries = []
|
||||
|
||||
mu, S = vp.mean, vp.variance
|
||||
|
||||
N, Q = mu.shape
|
||||
|
|
@ -265,8 +276,7 @@ class RBF(Stationary):
|
|||
variance_sq = float(np.square(self.variance))
|
||||
code = """
|
||||
double tmp, exponent_tmp;
|
||||
|
||||
#pragma omp parallel for private(tmp, exponent_tmp)
|
||||
%s
|
||||
for (int n=0; n<N; n++)
|
||||
{
|
||||
for (int m=0; m<M; m++)
|
||||
|
|
@ -290,20 +300,20 @@ class RBF(Stationary):
|
|||
tmp = -Zdist_sq(m,mm,q) - tmp - half_log_denom(n,q);
|
||||
exponent_tmp += tmp;
|
||||
}
|
||||
//compute psi2 by exponontiating
|
||||
//compute psi2 by exponentiating
|
||||
psi2(n,m,mm) = variance_sq * exp(exponent_tmp);
|
||||
psi2(n,mm,m) = psi2(n,m,mm);
|
||||
}
|
||||
}
|
||||
}
|
||||
"""
|
||||
""" % pragma_string
|
||||
|
||||
support_code = """
|
||||
#include <omp.h>
|
||||
%s
|
||||
#include <math.h>
|
||||
"""
|
||||
""" % header_string
|
||||
mu = param_to_array(mu)
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
weave.inline(code, support_code=support_code, libraries=libraries,
|
||||
arg_names=['N', 'M', 'Q', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'denom_l2', 'Zdist_sq', 'half_log_denom', 'psi2', 'variance_sq'],
|
||||
type_converters=weave.converters.blitz, **self.weave_options)
|
||||
|
||||
|
|
@ -315,12 +325,20 @@ class RBF(Stationary):
|
|||
#return 2.*np.einsum( 'ijk,ijk,ijkl,il->l', dL_dpsi2, psi2, Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2, 1./(2.*S + l2))*self.lengthscale
|
||||
|
||||
result = np.zeros(self.input_dim)
|
||||
if config.getboolean('parallel', 'openmp'):
|
||||
pragma_string = '#pragma omp parallel for reduction(+:tmp)'
|
||||
header_string = '#include <omp.h>'
|
||||
libraries = ['gomp']
|
||||
else:
|
||||
pragma_string = ''
|
||||
header_string = ''
|
||||
libraries = []
|
||||
code = """
|
||||
double tmp;
|
||||
for(int q=0; q<Q; q++)
|
||||
{
|
||||
tmp = 0.0;
|
||||
#pragma omp parallel for reduction(+:tmp)
|
||||
%s
|
||||
for(int n=0; n<N; n++)
|
||||
{
|
||||
for(int m=0; m<M; m++)
|
||||
|
|
@ -338,16 +356,16 @@ class RBF(Stationary):
|
|||
result(q) = tmp;
|
||||
}
|
||||
|
||||
"""
|
||||
""" % pragma_string
|
||||
support_code = """
|
||||
#include <omp.h>
|
||||
%s
|
||||
#include <math.h>
|
||||
"""
|
||||
""" % header_string
|
||||
N,Q = S.shape
|
||||
M = psi2.shape[-1]
|
||||
|
||||
S = param_to_array(S)
|
||||
weave.inline(code, support_code=support_code, libraries=['gomp'],
|
||||
weave.inline(code, support_code=support_code, libraries=libraries,
|
||||
arg_names=['psi2', 'dL_dpsi2', 'N', 'M', 'Q', 'mudist_sq', 'l2', 'Zdist_sq', 'S', 'result'],
|
||||
type_converters=weave.converters.blitz, **self.weave_options)
|
||||
|
||||
|
|
|
|||
|
|
@ -192,6 +192,27 @@ class Exponential(Stationary):
|
|||
def dK_dr(self, r):
|
||||
return -0.5*self.K_of_r(r)
|
||||
|
||||
|
||||
class OU(Stationary):
|
||||
"""
|
||||
OU kernel:
|
||||
|
||||
.. math::
|
||||
|
||||
k(r) = \\sigma^2 \exp(- r) \\ \\ \\ \\ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} }
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='OU'):
|
||||
super(OU, self).__init__(input_dim, variance, lengthscale, ARD, active_dims, name)
|
||||
|
||||
def K_of_r(self, r):
|
||||
return self.variance * np.exp(-r)
|
||||
|
||||
def dK_dr(self,r):
|
||||
return -1.*self.variance*np.exp(-r)
|
||||
|
||||
|
||||
class Matern32(Stationary):
|
||||
"""
|
||||
Matern 3/2 kernel:
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue