From a9d4006488c4d476face17ba231ad916b672ab5b Mon Sep 17 00:00:00 2001 From: Nicolas Date: Wed, 5 Jun 2013 11:44:48 +0100 Subject: [PATCH 1/5] Constrain fixed now unconstrains first --- GPy/core/parameterised.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py index d1abb9c3..7afeb1af 100644 --- a/GPy/core/parameterised.py +++ b/GPy/core/parameterised.py @@ -223,7 +223,15 @@ class parameterised(object): To fix multiple parameters to the same value, simply pass a regular expression which matches both parameter names, or pass both of the indexes """ matches = self.grep_param_names(regexp) - assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained" + + overlap = set(matches).intersection(set(self.all_constrained_indices())) + if overlap: + self.unconstrain(np.asarray(list(overlap))) + print 'Warning: re-constraining these parameters' + pn = self._get_param_names() + for i in overlap: + print pn[i] + self.fixed_indices.append(matches) if value != None: self.fixed_values.append(value) From dea4359b4e97fac39701a8030e4cbc4928ce2180 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 5 Jun 2013 12:08:33 +0100 Subject: [PATCH 2/5] changed the behaviour of checkgrad checkgrad usd to check the passed string (for name matching) against the list of _get_param_names(). Then it would index along _get_param_names_transformed()! this led to inconsistensies when fixed or tied variables were used, which screwed up the ordering of the variable names. We now match against _get_param_names_transformed(). --- GPy/core/model.py | 6 ++++- GPy/core/parameterised.py | 23 ++++++++++++++---- GPy/examples/regression.py | 48 +++++++++++++++++++------------------- 3 files changed, 48 insertions(+), 29 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 2acb9963..f9dcaae4 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -392,7 +392,11 @@ class model(parameterised): if target_param is None: param_list = range(len(x)) else: - param_list = self.grep_param_names(target_param) + param_list = self.grep_param_names(target_param, transformed=True, search=True) + if not param_list: + print "No free parameters to check" + return + for i in param_list: xx = x.copy() diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py index d1abb9c3..c8d8ce4d 100644 --- a/GPy/core/parameterised.py +++ b/GPy/core/parameterised.py @@ -119,7 +119,7 @@ class parameterised(object): """Unties all parameters by setting tied_indices to an empty list.""" self.tied_indices = [] - def grep_param_names(self, regexp): + def grep_param_names(self, regexp, transformed=False, search=False): """ :param regexp: regular expression to select parameter names :type regexp: re | str | int @@ -129,13 +129,21 @@ class parameterised(object): Other objects are passed through - i.e. integers which weren't meant for grepping """ + if transformed: + names = self._get_param_names_transformed() + else: + names = self._get_param_names() + if type(regexp) in [str, np.string_, np.str]: regexp = re.compile(regexp) - return np.nonzero([regexp.match(name) for name in self._get_param_names()])[0] elif type(regexp) is re._pattern_type: - return np.nonzero([regexp.match(name) for name in self._get_param_names()])[0] + pass else: return regexp + if search: + return np.nonzero([regexp.search(name) for name in names])[0] + else: + return np.nonzero([regexp.match(name) for name in names])[0] def Nparam_transformed(self): removed = 0 @@ -223,7 +231,14 @@ class parameterised(object): To fix multiple parameters to the same value, simply pass a regular expression which matches both parameter names, or pass both of the indexes """ matches = self.grep_param_names(regexp) - assert not np.any(matches[:, None] == self.all_constrained_indices()), "Some indices are already constrained" + overlap = set(matches).intersection(set(self.all_constrained_indices())) + if overlap: + self.unconstrain(np.asarray(list(overlap))) + print 'Warning: re-constraining these parameters' + pn = self._get_param_names() + for i in overlap: + print pn[i] + self.fixed_indices.append(matches) if value != None: self.fixed_values.append(value) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index fd2e85d4..880ea191 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -10,7 +10,7 @@ import numpy as np import GPy -def toy_rbf_1d(max_nb_eval_optim=100): +def toy_rbf_1d(optim_iters=100): """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" data = GPy.util.datasets.toy_rbf_1d() @@ -19,13 +19,13 @@ def toy_rbf_1d(max_nb_eval_optim=100): # optimize m.ensure_default_constraints() - m.optimize(max_f_eval=max_nb_eval_optim) + m.optimize(max_f_eval=optim_iters) # plot m.plot() print(m) return m -def rogers_girolami_olympics(max_nb_eval_optim=100): +def rogers_girolami_olympics(optim_iters=100): """Run a standard Gaussian process regression on the Rogers and Girolami olympics data.""" data = GPy.util.datasets.rogers_girolami_olympics() @@ -37,14 +37,14 @@ def rogers_girolami_olympics(max_nb_eval_optim=100): # optimize m.ensure_default_constraints() - m.optimize(max_f_eval=max_nb_eval_optim) + m.optimize(max_f_eval=optim_iters) # plot m.plot(plot_limits = (1850, 2050)) print(m) return m -def toy_rbf_1d_50(max_nb_eval_optim=100): +def toy_rbf_1d_50(optim_iters=100): """Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance.""" data = GPy.util.datasets.toy_rbf_1d_50() @@ -53,14 +53,14 @@ def toy_rbf_1d_50(max_nb_eval_optim=100): # optimize m.ensure_default_constraints() - m.optimize(max_f_eval=max_nb_eval_optim) + m.optimize(max_f_eval=optim_iters) # plot m.plot() print(m) return m -def silhouette(max_nb_eval_optim=100): +def silhouette(optim_iters=100): """Predict the pose of a figure given a silhouette. This is a task from Agarwal and Triggs 2004 ICML paper.""" data = GPy.util.datasets.silhouette() @@ -69,12 +69,12 @@ def silhouette(max_nb_eval_optim=100): # optimize m.ensure_default_constraints() - m.optimize(messages=True,max_f_eval=max_nb_eval_optim) + m.optimize(messages=True,max_f_eval=optim_iters) print(m) return m -def coregionalisation_toy2(max_nb_eval_optim=100): +def coregionalisation_toy2(optim_iters=100): """ A simple demonstration of coregionalisation on two sinusoidal functions. """ @@ -93,7 +93,7 @@ def coregionalisation_toy2(max_nb_eval_optim=100): m.constrain_fixed('.*rbf_var',1.) #m.constrain_positive('.*kappa') m.ensure_default_constraints() - m.optimize('sim',messages=1,max_f_eval=max_nb_eval_optim) + m.optimize('sim',messages=1,max_f_eval=optim_iters) pb.figure() Xtest1 = np.hstack((np.linspace(0,9,100)[:,None],np.zeros((100,1)))) @@ -106,7 +106,7 @@ def coregionalisation_toy2(max_nb_eval_optim=100): pb.plot(X2[:,0],Y2[:,0],'gx',mew=2) return m -def coregionalisation_toy(max_nb_eval_optim=100): +def coregionalisation_toy(optim_iters=100): """ A simple demonstration of coregionalisation on two sinusoidal functions. """ @@ -125,7 +125,7 @@ def coregionalisation_toy(max_nb_eval_optim=100): m.constrain_fixed('.*rbf_var',1.) #m.constrain_positive('kappa') m.ensure_default_constraints() - m.optimize(max_f_eval=max_nb_eval_optim) + m.optimize(max_f_eval=optim_iters) pb.figure() Xtest1 = np.hstack((np.linspace(0,9,100)[:,None],np.zeros((100,1)))) @@ -139,7 +139,7 @@ def coregionalisation_toy(max_nb_eval_optim=100): return m -def coregionalisation_sparse(max_nb_eval_optim=100): +def coregionalisation_sparse(optim_iters=100): """ A simple demonstration of coregionalisation on two sinusoidal functions using sparse approximations. """ @@ -162,9 +162,9 @@ def coregionalisation_sparse(max_nb_eval_optim=100): m.scale_factor = 10000. m.constrain_fixed('.*rbf_var',1.) #m.constrain_positive('kappa') - m.constrain_fixed('iip') + m.constrain_fixed('Iip') m.ensure_default_constraints() - m.optimize_restarts(5, robust=True, messages=1, max_f_eval=max_nb_eval_optim) + m.optimize_restarts(5, robust=True, messages=1, max_f_eval=optim_iters) pb.figure() Xtest1 = np.hstack((np.linspace(0,9,100)[:,None],np.zeros((100,1)))) @@ -181,7 +181,7 @@ def coregionalisation_sparse(max_nb_eval_optim=100): return m -def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000, max_nb_eval_optim=100): +def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000, optim_iters=100): """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher.""" # Contour over a range of length scales and signal/noise ratios. @@ -219,7 +219,7 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000 # optimize m.ensure_default_constraints() - m.optimize(xtol=1e-6, ftol=1e-6, max_f_eval=max_nb_eval_optim) + m.optimize(xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters) optim_point_x[1] = m.get('rbf_lengthscale') optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); @@ -266,7 +266,7 @@ def _contour_data(data, length_scales, log_SNRs, signal_kernel_call=GPy.kern.rbf lls.append(length_scale_lls) return np.array(lls) -def sparse_GP_regression_1D(N = 400, M = 5, max_nb_eval_optim=100): +def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100): """Run a 1D example of a sparse GP regression.""" # sample inputs and outputs X = np.random.uniform(-3.,3.,(N,1)) @@ -281,11 +281,11 @@ def sparse_GP_regression_1D(N = 400, M = 5, max_nb_eval_optim=100): m.ensure_default_constraints() m.checkgrad(verbose=1) - m.optimize('tnc', messages = 1, max_f_eval=max_nb_eval_optim) + m.optimize('tnc', messages = 1, max_f_eval=optim_iters) m.plot() return m -def sparse_GP_regression_2D(N = 400, M = 50, max_nb_eval_optim=100): +def sparse_GP_regression_2D(N = 400, M = 50, optim_iters=100): """Run a 2D example of a sparse GP regression.""" X = np.random.uniform(-3.,3.,(N,2)) Y = np.sin(X[:,0:1]) * np.sin(X[:,1:2])+np.random.randn(N,1)*0.05 @@ -306,12 +306,12 @@ def sparse_GP_regression_2D(N = 400, M = 50, max_nb_eval_optim=100): # optimize and plot pb.figure() - m.optimize('tnc', messages = 1, max_f_eval=max_nb_eval_optim) + m.optimize('tnc', messages = 1, max_f_eval=optim_iters) m.plot() print(m) return m -def uncertain_inputs_sparse_regression(max_nb_eval_optim=100): +def uncertain_inputs_sparse_regression(optim_iters=100): """Run a 1D example of a sparse GP regression with uncertain inputs.""" fig, axes = pb.subplots(1,2,figsize=(12,5)) @@ -327,7 +327,7 @@ def uncertain_inputs_sparse_regression(max_nb_eval_optim=100): # create simple GP model - no input uncertainty on this one m = GPy.models.sparse_GP_regression(X, Y, kernel=k, Z=Z) m.ensure_default_constraints() - m.optimize('scg', messages=1, max_f_eval=max_nb_eval_optim) + m.optimize('scg', messages=1, max_f_eval=optim_iters) m.plot(ax=axes[0]) axes[0].set_title('no input uncertainty') @@ -335,7 +335,7 @@ def uncertain_inputs_sparse_regression(max_nb_eval_optim=100): #the same model with uncertainty m = GPy.models.sparse_GP_regression(X, Y, kernel=k, Z=Z, X_variance=S) m.ensure_default_constraints() - m.optimize('scg', messages=1, max_f_eval=max_nb_eval_optim) + m.optimize('scg', messages=1, max_f_eval=optim_iters) m.plot(ax=axes[1]) axes[1].set_title('with input uncertainty') print(m) From aef3ba72abfde9344e7cfa2258f7af229440bf6d Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 5 Jun 2013 12:18:23 +0100 Subject: [PATCH 3/5] how did that happen? --- GPy/examples/regression.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 880ea191..310a0691 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -162,7 +162,7 @@ def coregionalisation_sparse(optim_iters=100): m.scale_factor = 10000. m.constrain_fixed('.*rbf_var',1.) #m.constrain_positive('kappa') - m.constrain_fixed('Iip') + m.constrain_fixed('iip') m.ensure_default_constraints() m.optimize_restarts(5, robust=True, messages=1, max_f_eval=optim_iters) From adb8a98cb36e2718a1b4eb5e2ae6f8fe1d8d1fe2 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 5 Jun 2013 13:08:51 +0100 Subject: [PATCH 4/5] fixed the multiple optima demo --- GPy/core/parameterised.py | 2 +- GPy/examples/regression.py | 51 +++++++++++++++++--------------------- 2 files changed, 24 insertions(+), 29 deletions(-) diff --git a/GPy/core/parameterised.py b/GPy/core/parameterised.py index c8d8ce4d..e091e820 100644 --- a/GPy/core/parameterised.py +++ b/GPy/core/parameterised.py @@ -85,7 +85,7 @@ class parameterised(object): else: return self._get_params()[matches] else: - raise AttributeError, "no parameter matches %s" % name + raise AttributeError, "no parameter matches %s" % regexp def __setitem__(self, name, val): """ diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index 310a0691..7989b3c1 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -181,7 +181,7 @@ def coregionalisation_sparse(optim_iters=100): return m -def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000, optim_iters=100): +def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000, optim_iters=300): """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher.""" # Contour over a range of length scales and signal/noise ratios. @@ -197,7 +197,7 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000 data['Y'] = data['Y'] - np.mean(data['Y']) lls = GPy.examples.regression._contour_data(data, length_scales, log_SNRs, GPy.kern.rbf) - pb.contour(length_scales, log_SNRs, np.exp(lls), 20) + pb.contour(length_scales, log_SNRs, np.exp(lls), 20, cmap=pb.cm.jet) ax = pb.gca() pb.xlabel('length scale') pb.ylabel('log_10 SNR') @@ -211,18 +211,20 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000 optim_point_y = np.empty(2) np.random.seed(seed=seed) for i in range(0, model_restarts): - kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + GPy.kern.white(1,variance=np.random.exponential(1.)) + #kern = GPy.kern.rbf(1, variance=np.random.exponential(1.), lengthscale=np.random.exponential(50.)) + kern = GPy.kern.rbf(1, variance=np.random.uniform(1e-3,1), lengthscale=np.random.uniform(5,50)) m = GPy.models.GP_regression(data['X'],data['Y'], kernel=kern) - optim_point_x[0] = m.get('rbf_lengthscale') - optim_point_y[0] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); + m['noise_variance'] = np.random.uniform(1e-3,1) + optim_point_x[0] = m['rbf_lengthscale'] + optim_point_y[0] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']); # optimize m.ensure_default_constraints() - m.optimize(xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters) + m.optimize('scg', xtol=1e-6, ftol=1e-6, max_f_eval=optim_iters) - optim_point_x[1] = m.get('rbf_lengthscale') - optim_point_y[1] = np.log10(m.get('rbf_variance')) - np.log10(m.get('white_variance')); + optim_point_x[1] = m['rbf_lengthscale'] + optim_point_y[1] = np.log10(m['rbf_variance']) - np.log10(m['noise_variance']); pb.arrow(optim_point_x[0], optim_point_y[0], optim_point_x[1]-optim_point_x[0], optim_point_y[1]-optim_point_y[0], label=str(i), head_length=1, head_width=0.5, fc='k', ec='k') models.append(m) @@ -231,39 +233,32 @@ def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000 ax.set_ylim(ylim) return (models, lls) -def _contour_data(data, length_scales, log_SNRs, signal_kernel_call=GPy.kern.rbf): +def _contour_data(data, length_scales, log_SNRs, kernel_call=GPy.kern.rbf): """Evaluate the GP objective function for a given data set for a range of signal to noise ratios and a range of lengthscales. :data_set: A data set from the utils.datasets director. :length_scales: a list of length scales to explore for the contour plot. :log_SNRs: a list of base 10 logarithm signal to noise ratios to explore for the contour plot. - :signal_kernel: a kernel to use for the 'signal' portion of the data.""" + :kernel: a kernel to use for the 'signal' portion of the data.""" lls = [] total_var = np.var(data['Y']) + kernel = kernel_call(1, variance=1., lengthscale=1.) + model = GPy.models.GP_regression(data['X'], data['Y'], kernel=kernel) for log_SNR in log_SNRs: - SNR = 10**log_SNR + SNR = 10.**log_SNR + noise_var = total_var/(1.+SNR) + signal_var = total_var - noise_var + model.kern['.*variance'] = signal_var + model['noise_variance'] = noise_var length_scale_lls = [] + for length_scale in length_scales: - noise_var = 1. - signal_var = SNR - noise_var = noise_var/(noise_var + signal_var)*total_var - signal_var = signal_var/(noise_var + signal_var)*total_var - - signal_kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale) - noise_kernel = GPy.kern.white(1, variance=noise_var) - kernel = signal_kernel + noise_kernel - K = kernel.K(data['X']) - total_var = (np.dot(np.dot(data['Y'].T,GPy.util.linalg.pdinv(K)[0]), data['Y'])/data['Y'].shape[0])[0,0] - noise_var *= total_var - signal_var *= total_var - - kernel = signal_kernel_call(1, variance=signal_var, lengthscale=length_scale) + GPy.kern.white(1, variance=noise_var) - - model = GPy.models.GP_regression(data['X'], data['Y'], kernel=kernel) - model.constrain_positive('') + model['.*lengthscale'] = length_scale length_scale_lls.append(model.log_likelihood()) + lls.append(length_scale_lls) + return np.array(lls) def sparse_GP_regression_1D(N = 400, M = 5, optim_iters=100): From d43f4c911cb0904c1fe1719f6cf99b0d513e2609 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 5 Jun 2013 13:58:04 +0100 Subject: [PATCH 5/5] preferred optimiser is now scg bydefault --- GPy/core/model.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index f9dcaae4..19e38080 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -23,7 +23,7 @@ class model(parameterised): self.priors = None self.optimization_runs = [] self.sampling_runs = [] - self.preferred_optimizer = 'tnc' + self.preferred_optimizer = 'scg' #self._set_params(self._get_params()) has been taken out as it should only be called on leaf nodes def _get_params(self): raise NotImplementedError, "this needs to be implemented to use the model class"