diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index d317e4b7..1ce5db1a 100644 --- a/GPy/kern/parts/rbf.py +++ b/GPy/kern/parts/rbf.py @@ -4,7 +4,6 @@ from kernpart import Kernpart import numpy as np -import hashlib from scipy import weave from ...util.linalg import tdot from ...util.misc import fast_array_equal @@ -112,7 +111,7 @@ class RBF(Kernpart): } """ num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim - weave.inline(code, arg_names=['num_data','num_inducing','input_dim','X','X2','target','dvardLdK','var_len3'], type_converters=weave.converters.blitz, **self.weave_options) + weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options) else: code = """ int q,i,j; @@ -128,8 +127,8 @@ class RBF(Kernpart): } """ num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim - #[np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)] - weave.inline(code, arg_names=['num_data','num_inducing','input_dim','X','X2','target','dvardLdK','var_len3'], type_converters=weave.converters.blitz, **self.weave_options) + # [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)] + weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options) else: target[1] += (self.variance / self.lengthscale) * np.sum(self._K_dvar * self._K_dist2 * dL_dK) @@ -225,7 +224,7 @@ class RBF(Kernpart): def _K_computations(self, X, X2): if not (fast_array_equal(X, self._X) and fast_array_equal(X2, self._X2) and fast_array_equal(self._params , self._get_params())): self._X = X.copy() - self._params == self._get_params().copy() + self._params = self._get_params().copy() if X2 is None: self._X2 = None X = X / self.lengthscale @@ -241,41 +240,40 @@ class RBF(Kernpart): def _psi_computations(self, Z, mu, S): # here are the "statistics" for psi1 and psi2 if not fast_array_equal(Z, self._Z): - #Z has changed, compute Z specific stuff - self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,Q - self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,Q - self._psi2_Zdist_sq = np.square(self._psi2_Zdist/self.lengthscale) # M,M,Q - self._Z = Z + # Z has changed, compute Z specific stuff + self._psi2_Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q + self._psi2_Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q + self._psi2_Zdist_sq = np.square(self._psi2_Zdist / self.lengthscale) # M,M,Q if not (fast_array_equal(Z, self._Z) and fast_array_equal(mu, self._mu) and fast_array_equal(S, self._S)): - #something's changed. recompute EVERYTHING + # something's changed. recompute EVERYTHING - #psi1 - self._psi1_denom = S[:,None,:]/self.lengthscale2 + 1. - self._psi1_dist = Z[None,:,:]-mu[:,None,:] - self._psi1_dist_sq = np.square(self._psi1_dist)/self.lengthscale2/self._psi1_denom - self._psi1_exponent = -0.5*np.sum(self._psi1_dist_sq+np.log(self._psi1_denom),-1) - self._psi1 = self.variance*np.exp(self._psi1_exponent) + # psi1 + self._psi1_denom = S[:, None, :] / self.lengthscale2 + 1. + self._psi1_dist = Z[None, :, :] - mu[:, None, :] + self._psi1_dist_sq = np.square(self._psi1_dist) / self.lengthscale2 / self._psi1_denom + self._psi1_exponent = -0.5 * np.sum(self._psi1_dist_sq + np.log(self._psi1_denom), -1) + self._psi1 = self.variance * np.exp(self._psi1_exponent) - #psi2 - self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,Q - self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu,self._psi2_Zhat) - #self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q - #self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) - #self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M,Q - self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M,Q + # psi2 + self._psi2_denom = 2.*S[:, None, None, :] / self.lengthscale2 + 1. # N,M,M,Q + self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu, self._psi2_Zhat) + # self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q + # self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) + # self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M,Q + self._psi2 = np.square(self.variance) * np.exp(self._psi2_exponent) # N,M,M,Q - #store matrices for caching - self._Z, self._mu, self._S = Z, mu,S + # store matrices for caching + self._Z, self._mu, self._S = Z, mu, S - def weave_psi2(self,mu,Zhat): - N,input_dim = mu.shape + def weave_psi2(self, mu, Zhat): + N, input_dim = mu.shape num_inducing = Zhat.shape[0] - mudist = np.empty((N,num_inducing,num_inducing,input_dim)) - mudist_sq = np.empty((N,num_inducing,num_inducing,input_dim)) - psi2_exponent = np.zeros((N,num_inducing,num_inducing)) - psi2 = np.empty((N,num_inducing,num_inducing)) + mudist = np.empty((N, num_inducing, num_inducing, input_dim)) + mudist_sq = np.empty((N, num_inducing, num_inducing, input_dim)) + psi2_exponent = np.zeros((N, num_inducing, num_inducing)) + psi2 = np.empty((N, num_inducing, num_inducing)) psi2_Zdist_sq = self._psi2_Zdist_sq _psi2_denom = self._psi2_denom.squeeze().reshape(N, self.input_dim) @@ -325,7 +323,7 @@ class RBF(Kernpart): #include """ weave.inline(code, support_code=support_code, libraries=['gomp'], - arg_names=['N','num_inducing','input_dim','mu','Zhat','mudist_sq','mudist','lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'], + arg_names=['N', 'num_inducing', 'input_dim', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'lengthscale2', '_psi2_denom', 'psi2_Zdist_sq', 'psi2_exponent', 'half_log_psi2_denom', 'psi2', 'variance_sq'], type_converters=weave.converters.blitz, **self.weave_options) return mudist, mudist_sq, psi2_exponent, psi2 diff --git a/GPy/kern/parts/rbf_inv.py b/GPy/kern/parts/rbf_inv.py index f933d395..05605056 100644 --- a/GPy/kern/parts/rbf_inv.py +++ b/GPy/kern/parts/rbf_inv.py @@ -73,7 +73,7 @@ class RBFInv(RBF): self.inv_lengthscale = x[1:] self.inv_lengthscale2 = np.square(self.inv_lengthscale) # TODO: We can rewrite everything with inv_lengthscale and never need to do the below - self.lengthscale = 1./self.inv_lengthscale + self.lengthscale = 1. / self.inv_lengthscale self.lengthscale2 = np.square(self.lengthscale) # reset cached results self._X, self._X2, self._params = np.empty(shape=(3, 1)) @@ -110,7 +110,7 @@ class RBFInv(RBF): } """ num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim - weave.inline(code, arg_names=['num_data','num_inducing','input_dim','X','X2','target','dvardLdK','var_len3', 'len2'], type_converters=weave.converters.blitz, **self.weave_options) + weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3', 'len2'], type_converters=weave.converters.blitz, **self.weave_options) else: code = """ int q,i,j; @@ -126,10 +126,10 @@ class RBFInv(RBF): } """ num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim - #[np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)] - weave.inline(code, arg_names=['num_data','num_inducing','input_dim','X','X2','target','dvardLdK','var_len3', 'len2'], type_converters=weave.converters.blitz, **self.weave_options) + # [np.add(target[1+q:2+q],var_len3[q]*np.sum(dvardLdK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.input_dim)] + weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3', 'len2'], type_converters=weave.converters.blitz, **self.weave_options) else: - target[1] += (self.variance / self.lengthscale) * np.sum(self._K_dvar * self._K_dist2 * dL_dK)*(-self.lengthscale2) + target[1] += (self.variance / self.lengthscale) * np.sum(self._K_dvar * self._K_dist2 * dL_dK) * (-self.lengthscale2) def dK_dX(self, dL_dK, X, X2, target): self._K_computations(X, X2) @@ -159,21 +159,21 @@ class RBFInv(RBF): def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S, target): self._psi_computations(Z, mu, S) - ##d_length = self._psi1[:, :, None] * (-0.5 * ((np.square((self._psi1_dist)/(self.inv_lengthscale * S[:,None,:] + 1))) + ((S[:, None, :])/(self.inv_lengthscale * S[:, None, :] + 1)))) - tmp = 1 + S[:,None,:]*self.inv_lengthscale2 - #inv_len3 = np.power(self.inv_lengthscale,3) - d_length = -(self._psi1[:, :, None] * ((np.square(self._psi1_dist) * self.inv_lengthscale)/(tmp**2) + (S[:,None,:]*self.inv_lengthscale)/(tmp))) + # #d_length = self._psi1[:, :, None] * (-0.5 * ((np.square((self._psi1_dist)/(self.inv_lengthscale * S[:,None,:] + 1))) + ((S[:, None, :])/(self.inv_lengthscale * S[:, None, :] + 1)))) + tmp = 1 + S[:, None, :] * self.inv_lengthscale2 + # inv_len3 = np.power(self.inv_lengthscale,3) + d_length = -(self._psi1[:, :, None] * ((np.square(self._psi1_dist) * self.inv_lengthscale) / (tmp ** 2) + (S[:, None, :] * self.inv_lengthscale) / (tmp))) target[0] += np.sum(dL_dpsi1 * self._psi1 / self.variance) dpsi1_dlength = d_length * dL_dpsi1[:, :, None] if not self.ARD: - target[1] += dpsi1_dlength.sum()#*(-self.lengthscale2) + target[1] += dpsi1_dlength.sum() # *(-self.lengthscale2) else: - target[1:] += dpsi1_dlength.sum(0).sum(0)#*(-self.lengthscale2) - #target[1:] = target[1:]*(-self.lengthscale2) + target[1:] += dpsi1_dlength.sum(0).sum(0) # *(-self.lengthscale2) + # target[1:] = target[1:]*(-self.lengthscale2) def dpsi1_dZ(self, dL_dpsi1, Z, mu, S, target): self._psi_computations(Z, mu, S) - dpsi1_dZ = -self._psi1[:, :, None] * ((self.inv_lengthscale2*self._psi1_dist)/self._psi1_denom) + dpsi1_dZ = -self._psi1[:, :, None] * ((self.inv_lengthscale2 * self._psi1_dist) / self._psi1_denom) target += np.sum(dL_dpsi1[:, :, None] * dpsi1_dZ, 0) def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S, target_mu, target_S): @@ -186,15 +186,15 @@ class RBFInv(RBF): """Shape N,num_inducing,num_inducing,Ntheta""" self._psi_computations(Z, mu, S) d_var = 2.*self._psi2 / self.variance - #d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom) + # d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / self.lengthscale2) / (self.lengthscale * self._psi2_denom) d_length = -2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] * self.inv_lengthscale2) / (self.inv_lengthscale * self._psi2_denom) target[0] += np.sum(dL_dpsi2 * d_var) dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] if not self.ARD: - target[1] += dpsi2_dlength.sum()#*(-self.lengthscale2) + target[1] += dpsi2_dlength.sum() # *(-self.lengthscale2) else: - target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0)#*(-self.lengthscale2) - #target[1:] = target[1:]*(-self.lengthscale2) + target[1:] += dpsi2_dlength.sum(0).sum(0).sum(0) # *(-self.lengthscale2) + # target[1:] = target[1:]*(-self.lengthscale2) def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, target): self._psi_computations(Z, mu, S) @@ -217,7 +217,7 @@ class RBFInv(RBF): def _K_computations(self, X, X2): if not (np.array_equal(X, self._X) and np.array_equal(X2, self._X2) and np.array_equal(self._params , self._get_params())): self._X = X.copy() - self._params == self._get_params().copy() + self._params = self._get_params().copy() if X2 is None: self._X2 = None X = X * self.inv_lengthscale @@ -233,41 +233,40 @@ class RBFInv(RBF): def _psi_computations(self, Z, mu, S): # here are the "statistics" for psi1 and psi2 if not np.array_equal(Z, self._Z): - #Z has changed, compute Z specific stuff - self._psi2_Zhat = 0.5*(Z[:,None,:] +Z[None,:,:]) # M,M,Q - self._psi2_Zdist = 0.5*(Z[:,None,:]-Z[None,:,:]) # M,M,Q - self._psi2_Zdist_sq = np.square(self._psi2_Zdist*self.inv_lengthscale) # M,M,Q - self._Z = Z + # Z has changed, compute Z specific stuff + self._psi2_Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q + self._psi2_Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q + self._psi2_Zdist_sq = np.square(self._psi2_Zdist * self.inv_lengthscale) # M,M,Q if not (np.array_equal(Z, self._Z) and np.array_equal(mu, self._mu) and np.array_equal(S, self._S)): - #something's changed. recompute EVERYTHING + # something's changed. recompute EVERYTHING - #psi1 - self._psi1_denom = S[:,None,:]*self.inv_lengthscale2 + 1. - self._psi1_dist = Z[None,:,:]-mu[:,None,:] - self._psi1_dist_sq = (np.square(self._psi1_dist)*self.inv_lengthscale2)/self._psi1_denom - self._psi1_exponent = -0.5*np.sum(self._psi1_dist_sq+np.log(self._psi1_denom),-1) - self._psi1 = self.variance*np.exp(self._psi1_exponent) + # psi1 + self._psi1_denom = S[:, None, :] * self.inv_lengthscale2 + 1. + self._psi1_dist = Z[None, :, :] - mu[:, None, :] + self._psi1_dist_sq = (np.square(self._psi1_dist) * self.inv_lengthscale2) / self._psi1_denom + self._psi1_exponent = -0.5 * np.sum(self._psi1_dist_sq + np.log(self._psi1_denom), -1) + self._psi1 = self.variance * np.exp(self._psi1_exponent) - #psi2 - self._psi2_denom = 2.*S[:,None,None,:]*self.inv_lengthscale2+1. # N,M,M,Q - self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu,self._psi2_Zhat) - #self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q - #self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) - #self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M,Q - self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M,Q + # psi2 + self._psi2_denom = 2.*S[:, None, None, :] * self.inv_lengthscale2 + 1. # N,M,M,Q + self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu, self._psi2_Zhat) + # self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q + # self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) + # self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M,Q + self._psi2 = np.square(self.variance) * np.exp(self._psi2_exponent) # N,M,M,Q - #store matrices for caching - self._Z, self._mu, self._S = Z, mu,S + # store matrices for caching + self._Z, self._mu, self._S = Z, mu, S - def weave_psi2(self,mu,Zhat): - N,input_dim = mu.shape + def weave_psi2(self, mu, Zhat): + N, input_dim = mu.shape num_inducing = Zhat.shape[0] - mudist = np.empty((N,num_inducing,num_inducing,input_dim)) - mudist_sq = np.empty((N,num_inducing,num_inducing,input_dim)) - psi2_exponent = np.zeros((N,num_inducing,num_inducing)) - psi2 = np.empty((N,num_inducing,num_inducing)) + mudist = np.empty((N, num_inducing, num_inducing, input_dim)) + mudist_sq = np.empty((N, num_inducing, num_inducing, input_dim)) + psi2_exponent = np.zeros((N, num_inducing, num_inducing)) + psi2 = np.empty((N, num_inducing, num_inducing)) psi2_Zdist_sq = self._psi2_Zdist_sq _psi2_denom = self._psi2_denom.squeeze().reshape(N, self.input_dim) @@ -317,7 +316,7 @@ class RBFInv(RBF): #include """ weave.inline(code, support_code=support_code, libraries=['gomp'], - arg_names=['N','num_inducing','input_dim','mu','Zhat','mudist_sq','mudist','inv_lengthscale2','_psi2_denom','psi2_Zdist_sq','psi2_exponent','half_log_psi2_denom','psi2','variance_sq'], + arg_names=['N', 'num_inducing', 'input_dim', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'inv_lengthscale2', '_psi2_denom', 'psi2_Zdist_sq', 'psi2_exponent', 'half_log_psi2_denom', 'psi2', 'variance_sq'], type_converters=weave.converters.blitz, **self.weave_options) return mudist, mudist_sq, psi2_exponent, psi2