From ba77b20a064ac14dffcd6c3e5e81e24fa6421a90 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Fri, 16 May 2014 15:41:02 +0100 Subject: [PATCH 01/14] Bug fixed --- .../latent_function_inference/expectation_propagation.py | 2 +- .../latent_function_inference/expectation_propagation_dtc.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/inference/latent_function_inference/expectation_propagation.py b/GPy/inference/latent_function_inference/expectation_propagation.py index de8b0931..1afc8100 100644 --- a/GPy/inference/latent_function_inference/expectation_propagation.py +++ b/GPy/inference/latent_function_inference/expectation_propagation.py @@ -32,7 +32,7 @@ class EP(LatentFunctionInference): pass def inference(self, kern, X, likelihood, Y, Y_metadata=None, Z=None): - num_data, output_dim = X.shape + num_data, output_dim = Y.shape assert output_dim ==1, "ep in 1D only (for now!)" K = kern.K(X) diff --git a/GPy/inference/latent_function_inference/expectation_propagation_dtc.py b/GPy/inference/latent_function_inference/expectation_propagation_dtc.py index 85d8cc89..3aeb4fbb 100644 --- a/GPy/inference/latent_function_inference/expectation_propagation_dtc.py +++ b/GPy/inference/latent_function_inference/expectation_propagation_dtc.py @@ -56,7 +56,7 @@ class EPDTC(LatentFunctionInference): self._ep_approximation = None def inference(self, kern, X, Z, likelihood, Y, Y_metadata=None): - num_data, output_dim = X.shape + num_data, output_dim = Y.shape assert output_dim ==1, "ep in 1D only (for now!)" Kmm = kern.K(Z) From e17e539bceee9f39caff5dd2260883facd7e2e1e Mon Sep 17 00:00:00 2001 From: Niu Date: Sat, 2 Aug 2014 14:06:24 +0100 Subject: [PATCH 02/14] initial implementation of hmc --- GPy/inference/optimization/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/inference/optimization/__init__.py b/GPy/inference/optimization/__init__.py index 1a8f043b..04e245e3 100644 --- a/GPy/inference/optimization/__init__.py +++ b/GPy/inference/optimization/__init__.py @@ -1,2 +1,3 @@ from scg import SCG from optimization import * +from hmc import HMC,Gmodel From 93b92111f821d611b08089a13c7a64ae612eb9ae Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Tue, 5 Aug 2014 08:27:40 -0700 Subject: [PATCH 03/14] [minor] minor changes --- GPy/core/model.py | 2 +- GPy/core/parameterization/parameter_core.py | 28 +-------------------- GPy/kern/_src/add.py | 4 +-- 3 files changed, 4 insertions(+), 30 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 1595e347..7feb72b2 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -349,7 +349,7 @@ class Model(Parameterized): numerical_gradient = (f1 - f2) / (2 * step) if np.all(gradient[xind] == 0): ratio = (f1 - f2) == gradient[xind] else: ratio = (f1 - f2) / (2 * step * gradient[xind]) - difference = np.abs((f1 - f2) / 2 / step - gradient[xind]) + difference = np.abs(numerical_gradient - gradient[xind]) if (np.abs(1. - ratio) < tolerance) or np.abs(difference) < tolerance: formatted_name = "\033[92m {0} \033[0m".format(names[nind]) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index e359409e..2a036378 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -699,36 +699,10 @@ class OptimizationHandlable(Indexable): def _get_params_transformed(self): raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!" -# # transformed parameters (apply un-transformation rules) -# p = self.param_array.copy() -# [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] -# if self.has_parent() and self.constraints[__fixed__].size != 0: -# fixes = np.ones(self.size).astype(bool) -# fixes[self.constraints[__fixed__]] = FIXED -# return p[fixes] -# elif self._has_fixes(): -# return p[self._fixes_] -# return p # def _set_params_transformed(self, p): raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!" -# """ -# Set parameters p, but make sure they get transformed before setting. -# This means, the optimizer sees p, whereas the model sees transformed(p), -# such that, the parameters the model sees are in the right domain. -# """ -# if not(p is self.param_array): -# if self.has_parent() and self.constraints[__fixed__].size != 0: -# fixes = np.ones(self.size).astype(bool) -# fixes[self.constraints[__fixed__]] = FIXED -# self.param_array.flat[fixes] = p -# elif self._has_fixes(): self.param_array.flat[self._fixes_] = p -# else: self.param_array.flat = p -# [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) -# for c, ind in self.constraints.iteritems() if c != __fixed__] -# self._trigger_params_changed() - def _trigger_params_changed(self, trigger_parent=True): """ First tell all children to update, @@ -736,7 +710,7 @@ class OptimizationHandlable(Indexable): If trigger_parent is True, we will tell the parent, otherwise not. """ - [p._trigger_params_changed(trigger_parent=False) for p in self.parameters] + [p._trigger_params_changed(trigger_parent=False) for p in self.parameters if not p.is_fixed] self.notify_observers(None, None if trigger_parent else -np.inf) def _size_transformed(self): diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 12f5d444..ee743f8b 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -10,7 +10,7 @@ class Add(CombinationKernel): """ Add given list of kernels together. propagates gradients through. - + This kernel will take over the active dims of it's subkernels passed in. """ def __init__(self, subkerns, name='add'): @@ -40,7 +40,7 @@ class Add(CombinationKernel): return reduce(np.add, (p.Kdiag(X) for p in which_parts)) def update_gradients_full(self, dL_dK, X, X2=None): - [p.update_gradients_full(dL_dK, X, X2) for p in self.parts] + [p.update_gradients_full(dL_dK, X, X2) for p in self.parts if not p.is_fixed] def update_gradients_diag(self, dL_dK, X): [p.update_gradients_diag(dL_dK, X) for p in self.parts] From 18a1513edb6b3ad9fe23d7e039eb755d19a764cb Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Tue, 5 Aug 2014 16:48:35 +0100 Subject: [PATCH 04/14] add hmc.py --- GPy/inference/optimization/hmc.py | 63 +++++++++++++++++++++++++++++++ 1 file changed, 63 insertions(+) create mode 100644 GPy/inference/optimization/hmc.py diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py new file mode 100644 index 00000000..a8bb0405 --- /dev/null +++ b/GPy/inference/optimization/hmc.py @@ -0,0 +1,63 @@ +"""HMC implementation""" + +import numpy as np + + +class HMC: + def __init__(self,model,M=None,stepsize=1e-1): + self.model = model + self.stepsize = stepsize + self.p = np.empty_like(model.optimizer_array.copy()) + if M is None: + self.M = np.eye(self.p.size) + else: + self.M = M + self.Minv = np.linalg.inv(self.M) + + def sample(self, m_iters=1000, hmc_iters=20): + thetas = np.empty((m_iters,self.p.size)) + ps = np.empty((m_iters,self.p.size)) + for i in xrange(m_iters): + #Gibbs + self.p[:] = np.random.multivariate_normal(np.ones(self.p.size),self.M) + H_old = self._computeH() + p_old = self.p.copy() + theta_old = self.model.optimizer_array.copy() + #Matropolis + self._update(hmc_iters) + H_new = self._computeH() + + k = np.exp(H_old-H_new) + print k + if np.random.rand() Date: Wed, 6 Aug 2014 11:01:08 +0100 Subject: [PATCH 05/14] debug hmc code --- GPy/inference/optimization/hmc.py | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py index a8bb0405..783aadc7 100644 --- a/GPy/inference/optimization/hmc.py +++ b/GPy/inference/optimization/hmc.py @@ -28,11 +28,14 @@ class HMC: H_new = self._computeH() k = np.exp(H_old-H_new) + print H_old,H_new print k - if np.random.rand() Date: Wed, 6 Aug 2014 11:10:34 +0100 Subject: [PATCH 06/14] remove prints --- GPy/inference/optimization/hmc.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py index 783aadc7..859877f2 100644 --- a/GPy/inference/optimization/hmc.py +++ b/GPy/inference/optimization/hmc.py @@ -28,14 +28,11 @@ class HMC: H_new = self._computeH() k = np.exp(H_old-H_new) - print H_old,H_new - print k a = np.random.rand() if a Date: Wed, 6 Aug 2014 14:44:56 +0100 Subject: [PATCH 07/14] coregionalized 2D plotting fixed --- GPy/plotting/matplot_dep/models_plots.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 7926410e..46a79ad8 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -150,7 +150,11 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if plot_raw: m, _ = model._raw_predict(Xgrid) else: - m, _ = model.predict(Xgrid) + if isinstance(model,GPCoregionalizedRegression) or isinstance(model,SparseGPCoregionalizedRegression): + meta = {'output_index': Xgrid[:,-1:].astype(np.int)} + else: + meta = None + m, v = model.predict(Xgrid, full_cov=False, Y_metadata=meta) for d in which_data_ycols: m_d = m[:,d].reshape(resolution, resolution).T plots['contour'] = ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) From d11a2c482c3547096a7ded079357ec3b20a173f0 Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Wed, 6 Aug 2014 16:33:35 +0100 Subject: [PATCH 08/14] correct the initial distribution of p --- GPy/inference/optimization/hmc.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py index 859877f2..a5b6fe19 100644 --- a/GPy/inference/optimization/hmc.py +++ b/GPy/inference/optimization/hmc.py @@ -18,8 +18,7 @@ class HMC: thetas = np.empty((m_iters,self.p.size)) ps = np.empty((m_iters,self.p.size)) for i in xrange(m_iters): - #Gibbs - self.p[:] = np.random.multivariate_normal(np.ones(self.p.size),self.M) + self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M) H_old = self._computeH() p_old = self.p.copy() theta_old = self.model.optimizer_array.copy() @@ -27,9 +26,11 @@ class HMC: self._update(hmc_iters) H_new = self._computeH() - k = np.exp(H_old-H_new) - a = np.random.rand() - if aH_new: + k = 1. + else: + k = np.exp(H_old-H_new) + if np.random.rand() Date: Wed, 6 Aug 2014 22:15:28 +0100 Subject: [PATCH 09/14] hmc shortcut --- GPy/inference/optimization/__init__.py | 2 +- GPy/inference/optimization/hmc.py | 106 +++++++++++++++++++++++-- 2 files changed, 99 insertions(+), 9 deletions(-) diff --git a/GPy/inference/optimization/__init__.py b/GPy/inference/optimization/__init__.py index 04e245e3..1590568f 100644 --- a/GPy/inference/optimization/__init__.py +++ b/GPy/inference/optimization/__init__.py @@ -1,3 +1,3 @@ from scg import SCG from optimization import * -from hmc import HMC,Gmodel +from hmc import HMC,HMC_shortcut diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py index a5b6fe19..93c0e5e3 100644 --- a/GPy/inference/optimization/hmc.py +++ b/GPy/inference/optimization/hmc.py @@ -48,13 +48,103 @@ class HMC: def _computeH(self,): return self.model.objective_function()+self.p.size*np.log(2*np.pi)/2.+np.log(np.linalg.det(self.M))/2.+np.dot(self.p, np.dot(self.Minv,self.p[:,None]))/2. -class Gmodel: - def __init__(self,): - self.cov = np.array([[1., 0.99],[0.99, 1.]]) - self.param_array = np.random.rand(2) +class HMC_shortcut: + def __init__(self,model,M=None,stepsize_range=[1e-6, 1e-1],groupsize=5, Hstd_th=[1e-3, 20.]): + self.model = model + self.stepsize_range = np.log10(stepsize_range) + self.p = np.empty_like(model.optimizer_array.copy()) + self.groupsize = groupsize + self.Hstd_th = Hstd_th + if M is None: + self.M = np.eye(self.p.size) + else: + self.M = M + self.Minv = np.linalg.inv(self.M) - def grad(self,): - return -np.dot(np.linalg.inv(self.cov),self.param_array[:,None]) + def sample(self, m_iters=1000, hmc_iters=20): + thetas = np.empty((m_iters,self.p.size)) + ps = np.empty((m_iters,self.p.size)) + for i in xrange(m_iters): + # sample a stepsize from the uniform distribution + stepsize = np.exp10(np.random.rand()*(self.stepsize_range[1]-self.stepsize_range[0])+self.stepsize_range[0]) + self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M) + H_old = self._computeH() + p_old = self.p.copy() + theta_old = self.model.optimizer_array.copy() + #Matropolis + self._update(hmc_iters, stepsize) + H_new = self._computeH() + + if H_old>H_new: + k = 1. + else: + k = np.exp(H_old-H_new) + if np.random.rand()(reversal[0]-pos): + pos_new = 2*reversal[0] - r - pos + else: + pos_new = 2*pos + r - reversal[0] + self.model.optimizer_array = theta_buf[hmc_iters+pos_new] + self.p[:] = p_buf[hmc_iters+pos_new] # the sign of momentum might be wrong! + break + + + def _testH(self, Hlist): + Hstd = np.std(Hlist) + if Hstdself.Hstd_th[1]: + return False + else: + return True + + def _computeH(self,): + return self.model.objective_function()+self.p.size*np.log(2*np.pi)/2.+np.log(np.linalg.det(self.M))/2.+np.dot(self.p, np.dot(self.Minv,self.p[:,None]))/2. - def objective_function(self,): - return np.dot(self.param_array, np.dot(np.linalg.inv(self.cov),self.param_array[:,None]))/2. From c3bb7b28a137ebb8c006c1b0f605fa663df5f890 Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Thu, 7 Aug 2014 14:28:46 +0100 Subject: [PATCH 10/14] debug HMC shortcut --- GPy/core/parameterization/parameter_core.py | 19 ++++++ GPy/inference/optimization/hmc.py | 69 ++++++++++++--------- 2 files changed, 57 insertions(+), 31 deletions(-) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index e359409e..e8ef186f 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -888,6 +888,25 @@ class Parameterizable(OptimizationHandlable): self._param_array_ = np.empty(self.size, dtype=np.float64) return self._param_array_ + @property + def unfixed_param_array(self): + """ + Array representing the parameters of this class. + There is only one copy of all parameters in memory, two during optimization. + + !WARNING!: setting the parameter array MUST always be done in memory: + m.param_array[:] = m_copy.param_array + """ + if self.__dict__.get('_param_array_', None) is None: + self._param_array_ = np.empty(self.size, dtype=np.float64) + + if self.constraints[__fixed__].size !=0: + fixes = np.ones(self.size).astype(bool) + fixes[self.constraints[__fixed__]] = FIXED + return self._param_array_[fixes] + else: + return self._param_array_ + @param_array.setter def param_array(self, arr): self._param_array_ = arr diff --git a/GPy/inference/optimization/hmc.py b/GPy/inference/optimization/hmc.py index 93c0e5e3..8c65cdf0 100644 --- a/GPy/inference/optimization/hmc.py +++ b/GPy/inference/optimization/hmc.py @@ -15,13 +15,12 @@ class HMC: self.Minv = np.linalg.inv(self.M) def sample(self, m_iters=1000, hmc_iters=20): - thetas = np.empty((m_iters,self.p.size)) - ps = np.empty((m_iters,self.p.size)) + params = np.empty((m_iters,self.p.size)) for i in xrange(m_iters): self.p[:] = np.random.multivariate_normal(np.zeros(self.p.size),self.M) H_old = self._computeH() - p_old = self.p.copy() theta_old = self.model.optimizer_array.copy() + params[i] = self.model.unfixed_param_array #Matropolis self._update(hmc_iters) H_new = self._computeH() @@ -31,13 +30,10 @@ class HMC: else: k = np.exp(H_old-H_new) if np.random.rand()pos: + pos = -1 + i += pos + self.model.optimizer_array = theta_buf[hmc_iters] + self.p[:] = -p_buf[hmc_iters] + else: + pos_new = pos-hmc_iters+i + self.model.optimizer_array = theta_buf[hmc_iters+pos_new] + self.p[:] = -p_buf[hmc_iters+pos_new] + break else: - Hlist = range(pos,pos+self.groupsize) + Hlist = range(hmc_iters+pos,hmc_iters+pos+self.groupsize) +# print Hlist +# print self._testH(H_buf[Hlist]) + if self._testH(H_buf[Hlist]): pos += -1 else: @@ -132,14 +136,17 @@ class HMC_shortcut: if r>(reversal[0]-pos): pos_new = 2*reversal[0] - r - pos else: - pos_new = 2*pos + r - reversal[0] + pos_new = pos + r self.model.optimizer_array = theta_buf[hmc_iters+pos_new] self.p[:] = p_buf[hmc_iters+pos_new] # the sign of momentum might be wrong! +# print reversal[0],pos,pos_new +# print H_buf break - def _testH(self, Hlist): Hstd = np.std(Hlist) +# print Hlist +# print Hstd if Hstdself.Hstd_th[1]: return False else: From 8273fd9e2a2b43630ca76e4395908a9576dd85bf Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Thu, 7 Aug 2014 08:24:41 -0700 Subject: [PATCH 11/14] [bgplvm] gradient settings --- GPy/models/bayesian_gplvm.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 6354f13d..a16c4ac1 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -84,6 +84,22 @@ class BayesianGPLVM(SparseGP): self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, dL_dpsi0=self.grad_dict['dL_dpsi0'], dL_dpsi1=self.grad_dict['dL_dpsi1'], dL_dpsi2=self.grad_dict['dL_dpsi2']) + # This is testing code ------------------------- +# i = np.random.randint(self.X.shape[0]) +# X_ = self.X.mean +# which = np.sqrt(((X_ - X_[i:i+1])**2).sum(1)).argsort()>(max(0, self.X.shape[0]-51)) +# _, _, grad_dict = self.inference_method.inference(self.kern, self.X[which], self.Z, self.likelihood, self.Y[which], self.Y_metadata) +# grad = self.kern.gradients_qX_expectations(variational_posterior=self.X[which], Z=self.Z, dL_dpsi0=grad_dict['dL_dpsi0'], dL_dpsi1=grad_dict['dL_dpsi1'], dL_dpsi2=grad_dict['dL_dpsi2']) +# +# self.X.mean.gradient[:] = 0 +# self.X.variance.gradient[:] = 0 +# self.X.mean.gradient[which] = grad[0] +# self.X.variance.gradient[which] = grad[1] + + # update for the KL divergence +# self.variational_prior.update_gradients_KL(self.X, which) + # ----------------------------------------------- + # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) From 6fc982faa089bb0411e8b6f6060b4d5dd593e4be Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Thu, 7 Aug 2014 08:28:16 -0700 Subject: [PATCH 12/14] [gradcheck] some performance enhancement --- GPy/core/parameterization/parameter_core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index f17f7d55..f6586d3c 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -281,7 +281,7 @@ class Gradcheckable(Pickleable, Parentable): """ if self.has_parent(): return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) - return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance) + return self._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) def _checkgrad(self, param, verbose=0, step=1e-6, tolerance=1e-3): """ From 2de4b92c38e2b60fad4882da473c8aacaa26ba6e Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Thu, 7 Aug 2014 08:39:53 -0700 Subject: [PATCH 13/14] [parameter core] offset for can be done without parameter slices --- GPy/core/parameterization/parameter_core.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index f6586d3c..0ecc1ebf 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -372,8 +372,9 @@ class Indexable(Nameable, Observable): basically just sums up the parameter sizes which come before param. """ if param.has_parent(): - if param._parent_._get_original(param) in self.parameters: - return self._param_slices_[param._parent_._get_original(param)._parent_index_].start + p = param._parent_._get_original(param) + if p in self.parameters: + return reduce(lambda a,b: a + b.size, self.parameters[:p._parent_index_], 0) return self._offset_for(param._parent_) + param._parent_._offset_for(param) return 0 From 479e59f3ef39124cd0b245d181442f1133632d60 Mon Sep 17 00:00:00 2001 From: mzwiessele Date: Thu, 7 Aug 2014 10:53:18 -0700 Subject: [PATCH 14/14] [setup] no import of os, nanana alan :) --- setup.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/setup.py b/setup.py index 40574a76..d9a6ab5e 100644 --- a/setup.py +++ b/setup.py @@ -1,13 +1,13 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- -#import os +import os from setuptools import setup # Version number version = '0.4.6' -from pkg_resources import Requirement, resource_string +from pkg_resources import Requirement def read(fname): return open(os.path.join(os.path.dirname(__file__), fname)).read()