From 78ceef01c3c54a0d6aaf8e383c9a8507729d8c18 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 21 Feb 2014 17:06:06 +0000 Subject: [PATCH 1/4] removed materns --- GPy/kern/__init__.py | 1 + GPy/kern/_src/Matern32.py | 139 ---------------------- GPy/kern/_src/Matern52.py | 145 ----------------------- GPy/kern/_src/exponential.py | 129 -------------------- GPy/kern/_src/stationary.py | 221 +++++++++++++++++++++++++++++++++++ GPy/util/__init__.py | 1 + GPy/util/diag.py | 40 ++++--- 7 files changed, 246 insertions(+), 430 deletions(-) delete mode 100644 GPy/kern/_src/Matern32.py delete mode 100644 GPy/kern/_src/Matern52.py delete mode 100644 GPy/kern/_src/exponential.py create mode 100644 GPy/kern/_src/stationary.py diff --git a/GPy/kern/__init__.py b/GPy/kern/__init__.py index 16c13066..e5dc6d35 100644 --- a/GPy/kern/__init__.py +++ b/GPy/kern/__init__.py @@ -3,6 +3,7 @@ from _src.white import White from _src.kern import Kern from _src.linear import Linear from _src.brownian import Brownian +from _src.stationary import Exponential, Matern32, Matern52, ExpQuad #from _src.bias import Bias #import coregionalize #import exponential diff --git a/GPy/kern/_src/Matern32.py b/GPy/kern/_src/Matern32.py deleted file mode 100644 index 08fa452c..00000000 --- a/GPy/kern/_src/Matern32.py +++ /dev/null @@ -1,139 +0,0 @@ -# Copyright (c) 2012, GPy authors (see AUTHORS.txt). -# Licensed under the BSD 3-clause license (see LICENSE.txt) - - -from kernpart import Kernpart -import numpy as np -from scipy import integrate - -class Matern32(Kernpart): - """ - Matern 3/2 kernel: - - .. math:: - - k(r) = \\sigma^2 (1 + \\sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } - - :param input_dim: the number of input dimensions - :type input_dim: int - :param variance: the variance :math:`\sigma^2` - :type variance: float - :param lengthscale: the vector of lengthscale :math:`\ell_i` - :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter) - :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension. - :type ARD: Boolean - :rtype: kernel object - - """ - - def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False): - self.input_dim = input_dim - self.ARD = ARD - if ARD == False: - self.num_params = 2 - self.name = 'Mat32' - if lengthscale is not None: - lengthscale = np.asarray(lengthscale) - assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel" - else: - lengthscale = np.ones(1) - else: - self.num_params = self.input_dim + 1 - self.name = 'Mat32' - if lengthscale is not None: - lengthscale = np.asarray(lengthscale) - assert lengthscale.size == self.input_dim, "bad number of lengthscales" - else: - lengthscale = np.ones(self.input_dim) - self._set_params(np.hstack((variance, lengthscale.flatten()))) - - def _get_params(self): - """return the value of the parameters.""" - return np.hstack((self.variance, self.lengthscale)) - - def _set_params(self, x): - """set the value of the parameters.""" - assert x.size == self.num_params - self.variance = x[0] - self.lengthscale = x[1:] - - def _get_param_names(self): - """return parameter names.""" - if self.num_params == 2: - return ['variance', 'lengthscale'] - else: - return ['variance'] + ['lengthscale_%i' % i for i in range(self.lengthscale.size)] - - def K(self, X, X2, target): - """Compute the covariance matrix between X and X2.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1)) - np.add(self.variance * (1 + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist), target, target) - - def Kdiag(self, X, target): - """Compute the diagonal of the covariance matrix associated to X.""" - np.add(target, self.variance, target) - - def _param_grad_helper(self, dL_dK, X, X2, target): - """derivative of the covariance matrix with respect to the parameters.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1)) - dvar = (1 + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist) - invdist = 1. / np.where(dist != 0., dist, np.inf) - dist2M = np.square(X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 3 - # dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - target[0] += np.sum(dvar * dL_dK) - if self.ARD == True: - dl = (self.variance * 3 * dist * np.exp(-np.sqrt(3.) * dist))[:, :, np.newaxis] * dist2M * invdist[:, :, np.newaxis] - # dl = self.variance*dvar[:,:,None]*dist2M*invdist[:,:,None] - target[1:] += (dl * dL_dK[:, :, None]).sum(0).sum(0) - else: - dl = (self.variance * 3 * dist * np.exp(-np.sqrt(3.) * dist)) * dist2M.sum(-1) * invdist - # dl = self.variance*dvar*dist2M.sum(-1)*invdist - target[1] += np.sum(dl * dL_dK) - - def dKdiag_dtheta(self, dL_dKdiag, X, target): - """derivative of the diagonal of the covariance matrix with respect to the parameters.""" - target[0] += np.sum(dL_dKdiag) - - def gradients_X(self, dL_dK, X, X2, target): - """derivative of the covariance matrix with respect to X.""" - if X2 is None: - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X[None, :, :]) / self.lengthscale), -1))[:, :, None] - ddist_dX = 2*(X[:, None, :] - X[None, :, :]) / self.lengthscale ** 2 / np.where(dist != 0., dist, np.inf) - - else: - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))[:, :, None] - ddist_dX = (X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 2 / np.where(dist != 0., dist, np.inf) - gradients_X = -np.transpose(3 * self.variance * dist * np.exp(-np.sqrt(3) * dist) * ddist_dX, (1, 0, 2)) - target += np.sum(gradients_X * dL_dK.T[:, :, None], 0) - - def dKdiag_dX(self, dL_dKdiag, X, target): - pass - - def Gram_matrix(self, F, F1, F2, lower, upper): - """ - Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1. - - :param F: vector of functions - :type F: np.array - :param F1: vector of derivatives of F - :type F1: np.array - :param F2: vector of second derivatives of F - :type F2: np.array - :param lower,upper: boundaries of the input domain - :type lower,upper: floats - """ - assert self.input_dim == 1 - def L(x, i): - return(3. / self.lengthscale ** 2 * F[i](x) + 2 * np.sqrt(3) / self.lengthscale * F1[i](x) + F2[i](x)) - n = F.shape[0] - G = np.zeros((n, n)) - for i in range(n): - for j in range(i, n): - G[i, j] = G[j, i] = integrate.quad(lambda x : L(x, i) * L(x, j), lower, upper)[0] - Flower = np.array([f(lower) for f in F])[:, None] - F1lower = np.array([f(lower) for f in F1])[:, None] - # print "OLD \n", np.dot(F1lower,F1lower.T), "\n \n" - # return(G) - return(self.lengthscale ** 3 / (12.*np.sqrt(3) * self.variance) * G + 1. / self.variance * np.dot(Flower, Flower.T) + self.lengthscale ** 2 / (3.*self.variance) * np.dot(F1lower, F1lower.T)) diff --git a/GPy/kern/_src/Matern52.py b/GPy/kern/_src/Matern52.py deleted file mode 100644 index 7d36254c..00000000 --- a/GPy/kern/_src/Matern52.py +++ /dev/null @@ -1,145 +0,0 @@ -# Copyright (c) 2012, GPy authors (see AUTHORS.txt). -# Licensed under the BSD 3-clause license (see LICENSE.txt) - - -from kernpart import Kernpart -import numpy as np -import hashlib -from scipy import integrate - -class Matern52(Kernpart): - """ - Matern 5/2 kernel: - - .. math:: - - k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } - - :param input_dim: the number of input dimensions - :type input_dim: int - :param variance: the variance :math:`\sigma^2` - :type variance: float - :param lengthscale: the vector of lengthscale :math:`\ell_i` - :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter) - :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension. - :type ARD: Boolean - :rtype: kernel object - - """ - def __init__(self,input_dim,variance=1.,lengthscale=None,ARD=False): - self.input_dim = input_dim - self.ARD = ARD - if ARD == False: - self.num_params = 2 - self.name = 'Mat52' - if lengthscale is not None: - lengthscale = np.asarray(lengthscale) - assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel" - else: - lengthscale = np.ones(1) - else: - self.num_params = self.input_dim + 1 - self.name = 'Mat52' - if lengthscale is not None: - lengthscale = np.asarray(lengthscale) - assert lengthscale.size == self.input_dim, "bad number of lengthscales" - else: - lengthscale = np.ones(self.input_dim) - self._set_params(np.hstack((variance,lengthscale.flatten()))) - - def _get_params(self): - """return the value of the parameters.""" - return np.hstack((self.variance,self.lengthscale)) - - def _set_params(self,x): - """set the value of the parameters.""" - assert x.size == self.num_params - self.variance = x[0] - self.lengthscale = x[1:] - - def _get_param_names(self): - """return parameter names.""" - if self.num_params == 2: - return ['variance','lengthscale'] - else: - return ['variance']+['lengthscale_%i'%i for i in range(self.lengthscale.size)] - - def K(self,X,X2,target): - """Compute the covariance matrix between X and X2.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) - np.add(self.variance*(1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist), target,target) - - def Kdiag(self,X,target): - """Compute the diagonal of the covariance matrix associated to X.""" - np.add(target,self.variance,target) - - def _param_grad_helper(self,dL_dK,X,X2,target): - """derivative of the covariance matrix with respect to the parameters.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1)) - invdist = 1./np.where(dist!=0.,dist,np.inf) - dist2M = np.square(X[:,None,:]-X2[None,:,:])/self.lengthscale**3 - dvar = (1+np.sqrt(5.)*dist+5./3*dist**2)*np.exp(-np.sqrt(5.)*dist) - dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - target[0] += np.sum(dvar*dL_dK) - if self.ARD: - dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist))[:,:,np.newaxis] * dist2M*invdist[:,:,np.newaxis] - target[1:] += (dl*dL_dK[:,:,None]).sum(0).sum(0) - else: - dl = (self.variance * 5./3 * dist * (1 + np.sqrt(5.)*dist ) * np.exp(-np.sqrt(5.)*dist)) * dist2M.sum(-1)*invdist - #dl = (self.variance* 3 * dist * np.exp(-np.sqrt(3.)*dist)) * dist2M.sum(-1)*invdist - target[1] += np.sum(dl*dL_dK) - - def dKdiag_dtheta(self,dL_dKdiag,X,target): - """derivative of the diagonal of the covariance matrix with respect to the parameters.""" - target[0] += np.sum(dL_dKdiag) - - def gradients_X(self,dL_dK,X,X2,target): - """derivative of the covariance matrix with respect to X.""" - if X2 is None: - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X[None,:,:])/self.lengthscale),-1))[:,:,None] - ddist_dX = 2*(X[:,None,:]-X[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf) - else: - dist = np.sqrt(np.sum(np.square((X[:,None,:]-X2[None,:,:])/self.lengthscale),-1))[:,:,None] - ddist_dX = (X[:,None,:]-X2[None,:,:])/self.lengthscale**2/np.where(dist!=0.,dist,np.inf) - gradients_X = - np.transpose(self.variance*5./3*dist*(1+np.sqrt(5)*dist)*np.exp(-np.sqrt(5)*dist)*ddist_dX,(1,0,2)) - target += np.sum(gradients_X*dL_dK.T[:,:,None],0) - - def dKdiag_dX(self,dL_dKdiag,X,target): - pass - - def Gram_matrix(self,F,F1,F2,F3,lower,upper): - """ - Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1. - - :param F: vector of functions - :type F: np.array - :param F1: vector of derivatives of F - :type F1: np.array - :param F2: vector of second derivatives of F - :type F2: np.array - :param F3: vector of third derivatives of F - :type F3: np.array - :param lower,upper: boundaries of the input domain - :type lower,upper: floats - """ - assert self.input_dim == 1 - def L(x,i): - return(5*np.sqrt(5)/self.lengthscale**3*F[i](x) + 15./self.lengthscale**2*F1[i](x)+ 3*np.sqrt(5)/self.lengthscale*F2[i](x) + F3[i](x)) - n = F.shape[0] - G = np.zeros((n,n)) - for i in range(n): - for j in range(i,n): - G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0] - G_coef = 3.*self.lengthscale**5/(400*np.sqrt(5)) - Flower = np.array([f(lower) for f in F])[:,None] - F1lower = np.array([f(lower) for f in F1])[:,None] - F2lower = np.array([f(lower) for f in F2])[:,None] - orig = 9./8*np.dot(Flower,Flower.T) + 9.*self.lengthscale**4/200*np.dot(F2lower,F2lower.T) - orig2 = 3./5*self.lengthscale**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T)) - return(1./self.variance* (G_coef*G + orig + orig2)) - - - diff --git a/GPy/kern/_src/exponential.py b/GPy/kern/_src/exponential.py deleted file mode 100644 index 372d4d9b..00000000 --- a/GPy/kern/_src/exponential.py +++ /dev/null @@ -1,129 +0,0 @@ -# Copyright (c) 2012, GPy authors (see AUTHORS.txt). -# Licensed under the BSD 3-clause license (see LICENSE.txt) - - -from kernpart import Kernpart -import numpy as np -from scipy import integrate - -class Exponential(Kernpart): - """ - Exponential kernel (aka Ornstein-Uhlenbeck or Matern 1/2) - - .. 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} } - - :param input_dim: the number of input dimensions - :type input_dim: int - :param variance: the variance :math:`\sigma^2` - :type variance: float - :param lengthscale: the vector of lengthscale :math:`\ell_i` - :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter) - :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension. - :type ARD: Boolean - :param name: the name of the kernel - :rtype: kernel object - - """ - def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='exp'): - self.input_dim = input_dim - self.ARD = ARD - self.variance = variance - self.name = name - if ARD == False: - self.num_params = 2 - if lengthscale is not None: - lengthscale = np.asarray(lengthscale) - assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel" - else: - lengthscale = np.ones(1) - else: - self.num_params = self.input_dim + 1 - if lengthscale is not None: - lengthscale = np.asarray(lengthscale) - assert lengthscale.size == self.input_dim, "bad number of lengthscales" - else: - lengthscale = np.ones(self.input_dim) - #self._set_params(np.hstack((variance, lengthscale.flatten()))) - self.set_as_parameter('variance', 'lengthscale') - -# def _get_params(self): -# """return the value of the parameters.""" -# return np.hstack((self.variance, self.lengthscale)) -# -# def _set_params(self, x): -# """set the value of the parameters.""" -# assert x.size == self.num_params -# self.variance = x[0] -# self.lengthscale = x[1:] -# -# def _get_param_names(self): -# """return parameter names.""" -# if self.num_params == 2: -# return ['variance', 'lengthscale'] -# else: -# return ['variance'] + ['lengthscale_%i' % i for i in range(self.lengthscale.size)] - - def K(self, X, X2, target): - """Compute the covariance matrix between X and X2.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1)) - np.add(self.variance * np.exp(-dist), target, target) - - def Kdiag(self, X, target): - """Compute the diagonal of the covariance matrix associated to X.""" - np.add(target, self.variance, target) - - def _param_grad_helper(self, dL_dK, X, X2, target): - """derivative of the covariance matrix with respect to the parameters.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1)) - invdist = 1. / np.where(dist != 0., dist, np.inf) - dist2M = np.square(X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 3 - dvar = np.exp(-dist) - target[0] += np.sum(dvar * dL_dK) - if self.ARD == True: - dl = self.variance * dvar[:, :, None] * dist2M * invdist[:, :, None] - target[1:] += (dl * dL_dK[:, :, None]).sum(0).sum(0) - else: - dl = self.variance * dvar * dist2M.sum(-1) * invdist - target[1] += np.sum(dl * dL_dK) - - def dKdiag_dtheta(self, dL_dKdiag, X, target): - """derivative of the diagonal of the covariance matrix with respect to the parameters.""" - # NB: derivative of diagonal elements wrt lengthscale is 0 - target[0] += np.sum(dL_dKdiag) - - def gradients_X(self, dL_dK, X, X2, target): - """derivative of the covariance matrix with respect to X.""" - if X2 is None: X2 = X - dist = np.sqrt(np.sum(np.square((X[:, None, :] - X2[None, :, :]) / self.lengthscale), -1))[:, :, None] - ddist_dX = (X[:, None, :] - X2[None, :, :]) / self.lengthscale ** 2 / np.where(dist != 0., dist, np.inf) - gradients_X = -np.transpose(self.variance * np.exp(-dist) * ddist_dX, (1, 0, 2)) - target += np.sum(gradients_X * dL_dK.T[:, :, None], 0) - - def dKdiag_dX(self, dL_dKdiag, X, target): - pass - - def Gram_matrix(self, F, F1, lower, upper): - """ - Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1. - - :param F: vector of functions - :type F: np.array - :param F1: vector of derivatives of F - :type F1: np.array - :param lower,upper: boundaries of the input domain - :type lower,upper: floats - """ - assert self.input_dim == 1 - def L(x, i): - return(1. / self.lengthscale * F[i](x) + F1[i](x)) - n = F.shape[0] - G = np.zeros((n, n)) - for i in range(n): - for j in range(i, n): - G[i, j] = G[j, i] = integrate.quad(lambda x : L(x, i) * L(x, j), lower, upper)[0] - Flower = np.array([f(lower) for f in F])[:, None] - return(self.lengthscale / 2. / self.variance * G + 1. / self.variance * np.dot(Flower, Flower.T)) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py new file mode 100644 index 00000000..aaa534ac --- /dev/null +++ b/GPy/kern/_src/stationary.py @@ -0,0 +1,221 @@ +# 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 +from ... import util +import numpy as np +from scipy import integrate + +class Stationary(Kern): + def __init__(self, input_dim, variance, lengthscale, ARD, name): + super(Stationary, self).__init__(input_dim, name) + self.ARD = ARD + if not ARD: + if lengthscale is None: + lengthscale = np.ones(1) + else: + lengthscale = np.asarray(lengthscale) + assert lengthscale.size == 1 "Only lengthscale needed for non-ARD kernel" + else: + if lengthscale is not None: + lengthscale = np.asarray(lengthscale) + assert lengthscale.size in [1, input_dim], "Bad lengthscales" + if lengthscale.size != input_dim: + lengthscale = np.ones(input_dim)*lengthscale + else: + lengthscale = np.ones(self.input_dim) + self.lengthscale = Param('lengthscale', lengthscale, Logexp()) + self.variance = Param('variance', variance, Logexp()) + assert self.variance.size==1 + self.add_parameters(self.variance, self.lengthscale) + + def _dist(self, X, X2): + if X2 is None: + X2 = X + return X[:, None, :] - X2[None, :, :] + + def _scaled_dist(self, X, X2=None): + return np.sqrt(np.sum(np.square(self._dist(X, X2) / self.lengthscale), -1)) + + def Kdiag(self, X): + ret = np.empty(X.shape[0]) + ret[:] = self.variance + return ret + + def update_gradients_diag(self, dL_dKdiag, X): + self.variance.gradient = np.sum(dL_dKdiag) + self.lengthscale.gradient = 0. + + def gradients_X_diag(self, dL_dKdiag, X): + return np.zeros(X.shape) + + def update_gradients_full(self, dL_dK, X, X2=None): + K = self.K(X, X2) + self.variance.gradient = np.sum(K * dL_dK)/self.variance + + rinv = self._inv_dist(X, X2) + dL_dr = self.dK_dr(X, X2) * dL_dK + x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 + + if self.ARD: + self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0) + else: + self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum() + + def _inv_dist(self, X, X2=None): + dist = self._scaled_dist(X, X2) + if X2 is None: + nondiag = util.diag.offdiag_view(dist) + nondiag[:] = 1./nondiag + return dist + else: + return 1./np.where(dist != 0., dist, np.inf) + + def gradients_X(self, dL_dK, X, X2=None): + dL_dr = self.dK_dr(X, X2) * dL_dK + invdist = self._inv_dist(X, X2) + ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2 + if X2 is None: + ret *= 2. + return ret + + + + +class Exponential(Stationary): + def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Exponential'): + super(Exponential, self).__init__(input_dim, variance, lengthscale, ARD, name) + + def K(self, X, X2=None): + dist = self._scaled_dist(X, X2) + return self.variance * np.exp(-0.5 * dist) + + def dK_dr(self, X, X2): + return -0.5*self.K(X, X2) + +class Matern32(Stationary): + """ + Matern 3/2 kernel: + + .. math:: + + k(r) = \\sigma^2 (1 + \\sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } + + :param input_dim: the number of input dimensions + :type input_dim: int + :param variance: the variance :math:`\sigma^2` + :type variance: float + :param lengthscale: the vector of lengthscale :math:`\ell_i` + :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter) + :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension. + :type ARD: Boolean + :rtype: kernel object + + """ + + def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat32'): + super(Matern32, self).__init__(input_dim, variance, lengthscale, ARD, name) + + def K(self, X, X2=None): + dist = self._scaled_dist(X, X2) + return self.variance * (1. + np.sqrt(3.) * dist) * np.exp(-np.sqrt(3.) * dist) + + def dK_dr(self, X, X2): + dist = self._scaled_dist(X, X2) + return -3.*self.variance*dist*np.exp(-np.sqrt(3.)*dist) + + def Gram_matrix(self, F, F1, F2, lower, upper): + """ + Return the Gram matrix of the vector of functions F with respect to the + RKHS norm. The use of this function is limited to input_dim=1. + + :param F: vector of functions + :type F: np.array + :param F1: vector of derivatives of F + :type F1: np.array + :param F2: vector of second derivatives of F + :type F2: np.array + :param lower,upper: boundaries of the input domain + :type lower,upper: floats + """ + assert self.input_dim == 1 + def L(x, i): + return(3. / self.lengthscale ** 2 * F[i](x) + 2 * np.sqrt(3) / self.lengthscale * F1[i](x) + F2[i](x)) + n = F.shape[0] + G = np.zeros((n, n)) + for i in range(n): + for j in range(i, n): + G[i, j] = G[j, i] = integrate.quad(lambda x : L(x, i) * L(x, j), lower, upper)[0] + Flower = np.array([f(lower) for f in F])[:, None] + F1lower = np.array([f(lower) for f in F1])[:, None] + return(self.lengthscale ** 3 / (12.*np.sqrt(3) * self.variance) * G + 1. / self.variance * np.dot(Flower, Flower.T) + self.lengthscale ** 2 / (3.*self.variance) * np.dot(F1lower, F1lower.T)) + + +class Matern52(Stationary): + """ + Matern 5/2 kernel: + + .. math:: + + k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } + """ + + def K(self, X, X2=None): + r = self._scaled_dist(X, X2) + return self.variance*(1+np.sqrt(5.)*r+5./3*r**2)*np.exp(-np.sqrt(5.)*r) + + def dK_dr(self, X, X2): + r = self._scaled_dist(X, X2) + return self.variance*(10./3*r -5.*r -5.*np.sqrt(5.)/3*r**2)*np.exp(-np.sqrt(5.)*r) + + def Gram_matrix(self,F,F1,F2,F3,lower,upper): + """ + Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1. + + :param F: vector of functions + :type F: np.array + :param F1: vector of derivatives of F + :type F1: np.array + :param F2: vector of second derivatives of F + :type F2: np.array + :param F3: vector of third derivatives of F + :type F3: np.array + :param lower,upper: boundaries of the input domain + :type lower,upper: floats + """ + assert self.input_dim == 1 + def L(x,i): + return(5*np.sqrt(5)/self.lengthscale**3*F[i](x) + 15./self.lengthscale**2*F1[i](x)+ 3*np.sqrt(5)/self.lengthscale*F2[i](x) + F3[i](x)) + n = F.shape[0] + G = np.zeros((n,n)) + for i in range(n): + for j in range(i,n): + G[i,j] = G[j,i] = integrate.quad(lambda x : L(x,i)*L(x,j),lower,upper)[0] + G_coef = 3.*self.lengthscale**5/(400*np.sqrt(5)) + Flower = np.array([f(lower) for f in F])[:,None] + F1lower = np.array([f(lower) for f in F1])[:,None] + F2lower = np.array([f(lower) for f in F2])[:,None] + orig = 9./8*np.dot(Flower,Flower.T) + 9.*self.lengthscale**4/200*np.dot(F2lower,F2lower.T) + orig2 = 3./5*self.lengthscale**2 * ( np.dot(F1lower,F1lower.T) + 1./8*np.dot(Flower,F2lower.T) + 1./8*np.dot(F2lower,Flower.T)) + return(1./self.variance* (G_coef*G + orig + orig2)) + + + + +class ExpQuad(Stationary): + def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='ExpQuad'): + super(ExpQuad, self).__init__(input_dim, variance, lengthscale, ARD, name) + + def K(self, X, X2=None): + r = self._scaled_dist(X, X2) + return self.variance * np.exp(-0.5 * r**2) + + def dK_dr(self, X, X2): + dist = self._scaled_dist(X, X2) + return -dist*self.K(X, X2) + + + diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py index c10fea4c..f93bb0ec 100644 --- a/GPy/util/__init__.py +++ b/GPy/util/__init__.py @@ -12,6 +12,7 @@ import decorators import classification import subarray_and_sorting import caching +import diag try: import sympy diff --git a/GPy/util/diag.py b/GPy/util/diag.py index 3d6b4dc9..3044ed54 100644 --- a/GPy/util/diag.py +++ b/GPy/util/diag.py @@ -11,14 +11,14 @@ import numpy as np def view(A, offset=0): """ Get a view on the diagonal elements of a 2D array. - - This is actually a view (!) on the diagonal of the array, so you can + + This is actually a view (!) on the diagonal of the array, so you can in-place adjust the view. - + :param :class:`ndarray` A: 2 dimensional numpy array :param int offset: view offset to give back (negative entries allowed) :rtype: :class:`ndarray` view of diag(A) - + >>> import numpy as np >>> X = np.arange(9).reshape(3,3) >>> view(X) @@ -36,7 +36,7 @@ def view(A, offset=0): """ from numpy.lib.stride_tricks import as_strided assert A.ndim == 2, "only implemented for 2 dimensions" - assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!" + assert A.shape[0] == A.shape[1], "attempting to get the view of non-square matrix?!" if offset > 0: return as_strided(A[0, offset:], shape=(A.shape[0] - offset, ), strides=((A.shape[0]+1)*A.itemsize, )) elif offset < 0: @@ -44,6 +44,12 @@ def view(A, offset=0): else: return as_strided(A, shape=(A.shape[0], ), strides=((A.shape[0]+1)*A.itemsize, )) +def offdiag_view(A, offset=0): + from numpy.lib.stride_tricks import as_strided + assert A.ndim == 2, "only implemented for 2 dimensions" + Af = as_strided(A, shape=(A.size,), strides=(A.itemsize,)) + return as_strided(Af[(1+offset):], shape=(A.shape[0]-1, A.shape[1]), strides=(A.strides[0] + A.itemsize, A.strides[1])) + def _diag_ufunc(A,b,offset,func): dA = view(A, offset); func(dA,b,dA) return A @@ -51,11 +57,11 @@ def _diag_ufunc(A,b,offset,func): def times(A, b, offset=0): """ Times the view of A with b in place (!). - Returns modified A + Returns modified A Broadcasting is allowed, thus b can be scalar. - + if offset is not zero, make sure b is of right shape! - + :param ndarray A: 2 dimensional array :param ndarray-like b: either one dimensional or scalar :param int offset: same as in view. @@ -67,11 +73,11 @@ multiply = times def divide(A, b, offset=0): """ Divide the view of A by b in place (!). - Returns modified A + Returns modified A Broadcasting is allowed, thus b can be scalar. - + if offset is not zero, make sure b is of right shape! - + :param ndarray A: 2 dimensional array :param ndarray-like b: either one dimensional or scalar :param int offset: same as in view. @@ -84,9 +90,9 @@ def add(A, b, offset=0): Add b to the view of A in place (!). Returns modified A. Broadcasting is allowed, thus b can be scalar. - + if offset is not zero, make sure b is of right shape! - + :param ndarray A: 2 dimensional array :param ndarray-like b: either one dimensional or scalar :param int offset: same as in view. @@ -99,16 +105,16 @@ def subtract(A, b, offset=0): Subtract b from the view of A in place (!). Returns modified A. Broadcasting is allowed, thus b can be scalar. - + if offset is not zero, make sure b is of right shape! - + :param ndarray A: 2 dimensional array :param ndarray-like b: either one dimensional or scalar :param int offset: same as in view. :rtype: view of A, which is adjusted inplace """ return _diag_ufunc(A, b, offset, np.subtract) - + if __name__ == '__main__': import doctest - doctest.testmod() \ No newline at end of file + doctest.testmod() From fddc663f286e94feef218c74a4f555903e097bee Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 21 Feb 2014 17:32:40 +0000 Subject: [PATCH 2/4] working on coregionalize --- GPy/kern/_src/coregionalize.py | 91 ++++++++++++++++------------------ 1 file changed, 44 insertions(+), 47 deletions(-) diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py index 69fc27ef..0d99ce21 100644 --- a/GPy/kern/_src/coregionalize.py +++ b/GPy/kern/_src/coregionalize.py @@ -5,6 +5,7 @@ from kern import Kern import numpy as np from scipy import weave from ...core.parameterization import Param +from ...core.parameterization.transformations import Logexp class Coregionalize(Kern): """ @@ -20,7 +21,7 @@ class Coregionalize(Kern): k_2(x, y)=\mathbf{B} k(x, y) it is obtained as the tensor product between a covariance function - k(x,y) and B. + k(x, y) and B. :param output_dim: number of outputs to coregionalize :type output_dim: int @@ -29,7 +30,7 @@ class Coregionalize(Kern): :param W: a low rank matrix that determines the correlations between the different outputs, together with kappa it forms the coregionalization matrix B :type W: numpy array of dimensionality (num_outpus, W_columns) :param kappa: a vector which allows the outputs to behave independently - :type kappa: numpy array of dimensionality (output_dim,) + :type kappa: numpy array of dimensionality (output_dim, ) .. note: see coregionalization examples in GPy.examples.regression for some usage. """ @@ -37,18 +38,18 @@ class Coregionalize(Kern): super(Coregionalize, self).__init__(input_dim=1, name=name) self.output_dim = output_dim self.rank = rank - if self.rank>output_dim-1: + if self.rank>output_dim: print("Warning: Unusual choice of rank, it should normally be less than the output_dim.") if W is None: - W = 0.5*np.random.randn(self.output_dim,self.rank)/np.sqrt(self.rank) + W = 0.5*np.random.randn(self.output_dim, self.rank)/np.sqrt(self.rank) else: - assert W.shape==(self.output_dim,self.rank) - self.W = Param('W',W) + assert W.shape==(self.output_dim, self.rank) + self.W = Param('W', W) if kappa is None: kappa = 0.5*np.ones(self.output_dim) else: - assert kappa.shape==(self.output_dim,) - self.kappa = Param('kappa', kappa) + assert kappa.shape==(self.output_dim, ) + self.kappa = Param('kappa', kappa, Logexp()) self.add_parameters(self.W, self.kappa) self.parameters_changed() @@ -56,54 +57,58 @@ class Coregionalize(Kern): def parameters_changed(self): self.B = np.dot(self.W, self.W.T) + np.diag(self.kappa) - def K(self,index,index2,target): - index = np.asarray(index,dtype=np.int) + def K(self, X, X2=None): + index = np.asarray(X, dtype=np.int) #here's the old code (numpy) #if index2 is None: #index2 = index #else: - #index2 = np.asarray(index2,dtype=np.int) + #index2 = np.asarray(index2, dtype=np.int) #false_target = target.copy() - #ii,jj = np.meshgrid(index,index2) - #ii,jj = ii.T, jj.T - #false_target += self.B[ii,jj] + #ii, jj = np.meshgrid(index, index2) + #ii, jj = ii.T, jj.T + #false_target += self.B[ii, jj] - if index2 is None: + + if X2 is None: + target = np.empty((X.shape[0], X.shape[0]), dtype=np.float64) code=""" for(int i=0;i Date: Fri, 21 Feb 2014 17:39:02 +0000 Subject: [PATCH 3/4] tidying --- GPy/kern/_src/coregionalize.py | 4 ++++ GPy/kern/_src/stationary.py | 16 +++------------- 2 files changed, 7 insertions(+), 13 deletions(-) diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py index 0d99ce21..74cd2a1d 100644 --- a/GPy/kern/_src/coregionalize.py +++ b/GPy/kern/_src/coregionalize.py @@ -135,3 +135,7 @@ class Coregionalize(Kern): def gradients_X(self, dL_dK, X, X2=None): return np.zeros(X.shape) + + def gradients_X_diag(self, dL_dKdiag, X): + return np.zeros(X.shape) + diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index aaa534ac..7cc2e695 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -49,9 +49,6 @@ class Stationary(Kern): self.variance.gradient = np.sum(dL_dKdiag) self.lengthscale.gradient = 0. - def gradients_X_diag(self, dL_dKdiag, X): - return np.zeros(X.shape) - def update_gradients_full(self, dL_dK, X, X2=None): K = self.K(X, X2) self.variance.gradient = np.sum(K * dL_dK)/self.variance @@ -82,6 +79,9 @@ class Stationary(Kern): ret *= 2. return ret + def gradients_X_diag(self, dL_dKdiag, X): + return np.zeros(X.shape) + @@ -104,16 +104,6 @@ class Matern32(Stationary): k(r) = \\sigma^2 (1 + \\sqrt{3} r) \exp(- \sqrt{3} r) \\ \\ \\ \\ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } - :param input_dim: the number of input dimensions - :type input_dim: int - :param variance: the variance :math:`\sigma^2` - :type variance: float - :param lengthscale: the vector of lengthscale :math:`\ell_i` - :type lengthscale: array or list of the appropriate size (or float if there is only one lengthscale parameter) - :param ARD: Auto Relevance Determination. If equal to "False", the kernel is isotropic (ie. one single lengthscale parameter \ell), otherwise there is one lengthscale parameter per dimension. - :type ARD: Boolean - :rtype: kernel object - """ def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat32'): From 659643038fe0c6937c69e48cf12c4efd32e41edf Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 21 Feb 2014 17:53:44 +0000 Subject: [PATCH 4/4] parameterized now supports deleting of parameters --- GPy/core/model.py | 17 ++++------ GPy/core/parameterization/index_operations.py | 18 +++++++--- GPy/core/parameterization/parameter_core.py | 12 ++++--- GPy/core/parameterization/parameterized.py | 18 ++++++---- GPy/examples/dimensionality_reduction.py | 19 +++++++---- .../latent_function_inference/var_dtc.py | 10 +++--- GPy/kern/_src/kern.py | 2 +- .../matplot_dep/dim_reduction_plots.py | 20 +++++------ GPy/plotting/matplot_dep/kernel_plots.py | 34 +++++++++---------- GPy/plotting/matplot_dep/models_plots.py | 23 +++++++------ GPy/testing/index_operations_tests.py | 7 ++++ GPy/testing/parameterized_tests.py | 16 ++++----- 12 files changed, 113 insertions(+), 83 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index c067d51d..21bcf0c7 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -485,20 +485,17 @@ class Model(Parameterized): if not hasattr(self, 'kern'): raise ValueError, "this model has no kernel" - k = [p for p in self.kern._parameters_ if hasattr(p, "ARD") and p.ARD] - if (not len(k) == 1): - raise ValueError, "cannot determine sensitivity for this kernel" - k = k[0] - from ..kern.parts.rbf import RBF - from ..kern.parts.rbf_inv import RBFInv - from ..kern.parts.linear import Linear + k = self.kern#[p for p in self.kern._parameters_ if hasattr(p, "ARD") and p.ARD] + from ..kern import RBF, Linear#, RBFInv + if isinstance(k, RBF): return 1. / k.lengthscale - elif isinstance(k, RBFInv): - return k.inv_lengthscale + #elif isinstance(k, RBFInv): + # return k.inv_lengthscale elif isinstance(k, Linear): return k.variances - + else: + raise ValueError, "cannot determine sensitivity for this kernel" def pseudo_EM(self, stop_crit=.1, **kwargs): """ diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index bfd0bf21..b5399741 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -83,11 +83,21 @@ class ParameterIndexOperations(object): def iterproperties(self): return self._properties.iterkeys() - def shift(self, start, size): + def shift_right(self, start, size): for ind in self.iterindices(): toshift = ind>=start - if toshift.size > 0: - ind[toshift] += size + ind[toshift] += size + + def shift_left(self, start, size): + for v, ind in self.items(): + todelete = (ind>=start) * (ind=start + if toshift.size != 0: + ind[toshift] -= size + if ind.size != 0: self._properties[v] = ind + else: del self._properties[v] def clear(self): self._properties.clear() @@ -183,7 +193,7 @@ class ParameterIndexOperationsView(object): yield i - def shift(self, start, size): + def shift_right(self, start, size): raise NotImplementedError, 'Shifting only supported in original ParamIndexOperations' diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 45b57eab..c2c8a05a 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -390,6 +390,7 @@ class Parameterizable(Constrainable): import copy from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView from .array_core import ParamList + dc = dict() for k, v in self.__dict__.iteritems(): if k not in ['_direct_parent_', '_parameters_', '_parent_index_'] + self.parameter_names(): @@ -399,18 +400,21 @@ class Parameterizable(Constrainable): dc[k] = copy.deepcopy(v) if k == '_parameters_': params = [p.copy() for p in v] - # dc = copy.deepcopy(self.__dict__) + dc['_direct_parent_'] = None dc['_parent_index_'] = None dc['_parameters_'] = ParamList() + dc['constraints'].clear() + dc['priors'].clear() + dc['size'] = 0 + s = self.__new__(self.__class__) s.__dict__ = dc - # import ipdb;ipdb.set_trace() + for p in params: s.add_parameter(p) - # dc._notify_parent_change() + return s - # return copy.deepcopy(self) def _notify_parameters_changed(self): self.parameters_changed() diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 177cc217..d463ed43 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -87,8 +87,8 @@ class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable): self._parameters_.append(param) else: start = sum(p.size for p in self._parameters_[:index]) - self.constraints.shift(start, param.size) - self.priors.shift(start, param.size) + self.constraints.shift_right(start, param.size) + self.priors.shift_right(start, param.size) self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) @@ -113,15 +113,19 @@ class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable): """ if not param in self._parameters_: raise RuntimeError, "Parameter {} does not belong to this object, remove parameters directly from their respective parents".format(param._short()) - del self._parameters_[param._parent_index_] + + start = sum([p.size for p in self._parameters_[:param._parent_index_]]) + self._remove_parameter_name(param) self.size -= param.size + del self._parameters_[param._parent_index_] param._disconnect_parent() - self._remove_parameter_name(param) - - #self._notify_parent_change() + self.constraints.shift_left(start, param.size) self._connect_fixes() - + self._connect_parameters() + self._notify_parent_change() + + def _connect_parameters(self): # connect parameterlist to this parameterized object # This just sets up the right connection for the params objects diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index c8e79e6c..3ba54d34 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -74,7 +74,7 @@ def gplvm_oil_100(optimize=True, verbose=1, plot=True): data = GPy.util.datasets.oil_100() Y = data['X'] # create simple GP model - kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.bias(6) + kernel = GPy.kern.RBF(6, ARD=True) + GPy.kern.Bias(6) m = GPy.models.GPLVM(Y, 6, kernel=kernel) m.data_labels = data['Y'].argmax(axis=1) if optimize: m.optimize('scg', messages=verbose) @@ -190,17 +190,22 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False): _np.random.seed(1234) x = _np.linspace(0, 4 * _np.pi, N)[:, None] - s1 = _np.vectorize(lambda x: _np.sin(x)) + s1 = _np.vectorize(lambda x: -_np.sin(x)) s2 = _np.vectorize(lambda x: _np.cos(x)) s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x))) - sS = _np.vectorize(lambda x: _np.sin(2 * x)) + sS = _np.vectorize(lambda x: x*_np.sin(x)) s1 = s1(x) s2 = s2(x) s3 = s3(x) sS = sS(x) - S1 = _np.hstack([s1, sS]) + s1 -= s1.mean(); s1 /= s1.std(0) + s2 -= s2.mean(); s2 /= s2.std(0) + s3 -= s3.mean(); s3 /= s3.std(0) + sS -= sS.mean(); sS /= sS.std(0) + + S1 = _np.hstack([s1, s2, sS]) S2 = _np.hstack([s2, s3, sS]) S3 = _np.hstack([s3, sS]) @@ -271,7 +276,7 @@ def bgplvm_simulation(optimize=True, verbose=1, D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 3, 10 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] - k = kern.linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) + k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k) if optimize: @@ -291,10 +296,10 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, from GPy.models import BayesianGPLVM from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData - D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 3, 10 + D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 5, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] - k = kern.linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) + k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) inan = _np.random.binomial(1, .6, size=Y.shape).astype(bool) m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing, kernel=k) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index a81bb711..5e88569c 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -308,14 +308,14 @@ class VarDTCMissingData(object): # gradients: if uncertain_inputs: grad_dict = {'dL_dKmm': dL_dKmm, - 'dL_dpsi0':dL_dpsi0, - 'dL_dpsi1':dL_dpsi1, - 'dL_dpsi2':dL_dpsi2, + 'dL_dpsi0':dL_dpsi0_all, + 'dL_dpsi1':dL_dpsi1_all, + 'dL_dpsi2':dL_dpsi2_all, 'partial_for_likelihood':partial_for_likelihood} else: grad_dict = {'dL_dKmm': dL_dKmm, - 'dL_dKdiag':dL_dpsi0, - 'dL_dKnm':dL_dpsi1, + 'dL_dKdiag':dL_dpsi0_all, + 'dL_dKnm':dL_dpsi1_all, 'partial_for_likelihood':partial_for_likelihood} #get sufficient things for posterior prediction diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index b3ee57cd..f436d322 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -67,7 +67,7 @@ class Kern(Parameterized): See GPy.plotting.matplot_dep.plot_ARD """ assert "matplotlib" in sys.modules, "matplotlib package has not been imported." - from ..plotting.matplot_dep import kernel_plots + from ...plotting.matplot_dep import kernel_plots return kernel_plots.plot_ARD(self,*args) diff --git a/GPy/plotting/matplot_dep/dim_reduction_plots.py b/GPy/plotting/matplot_dep/dim_reduction_plots.py index 74292c05..3f4ea9b0 100644 --- a/GPy/plotting/matplot_dep/dim_reduction_plots.py +++ b/GPy/plotting/matplot_dep/dim_reduction_plots.py @@ -1,8 +1,8 @@ import pylab as pb import numpy as np -from ... import util from latent_space_visualizations.controllers.imshow_controller import ImshowController,ImAnnotateController -from GPy.util.misc import param_to_array +from ...util.misc import param_to_array +from .base_plots import x_frame2D import itertools import Tango from matplotlib.cm import get_cmap @@ -37,7 +37,7 @@ def plot_latent(model, labels=None, which_indices=None, if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - util.plot.Tango.reset() + Tango.reset() if labels is None: labels = np.ones(model.num_data) @@ -46,7 +46,7 @@ def plot_latent(model, labels=None, which_indices=None, X = param_to_array(model.X) # first, plot the output variance as a function of the latent space - Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(X[:, [input_1, input_2]], resolution=resolution) + Xtest, xx, yy, xmin, xmax = x_frame2D(X[:, [input_1, input_2]], resolution=resolution) Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1])) def plot_function(x): @@ -87,7 +87,7 @@ def plot_latent(model, labels=None, which_indices=None, else: x = X[index, input_1] y = X[index, input_2] - ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label) + ax.scatter(x, y, marker=m, s=s, color=Tango.nextMedium(), label=this_label) ax.set_xlabel('latent dimension %i' % input_1) ax.set_ylabel('latent dimension %i' % input_2) @@ -120,7 +120,7 @@ def plot_magnification(model, labels=None, which_indices=None, if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - util.plot.Tango.reset() + Tango.reset() if labels is None: labels = np.ones(model.num_data) @@ -128,7 +128,7 @@ def plot_magnification(model, labels=None, which_indices=None, input_1, input_2 = most_significant_input_dimensions(model, which_indices) # first, plot the output variance as a function of the latent space - Xtest, xx, yy, xmin, xmax = util.plot.x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution) + Xtest, xx, yy, xmin, xmax = x_frame2D(model.X[:, [input_1, input_2]], resolution=resolution) Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1])) def plot_function(x): @@ -165,7 +165,7 @@ def plot_magnification(model, labels=None, which_indices=None, else: x = model.X[index, input_1] y = model.X[index, input_2] - ax.scatter(x, y, marker=m, s=s, color=util.plot.Tango.nextMedium(), label=this_label) + ax.scatter(x, y, marker=m, s=s, color=Tango.nextMedium(), label=this_label) ax.set_xlabel('latent dimension %i' % input_1) ax.set_ylabel('latent dimension %i' % input_2) @@ -205,7 +205,7 @@ def plot_steepest_gradient_map(model, fignum=None, ax=None, which_indices=None, return dmu_dX[indices, argmax], np.array(labels)[argmax] if ax is None: - fig = pyplot.figure(num=fignum) + fig = pb.figure(num=fignum) ax = fig.add_subplot(111) if data_labels is None: @@ -241,7 +241,7 @@ def plot_steepest_gradient_map(model, fignum=None, ax=None, which_indices=None, ax.legend() ax.figure.tight_layout() if updates: - pyplot.show() + pb.show() clear = raw_input('Enter to continue') if clear.lower() in 'yes' or clear == '': controller.deactivate() diff --git a/GPy/plotting/matplot_dep/kernel_plots.py b/GPy/plotting/matplot_dep/kernel_plots.py index 30157294..3436c4ff 100644 --- a/GPy/plotting/matplot_dep/kernel_plots.py +++ b/GPy/plotting/matplot_dep/kernel_plots.py @@ -1,7 +1,6 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -import sys import numpy as np import pylab as pb import Tango @@ -29,22 +28,23 @@ def plot_ARD(kernel, fignum=None, ax=None, title='', legend=False): xticklabels = [] bars = [] x0 = 0 - for p in kernel._parameters_: - c = Tango.nextMedium() - if hasattr(p, 'ARD') and p.ARD: - if title is None: - ax.set_title('ARD parameters, %s kernel' % p.name) - else: - ax.set_title(title) - if isinstance(p, Linear): - ard_params = p.variances - else: - ard_params = 1. / p.lengthscale - - x = np.arange(x0, x0 + len(ard_params)) - bars.append(ax.bar(x, ard_params, align='center', color=c, edgecolor='k', linewidth=1.2, label=p.name.replace("_"," "))) - xticklabels.extend([r"$\mathrm{{{name}}}\ {x}$".format(name=p.name, x=i) for i in np.arange(len(ard_params))]) - x0 += len(ard_params) + #for p in kernel._parameters_: + p = kernel + c = Tango.nextMedium() + if hasattr(p, 'ARD') and p.ARD: + if title is None: + ax.set_title('ARD parameters, %s kernel' % p.name) + else: + ax.set_title(title) + if isinstance(p, Linear): + ard_params = p.variances + else: + ard_params = 1. / p.lengthscale + x = np.arange(x0, x0 + len(ard_params)) + from ...util.misc import param_to_array + bars.append(ax.bar(x, param_to_array(ard_params), align='center', color=c, edgecolor='k', linewidth=1.2, label=p.name.replace("_"," "))) + xticklabels.extend([r"$\mathrm{{{name}}}\ {x}$".format(name=p.name, x=i) for i in np.arange(len(ard_params))]) + x0 += len(ard_params) x = np.arange(x0) transOffset = offset_copy(ax.transData, fig=fig, x=0., y= -2., units='points') diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 47c8642e..59c32775 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -56,7 +56,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - + + X, Y = param_to_array(model.X, model.Y) + if model.has_uncertain_inputs(): X_variance = model.X_variance + #work out what the inputs are for plotting (1D or 2D) fixed_dims = np.array([i for i,v in fixed_inputs]) free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims) @@ -66,7 +69,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #define the frame on which to plot resolution = resolution or 200 - Xnew, xmin, xmax = x_frame1D(model.X[:,free_dims], plot_limits=plot_limits) + Xnew, xmin, xmax = x_frame1D(X[:,free_dims], plot_limits=plot_limits) Xgrid = np.empty((Xnew.shape[0],model.input_dim)) Xgrid[:,free_dims] = Xnew for i,v in fixed_inputs: @@ -77,13 +80,13 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', m, v = model._raw_predict(Xgrid) lower = m - 2*np.sqrt(v) upper = m + 2*np.sqrt(v) - Y = model.Y + Y = Y else: m, v, lower, upper = model.predict(Xgrid) - Y = model.Y + Y = Y for d in which_data_ycols: gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol) - ax.plot(model.X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5) + ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5) #optionally plot some samples if samples: #NOTE not tested with fixed_inputs @@ -95,8 +98,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #add error bars for uncertain (if input uncertainty is being modelled) if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): - ax.errorbar(model.X[which_data_rows, free_dims], model.Y[which_data_rows, which_data_ycols], - xerr=2 * np.sqrt(model.X_variance[which_data_rows, free_dims]), + ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), + xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()), ecolor='k', fmt=None, elinewidth=.5, alpha=.5) @@ -120,7 +123,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #define the frame for plotting on resolution = resolution or 50 - Xnew, _, _, xmin, xmax = x_frame2D(model.X[:,free_dims], plot_limits, resolution) + Xnew, _, _, xmin, xmax = x_frame2D(X[:,free_dims], plot_limits, resolution) Xgrid = np.empty((Xnew.shape[0],model.input_dim)) Xgrid[:,free_dims] = Xnew for i,v in fixed_inputs: @@ -130,14 +133,14 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #predict on the frame and plot if plot_raw: m, _ = model._raw_predict(Xgrid) - Y = model.Y + Y = Y else: m, _, _, _ = model.predict(Xgrid) Y = model.data for d in which_data_ycols: m_d = m[:,d].reshape(resolution, resolution).T ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) - ax.scatter(model.X[which_data_rows, free_dims[0]], model.X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) + ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) #set the limits of the plot to some sensible values ax.set_xlim(xmin[0], xmax[0]) diff --git a/GPy/testing/index_operations_tests.py b/GPy/testing/index_operations_tests.py index d5ef7007..171db5cc 100644 --- a/GPy/testing/index_operations_tests.py +++ b/GPy/testing/index_operations_tests.py @@ -24,6 +24,13 @@ class Test(unittest.TestCase): self.param_index.remove(one, [1]) self.assertListEqual(self.param_index[one].tolist(), [3]) + def test_shift_left(self): + self.param_index.shift_left(1, 2) + self.assertListEqual(self.param_index[three].tolist(), [2,5]) + self.assertListEqual(self.param_index[two].tolist(), [0,3]) + self.assertListEqual(self.param_index[one].tolist(), [1]) + + def test_index_view(self): #======================================================================= # 0 1 2 3 4 5 6 7 8 9 diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index ff57606a..6f13d294 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -10,8 +10,8 @@ import numpy as np class Test(unittest.TestCase): def setUp(self): - self.rbf = GPy.kern.rbf(1) - self.white = GPy.kern.white(1) + self.rbf = GPy.kern.RBF(1) + self.white = GPy.kern.White(1) from GPy.core.parameterization import Param from GPy.core.parameterization.transformations import Logistic self.param = Param('param', np.random.rand(25,2), Logistic(0, 1)) @@ -39,14 +39,13 @@ class Test(unittest.TestCase): def test_remove_parameter(self): - from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__ + from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__, Logexp self.white.fix() self.test1.remove_parameter(self.white) self.assertIs(self.test1._fixes_,None) self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) - self.assertIs(self.white.constraints,self.white.white.constraints._param_index_ops) - self.assertEquals(self.white.white.constraints._offset, 0) + self.assertEquals(self.white.constraints._offset, 0) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) @@ -57,18 +56,19 @@ class Test(unittest.TestCase): self.assertListEqual(self.test1.constraints[__fixed__].tolist(), [0]) self.assertIs(self.white._fixes_,None) self.assertListEqual(self.test1._fixes_.tolist(),[FIXED] + [UNFIXED] * 52) + self.test1.remove_parameter(self.white) self.assertIs(self.test1._fixes_,None) self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) - self.assertIs(self.white.constraints,self.white.white.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) - self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) + self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1]) def test_add_parameter_already_in_hirarchy(self): self.test1.add_parameter(self.white._parameters_[0]) def test_default_constraints(self): - self.assertIs(self.rbf.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops) + self.assertIs(self.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) self.assertListEqual(self.rbf.constraints.indices()[0].tolist(), range(2)) from GPy.core.parameterization.transformations import Logexp