From 4aca883df35c670c1684200269c79e38f90452bc Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 10 Apr 2013 15:50:31 +0100 Subject: [PATCH 1/6] weaved some rbf code --- GPy/core/model.py | 11 ++++++++--- GPy/kern/rbf.py | 48 +++++++++++++++++++++++++++++++++++++++++++++-- 2 files changed, 54 insertions(+), 5 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 645f7228..c7e3c4b8 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -229,13 +229,18 @@ class model(parameterised): Ensure that any variables which should clearly be positive have been constrained somehow. """ positive_strings = ['variance','lengthscale', 'precision'] + param_names = self._get_param_names() + currently_constrained = self.all_constrained_indices() + to_make_positive = [] for s in positive_strings: for i in self.grep_param_names(s): - if not (i in self.all_constrained_indices()): - name = self._get_param_names()[i] - self.constrain_positive(name) + if not (i in currently_constrained): + to_make_positive.append(param_names[i]) if warn: print "Warning! constraining %s postive"%name + if len(to_make_positive): + self.constrain_positive('('+'|'.join(to_make_positive)+')') + def objective_function(self, x): diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index 84f7d68d..d473a5cc 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -5,6 +5,7 @@ from kernpart import kernpart import numpy as np import hashlib +from scipy import weave class rbf(kernpart): """ @@ -217,8 +218,51 @@ class rbf(kernpart): #psi2 self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,Q self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q - self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) - self._psi2_exponent = np.sum(-self._psi2_Zdist_sq/4. -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M + self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_stuff() + #self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) + #self._psi2_exponent = np.sum(-self._psi2_Zdist_sq/4. -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M + #store matrices for caching self._Z, self._mu, self._S = Z, mu,S + + def weave_psi2(self): + weave_options = {'extra_compile_args': ['-O3']} + N,M,M,Q = self._psi2_mudist.shape + mudist = self._psi2_mudist + psi2_Zdist_sq = self._psi2_Zdist_sq + mudist_sq = np.empty((N,M,M,Q)) + psi2_exponent = np.zeros((N,M,M)) + psi2 = np.empty((N,M,M)) + half_log_psi2_denom = 0.5*np.log(self._psi2_denom).squeeze() + variance_sq = float(np.square(self.variance)) + if self.ARD: + lengthscale2 = self.lengthscale2 + else: + lengthscale2 = np.ones(Q)*self.lengthscale2 + _psi2_denom = self._psi2_denom.squeeze() + code = """ + double tmp; + for (int n=0; n Date: Wed, 10 Apr 2013 16:12:09 +0100 Subject: [PATCH 2/6] improved weaving --- GPy/kern/rbf.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index d473a5cc..b774db1c 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -217,8 +217,8 @@ class rbf(kernpart): #psi2 self._psi2_denom = 2.*S[:,None,None,:]/self.lengthscale2+1. # N,M,M,Q - self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q - self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_stuff() + self._psi2_mudist, self._psi2_mudist_sq, self._psi2_exponent, _ = self.weave_psi2(mu,self._psi2_Zhat) + #self._psi2_mudist = mu[:,None,None,:]-self._psi2_Zhat #N,M,M,Q #self._psi2_mudist_sq = np.square(self._psi2_mudist)/(self.lengthscale2*self._psi2_denom) #self._psi2_exponent = np.sum(-self._psi2_Zdist_sq/4. -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M self._psi2 = np.square(self.variance)*np.exp(self._psi2_exponent) # N,M,M @@ -226,14 +226,18 @@ class rbf(kernpart): #store matrices for caching self._Z, self._mu, self._S = Z, mu,S - def weave_psi2(self): + def weave_psi2(self,mu,Zhat): weave_options = {'extra_compile_args': ['-O3']} - N,M,M,Q = self._psi2_mudist.shape - mudist = self._psi2_mudist - psi2_Zdist_sq = self._psi2_Zdist_sq + + N,Q = mu.shape + M = Zhat.shape[0] + + mudist = np.empty((N,M,M,Q)) mudist_sq = np.empty((N,M,M,Q)) psi2_exponent = np.zeros((N,M,M)) psi2 = np.empty((N,M,M)) + + psi2_Zdist_sq = self._psi2_Zdist_sq half_log_psi2_denom = 0.5*np.log(self._psi2_denom).squeeze() variance_sq = float(np.square(self.variance)) if self.ARD: @@ -247,14 +251,23 @@ class rbf(kernpart): for (int m=0; m Date: Wed, 10 Apr 2013 16:50:02 +0100 Subject: [PATCH 3/6] OMP for psi2 computations in RBF --- GPy/kern/rbf.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index b774db1c..a26bb79c 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -227,7 +227,10 @@ class rbf(kernpart): self._Z, self._mu, self._S = Z, mu,S def weave_psi2(self,mu,Zhat): - weave_options = {'extra_compile_args': ['-O3']} + weave_options = {'headers' : [''], + 'extra_compile_args': ['-fopenmp -march=native'], + 'extra_link_args' : ['-lgomp'], + 'compiler' : 'gcc'} N,Q = mu.shape M = Zhat.shape[0] @@ -247,6 +250,8 @@ class rbf(kernpart): _psi2_denom = self._psi2_denom.squeeze() code = """ double tmp; + + #pragma omp parallel for private(tmp) for (int n=0; n + #include + """ + 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) + + return mudist,mudist_sq, psi2_exponent, psi2 From 6499a76e246ed81ccb6005065d9a3a6a4243ade5 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 10 Apr 2013 16:50:34 +0100 Subject: [PATCH 4/6] pdinv now uses dpotri instead of dtrtri and dot --- GPy/util/linalg.py | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 59f598f9..d82bb50f 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -51,6 +51,26 @@ def _mdot_r(a,b): return np.dot(a,b) def jitchol(A,maxtries=5): + A = np.asfortranarray(A) + L,info = linalg.lapack.flapack.dpotrf(A,lower=1) + if info ==0: + return L + else: + diagA = np.diag(A) + if np.any(diagA<0.): + raise linalg.LinAlgError, "not pd: negative diagonal elements" + jitter= diagA.mean()*1e-6 + for i in range(1,maxtries+1): + print 'Warning: adding jitter of '+str(jitter) + try: + return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True) + except: + jitter *= 10 + raise linalg.LinAlgError,"not positive definite, even with jitter." + + + +def jitchol_old(A,maxtries=5): """ :param A : An almost pd square matrix @@ -71,7 +91,7 @@ def jitchol(A,maxtries=5): for i in range(1,maxtries+1): print 'Warning: adding jitter of '+str(jitter) try: - return linalg.cholesky(A+np.eye(A.shape[0])*jitter, lower = True) + return linalg.cholesky(A+np.eye(A.shape[0]).T*jitter, lower = True) except: jitter *= 10 @@ -93,7 +113,8 @@ def pdinv(A): L = jitchol(A) logdet = 2.*np.sum(np.log(np.diag(L))) Li = chol_inv(L) - Ai = np.dot(Li.T,Li) #TODO: get the flapack routine form multiplying triangular matrices + Ai = linalg.lapack.flapack.dpotri(L)[0] + Ai = np.tril(Ai) + np.tril(Ai,-1).T return Ai, L, Li, logdet From 87304a0778f9b46800ffa744e1e762017120085d Mon Sep 17 00:00:00 2001 From: Nicolo Fusi Date: Wed, 10 Apr 2013 17:04:44 +0100 Subject: [PATCH 5/6] merged master back into devel (to sync bugfixes) --- GPy/examples/dimensionality_reduction.py | 29 ------------------------ GPy/models/GPLVM.py | 5 ++-- 2 files changed, 3 insertions(+), 31 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index f1fdaaf1..61a4abd8 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -133,32 +133,3 @@ def stick(): plt.close('all') return m - - -def BGPLVM_oil(): - data = GPy.util.datasets.oil() - Y, X = data['Y'], data['X'] - X -= X.mean(axis=0) - X /= X.std(axis=0) - - Q = 10 - M = 30 - - kernel = GPy.kern.rbf(Q, ARD = True) + GPy.kern.bias(Q) + GPy.kern.white(Q) - m = GPy.models.Bayesian_GPLVM(X, Q, kernel=kernel, M=M) - # m.scale_factor = 100.0 - m.constrain_positive('(white|noise|bias|X_variance|rbf_variance|rbf_length)') - from sklearn import cluster - km = cluster.KMeans(M, verbose=10) - Z = km.fit(m.X).cluster_centers_ - # Z = GPy.util.misc.kmm_init(m.X, M) - m.set('iip', Z) - m.set('bias', 1e-4) - # optimize - # m.ensure_default_constraints() - - import pdb; pdb.set_trace() - m.optimize('tnc', messages=1) - print m - m.plot_latent(labels=data['Y'].argmax(axis=1)) - return m diff --git a/GPy/models/GPLVM.py b/GPy/models/GPLVM.py index 470aff96..2ce55dda 100644 --- a/GPy/models/GPLVM.py +++ b/GPy/models/GPLVM.py @@ -89,7 +89,7 @@ class GPLVM(GP): Xtest_full = np.zeros((Xtest.shape[0], self.X.shape[1])) Xtest_full[:, :2] = Xtest mu, var, low, up = self.predict(Xtest_full) - var = var[:, :1] # FIXME: this was a :2 + var = var[:, :1] pb.imshow(var.reshape(resolution,resolution).T[::-1,:], extent=[xmin[0], xmax[0], xmin[1], xmax[1]], cmap=pb.cm.binary,interpolation='bilinear') @@ -119,5 +119,6 @@ class GPLVM(GP): pb.xlim(xmin[0],xmax[0]) pb.ylim(xmin[1],xmax[1]) pb.grid(b=False) # remove the grid if present, it doesn't look good + ax = pb.gca() ax.set_aspect('auto') # set a nice aspect ratio - return pb.gca() #input_1, input_2 temporary removal, to return axes. + return ax From 48b0ac6399d8c9e0df0114daaf4c21d9ca1cb879 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 10 Apr 2013 20:02:22 +0100 Subject: [PATCH 6/6] some minor improvements in visualize --- GPy/examples/dimensionality_reduction.py | 8 ++++---- GPy/kern/rbf.py | 9 ++++----- GPy/models/Bayesian_GPLVM.py | 1 + GPy/util/visualize.py | 17 ++++++++++------- 4 files changed, 19 insertions(+), 16 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 61a4abd8..f2558341 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -60,7 +60,7 @@ def GPLVM_oil_100(optimize=True,M=15): m.plot_latent(labels=m.data_labels) return m -def BGPLVM_oil(optimize=True,N=100,Q=10,M=15): +def BGPLVM_oil(optimize=True,N=100,Q=10,M=15,max_f_eval=300): data = GPy.util.datasets.oil() # create simple GP model @@ -72,10 +72,10 @@ def BGPLVM_oil(optimize=True,N=100,Q=10,M=15): if optimize: m.constrain_fixed('noise',0.05) m.ensure_default_constraints() - m.optimize('scg',messages=1) + m.optimize('scg',messages=1,max_f_eval=max(80,max_f_eval)) m.unconstrain('noise') m.constrain_positive('noise') - m.optimize('scg',messages=1) + m.optimize('scg',messages=1,max_f_eval=max(0,max_f_eval-80)) else: m.ensure_default_constraints() @@ -120,7 +120,7 @@ def brendan_faces(): def stick(): data = GPy.util.datasets.stick() m = GPy.models.GPLVM(data['Y'], 2) - + # optimize m.ensure_default_constraints() m.optimize(messages=1, max_f_eval=10000) diff --git a/GPy/kern/rbf.py b/GPy/kern/rbf.py index a26bb79c..84521cf9 100644 --- a/GPy/kern/rbf.py +++ b/GPy/kern/rbf.py @@ -173,7 +173,7 @@ class rbf(kernpart): """Think N,M,M,Q """ self._psi_computations(Z,mu,S) tmp = self._psi2[:,:,:,None]/self.lengthscale2/self._psi2_denom - target_mu += (dL_dpsi2[:,:,:,None]*-tmp*2.*self._psi2_mudist).sum(1).sum(1) + target_mu += -2.*(dL_dpsi2[:,:,:,None]*tmp*self._psi2_mudist).sum(1).sum(1) target_S += (dL_dpsi2[:,:,:,None]*tmp*(2.*self._psi2_mudist_sq-1)).sum(1).sum(1) @@ -207,7 +207,6 @@ class rbf(kernpart): if not (np.all(Z==self._Z) and np.all(mu==self._mu) and np.all(S==self._S)): #something's changed. recompute EVERYTHING - #TODO: make more efficient for large Q (using NDL's dot product trick) #psi1 self._psi1_denom = S[:,None,:]/self.lengthscale2 + 1. self._psi1_dist = Z[None,:,:]-mu[:,None,:] @@ -250,7 +249,7 @@ class rbf(kernpart): _psi2_denom = self._psi2_denom.squeeze() code = """ double tmp; - + #pragma omp parallel for private(tmp) for (int n=0; n #include """ - weave.inline(code, support_code=support_code, libraries=['gomp'], + 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) - + return mudist,mudist_sq, psi2_exponent, psi2 diff --git a/GPy/models/Bayesian_GPLVM.py b/GPy/models/Bayesian_GPLVM.py index ba9603bb..ee485e76 100644 --- a/GPy/models/Bayesian_GPLVM.py +++ b/GPy/models/Bayesian_GPLVM.py @@ -95,3 +95,4 @@ class Bayesian_GPLVM(sparse_GP, GPLVM): input_1, input_2 = which_indices ax = GPLVM.plot_latent(self, which_indices=[input_1, input_2],*args, **kwargs) ax.plot(self.Z[:, input_1], self.Z[:, input_2], '^w') + return ax diff --git a/GPy/util/visualize.py b/GPy/util/visualize.py index dde9cd32..f9538942 100644 --- a/GPy/util/visualize.py +++ b/GPy/util/visualize.py @@ -4,7 +4,7 @@ import GPy import numpy as np class lvm: - def __init__(self, model, data_visualize, latent_axis): + def __init__(self, model, data_visualize, latent_axis, latent_index=[0,1], latent_dim=2): self.cid = latent_axis.figure.canvas.mpl_connect('button_press_event', self.on_click) self.cid = latent_axis.figure.canvas.mpl_connect('motion_notify_event', self.on_move) self.data_visualize = data_visualize @@ -12,6 +12,8 @@ class lvm: self.latent_axis = latent_axis self.called = False self.move_on = False + self.latent_index = latent_index + self.latent_dim = latent_dim def on_click(self, event): #print 'click', event.xdata, event.ydata @@ -32,7 +34,8 @@ class lvm: if self.called and self.move_on: # Call modify code on move #print 'move', event.xdata, event.ydata - latent_values = np.array((event.xdata, event.ydata)) + latent_values = np.zeros((1,self.latent_dim)) + latent_values[0,self.latent_index] = np.array([event.xdata, event.ydata]) y = self.model.predict(latent_values)[0] self.data_visualize.modify(y) #print 'y', y @@ -45,7 +48,7 @@ class data_show: # If no axes are defined, create some. if axis==None: fig = plt.figure() - self.axis = fig.add_subplot(111) + self.axis = fig.add_subplot(111) else: self.axis = axis @@ -57,7 +60,7 @@ class vector_show(data_show): def __init__(self, vals, axis=None): data_show.__init__(self, vals, axis) self.vals = vals.T - self.handle = plt.plot(np.arange(0, len(vals))[:, None], self.vals)[0] + self.handle = self.axis.plot(np.arange(0, len(vals))[:, None], self.vals)[0] def modify(self, vals): xdata, ydata = self.handle.get_data() @@ -84,7 +87,7 @@ class image_show(data_show): self.handle.set_array(self.vals) #self.axis.figure.canvas.draw() plt.show() - + def set_image(self, vals): self.vals = np.reshape(vals, self.dimensions, order='F') if self.transpose: @@ -94,7 +97,7 @@ class image_show(data_show): #if self.invert: # self.vals = -self.vals -class stick_show(data_show): +class stick_show(data_show): """Show a three dimensional point cloud as a figure. Connect elements of the figure together using the matrix connect.""" def __init__(self, vals, axis=None, connect=None): @@ -159,6 +162,6 @@ class stick_show(data_show): self.line_handle = self.axis.plot(np.array(x), np.array(y), np.array(z), 'b-') self.axis.figure.canvas.draw() - +