diff --git a/GPy/core/model.py b/GPy/core/model.py index 2cdecdf9..4b569f98 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -239,7 +239,7 @@ class Model(Parameterized): print 'nothing to optimize' if not self.update_model(): - print "setting updates on again" + print "Updates were off, setting updates on again" self.update_model(True) if start == None: @@ -400,7 +400,11 @@ class Model(Parameterized): ['Log-likelihood', '{}
'.format(float(self.log_likelihood()))], ["Number of Parameters", '{}
'.format(self.size)]] from operator import itemgetter - to_print = [""] + ["{}: {}".format(name, detail) for name, detail in model_details] + ["
Parameters:"] + to_print = ["""\n"""] + ["

"] + ["{}: {}".format(name, detail) for name, detail in model_details] + ["

"] to_print.append(super(Model, self)._repr_html_()) return "\n".join(to_print) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index c7d6be5d..78bc4fa2 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -264,15 +264,21 @@ class Param(Parameterizable, ObsAr): ties = [' '.join(map(lambda x: x, t)) for t in ties] header_format = """ - {i} - {x} - {c} - {p} - {t} + {i} + {x} + {c} + {p} + {t} """ header = header_format.format(x=self.hierarchy_name(), c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # nice header for printing if not ties: ties = itertools.cycle(['']) - return "\n".join([''] + [header] + ["".format(x=x, c=" ".join(map(str, c)), p=" ".join(map(str, p)), t=(t or ''), i=i) for i, x, c, t, p in itertools.izip(indices, vals, constr_matrix, ties, prirs)] + ["
{i}{x}{c}{p}{t}
"]) + return "\n".join([""""""] + [''] + [header] + ["".format(x=x, c=" ".join(map(str, c)), p=" ".join(map(str, p)), t=(t or ''), i=i) for i, x, c, t, p in itertools.izip(indices, vals, constr_matrix, ties, prirs)] + ["
{i}{x}{c}{p}{t}
"]) def __str__(self, constr_matrix=None, indices=None, prirs=None, ties=None, lc=None, lx=None, li=None, lp=None, lt=None, only_name=False): filter_ = self._current_slice_ diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 897c53e3..ef44cfa9 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -377,7 +377,7 @@ class Parameterized(Parameterizable): cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]]) tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]]) pl = max([len(str(x)) if x else 0 for x in prirs + ["Prior"]]) - format_spec = "{{name:<{0}s}}{{desc:>{1}s}}{{const:^{2}s}}{{pri:^{3}s}}{{t:^{4}s}}".format(nl, sl, cl, pl, tl) + format_spec = "{{name:<{0}s}}{{desc:>{1}s}}{{const:^{2}s}}{{pri:^{3}s}}{{t:^{4}s}}".format(nl, sl, cl, pl, tl) to_print = [] for n, d, c, t, p in itertools.izip(names, desc, constrs, ts, prirs): to_print.append(format_spec.format(name=n, desc=d, const=c, t=t, pri=p)) @@ -385,13 +385,21 @@ class Parameterized(Parameterizable): if header: header = """ - {name} - Value - Constraint - Prior - Tied to""".format(name=name) + {name} + Value + Constraint + Prior + Tied to +""".format(name=name) to_print.insert(0, header) - return '' + '\n'.format(sep).join(to_print) + '\n
' + style = """""" + return style + '\n' + '' + '\n'.format(sep).join(to_print) + '\n
' def __str__(self, header=True): name = adjust_name_for_printing(self.name) + "." diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index 291076a1..235b9d1c 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -243,6 +243,41 @@ class NormalNaturalThroughTheta(NormalTheta): dmuvar[self.mu_indices] -= 2*mu*dmuvar[self.var_indices] #======================================================================= + #======================================================================= + # This is by going through theta fully and then going into eta direction: + #dmu = dmuvar[self.mu_indices] + #dmuvar[self.var_indices] += dmu*mu*(var + 4/var) + #======================================================================= + return dmuvar # which is now the gradient multiplicator + + def __str__(self): + return "natgrad" + + +class NormalNaturalWhooot(NormalTheta): + _instances = [] + def __new__(cls, mu_indices, var_indices): + if cls._instances: + cls._instances[:] = [instance for instance in cls._instances if instance()] + for instance in cls._instances: + if np.all(instance().mu_indices==mu_indices, keepdims=False) and np.all(instance().var_indices==var_indices, keepdims=False): + return instance() + o = super(Transformation, cls).__new__(cls, mu_indices, var_indices) + cls._instances.append(weakref.ref(o)) + return cls._instances[-1]() + + def __init__(self, mu_indices, var_indices): + self.mu_indices = mu_indices + self.var_indices = var_indices + + def gradfactor(self, muvar, dmuvar): + #mu = muvar[self.mu_indices] + #var = muvar[self.var_indices] + + #======================================================================= + # This is just eta direction: + #dmuvar[self.mu_indices] -= 2*mu*dmuvar[self.var_indices] + #======================================================================= #======================================================================= # This is by going through theta fully and then going into eta direction: diff --git a/GPy/core/parameterization/updateable.py b/GPy/core/parameterization/updateable.py index daf07d3c..593f3c05 100644 --- a/GPy/core/parameterization/updateable.py +++ b/GPy/core/parameterization/updateable.py @@ -15,14 +15,6 @@ class Updateable(Observable): def __init__(self, *args, **kwargs): super(Updateable, self).__init__(*args, **kwargs) - @property - def updates(self): - raise DeprecationWarning("updates is now a function, see update(True|False|None)") - - @updates.setter - def updates(self, ups): - raise DeprecationWarning("updates is now a function, see update(True|False|None)") - def update_model(self, updates=None): """ Get or set, whether automatic updates are performed. When updates are diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index beb69138..51dbd5db 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -121,12 +121,15 @@ class SparseGP(GP): Kxx = kern.Kdiag(Xnew) var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T else: - Kx = kern.psi1(self.Z, Xnew) - mu = np.dot(Kx, self.posterior.woodbury_vector) + Kx = kern.psi1(self.Z, Xnew).T + mu = np.dot(Kx.T, self.posterior.woodbury_vector) if full_cov: - raise NotImplementedError, "TODO" + Kxx = kern.K(Xnew.mean) + if self.posterior.woodbury_inv.ndim == 2: + var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx)) + elif self.posterior.woodbury_inv.ndim == 3: + var = Kxx[:,:,None] - np.tensordot(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx).T, Kx, [1,0]).swapaxes(1,2) else: Kxx = kern.psi0(self.Z, Xnew) - psi2 = kern.psi2(self.Z, Xnew) - var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1) + var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T return mu, var diff --git a/GPy/kern/_src/mlp.py b/GPy/kern/_src/mlp.py index badbd60d..16e84363 100644 --- a/GPy/kern/_src/mlp.py +++ b/GPy/kern/_src/mlp.py @@ -79,8 +79,14 @@ class MLP(Kern): + 2*self.bias_variance + 2.))*base_cov_grad).sum() def update_gradients_diag(self, X): - raise NotImplementedError, "TODO" - + self._K_diag_computations(X) + self.variance.gradient = np.sum(self._K_diag_dvar*dL_dKdiag) + + base = four_over_tau*self.variance/np.sqrt(1-self._K_diag_asin_arg*self._K_diag_asin_arg) + base_cov_grad = base*dL_dKdiag/np.square(self._K_diag_denom) + + self.weight_variance.gradient = (base_cov_grad*np.square(X).sum(axis=1)).sum() + self.bias_variance.gradient = base_cov_grad.sum() def gradients_X(self, dL_dK, X, X2): """Derivative of the covariance matrix with respect to X""" diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 443871af..06671b23 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -159,7 +159,7 @@ class Stationary(Kern): #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 tmp = dL_dr*self._inv_dist(X, X2) if X2 is None: X2 = X - + if config.getboolean('weave', 'working'): try: @@ -261,7 +261,7 @@ class Stationary(Kern): ret(n,d) = retnd; } } - + """ if hasattr(X, 'values'):X = X.values #remove the GPy wrapping to make passing into weave safe if hasattr(X2, 'values'):X2 = X2.values @@ -278,12 +278,12 @@ class Stationary(Kern): 'extra_link_args' : ['-lgomp']} weave.inline(code, ['ret', 'N', 'D', 'M', 'tmp', 'X', 'X2'], type_converters=weave.converters.blitz, support_code=support_code, **weave_options) return ret/self.lengthscale**2 - + def gradients_X_diag(self, dL_dKdiag, X): return np.zeros(X.shape) def input_sensitivity(self, summarize=True): - return np.ones(self.input_dim)/self.lengthscale**2 + return self.variance*np.ones(self.input_dim)/self.lengthscale**2 class Exponential(Stationary): def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, active_dims=None, name='Exponential'): diff --git a/GPy/models/bayesian_gplvm_minibatch.py b/GPy/models/bayesian_gplvm_minibatch.py index f164b466..64aed246 100644 --- a/GPy/models/bayesian_gplvm_minibatch.py +++ b/GPy/models/bayesian_gplvm_minibatch.py @@ -8,6 +8,7 @@ from ..core.parameterization.variational import NormalPosterior, NormalPrior from ..inference.latent_function_inference.var_dtc_parallel import VarDTC_minibatch import logging from GPy.models.sparse_gp_minibatch import SparseGPMiniBatch +from GPy.core.parameterization.param import Param class BayesianGPLVMMiniBatch(SparseGPMiniBatch): """ @@ -35,15 +36,20 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): self.init = init - if X_variance is None: - self.logger.info("initializing latent space variance ~ uniform(0,.1)") - X_variance = np.random.uniform(0,.1,X.shape) - if Z is None: self.logger.info("initializing inducing inputs") Z = np.random.permutation(X.copy())[:num_inducing] assert Z.shape[1] == X.shape[1] + if X_variance == False: + self.logger.info('no variance on X, activating sparse GPLVM') + X = Param("latent space", X) + elif X_variance is None: + self.logger.info("initializing latent space variance ~ uniform(0,.1)") + X_variance = np.random.uniform(0,.1,X.shape) + self.variational_prior = NormalPrior() + X = NormalPosterior(X, X_variance) + if kernel is None: self.logger.info("initializing kernel RBF") kernel = kern.RBF(input_dim, lengthscale=1./fracs, ARD=True) #+ kern.Bias(input_dim) + kern.White(input_dim) @@ -51,9 +57,6 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): if likelihood is None: likelihood = Gaussian() - self.variational_prior = NormalPrior() - X = NormalPosterior(X, X_variance) - self.kl_factr = 1. if inference_method is None: @@ -83,36 +86,42 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): def _inner_parameters_changed(self, kern, X, Z, likelihood, Y, Y_metadata, Lm=None, dL_dKmm=None, subset_indices=None): posterior, log_marginal_likelihood, grad_dict, current_values, value_indices = super(BayesianGPLVMMiniBatch, self)._inner_parameters_changed(kern, X, Z, likelihood, Y, Y_metadata, Lm=Lm, dL_dKmm=dL_dKmm, subset_indices=subset_indices) - current_values['meangrad'], current_values['vargrad'] = self.kern.gradients_qX_expectations( - variational_posterior=X, - Z=Z, dL_dpsi0=grad_dict['dL_dpsi0'], - dL_dpsi1=grad_dict['dL_dpsi1'], - dL_dpsi2=grad_dict['dL_dpsi2']) + if self.has_uncertain_inputs(): + current_values['meangrad'], current_values['vargrad'] = self.kern.gradients_qX_expectations( + variational_posterior=X, + Z=Z, dL_dpsi0=grad_dict['dL_dpsi0'], + dL_dpsi1=grad_dict['dL_dpsi1'], + dL_dpsi2=grad_dict['dL_dpsi2']) + else: + current_values['Xgrad'] = self.kern.gradients_X(grad_dict['dL_dKnm'], X, Z) + current_values['Xgrad'] += self.kern.gradients_X_diag(grad_dict['dL_dKdiag'], X) + if subset_indices is not None: + value_indices['Xgrad'] = subset_indices['samples'] kl_fctr = self.kl_factr - if self.missing_data: - d = self.output_dim - log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X)/d - else: - log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X) + if self.has_uncertain_inputs(): + if self.missing_data: + d = self.output_dim + log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X)/d + else: + log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(X) + # Subsetting Variational Posterior objects, makes the gradients + # empty. We need them to be 0 though: + X.mean.gradient[:] = 0 + X.variance.gradient[:] = 0 - # Subsetting Variational Posterior objects, makes the gradients - # empty. We need them to be 0 though: - X.mean.gradient[:] = 0 - X.variance.gradient[:] = 0 + self.variational_prior.update_gradients_KL(X) + if self.missing_data: + current_values['meangrad'] += kl_fctr*X.mean.gradient/d + current_values['vargrad'] += kl_fctr*X.variance.gradient/d + else: + current_values['meangrad'] += kl_fctr*X.mean.gradient + current_values['vargrad'] += kl_fctr*X.variance.gradient - self.variational_prior.update_gradients_KL(X) - if self.missing_data: - current_values['meangrad'] += kl_fctr*X.mean.gradient/d - current_values['vargrad'] += kl_fctr*X.variance.gradient/d - else: - current_values['meangrad'] += kl_fctr*X.mean.gradient - current_values['vargrad'] += kl_fctr*X.variance.gradient - - if subset_indices is not None: - value_indices['meangrad'] = subset_indices['samples'] - value_indices['vargrad'] = subset_indices['samples'] + if subset_indices is not None: + value_indices['meangrad'] = subset_indices['samples'] + value_indices['vargrad'] = subset_indices['samples'] return posterior, log_marginal_likelihood, grad_dict, current_values, value_indices def _outer_values_update(self, full_values): @@ -121,42 +130,24 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch): E.g. set the gradients of parameters, etc. """ super(BayesianGPLVMMiniBatch, self)._outer_values_update(full_values) - self.X.mean.gradient = full_values['meangrad'] - self.X.variance.gradient = full_values['vargrad'] + if self.has_uncertain_inputs(): + self.X.mean.gradient = full_values['meangrad'] + self.X.variance.gradient = full_values['vargrad'] + else: + self.X.gradient = full_values['Xgrad'] def _outer_init_full_values(self): - return dict(meangrad=np.zeros(self.X.mean.shape), - vargrad=np.zeros(self.X.variance.shape)) + if self.has_uncertain_inputs(): + return dict(meangrad=np.zeros(self.X.mean.shape), + vargrad=np.zeros(self.X.variance.shape)) + else: + return dict(Xgrad=np.zeros(self.X.shape)) def parameters_changed(self): super(BayesianGPLVMMiniBatch,self).parameters_changed() if isinstance(self.inference_method, VarDTC_minibatch): return - #super(BayesianGPLVM, self).parameters_changed() - #self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) - - #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) - def plot_latent(self, labels=None, which_indices=None, resolution=50, ax=None, marker='o', s=40, fignum=None, plot_inducing=True, legend=True, diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 645cdf88..f3e643c9 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -111,9 +111,6 @@ class MRD(BayesianGPLVMMiniBatch): assert all([isinstance(k, Kern) for k in kernel]), "invalid kernel object detected!" kernels = kernel - if X_variance is None: - X_variance = np.random.uniform(0.1, 0.2, X.shape) - self.variational_prior = NormalPrior() #self.X = NormalPosterior(X, X_variance) @@ -174,10 +171,13 @@ class MRD(BayesianGPLVMMiniBatch): self.Z.gradient[:] += b.full_values['Zgrad'] grad_dict = b.full_values - self.X.mean.gradient += grad_dict['meangrad'] - self.X.variance.gradient += grad_dict['vargrad'] + if self.has_uncertain_inputs(): + self.X.mean.gradient += grad_dict['meangrad'] + self.X.variance.gradient += grad_dict['vargrad'] + else: + self.X.gradient += grad_dict['Xgrad'] - if isinstance(self.X, VariationalPosterior): + if self.has_uncertain_inputs(): # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) diff --git a/setup.py b/setup.py index c4963bcc..0562c9d8 100644 --- a/setup.py +++ b/setup.py @@ -5,7 +5,7 @@ import os from setuptools import setup # Version number -version = '0.6.0' +version = '0.6.1' def read(fname): return open(os.path.join(os.path.dirname(__file__), fname)).read()