From 52ba8e4ba36fdfbcb0f0e643c7e1a366065fe250 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Sun, 28 Apr 2013 17:22:04 +0100 Subject: [PATCH] remo0ved slices from models slices are now handles by special indexing kern parts, such as coregionalisation, independent_outputs. The old slicing functionality has been removed simply to clean up the code a little. Now that input_slices still exist (and will continue to be useful) in kern.py. They do need a little work though, for the psi-statistics --- GPy/kern/kern.py | 152 ++++++++++------------------- GPy/models/GP.py | 63 +++++------- GPy/models/GP_regression.py | 10 +- GPy/models/generalized_FITC.py | 17 ++-- GPy/models/sparse_GP.py | 20 ++-- GPy/models/sparse_GP_regression.py | 12 +-- GPy/models/warped_GP.py | 4 +- 7 files changed, 103 insertions(+), 175 deletions(-) diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index a6551e11..4547fadc 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -13,15 +13,9 @@ from prod import prod class kern(parameterised): def __init__(self, D, parts=[], input_slices=None): """ - This kernel does 'compound' structures. + This is the main kernel class for GPy. It handles multiple (additive) kernel functions, and keeps track of variaous things like which parameters live where. - The compund structure enables many features of GPy, including - - Hierarchical models - - Correleated output models - - multi-view learning - - Hadamard product and outer-product kernels will require a new class. - This feature is currently WONTFIX. for small number sof inputs, you can use the sympy kernel for this. + The technical code for kernels is divided into _parts_ (see e.g. rbf.py). This obnject contains a list of parts, which are computed additively. For multiplication, special _prod_ parts are used. :param D: The dimensioality of the kernel's input space :type D: int @@ -94,34 +88,6 @@ class kern(parameterised): self.param_slices.append(slice(count, count + p.Nparam)) count += p.Nparam - def _process_slices(self, slices1=None, slices2=None): - """ - Format the slices so that they can easily be used. - Both slices can be any of three things: - - If None, the new points covary through every kernel part (default) - - If a list of slices, the i^th slice specifies which data are affected by the i^th kernel part - - If a list of booleans, specifying which kernel parts are active - - if the second arg is False, return only slices1 - - returns actual lists of slice objects - """ - if slices1 is None: - slices1 = [slice(None)] * self.Nparts - elif all([type(s_i) is bool for s_i in slices1]): - slices1 = [slice(None) if s_i else slice(0) for s_i in slices1] - else: - assert all([type(s_i) is slice for s_i in slices1]), "invalid slice objects" - if slices2 is None: - slices2 = [slice(None)] * self.Nparts - elif slices2 is False: - return slices1 - elif all([type(s_i) is bool for s_i in slices2]): - slices2 = [slice(None) if s_i else slice(0) for s_i in slices2] - else: - assert all([type(s_i) is slice for s_i in slices2]), "invalid slice objects" - return slices1, slices2 - def __add__(self, other): assert self.D == other.D newkern = kern(self.D, self.parts + other.parts, self.input_slices + other.input_slices) @@ -142,7 +108,7 @@ class kern(parameterised): :param other: the other kernel to be added :type other: GPy.kern """ - return self +other + return self + other def add_orthogonal(self, other): """ @@ -285,18 +251,19 @@ class kern(parameterised): return sum([[name + '_' + n for n in k._get_param_names()] for name, k in zip(names, self.parts)], []) - def K(self, X, X2=None, slices1=None, slices2=None): + def K(self, X, X2=None, which_parts='all'): + if which_parts=='all': + which_parts = [True]*self.Nparts assert X.shape[1] == self.D - slices1, slices2 = self._process_slices(slices1, slices2) if X2 is None: target = np.zeros((X.shape[0], X.shape[0])) - [p.K(X[s1, i_s], None, target=target[s1, s2]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.K(X[:, i_s], None, target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used] else: target = np.zeros((X.shape[0], X2.shape[0])) - [p.K(X[s1, i_s], X2[s2, i_s], target=target[s1, s2]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.K(X[:, i_s], X2[:,i_s], target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used] return target - def dK_dtheta(self, dL_dK, X, X2=None, slices1=None, slices2=None): + def dK_dtheta(self, dL_dK, X, X2=None): """ :param dL_dK: An array of dL_dK derivaties, dL_dK :type dL_dK: Np.ndarray (N x M) @@ -304,109 +271,94 @@ class kern(parameterised): :type X: np.ndarray (N x D) :param X2: Observed dara inputs (optional, defaults to X) :type X2: np.ndarray (M x D) - :param slices1: a slice object for each kernel part, describing which data are affected by each kernel part - :type slices1: list of slice objects, or list of booleans - :param slices2: slices for X2 """ assert X.shape[1] == self.D - slices1, slices2 = self._process_slices(slices1, slices2) target = np.zeros(self.Nparam) if X2 is None: - [p.dK_dtheta(dL_dK[s1, s2], X[s1, i_s], None, target[ps]) for p, i_s, ps, s1, s2 in zip(self.parts, self.input_slices,self.param_slices, slices1, slices2)] + [p.dK_dtheta(dL_dK, X[:, i_s], None, target[ps]) for p, i_s, ps, in zip(self.parts, self.input_slices, self.param_slices)] else: - [p.dK_dtheta(dL_dK[s1, s2], X[s1, i_s], X2[s2, i_s], target[ps]) for p, i_s, ps, s1, s2 in zip(self.parts, self.input_slices,self.param_slices, slices1, slices2)] - + [p.dK_dtheta(dL_dK, X[:, i_s], X2[:, i_s], target[ps]) for p, i_s, ps, in zip(self.parts, self.input_slices, self.param_slices)] return self._transform_gradients(target) - def dK_dX(self, dL_dK, X, X2=None, slices1=None, slices2=None): + def dK_dX(self, dL_dK, X, X2=None): if X2 is None: X2 = X - slices1, slices2 = self._process_slices(slices1, slices2) target = np.zeros_like(X) - [p.dK_dX(dL_dK[s1, s2], X[s1, i_s], X2[s2, i_s], target[s1, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + if X2 is None: + [p.dK_dX(dL_dK, X[:, i_s], None, target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] + else: + [p.dK_dX(dL_dK, X[:, i_s], X2[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] return target - def Kdiag(self, X, slices=None): + def Kdiag(self, X, which_parts='all'): + if which_parts=='all': + which_parts = [True]*self.Nparts assert X.shape[1] == self.D - slices = self._process_slices(slices, False) target = np.zeros(X.shape[0]) - [p.Kdiag(X[s, i_s], target=target[s]) for p, i_s, s in zip(self.parts, self.input_slices, slices)] + [p.Kdiag(X[:, i_s], target=target) for p, i_s in zip(self.parts, self.input_slices)] return target - def dKdiag_dtheta(self, dL_dKdiag, X, slices=None): + def dKdiag_dtheta(self, dL_dKdiag, X): assert X.shape[1] == self.D - assert len(dL_dKdiag.shape) == 1 assert dL_dKdiag.size == X.shape[0] - slices = self._process_slices(slices, False) target = np.zeros(self.Nparam) - [p.dKdiag_dtheta(dL_dKdiag[s], X[s, i_s], target[ps]) for p, i_s, s, ps in zip(self.parts, self.input_slices, slices, self.param_slices)] + [p.dKdiag_dtheta(dL_dKdiag, X[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)] return self._transform_gradients(target) - def dKdiag_dX(self, dL_dKdiag, X, slices=None): + def dKdiag_dX(self, dL_dKdiag, X): assert X.shape[1] == self.D - slices = self._process_slices(slices, False) target = np.zeros_like(X) - [p.dKdiag_dX(dL_dKdiag[s], X[s, i_s], target[s, i_s]) for p, i_s, s in zip(self.parts, self.input_slices, slices)] + [p.dKdiag_dX(dL_dKdiag, X[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] return target - def psi0(self, Z, mu, S, slices=None): - slices = self._process_slices(slices, False) + def psi0(self, Z, mu, S): target = np.zeros(mu.shape[0]) - [p.psi0(Z, mu[s], S[s], target[s]) for p, s in zip(self.parts, slices)] + [p.psi0(Z[:,i_s], mu[:,i_s], S[:,i_s], target) for p, i_s in zip(self.parts, self.input_slices)] return target - def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S, slices=None): - slices = self._process_slices(slices, False) + def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S): target = np.zeros(self.Nparam) - [p.dpsi0_dtheta(dL_dpsi0[s], Z, mu[s], S[s], target[ps]) for p, ps, s in zip(self.parts, self.param_slices, slices)] + [p.dpsi0_dtheta(dL_dpsi0, Z[:,i_s], mu[:,i_s], S[:,i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)] return self._transform_gradients(target) - def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S, slices=None): - slices = self._process_slices(slices, False) + def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S): target_mu, target_S = np.zeros_like(mu), np.zeros_like(S) - [p.dpsi0_dmuS(dL_dpsi0, Z, mu[s], S[s], target_mu[s], target_S[s]) for p, s in zip(self.parts, slices)] + [p.dpsi0_dmuS(dL_dpsi0, Z[:,i_s], mu[:,i_s], S[:,i_s], target_mu[:,i_s], target_S[:,i_s]) for p, i_s in zip(self.parts, self.input_slices)] return target_mu, target_S - def psi1(self, Z, mu, S, slices1=None, slices2=None): - """Think N,M,Q """ - slices1, slices2 = self._process_slices(slices1, slices2) + def psi1(self, Z, mu, S): target = np.zeros((mu.shape[0], Z.shape[0])) - [p.psi1(Z[s2], mu[s1], S[s1], target[s1, s2]) for p, s1, s2 in zip(self.parts, slices1, slices2)] + [p.psi1(Z[:,i_s], mu[:,i_s], S[:,i_s], target) for p, i_s in zip(self.parts, self.input_slices)] return target - def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S, slices1=None, slices2=None): - """N,M,(Ntheta)""" - slices1, slices2 = self._process_slices(slices1, slices2) + def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S): target = np.zeros((self.Nparam)) - [p.dpsi1_dtheta(dL_dpsi1[s2, s1], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[ps]) for p, ps, s1, s2, i_s in zip(self.parts, self.param_slices, slices1, slices2, self.input_slices)] + [p.dpsi1_dtheta(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)] return self._transform_gradients(target) - def dpsi1_dZ(self, dL_dpsi1, Z, mu, S, slices1=None, slices2=None): - """N,M,Q""" - slices1, slices2 = self._process_slices(slices1, slices2) + def dpsi1_dZ(self, dL_dpsi1, Z, mu, S): target = np.zeros_like(Z) - [p.dpsi1_dZ(dL_dpsi1[s2, s1], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[s2, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.dpsi1_dZ(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] return target - def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S, slices1=None, slices2=None): + def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S): """return shapes are N,M,Q""" - slices1, slices2 = self._process_slices(slices1, slices2) target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1])) - [p.dpsi1_dmuS(dL_dpsi1[s2, s1], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target_mu[s1, i_s], target_S[s1, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] return target_mu, target_S - def psi2(self, Z, mu, S, slices1=None, slices2=None): + def psi2(self, Z, mu, S): """ :param Z: np.ndarray of inducing inputs (M x Q) :param mu, S: np.ndarrays of means and variances (each N x Q) :returns psi2: np.ndarray (N,M,M) """ target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0])) - slices1, slices2 = self._process_slices(slices1, slices2) - [p.psi2(Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[s1, s2, s2]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)] # compute the "cross" terms + #TODO: input_slices needed for p1, p2 in itertools.combinations(self.parts, 2): # white doesn;t combine with anything if p1.name == 'white' or p2.name == 'white': @@ -434,14 +386,12 @@ class kern(parameterised): raise NotImplementedError, "psi2 cannot be computed for this kernel" return target - def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None): - """Returns shape (N,M,M,Ntheta)""" - slices1, slices2 = self._process_slices(slices1, slices2) + def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S): target = np.zeros(self.Nparam) - [p.dpsi2_dtheta(dL_dpsi2[s1, s2, s2], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[ps]) for p, i_s, s1, s2, ps in zip(self.parts, self.input_slices, slices1, slices2, self.param_slices)] + [p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)] # compute the "cross" terms - # TODO: better looping + # TODO: better looping, input_slices for i1, i2 in itertools.combinations(range(len(self.parts)), 2): p1, p2 = self.parts[i1], self.parts[i2] # ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2] @@ -478,12 +428,12 @@ class kern(parameterised): return self._transform_gradients(target) - def dpsi2_dZ(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None): - slices1, slices2 = self._process_slices(slices1, slices2) + def dpsi2_dZ(self, dL_dpsi2, Z, mu, S): target = np.zeros_like(Z) - [p.dpsi2_dZ(dL_dpsi2[s1, s2, s2], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target[s2, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] # compute the "cross" terms + #TODO: we need input_slices here. for p1, p2 in itertools.combinations(self.parts, 2): # white doesn;t combine with anything if p1.name == 'white' or p2.name == 'white': @@ -506,16 +456,14 @@ class kern(parameterised): else: raise NotImplementedError, "psi2 cannot be computed for this kernel" - return target * 2. - def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S, slices1=None, slices2=None): - """return shapes are N,M,M,Q""" - slices1, slices2 = self._process_slices(slices1, slices2) + def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S): target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1])) - [p.dpsi2_dmuS(dL_dpsi2[s1, s2, s2], Z[s2, i_s], mu[s1, i_s], S[s1, i_s], target_mu[s1, i_s], target_S[s1, i_s]) for p, i_s, s1, s2 in zip(self.parts, self.input_slices, slices1, slices2)] + [p.dpsi2_dmuS(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)] # compute the "cross" terms + #TODO: we need input_slices here. for p1, p2 in itertools.combinations(self.parts, 2): # white doesn;t combine with anything if p1.name == 'white' or p2.name == 'white': diff --git a/GPy/models/GP.py b/GPy/models/GP.py index c6e46bea..45ed61ca 100644 --- a/GPy/models/GP.py +++ b/GPy/models/GP.py @@ -19,7 +19,6 @@ class GP(model): :parm likelihood: a GPy likelihood :param normalize_X: whether to normalize the input data before computing (predictions will be in original scales) :type normalize_X: False|True - :param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing) :rtype: model object :param epsilon_ep: convergence criterion for the Expectation Propagation algorithm, defaults to 0.1 :param powerep: power-EP parameters [$\eta$,$\delta$], defaults to [1.,1.] @@ -28,10 +27,9 @@ class GP(model): .. Note:: Multiple independent outputs are allowed using columns of Y """ - def __init__(self, X, likelihood, kernel, normalize_X=False, Xslices=None): + def __init__(self, X, likelihood, kernel, normalize_X=False): # parse arguments - self.Xslices = Xslices self.X = X assert len(self.X.shape) == 2 self.N, self.Q = self.X.shape @@ -64,12 +62,12 @@ class GP(model): return np.zeros_like(self.Z) def _set_params(self, p): - self.kern._set_params_transformed(p[:self.kern.Nparam]) + self.kern._set_params_transformed(p[:self.kern.Nparam_transformed()]) # self.likelihood._set_params(p[self.kern.Nparam:]) # test by Nicolas self.likelihood._set_params(p[self.kern.Nparam_transformed():]) # test by Nicolas - self.K = self.kern.K(self.X, slices1=self.Xslices, slices2=self.Xslices) + self.K = self.kern.K(self.X) self.K += self.likelihood.covariance_matrix self.Ki, self.L, self.Li, self.K_logdet = pdinv(self.K) @@ -92,7 +90,7 @@ class GP(model): """ Approximates a non-gaussian likelihood using Expectation Propagation - For a Gaussian (or direct: TODO) likelihood, no iteration is required: + For a Gaussian likelihood, no iteration is required: this function does nothing """ self.likelihood.fit_full(self.kern.K(self.X)) @@ -122,31 +120,33 @@ class GP(model): """ The gradient of all parameters. - For the kernel parameters, use the chain rule via dL_dK - - For the likelihood parameters, pass in alpha = K^-1 y + Note, we use the chain rule: dL_dtheta = dL_dK * d_K_dtheta """ - return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X, slices1=self.Xslices, slices2=self.Xslices), self.likelihood._gradients(partial=np.diag(self.dL_dK)))) + return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK)))) - def _raw_predict(self, _Xnew, slices=None, full_cov=False): + def _raw_predict(self, _Xnew, which_parts='all', full_cov=False): """ Internal helper function for making predictions, does not account for normalization or likelihood + + #TODO: which_parts does nothing + + """ - Kx = self.kern.K(self.X, _Xnew, slices1=self.Xslices, slices2=slices) + Kx = self.kern.K(self.X, _Xnew,which_parts=which_parts) mu = np.dot(np.dot(Kx.T, self.Ki), self.likelihood.Y) KiKx = np.dot(self.Ki, Kx) if full_cov: - Kxx = self.kern.K(_Xnew, slices1=slices, slices2=slices) + Kxx = self.kern.K(_Xnew, which_parts=which_parts) var = Kxx - np.dot(KiKx.T, Kx) else: - Kxx = self.kern.Kdiag(_Xnew, slices=slices) + Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts) var = Kxx - np.sum(np.multiply(KiKx, Kx), 0) var = var[:, None] return mu, var - def predict(self, Xnew, slices=None, full_cov=False): + def predict(self, Xnew, which_parts='all', full_cov=False): """ Predict the function(s) at the new point(s) Xnew. @@ -154,19 +154,14 @@ class GP(model): --------- :param Xnew: The points at which to make a prediction :type Xnew: np.ndarray, Nnew x self.Q - :param slices: specifies which outputs kernel(s) the Xnew correspond to (see below) - :type slices: (None, list of slice objects, list of ints) + :param which_parts: specifies which outputs kernel(s) to use in prediction + :type which_parts: ('all', list of bools) :param full_cov: whether to return the folll covariance matrix, or just the diagonal :type full_cov: bool :rtype: posterior mean, a Numpy array, Nnew x self.D :rtype: posterior variance, a Numpy array, Nnew x 1 if full_cov=False, Nnew x Nnew otherwise :rtype: lower and upper boundaries of the 95% confidence intervals, Numpy arrays, Nnew x self.D - .. Note:: "slices" specifies how the the points X_new co-vary wich the training points. - - - If None, the new points covary throigh every kernel part (default) - - If a list of slices, the i^th slice specifies which data are affected by the i^th kernel part - - If a list of booleans, specifying which kernel parts are active If full_cov and self.D > 1, the return shape of var is Nnew x Nnew x self.D. If self.D == 1, the return shape is Nnew x Nnew. This is to allow for different normalizations of the output dimensions. @@ -174,15 +169,15 @@ class GP(model): """ # normalize X values Xnew = (Xnew.copy() - self._Xmean) / self._Xstd - mu, var = self._raw_predict(Xnew, slices, full_cov) + mu, var = self._raw_predict(Xnew, which_parts, full_cov) - # now push through likelihood TODO + # now push through likelihood mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov) return mean, var, _025pm, _975pm - def plot_f(self, samples=0, plot_limits=None, which_data='all', which_functions='all', resolution=None, full_cov=False): + def plot_f(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False): """ Plot the GP's view of the world, where the data is normalized and the likelihood is Gaussian @@ -190,8 +185,8 @@ class GP(model): :param which_data: which if the training data to plot (default all) :type which_data: 'all' or a slice object to slice self.X, self.Y :param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits - :param which_functions: which of the kernel functions to plot (additively) - :type which_functions: list of bools + :param which_parts: which of the kernel functions to plot (additively) + :type which_parts: 'all', or list of bools :param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D Plot the posterior of the GP. @@ -202,19 +197,17 @@ class GP(model): Can plot only part of the data and part of the posterior functions using which_data and which_functions Plot the data's view of the world, with non-normalized values and GP predictions passed through the likelihood """ - if which_functions == 'all': - which_functions = [True] * self.kern.Nparts if which_data == 'all': which_data = slice(None) if self.X.shape[1] == 1: Xnew, xmin, xmax = x_frame1D(self.X, plot_limits=plot_limits) if samples == 0: - m, v = self._raw_predict(Xnew, slices=which_functions) + m, v = self._raw_predict(Xnew, which_parts=which_parts) gpplot(Xnew, m, m - 2 * np.sqrt(v), m + 2 * np.sqrt(v)) pb.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5) else: - m, v = self._raw_predict(Xnew, slices=which_functions, full_cov=True) + m, v = self._raw_predict(Xnew, which_parts=which_parts, full_cov=True) Ysim = np.random.multivariate_normal(m.flatten(), v, samples) gpplot(Xnew, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None]) for i in range(samples): @@ -230,7 +223,7 @@ class GP(model): elif self.X.shape[1] == 2: resolution = resolution or 50 Xnew, xmin, xmax, xx, yy = x_frame2D(self.X, plot_limits, resolution) - m, v = self._raw_predict(Xnew, slices=which_functions) + m, v = self._raw_predict(Xnew, which_parts=which_parts) m = m.reshape(resolution, resolution).T pb.contour(xx, yy, m, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) pb.scatter(Xorig[:, 0], Xorig[:, 1], 40, Yorig, linewidth=0, cmap=pb.cm.jet, vmin=m.min(), vmax=m.max()) @@ -246,8 +239,6 @@ class GP(model): """ # TODO include samples - if which_functions == 'all': - which_functions = [True] * self.kern.Nparts if which_data == 'all': which_data = slice(None) @@ -256,7 +247,7 @@ class GP(model): Xu = self.X * self._Xstd + self._Xmean # NOTE self.X are the normalized values now Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits) - m, var, lower, upper = self.predict(Xnew, slices=which_functions) + m, var, lower, upper = self.predict(Xnew, which_parts=which_parts) gpplot(Xnew, m, lower, upper) pb.plot(Xu[which_data], self.likelihood.data[which_data], 'kx', mew=1.5) if self.has_uncertain_inputs: @@ -277,7 +268,7 @@ class GP(model): resolution = resolution or 50 Xnew, xx, yy, xmin, xmax = x_frame2D(self.X, plot_limits, resolution) x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution) - m, var, lower, upper = self.predict(Xnew, slices=which_functions) + m, var, lower, upper = self.predict(Xnew, which_parts=which_parts) m = m.reshape(resolution, resolution).T pb.contour(x, y, m, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) Yf = self.likelihood.Y.flatten() diff --git a/GPy/models/GP_regression.py b/GPy/models/GP_regression.py index 5f9f9f3e..7f2673a6 100644 --- a/GPy/models/GP_regression.py +++ b/GPy/models/GP_regression.py @@ -11,26 +11,24 @@ class GP_regression(GP): """ Gaussian Process model for regression - This is a thin wrapper around the GP class, with a set of sensible defalts + This is a thin wrapper around the models.GP class, with a set of sensible defalts :param X: input observations :param Y: observed values - :param kernel: a GPy kernel, defaults to rbf+white + :param kernel: a GPy kernel, defaults to rbf :param normalize_X: whether to normalize the input data before computing (predictions will be in original scales) :type normalize_X: False|True :param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales) :type normalize_Y: False|True - :param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing) - :rtype: model object .. Note:: Multiple independent outputs are allowed using columns of Y """ - def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False, Xslices=None): + def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False): if kernel is None: kernel = kern.rbf(X.shape[1]) likelihood = likelihoods.Gaussian(Y,normalize=normalize_Y) - GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, Xslices=Xslices) + GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X) diff --git a/GPy/models/generalized_FITC.py b/GPy/models/generalized_FITC.py index 26875f64..25b6c18f 100644 --- a/GPy/models/generalized_FITC.py +++ b/GPy/models/generalized_FITC.py @@ -23,20 +23,19 @@ class generalized_FITC(sparse_GP): :type X_variance: np.ndarray (N x Q) | None :param Z: inducing inputs (optional, see note) :type Z: np.ndarray (M x Q) | None - :param Zslices: slices for the inducing inputs (see slicing TODO: link) :param M : Number of inducing points (optional, default 10. Ignored if Z is not None) :type M: int :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales) :type normalize_(X|Y): bool """ - def __init__(self, X, likelihood, kernel, Z, X_variance=None, Xslices=None,Zslices=None, normalize_X=False): + def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False): self.Z = Z self.M = self.Z.shape[0] self._precision = likelihood.precision - sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_variance=None, Xslices=None,Zslices=None, normalize_X=False) + sparse_GP.__init__(self, X, likelihood, kernel=kernel, Z=self.Z, X_variance=None, normalize_X=False) def _set_params(self, p): self.Z = p[:self.M*self.Q].reshape(self.M, self.Q) @@ -145,7 +144,7 @@ class generalized_FITC(sparse_GP): D = 0.5*np.trace(self.Cpsi1VVpsi1) return A+C+D - def _raw_predict(self, Xnew, slices, full_cov=False): + def _raw_predict(self, Xnew, which_parts, full_cov=False): if self.likelihood.is_heteroscedastic: """ Make a prediction for the generalized FITC model @@ -174,16 +173,16 @@ class generalized_FITC(sparse_GP): self.mu_H = mu_H Sigma_H = C + np.dot(mu_u,np.dot(self.Sigma,mu_u.T)) # q(f_star|y) = N(f_star|mu_star,sigma2_star) - Kx = self.kern.K(self.Z, Xnew) + Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts) KR0T = np.dot(Kx.T,self.Lmi.T) mu_star = np.dot(KR0T,mu_H) if full_cov: - Kxx = self.kern.K(Xnew) + Kxx = self.kern.K(Xnew,which_parts=which_parts) var = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) else: - Kxx = self.kern.Kdiag(Xnew) - Kxx_ = self.kern.K(Xnew) - var_ = Kxx_ + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) + Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts) + Kxx_ = self.kern.K(Xnew,which_parts=which_parts) # TODO: RA, is this line needed? + var_ = Kxx_ + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.M),KR0T.T)) # TODO: RA, is this line needed? var = (Kxx + np.sum(KR0T.T*np.dot(Sigma_H - np.eye(self.M),KR0T.T),0))[:,None] return mu_star[:,None],var else: diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 697a9978..20caa1a8 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -9,10 +9,6 @@ from .. import kern from GP import GP from scipy import linalg -#Still TODO: -# make use of slices properly (kernel can now do this) -# enable heteroscedatic noise (kernel will need to compute psi2 as a (NxMxM) array) - class sparse_GP(GP): """ Variational sparse GP model @@ -27,19 +23,16 @@ class sparse_GP(GP): :type X_variance: np.ndarray (N x Q) | None :param Z: inducing inputs (optional, see note) :type Z: np.ndarray (M x Q) | None - :param Zslices: slices for the inducing inputs (see slicing TODO: link) :param M : Number of inducing points (optional, default 10. Ignored if Z is not None) :type M: int :param normalize_(X|Y) : whether to normalize the data before computing (predictions will be in original scales) :type normalize_(X|Y): bool """ - def __init__(self, X, likelihood, kernel, Z, X_variance=None, Xslices=None,Zslices=None, normalize_X=False): + def __init__(self, X, likelihood, kernel, Z, X_variance=None, normalize_X=False): self.scale_factor = 100.0# a scaling factor to help keep the algorithm stable self.auto_scale_factor = False self.Z = Z - self.Zslices = Zslices - self.Xslices = Xslices self.M = Z.shape[0] self.likelihood = likelihood @@ -50,7 +43,7 @@ class sparse_GP(GP): self.has_uncertain_inputs=True self.X_variance = X_variance - GP.__init__(self, X, likelihood, kernel=kernel, normalize_X=normalize_X, Xslices=Xslices) + GP.__init__(self, X, likelihood, kernel=kernel, normalize_X=normalize_X) #normalize X uncertainty also if self.has_uncertain_inputs: @@ -65,13 +58,12 @@ class sparse_GP(GP): self.psi1 = self.kern.psi1(self.Z,self.X, self.X_variance).T self.psi2 = self.kern.psi2(self.Z,self.X, self.X_variance) else: - self.psi0 = self.kern.Kdiag(self.X,slices=self.Xslices) + self.psi0 = self.kern.Kdiag(self.X) self.psi1 = self.kern.K(self.Z,self.X) self.psi2 = None def _computations(self): #TODO: find routine to multiply triangular matrices - #TODO: slices for psi statistics (easy enough) sf = self.scale_factor sf2 = sf**2 @@ -252,16 +244,16 @@ class sparse_GP(GP): dL_dZ += self.kern.dK_dX(self.dL_dpsi1,self.Z,self.X) return dL_dZ - def _raw_predict(self, Xnew, slices, full_cov=False): + def _raw_predict(self, Xnew, which_parts='all', full_cov=False): """Internal helper function for making predictions, does not account for normalization""" Kx = self.kern.K(self.Z, Xnew) mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) if full_cov: - Kxx = self.kern.K(Xnew) + Kxx = self.kern.K(Xnew,which_parts=which_parts) var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx) #NOTE this won't work for plotting else: - Kxx = self.kern.Kdiag(Xnew) + Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts) var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0) return mu,var[:,None] diff --git a/GPy/models/sparse_GP_regression.py b/GPy/models/sparse_GP_regression.py index 0ef78c32..84a5d37c 100644 --- a/GPy/models/sparse_GP_regression.py +++ b/GPy/models/sparse_GP_regression.py @@ -13,7 +13,7 @@ class sparse_GP_regression(sparse_GP): """ Gaussian Process model for regression - This is a thin wrapper around the GP class, with a set of sensible defalts + This is a thin wrapper around the sparse_GP class, with a set of sensible defalts :param X: input observations :param Y: observed values @@ -22,25 +22,25 @@ class sparse_GP_regression(sparse_GP): :type normalize_X: False|True :param normalize_Y: whether to normalize the input data before computing (predictions will be in original scales) :type normalize_Y: False|True - :param Xslices: how the X,Y data co-vary in the kernel (i.e. which "outputs" they correspond to). See (link:slicing) :rtype: model object .. Note:: Multiple independent outputs are allowed using columns of Y """ - def __init__(self,X,Y,kernel=None,normalize_X=False,normalize_Y=False, Xslices=None,Z=None, M=10): - #kern defaults to rbf + def __init__(self, X, Y, kernel=None, normalize_X=False, normalize_Y=False, Z=None, M=10): + #kern defaults to rbf (plus white for stability) if kernel is None: kernel = kern.rbf(X.shape[1]) + kern.white(X.shape[1],1e-3) #Z defaults to a subset of the data if Z is None: - Z = np.random.permutation(X.copy())[:M] + i = np.random.permutation(X.shape[0])[:M] + Z = X[i].copy() else: assert Z.shape[1]==X.shape[1] #likelihood defaults to Gaussian likelihood = likelihoods.Gaussian(Y,normalize=normalize_Y) - sparse_GP.__init__(self, X, likelihood, kernel, Z, normalize_X=normalize_X, Xslices=Xslices) + sparse_GP.__init__(self, X, likelihood, kernel, Z, normalize_X=normalize_X) diff --git a/GPy/models/warped_GP.py b/GPy/models/warped_GP.py index 052f8d8e..9c3ce401 100644 --- a/GPy/models/warped_GP.py +++ b/GPy/models/warped_GP.py @@ -14,7 +14,7 @@ from .. import likelihoods from .. import kern class warpedGP(GP): - def __init__(self, X, Y, kernel=None, warping_function = None, warping_terms = 3, normalize_X=False, normalize_Y=False, Xslices=None): + def __init__(self, X, Y, kernel=None, warping_function = None, warping_terms = 3, normalize_X=False, normalize_Y=False): if kernel is None: kernel = kern.rbf(X.shape[1]) @@ -28,7 +28,7 @@ class warpedGP(GP): self.predict_in_warped_space = False likelihood = likelihoods.Gaussian(self.transform_data(), normalize=normalize_Y) - GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X, Xslices=Xslices) + GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X) def _set_params(self, x): self.warping_params = x[:self.warping_function.num_parameters]