apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-07 19:12:40 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

modifiying stationary.py

Browse source
This commit is contained in:
James Hensman 2015-04-28 07:40:09 +01:00
parent 59cf0b3c9a
commit 780cf85687
1 changed files with 21 additions and 98 deletions
Download patch file Download diff file Expand all files Collapse all files

119
GPy/kern/_src/stationary.py
View file

@ -9,14 +9,15 @@ from ...util.linalg import tdot
from ... import util from ... import util
import numpy as np import numpy as np
from scipy import integrate from scipy import integrate
from ...util.config import config # for assesing whether to use weave from ...util.config import config # for assesing whether to use cython
from ...util.caching import Cache_this from ...util.caching import Cache_this
import stationary_cython
try: try:
from scipy import weave import stationary_cython
except ImportError: except ImportError:
config.set('weave', 'working', 'False') print('warning: failed to import cython module: falling back to numpy')
config.set('cython', 'working', 'false')
class Stationary(Kern): class Stationary(Kern):
""" """
@ -159,21 +160,11 @@ class Stationary(Kern):
#now the lengthscale gradient(s) #now the lengthscale gradient(s)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
if self.ARD: if self.ARD:
#rinv = self._inv_dis# this is rather high memory? Should we loop instead?t(X, X2)
#d = X[:, None, :] - X2[None, :, :]
#x_xl3 = np.square(d)
#self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3
tmp = dL_dr*self._inv_dist(X, X2) tmp = dL_dr*self._inv_dist(X, X2)
if X2 is None: X2 = X if X2 is None: X2 = X
if config.getboolean('cython', 'working'):
self.lengthscale.gradient = self._lengthscale_grads_cython(tmp, X, X2)
if config.getboolean('weave', 'working'):
try:
self.lengthscale.gradient = self.weave_lengthscale_grads(tmp, X, X2)
except:
print("\n Weave compilation failed. Falling back to (slower) numpy implementation\n")
config.set('weave', 'working', 'False')
self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in range(self.input_dim)])
else: else:
self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in range(self.input_dim)]) self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in range(self.input_dim)])
else: else:
@ -190,29 +181,7 @@ class Stationary(Kern):
dist = self._scaled_dist(X, X2).copy() dist = self._scaled_dist(X, X2).copy()
return 1./np.where(dist != 0., dist, np.inf) return 1./np.where(dist != 0., dist, np.inf)
def weave_lengthscale_grads(self, tmp, X, X2): def _lengthscale_grads_cython(self, tmp, X, X2):
"""Use scipy.weave to compute derivatives wrt the lengthscales"""
N,M = tmp.shape
Q = X.shape[1]
if hasattr(X, 'values'):X = X.values
if hasattr(X2, 'values'):X2 = X2.values
grads = np.zeros(self.input_dim)
code = """
double gradq;
for(int q=0; q<Q; q++){
gradq = 0;
for(int n=0; n<N; n++){
for(int m=0; m<M; m++){
gradq += tmp(n,m)*(X(n,q)-X2(m,q))*(X(n,q)-X2(m,q));
}
}
grads(q) = gradq;
}
"""
weave.inline(code, ['tmp', 'X', 'X2', 'grads', 'N', 'M', 'Q'], type_converters=weave.converters.blitz, support_code="#include <math.h>")
return -grads/self.lengthscale**3
def cython_lengthscale_grads(self, tmp, X, X2):
N,M = tmp.shape N,M = tmp.shape
Q = X.shape[1] Q = X.shape[1]
if hasattr(X, 'values'):X = X.values if hasattr(X, 'values'):X = X.values
@ -221,23 +190,16 @@ class Stationary(Kern):
stationary_cython.lengthscale_grads(N, M, Q, tmp, X, X2, grads) stationary_cython.lengthscale_grads(N, M, Q, tmp, X, X2, grads)
return -grads/self.lengthscale**3 return -grads/self.lengthscale**3
def gradients_X(self, dL_dK, X, X2=None): def gradients_X(self, dL_dK, X, X2=None):
""" """
Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X
""" """
if config.getboolean('weave', 'working'): if config.getboolean('cython', 'working'):
try: return self._gradients_X_cython(dL_dK, X, X2)
return self.gradients_X_weave(dL_dK, X, X2)
except:
print("\n Weave compilation failed. Falling back to (slower) numpy implementation\n")
config.set('weave', 'working', 'False')
return self.gradients_X_(dL_dK, X, X2)
else: else:
return self.gradients_X_(dL_dK, X, X2) return self._gradients_X_pure(dL_dK, X, X2)
def gradients_X_(self, dL_dK, X, X2=None): def _gradients_X_pure(self, dL_dK, X, X2=None):
invdist = self._inv_dist(X, X2) invdist = self._inv_dist(X, X2)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
tmp = invdist*dL_dr tmp = invdist*dL_dr
@ -247,17 +209,15 @@ class Stationary(Kern):
#The high-memory numpy way: #The high-memory numpy way:
#d = X[:, None, :] - X2[None, :, :] #d = X[:, None, :] - X2[None, :, :]
#ret = np.sum(tmp[:,:,None]*d,1)/self.lengthscale**2 #grad = np.sum(tmp[:,:,None]*d,1)/self.lengthscale**2
#the lower memory way with a loop #the lower memory way with a loop
ret = np.empty(X.shape, dtype=np.float64) grad = np.empty(X.shape, dtype=np.float64)
for q in range(self.input_dim): for q in range(self.input_dim):
np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), axis=1, out=ret[:,q]) np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), axis=1, out=grad[:,q])
ret /= self.lengthscale**2 return grad/self.lengthscale**2
return ret def _gradients_X_cython(self, dL_dK, X, X2=None):
def gradients_X_cython(self, dL_dK, X, X2=None):
invdist = self._inv_dist(X, X2) invdist = self._inv_dist(X, X2)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
tmp = invdist*dL_dr tmp = invdist*dL_dr
@ -270,52 +230,15 @@ class Stationary(Kern):
stationary_cython.grad_X(X.shape[0], X.shape[1], X2.shape[0], X, X2, tmp, grad) stationary_cython.grad_X(X.shape[0], X.shape[1], X2.shape[0], X, X2, tmp, grad)
return grad/self.lengthscale**2 return grad/self.lengthscale**2
def gradients_X_weave(self, dL_dK, X, X2=None):
invdist = self._inv_dist(X, X2)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
tmp = invdist*dL_dr
if X2 is None:
tmp = tmp + tmp.T
X2 = X
code = """
int n,m,d;
double retnd;
#pragma omp parallel for private(n,d, retnd, m)
for(d=0;d<D;d++){
for(n=0;n<N;n++){
retnd = 0.0;
for(m=0;m<M;m++){
retnd += tmp(n,m)*(X(n,d)-X2(m,d));
}
ret(n,d) = retnd;
}
}
"""
if hasattr(X, 'values'):X = X.values #remove the GPy wrapping to make passing into weave safe
if hasattr(X2, 'values'):X2 = X2.values
ret = np.zeros(X.shape)
N,D = X.shape
N,M = tmp.shape
from scipy import weave
support_code = """
#include <omp.h>
#include <stdio.h>
"""
weave_options = {'headers' : ['<omp.h>'],
'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
'extra_link_args' : ['-lgomp']}
weave.inline(code, ['ret', 'N', 'D', 'M', 'tmp', 'X', 'X2'], type_converters=weave.converters.blitz, support_code=support_code, **weave_options)
return ret/self.lengthscale**2
def gradients_X_diag(self, dL_dKdiag, X): def gradients_X_diag(self, dL_dKdiag, X):
return np.zeros(X.shape) return np.zeros(X.shape)
def input_sensitivity(self, summarize=True): def input_sensitivity(self, summarize=True):
return self.variance*np.ones(self.input_dim)/self.lengthscale**2 return self.variance*np.ones(self.input_dim)/self.lengthscale**2
class Exponential(Stationary): class Exponential(Stationary):
def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='Exponential'): def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='Exponential'):
super(Exponential, self).__init__(input_dim, variance, lengthscale, ARD, active_dims, name) super(Exponential, self).__init__(input_dim, variance, lengthscale, ARD, active_dims, name)