From d3a4f99b898e6b66f78dcd73bdd3b6ba0daaa7ef Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 7 Jun 2013 18:49:19 +0100 Subject: [PATCH 1/5] pypi release update --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index a683928b..4e9873b8 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ import os from setuptools import setup # Version number -version = '0.4.3' +version = '0.4.5' def read(fname): return open(os.path.join(os.path.dirname(__file__), fname)).read() From 7071dff030f2ca540d5fb1483027879d4612f6bd Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Thu, 13 Jun 2013 06:38:11 +0100 Subject: [PATCH 2/5] Mods to visualize and dimensionality to make stick demos work for summer school. --- GPy/examples/dimensionality_reduction.py | 39 +++++++++++++++++++++--- GPy/util/visualize.py | 5 +-- 2 files changed, 37 insertions(+), 7 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 28ee2bde..429b38de 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -261,6 +261,7 @@ def bgplvm_simulation(optimize='scg', k = kern.linear(Q, ARD=True) + kern.bias(Q, np.exp(-2)) + kern.white(Q, np.exp(-2)) # + kern.bias(Q) m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k, _debug=True) + # m.constrain('variance|noise', logexp_clipped()) m.ensure_default_constraints() m['noise'] = Y.var() / 100. @@ -327,28 +328,56 @@ def brendan_faces(): data_show = GPy.util.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, invert=False, scale=False) lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) raw_input('Press enter to finish') - lvm_visualizer.close() return m +def stick_play(range=None, frame_rate=15): + data = GPy.util.datasets.stick() + # optimize + if range==None: + Y = data['Y'].copy() + else: + Y = data['Y'][range[0]:range[1], :].copy() + y = Y[0, :] + data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect']) + GPy.util.visualize.data_play(Y, data_show, frame_rate) + return Y def stick(): data = GPy.util.datasets.stick() - m = GPy.models.GPLVM(data['Y'], 2) - # optimize + m = GPy.models.GPLVM(data['Y'], 2) m.ensure_default_constraints() m.optimize(messages=1, max_f_eval=10000) m._set_params(m._get_params()) - + plt.clf ax = m.plot_latent() y = m.likelihood.Y[0, :] data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect']) lvm_visualizer = GPy.util.visualize.lvm(m.X[0, :].copy(), m, data_show, ax) raw_input('Press enter to finish') - lvm_visualizer.close() return m +def stick_bgplvm(model=None): + data = GPy.util.datasets.stick() + Q = 6 + kernel = GPy.kern.rbf(Q, ARD=True) + GPy.kern.bias(Q, np.exp(-2)) + GPy.kern.white(Q, np.exp(-2)) + m = BayesianGPLVM(data['Y'], Q, init="PCA", num_inducing=20,kernel=kernel) + # optimize + m.ensure_default_constraints() + m.optimize(messages=1, max_f_eval=3000,xtol=1e-300,ftol=1e-300) + m._set_params(m._get_params()) + plt.clf, (latent_axes, sense_axes) = plt.subplots(1, 2) + plt.sca(latent_axes) + m.plot_latent() + y = m.likelihood.Y[0, :].copy() + data_show = GPy.util.visualize.stick_show(y[None, :], connect=data['connect']) + lvm_visualizer = GPy.util.visualize.lvm_dimselect(m.X[0, :].copy(), m, data_show, latent_axes=latent_axes, sense_axes=sense_axes) + raw_input('Press enter to finish') + + return m + + def cmu_mocap(subject='35', motion=['01'], in_place=True): data = GPy.util.datasets.cmu_mocap(subject, motion) diff --git a/GPy/util/visualize.py b/GPy/util/visualize.py index 11730704..8803d6a3 100644 --- a/GPy/util/visualize.py +++ b/GPy/util/visualize.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib as mpl import time import Image -#import visual +import visual class data_show: """ @@ -203,6 +203,7 @@ class lvm_dimselect(lvm): self.sense_axes = sense_axes self.labels = labels lvm.__init__(self,vals,Model,data_visualize,latent_axes,sense_axes,latent_index) + self.show_sensitivities() print "use left and right mouse butons to select dimensions" @@ -506,5 +507,5 @@ def data_play(Y, visualizer, frame_rate=30): for y in Y: - visualizer.modify(y) + visualizer.modify(y[None, :]) time.sleep(1./float(frame_rate)) From d85aaaea962ad4bcbfe5faecfcf805fdc7039d23 Mon Sep 17 00:00:00 2001 From: Neil Lawrence Date: Thu, 13 Jun 2013 06:39:13 +0100 Subject: [PATCH 3/5] Comment import visual in visualize due to issues in Windows/OSX. --- GPy/util/visualize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/util/visualize.py b/GPy/util/visualize.py index 8803d6a3..529f0eff 100644 --- a/GPy/util/visualize.py +++ b/GPy/util/visualize.py @@ -5,7 +5,7 @@ import numpy as np import matplotlib as mpl import time import Image -import visual +#import visual class data_show: """ From 259c81273a4ea8cc1bbb42f8f458e0462b38ad71 Mon Sep 17 00:00:00 2001 From: Nicolas Date: Thu, 13 Jun 2013 18:14:20 +0100 Subject: [PATCH 4/5] corrected minor bug in Brownian kernel --- GPy/kern/Brownian.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/kern/Brownian.py b/GPy/kern/Brownian.py index 76e103af..bdfa0df5 100644 --- a/GPy/kern/Brownian.py +++ b/GPy/kern/Brownian.py @@ -21,7 +21,7 @@ class Brownian(Kernpart): def __init__(self,input_dim,variance=1.): self.input_dim = input_dim assert self.input_dim==1, "Brownian motion in 1D only" - self.num_params = 1. + self.num_params = 1 self.name = 'Brownian' self._set_params(np.array([variance]).flatten()) From 0c12641ecaf1e6ccf162fb265190f72d194ac8a6 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 14 Jun 2013 14:04:53 +0100 Subject: [PATCH 5/5] changes to the efficiency of the sparse GP when there are many outputs --- GPy/core/sparse_gp.py | 33 ++++++++++++++++++++++----------- GPy/likelihoods/ep.py | 6 ++++++ GPy/likelihoods/gaussian.py | 5 ++++- 3 files changed, 32 insertions(+), 12 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 3183cff0..04119071 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -63,6 +63,7 @@ class SparseGP(GPBase): def _computations(self): + # factor Kmm self.Lm = jitchol(self.Kmm) @@ -89,17 +90,18 @@ class SparseGP(GPBase): self.B = np.eye(self.num_inducing) + self.A self.LB = jitchol(self.B) - # TODO: make a switch for either first compute psi1V, or VV.T - self.psi1V = np.dot(self.psi1.T, self.likelihood.V) + #VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! + self.psi1Vf = np.dot(self.psi1.T, self.likelihood.VVT_factor) - # back substutue C into psi1V - tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1V), lower=1, trans=0) - self._LBi_Lmi_psi1V, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0) + # back substutue C into psi1Vf + tmp, info1 = dtrtrs(self.Lm, np.asfortranarray(self.psi1Vf), lower=1, trans=0) + self._LBi_Lmi_psi1Vf, _ = dtrtrs(self.LB, np.asfortranarray(tmp), lower=1, trans=0) tmp, info2 = dpotrs(self.LB, tmp, lower=1) - self.Cpsi1V, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1) + self.Cpsi1Vf, info3 = dtrtrs(self.Lm, tmp, lower=1, trans=1) # Compute dL_dKmm - tmp = tdot(self._LBi_Lmi_psi1V) + tmp = tdot(self._LBi_Lmi_psi1Vf) + self.data_fit = np.trace(tmp) self.DBi_plus_BiPBi = backsub_both_sides(self.LB, self.output_dim * np.eye(self.num_inducing) + tmp) tmp = -0.5 * self.DBi_plus_BiPBi tmp += -0.5 * self.B * self.output_dim @@ -108,7 +110,7 @@ class SparseGP(GPBase): # Compute dL_dpsi # FIXME: this is untested for the heterscedastic + uncertain inputs case self.dL_dpsi0 = -0.5 * self.output_dim * (self.likelihood.precision * np.ones([self.num_data, 1])).flatten() - self.dL_dpsi1 = np.dot(self.Cpsi1V, self.likelihood.V.T).T + self.dL_dpsi1 = np.dot(self.likelihood.VVT_factor, self.Cpsi1Vf.T) dL_dpsi2_beta = 0.5 * backsub_both_sides(self.Lm, self.output_dim * np.eye(self.num_inducing) - self.DBi_plus_BiPBi) if self.likelihood.is_heteroscedastic: @@ -138,18 +140,18 @@ class SparseGP(GPBase): # likelihood is not heterscedatic self.partial_for_likelihood = -0.5 * self.num_data * self.output_dim * self.likelihood.precision + 0.5 * self.likelihood.trYYT * self.likelihood.precision ** 2 self.partial_for_likelihood += 0.5 * self.output_dim * (self.psi0.sum() * self.likelihood.precision ** 2 - np.trace(self.A) * self.likelihood.precision) - self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - np.sum(np.square(self._LBi_Lmi_psi1V))) + self.partial_for_likelihood += self.likelihood.precision * (0.5 * np.sum(self.A * self.DBi_plus_BiPBi) - self.data_fit) def log_likelihood(self): """ Compute the (lower bound on the) log marginal likelihood """ if self.likelihood.is_heteroscedastic: - A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V * self.likelihood.Y) + A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.likelihood.precision)) - 0.5 * np.sum(self.likelihood.V*self.likelihood.Y) B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision.flatten() * self.psi0) - np.trace(self.A)) else: A = -0.5 * self.num_data * self.output_dim * (np.log(2.*np.pi) - np.log(self.likelihood.precision)) - 0.5 * self.likelihood.precision * self.likelihood.trYYT B = -0.5 * self.output_dim * (np.sum(self.likelihood.precision * self.psi0) - np.trace(self.A)) C = -self.output_dim * (np.sum(np.log(np.diag(self.LB)))) # + 0.5 * self.num_inducing * np.log(sf2)) - D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V)) + D = 0.5 * self.data_fit return A + B + C + D + self.likelihood.Z def _set_params(self, p): @@ -158,6 +160,7 @@ class SparseGP(GPBase): self.likelihood._set_params(p[self.Z.size + self.kern.num_params:]) self._compute_kernel_matrices() self._computations() + self.Cpsi1V = None def _get_params(self): return np.hstack([self.Z.flatten(), self.kern._get_params_transformed(), self.likelihood._get_params()]) @@ -224,6 +227,14 @@ class SparseGP(GPBase): symmetrify(Bi) Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi) + if self.Cpsi1V is None: + psi1V = np.dot(self.psi1.T,self.likelihood.V) + tmp, _ = dtrtrs(self.Lm, np.asfortranarray(psi1V), lower=1, trans=0) + tmp, _ = dpotrs(self.LB, tmp, lower=1) + self.Cpsi1V, _ = dtrtrs(self.Lm, tmp, lower=1, trans=1) + + + if X_variance_new is None: Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts) mu = np.dot(Kx.T, self.Cpsi1V) diff --git a/GPy/likelihoods/ep.py b/GPy/likelihoods/ep.py index fb9a55c7..94f760e9 100644 --- a/GPy/likelihoods/ep.py +++ b/GPy/likelihoods/ep.py @@ -34,6 +34,8 @@ class EP(likelihood): self.Z = 0 self.YYT = None self.V = self.precision * self.Y + self.VVT_factor = self.V + self.trYYT = 0. def restart(self): self.tau_tilde = np.zeros(self.N) @@ -44,6 +46,8 @@ class EP(likelihood): self.Z = 0 self.YYT = None self.V = self.precision * self.Y + self.VVT_factor = self.V + self.trYYT = 0. def predictive_values(self,mu,var,full_cov): if full_cov: @@ -71,6 +75,8 @@ class EP(likelihood): self.covariance_matrix = np.diag(1./self.tau_tilde) self.precision = self.tau_tilde[:,None] self.V = self.precision * self.Y + self.VVT_factor = self.V + self.trYYT = np.trace(self.YYT) def fit_full(self,K): """ diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index 59d8fe86..8dbd8863 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -40,9 +40,11 @@ class Gaussian(likelihood): if D > self.N: self.YYT = np.dot(self.Y, self.Y.T) self.trYYT = np.trace(self.YYT) + self.YYT_factor = jitchol(self.YYT) else: self.YYT = None self.trYYT = np.sum(np.square(self.Y)) + self.YYT_factor = self.Y def _get_params(self): return np.asarray(self._variance) @@ -53,12 +55,13 @@ class Gaussian(likelihood): def _set_params(self, x): x = np.float64(x) if np.all(self._variance != x): - if x == 0.: + if x == 0.:#special case of zero noise self.precision = np.inf self.V = None else: self.precision = 1. / x self.V = (self.precision) * self.Y + self.VVT_factor = self.precision * self.YYT_factor self.covariance_matrix = np.eye(self.N) * x self._variance = x