From 9b4bb41fe6f8072cfd3f5866196d4c6f3710871b Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 24 May 2013 13:08:17 +0100 Subject: [PATCH 1/5] added empty directoy for cmu mocap data --- GPy/util/datasets/mocap/cmu/README.txt | 1 + 1 file changed, 1 insertion(+) create mode 100644 GPy/util/datasets/mocap/cmu/README.txt diff --git a/GPy/util/datasets/mocap/cmu/README.txt b/GPy/util/datasets/mocap/cmu/README.txt new file mode 100644 index 00000000..83bec0b7 --- /dev/null +++ b/GPy/util/datasets/mocap/cmu/README.txt @@ -0,0 +1 @@ +this otherwise empty directory is for storing mnocap data, which we don;t distribute From a491b5237133a0d1c4c9e0b7eb837e4ed1fbb0a2 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 28 May 2013 11:11:49 +0200 Subject: [PATCH 2/5] added max to authors --- AUTHORS.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/AUTHORS.txt b/AUTHORS.txt index fc9e0fb7..f81db5ec 100644 --- a/AUTHORS.txt +++ b/AUTHORS.txt @@ -3,4 +3,5 @@ Nicolo Fusi Ricardo Andrade Nicolas Durrande Alan Saul +Max Zwiessele Neil D. Lawrence From 3853721355f792481890115d00c1b2c5c9a1cf6d Mon Sep 17 00:00:00 2001 From: Andreas Date: Tue, 28 May 2013 15:50:18 +0100 Subject: [PATCH 3/5] Fixed symmetrify for when C/F compiler doesn't work --- GPy/util/linalg.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index abfb1900..da23a0a8 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -313,7 +313,10 @@ def symmetrify(A,upper=False): elif A.flags['F_CONTIGUOUS'] and not upper: weave.inline(f_contig_code,['A','N'], extra_compile_args=['-O3']) else: - tmp = np.tril(A) + if upper: + tmp = np.tril(A.T) + else: + tmp = np.tril(A) A[:] = 0.0 A += tmp A += np.tril(tmp,-1).T From d461f5f9bfcea3986458fad43fde7d6c409eb53d Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 29 May 2013 13:22:25 +0100 Subject: [PATCH 4/5] sparse_GP now has a separate predict function GP and sparse_GP used t share a predict fumction. Since we'd like to propagate uncertainty in predictions, sparse_GP.predict needs to accept X_new_variance. --- GPy/core/model.py | 2 +- GPy/models/sparse_GP.py | 37 +++++++++++++++++++++++++++++++++++++ 2 files changed, 38 insertions(+), 1 deletion(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index f2b188d9..5dc6b254 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -447,7 +447,7 @@ class model(parameterised): assert isinstance(self.likelihood, likelihoods.EP), "EPEM is only available for EP likelihoods" ll_change = epsilon + 1. iteration = 0 - last_ll = -np.exp(1000) + last_ll = -np.inf convergence = False alpha = 0 diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py index 93877cbc..cc7105ae 100644 --- a/GPy/models/sparse_GP.py +++ b/GPy/models/sparse_GP.py @@ -245,3 +245,40 @@ class sparse_GP(GP): var = Kxx - np.sum(np.sum(psi2*Kmmi_LmiBLmi[None,:,:],1),1) return mu, var[:, None] + + def predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): + """ + Predict the function(s) at the new point(s) Xnew. + + Arguments + --------- + :param Xnew: The points at which to make a prediction + :type Xnew: np.ndarray, Nnew x self.Q + :param X_variance_new: The uncertainty in the prediction points + :type X_variance_new: np.ndarray, Nnew x self.Q + :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 + + + 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. + + """ + # normalize X values + Xnew = (Xnew.copy() - self._Xmean) / self._Xstd + if X_variance_new is not None: + X_variance_new = X_variance_new / self._Xstd**2 + + #here's the actual prediction by the GP model + mu, var = self._raw_predict(Xnew, X_variance_new, full_cov=full_cov, which_parts=which_parts) + + # now push through likelihood + mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov) + + return mean, var, _025pm, _975pm + From dc33aa1b8c334ed5cd05d574aa4690a636a0cc17 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 30 May 2013 09:29:26 +0100 Subject: [PATCH 5/5] used scipy.weave to improve the speed of rbf grads for large number of input dimensions with ARD, wqe have approx tenfold speedup. --- GPy/kern/rbf.py | 52 ++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 45 insertions(+), 7 deletions(-) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index c413b469..8b644241 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -56,6 +56,13 @@ class rbf(kernpart): self._Z, self._mu, self._S = np.empty(shape=(3,1)) self._X, self._X2, self._params = np.empty(shape=(3,1)) + #a set of optional args to pass to weave + self.weave_options = {'headers' : [''], + 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], + 'extra_link_args' : ['-lgomp']} + + + def _get_params(self): return np.hstack((self.variance,self.lengthscale)) @@ -85,8 +92,43 @@ class rbf(kernpart): self._K_computations(X,X2) target[0] += np.sum(self._K_dvar*dL_dK) if self.ARD: - if X2 is None: X2 = X - [np.add(target[1+q:2+q],(self.variance/self.lengthscale[q]**3)*np.sum(self._K_dvar*dL_dK*np.square(X[:,q][:,None]-X2[:,q][None,:])),target[1+q:2+q]) for q in range(self.D)] + dvardLdK = self._K_dvar*dL_dK + var_len3 = self.variance/np.power(self.lengthscale,3) + if X2 is None: + #save computation for the symmetrical case + dvardLdK += dvardLdK.T + code = """ + int q,i,j; + double tmp; + for(q=0; q'], - 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], - 'extra_link_args' : ['-lgomp']} - N,Q = mu.shape M = Zhat.shape[0] @@ -288,6 +326,6 @@ class rbf(kernpart): """ weave.inline(code, support_code=support_code, libraries=['gomp'], arg_names=['N','M','Q','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,**weave_options) + type_converters=weave.converters.blitz,**self.weave_options) return mudist,mudist_sq, psi2_exponent, psi2