From f50fac810418a8d56e2ecf1729f6eca14f346f0e Mon Sep 17 00:00:00 2001 From: Ricardo Date: Mon, 3 Feb 2014 15:29:06 +0000 Subject: [PATCH 01/51] minor changes --- GPy/kern/constructors.py | 8 ++++---- GPy/kern/kern.py | 4 ++-- GPy/kern/parts/prod.py | 25 +++++++++++++------------ 3 files changed, 19 insertions(+), 18 deletions(-) diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index 11da5e9b..a362aff8 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -411,7 +411,7 @@ def periodic_Matern52(input_dim, variance=1., lengthscale=None, period=2 * np.pi part = parts.periodic_Matern52.PeriodicMatern52(input_dim, variance, lengthscale, period, n_freq, lower, upper) return kern(input_dim, [part]) -def prod(k1,k2,tensor=False): +def prod(k1,k2,**kwargs): """ Construct a product kernel over input_dim from two kernels over input_dim @@ -422,7 +422,7 @@ def prod(k1,k2,tensor=False): :rtype: kernel object """ - part = parts.prod.Prod(k1, k2, tensor) + part = parts.prod.Prod(k1, k2, **kwargs) return kern(part.input_dim, [part]) def symmetric(k): @@ -433,7 +433,7 @@ def symmetric(k): k_.parts = [symmetric.Symmetric(p) for p in k.parts] return k_ -def coregionalize(output_dim,rank=1, W=None, kappa=None): +def coregionalize(output_dim,rank=1, W=None, kappa=None,name='coregion'): """ Coregionlization matrix B, of the form: @@ -458,7 +458,7 @@ def coregionalize(output_dim,rank=1, W=None, kappa=None): :rtype: kernel object """ - p = parts.coregionalize.Coregionalize(output_dim,rank,W,kappa) + p = parts.coregionalize.Coregionalize(output_dim,rank,W,kappa,name) return kern(1,[p]) diff --git a/GPy/kern/kern.py b/GPy/kern/kern.py index de87ff14..6c24ab7d 100644 --- a/GPy/kern/kern.py +++ b/GPy/kern/kern.py @@ -170,7 +170,7 @@ class kern(Parameterized): """ return self.prod(other, tensor=True) - def prod(self, other, tensor=False): + def prod(self, other, tensor=False, **kwargs): """ Multiply two kernels (either on the same space, or on the tensor product of the input space). @@ -189,7 +189,7 @@ class kern(Parameterized): s1[sl1], s2[sl2] = [True], [True] slices += [s1 + s2] - newkernparts = [prod(k1, k2, tensor) for k1, k2 in itertools.product(K1._parameters_, K2._parameters_)] + newkernparts = [prod(k1, k2, tensor, **kwargs) for k1, k2 in itertools.product(K1._parameters_, K2._parameters_)] if tensor: newkern = kern(K1.input_dim + K2.input_dim, newkernparts, slices) diff --git a/GPy/kern/parts/prod.py b/GPy/kern/parts/prod.py index 2569c51c..3e624756 100644 --- a/GPy/kern/parts/prod.py +++ b/GPy/kern/parts/prod.py @@ -17,21 +17,21 @@ class Prod(Kernpart): :rtype: kernel object """ - def __init__(self,k1,k2,tensor=False): + def __init__(self,k1,k2,tensor=False,name="product"): if tensor: - super(Prod, self).__init__(k1.input_dim + k2.input_dim, '['+k1.name + '**' + k2.name +']') + super(Prod, self).__init__(k1.input_dim + k2.input_dim, name) else: assert k1.input_dim == k2.input_dim, "Error: The input spaces of the kernels to sum don't have the same dimension." - super(Prod, self).__init__(k1.input_dim, '['+k1.name + '*' + k2.name +']') + super(Prod, self).__init__(k1.input_dim, name) #self.num_params = k1.num_params + k2.num_params self.k1 = k1 self.k2 = k2 -# if tensor: -# self.slice1 = slice(0,self.k1.input_dim) -# self.slice2 = slice(self.k1.input_dim,self.k1.input_dim+self.k2.input_dim) -# else: -# self.slice1 = slice(0,self.input_dim) -# self.slice2 = slice(0,self.input_dim) + if tensor: + self.slice1 = slice(0,self.k1.input_dim) + self.slice2 = slice(self.k1.input_dim,self.k1.input_dim+self.k2.input_dim) + else: + self.slice1 = slice(0,self.input_dim) + self.slice2 = slice(0,self.input_dim) self._X, self._X2 = np.empty(shape=(2,1)) self.add_parameters(self.k1, self.k2) @@ -117,18 +117,19 @@ class Prod(Kernpart): def _K_computations(self,X,X2): if not (np.array_equal(X,self._X) and np.array_equal(X2,self._X2) and np.array_equal(self._params , self._get_params())): - self._X = X.copy() - self._params == self._get_params().copy() + #self._X = X.copy() + #self._params == self._get_params().copy() if X2 is None: self._X2 = None self._K1 = np.zeros((X.shape[0],X.shape[0])) self._K2 = np.zeros((X.shape[0],X.shape[0])) self.k1.K(X[:,self.slice1],None,self._K1) self.k2.K(X[:,self.slice2],None,self._K2) + #self.k1.K(X[:,self.k1.input_slices],None,self._K1) + #self.k2.K(X[:,self.k2_input_slices],None,self._K2) else: self._X2 = X2.copy() self._K1 = np.zeros((X.shape[0],X2.shape[0])) self._K2 = np.zeros((X.shape[0],X2.shape[0])) self.k1.K(X[:,self.slice1],X2[:,self.slice1],self._K1) self.k2.K(X[:,self.slice2],X2[:,self.slice2],self._K2) - From 80acca640f80b6d596c13301bf529469dd527af1 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 25 Feb 2014 17:15:38 +0000 Subject: [PATCH 02/51] messing with kernels --- GPy/core/sparse_gp.py | 30 +++++++++++++++++++++++---- GPy/kern/_src/add.py | 12 +++++------ GPy/kern/_src/kern.py | 39 +++++++++++++++--------------------- GPy/kern/_src/linear.py | 13 ++++-------- GPy/kern/_src/rbf.py | 11 +++------- GPy/kern/_src/static.py | 12 +++++------ GPy/kern/_src/stationary.py | 4 ++++ GPy/models/bayesian_gplvm.py | 2 +- 8 files changed, 66 insertions(+), 57 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 00a80c7b..7677fea2 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -58,11 +58,33 @@ class SparseGP(GP): self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.Z, self.likelihood, self.Y) self.likelihood.update_gradients(self.grad_dict.pop('partial_for_likelihood')) if isinstance(self.X, VariationalPosterior): - self.kern.update_gradients_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) - self.Z.gradient = self.kern.gradients_Z_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) + #gradients wrt kernel + dL_dKmm = self.grad_dict.pop('dL_dKmm') + self.kern.update_gradients_full(dL_dKmm, self.Z, None) + target = np.zeros(self.kern.size) + self.kern._collect_gradient(target) + self.kern.update_gradients_expectations(variational_posterior=self.X, Z=self.Z, **self.grad_dict) + self.kern._collect_gradient(target) + self.kern._set_gradient(target) + + #gradients wrt Z + self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) + self.Z.gradient += self.kern.gradients_Z_expectations( + self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpis2'], Z=self.Z, variational_posterior=self.X) else: - self.kern.update_gradients_sparse(X=self.X, Z=self.Z, **self.grad_dict) - self.Z.gradient = self.kern.gradients_Z_sparse(X=self.X, Z=self.Z, **self.grad_dict) + #gradients wrt kernel + target = np.zeros(self.kern.size) + self.kern.update_gradients_diag(self.grad_dict['dL_dKdiag'], self.X) + self.kern._collect_gradient(target) + self.kern.update_gradients_full(self.grad_dict['dL_dKnm'], self.X, self.Z) + self.kern._collect_gradient(target) + self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None) + self.kern._collect_gradient(target) + self.kern._set_gradient(target) + + #gradients wrt Z + self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) + self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) def _raw_predict(self, Xnew, X_variance_new=None, full_cov=False): """ diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 1022124d..3e52d5af 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -101,7 +101,7 @@ class Add(Kern): raise NotImplementedError, "psi2 cannot be computed for this kernel" return psi2 - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, variational_posterior, Z): + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from white import White from rbf import RBF #from rbf_inv import RBFInv @@ -124,10 +124,10 @@ class Add(Kern): eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], mu[:,is2], S[:,is2]) * 2. - p1.update_gradients_variational(dL_dKmm, dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, mu[:,is1], S[:,is1], Z[:,is1]) + p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, mu[:,is1], S[:,is1], Z[:,is1]) - def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z): + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from white import White from rbf import RBF #from rbf_inv import rbfinv @@ -151,10 +151,10 @@ class Add(Kern): eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], mu[:,is2], S[:,is2]) * 2. - target += p1.gradients_z_variational(dL_dKmm, dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, mu[:,is1], S[:,is1], Z[:,is1]) + target += p1.gradients_z_variational(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, mu[:,is1], S[:,is1], Z[:,is1]) return target - def gradients_muS_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, mu, S, Z): + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): from white import white from rbf import rbf #from rbf_inv import rbfinv @@ -179,7 +179,7 @@ class Add(Kern): eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(z[:,is2], mu[:,is2], s[:,is2]) * 2. - a, b = p1.gradients_muS_variational(dL_dkmm, dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, mu[:,is1], s[:,is1], z[:,is1]) + a, b = p1.gradients_qX_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, mu[:,is1], s[:,is1], z[:,is1]) target_mu += a target_S += b return target_mu, target_S diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 6b23a69e..2e412688 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -39,28 +39,21 @@ class Kern(Parameterized): def update_gradients_full(self, dL_dK, X, X2): """Set the gradients of all parameters when doing full (N) inference.""" raise NotImplementedError - def update_gradients_sparse(self, dL_dKmm, dL_dKnm, dL_dKdiag, X, Z): - target = np.zeros(self.size) - self.update_gradients_diag(dL_dKdiag, X) - self._collect_gradient(target) - self.update_gradients_full(dL_dKnm, X, Z) - self._collect_gradient(target) - self.update_gradients_full(dL_dKmm, Z, None) - self._collect_gradient(target) - self._set_gradient(target) + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + """ + Set the gradients of all parameters when doing inference with + uncertain inputs, using expectations of the kernel. + """ + raise NotImplementedError + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + raise NotImplementedError + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + """ + Compute the gradients wrt the parameters of the variational + distruibution q(X), chain-ruling via the expectations of the kernel + """ + raise NotImplementedError - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - """Set the gradients of all parameters when doing variational (M) inference with uncertain inputs.""" - raise NotImplementedError - def gradients_Z_sparse(self, dL_dKmm, dL_dKnm, dL_dKdiag, X, Z): - grad = self.gradients_X(dL_dKmm, Z) - grad += self.gradients_X(dL_dKnm.T, Z, X) - return grad - def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - raise NotImplementedError - def gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - raise NotImplementedError - def plot_ARD(self, *args, **kw): if "matplotlib" in sys.modules: from ...plotting.matplot_dep import kernel_plots @@ -68,13 +61,13 @@ class Kern(Parameterized): assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ...plotting.matplot_dep import kernel_plots return kernel_plots.plot_ARD(self,*args,**kw) - + def input_sensitivity(self): """ Returns the sensitivity for each dimension of this kernel. """ return np.zeros(self.input_dim) - + def __add__(self, other): """ Overloading of the '+' operator. for more control, see self.add """ return self.add(other) diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py index 1d4f4611..e503180a 100644 --- a/GPy/kern/_src/linear.py +++ b/GPy/kern/_src/linear.py @@ -117,7 +117,7 @@ class Linear(Kern): ZAinner = self._ZAinner(variational_posterior, Z) return np.dot(ZAinner, ZA.T) - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, variational_posterior, Z): + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): mu, S = variational_posterior.mean, variational_posterior.variance # psi0: tmp = dL_dpsi0[:, None] * self._mu2S(variational_posterior) @@ -130,20 +130,15 @@ class Linear(Kern): tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(variational_posterior, Z)[:, :, None, :] * (2. * Z)[None, None, :, :]) if self.ARD: grad += tmp.sum(0).sum(0).sum(0) else: grad += tmp.sum() - #from Kmm - self.update_gradients_full(dL_dKmm, Z, None) - self.variances.gradient += grad - def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, variational_posterior, Z): - # Kmm - grad = self.gradients_X(dL_dKmm, Z, None) + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): #psi1 - grad += self.gradients_X(dL_dpsi1.T, Z, variational_posterior.mean) + grad = self.gradients_X(dL_dpsi1.T, Z, variational_posterior.mean) #psi2 self._weave_dpsi2_dZ(dL_dpsi2, Z, variational_posterior, grad) return grad - def gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, variational_posterior, Z): + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): grad_mu, grad_S = np.zeros(variational_posterior.mean.shape), np.zeros(variational_posterior.mean.shape) # psi0 grad_mu += dL_dpsi0[:, None] * (2.0 * variational_posterior.mean * self.variances) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index c80fb646..7c43b18d 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -55,10 +55,7 @@ class RBF(Stationary): self._psi_computations(Z, mu, S) return self._psi2 - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - #contributions from Kmm - sself.update_gradients_full(dL_dKmm, Z) - + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): mu = variational_posterior.mean S = variational_posterior.variance self._psi_computations(Z, mu, S) @@ -87,7 +84,7 @@ class RBF(Stationary): else: self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) - def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): mu = variational_posterior.mean S = variational_posterior.variance self._psi_computations(Z, mu, S) @@ -104,11 +101,9 @@ class RBF(Stationary): dZ = self._psi2[:, :, :, None] * (term1[None] + term2) grad += 2*(dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0) - grad += self.gradients_X(dL_dKmm, Z, None) - return grad - def gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): mu = variational_posterior.mean S = variational_posterior.variance self._psi_computations(Z, mu, S) diff --git a/GPy/kern/_src/static.py b/GPy/kern/_src/static.py index f4400ed7..135e3f9e 100644 --- a/GPy/kern/_src/static.py +++ b/GPy/kern/_src/static.py @@ -25,10 +25,10 @@ class Static(Kern): def gradients_X_diag(self, dL_dKdiag, X): return np.zeros(X.shape) - def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): return np.zeros(Z.shape) - def gradients_muS_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): return np.zeros(variational_posterior.shape), np.zeros(variational_posterior.shape) def psi0(self, Z, variational_posterior): @@ -61,8 +61,8 @@ class White(Static): def update_gradients_diag(self, dL_dKdiag, X): self.variance.gradient = dL_dKdiag.sum() - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - self.variance.gradient = np.trace(dL_dKmm) + dL_dpsi0.sum() + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + self.variance.gradient = dL_dpsi0.sum() class Bias(Static): @@ -86,6 +86,6 @@ class Bias(Static): ret[:] = self.variance**2 return ret - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - self.variance.gradient = dL_dKmm.sum() + dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum() + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + self.variance.gradient = dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum() diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index b998969c..2d0d284a 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -312,4 +312,8 @@ class RatQuad(Stationary): grad = np.sum(dL_dK*dK_dpow) self.power.gradient = grad + def update_gradients_diag(self, dL_dKdiag, X): + super(RatQuad, self).update_gradients_diag(dL_dKdiag, X) + self.power.gradient = 0. + diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 50fc2810..366995dc 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -66,7 +66,7 @@ class BayesianGPLVM(SparseGP): 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_q_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) + self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_q_variational(variational_posterior=self.X, Z=self.Z, **self.grad_dict) # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) From 51dca0fcbce9e9dec050ab6200d00ea2287525d0 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 08:21:14 +0000 Subject: [PATCH 03/51] ard plotting --- GPy/kern/_src/add.py | 4 +-- GPy/kern/_src/kern.py | 7 +++-- GPy/plotting/matplot_dep/base_plots.py | 37 +++++++++++++++++++----- GPy/plotting/matplot_dep/kernel_plots.py | 16 ++++------ 4 files changed, 41 insertions(+), 23 deletions(-) diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 45800dbf..a91a9c9e 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -196,9 +196,9 @@ class Add(Kern): kernel_plots.plot(self,*args) def input_sensitivity(self): - in_sen = np.zeros((self.input_dim, self.num_params)) + in_sen = np.zeros((self.num_params, self.input_dim)) for i, [p, i_s] in enumerate(zip(self._parameters_, self.input_slices)): - in_sen[i_s, i] = p.input_sensitivity() + in_sen[i, i_s] = p.input_sensitivity() return in_sen def _getstate(self): diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 1eec7af5..9f7dcb0a 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -62,9 +62,10 @@ class Kern(Parameterized): raise NotImplementedError def plot_ARD(self, *args, **kw): - if "matplotlib" in sys.modules: - from ...plotting.matplot_dep import kernel_plots - self.plot_ARD.__doc__ += kernel_plots.plot_ARD.__doc__ + """ + See :class:`~GPy.plotting.matplot_dep.kernel_plots` + """ + import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ...plotting.matplot_dep import kernel_plots return kernel_plots.plot_ARD(self,*args,**kw) diff --git a/GPy/plotting/matplot_dep/base_plots.py b/GPy/plotting/matplot_dep/base_plots.py index d5d4d6ee..a9d25223 100644 --- a/GPy/plotting/matplot_dep/base_plots.py +++ b/GPy/plotting/matplot_dep/base_plots.py @@ -6,27 +6,48 @@ import Tango import pylab as pb import numpy as np -def gpplot(x,mu,lower,upper,edgecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue'],axes=None,**kwargs): - if axes is None: - axes = pb.gca() +def ax_default(fignum, ax): + if ax is None: + fig = pb.figure(fignum) + ax = fig.add_subplot(111) + else: + fig = ax.figure + return fig, ax + +def meanplot(x, mu, color=Tango.colorsHex['darkBlue'], ax=None, fignum=None, linewidth=2,**kw): + _, axes = ax_default(fignum, ax) + #here's the mean + return axes.plot(x,mu,color=color,linewidth=linewidth,**kw) + +def gpplot(x,mu,lower,upper,edgecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue'],ax=None,fignum=None,xlabel='x',ylabel='y',**kwargs): + _, axes = ax_default(ax, fignum) + mu = mu.flatten() x = x.flatten() lower = lower.flatten() upper = upper.flatten() + plots = [] + #here's the mean - axes.plot(x,mu,color=edgecol,linewidth=2) + plots.append(meanplot(x, mu, edgecol, axes)) #here's the box kwargs['linewidth']=0.5 if not 'alpha' in kwargs.keys(): kwargs['alpha'] = 0.3 - axes.fill(np.hstack((x,x[::-1])),np.hstack((upper,lower[::-1])),color=fillcol,**kwargs) + plots.append(axes.fill(np.hstack((x,x[::-1])),np.hstack((upper,lower[::-1])),color=fillcol,**kwargs)) #this is the edge: - axes.plot(x,upper,color=edgecol,linewidth=0.2) - axes.plot(x,lower,color=edgecol,linewidth=0.2) - + plots.append(meanplot(x, upper,color=edgecol,linewidth=0.2,axes=axes)) + plots.append(meanplot(x, lower,color=edgecol,linewidth=0.2,axes=axes)) + + axes.set_xlabel(xlabel) + axes.set_ylabel(ylabel) + + return plots + + def removeRightTicks(ax=None): ax = ax or pb.gca() for i, line in enumerate(ax.get_yticklines()): diff --git a/GPy/plotting/matplot_dep/kernel_plots.py b/GPy/plotting/matplot_dep/kernel_plots.py index 6d4a7f0f..b55a0e53 100644 --- a/GPy/plotting/matplot_dep/kernel_plots.py +++ b/GPy/plotting/matplot_dep/kernel_plots.py @@ -6,7 +6,7 @@ import pylab as pb import Tango from matplotlib.textpath import TextPath from matplotlib.transforms import offset_copy -from ...kern import Linear +from .base_plots import ax_default @@ -52,11 +52,7 @@ def plot_ARD(kernel, fignum=None, ax=None, title='', legend=False): pass '' to not print a title pass None for a generic title """ - if ax is None: - fig = pb.figure(fignum) - ax = fig.add_subplot(111) - else: - fig = ax.figure + fig, ax = ax_default(fignum,ax) if title is None: ax.set_title('ARD parameters, %s kernel' % kernel.name) @@ -70,13 +66,13 @@ def plot_ARD(kernel, fignum=None, ax=None, title='', legend=False): bottom = 0 x = np.arange(kernel.input_dim) - for i in range(ard_params.shape[-1]): + for i in range(ard_params.shape[0]): c = Tango.nextMedium() - bars.append(plot_bars(fig, ax, x, ard_params[:,i], c, kernel._parameters_[i].name, bottom=bottom)) - bottom += ard_params[:,i] + bars.append(plot_bars(fig, ax, x, ard_params[i,:], c, kernel._parameters_[i].name, bottom=bottom)) + bottom += ard_params[i,:] ax.set_xlim(-.5, kernel.input_dim - .5) - add_bar_labels(fig, ax, [bars[-1]], bottom=bottom-ard_params[:,i]) + add_bar_labels(fig, ax, [bars[-1]], bottom=bottom-ard_params[i,:]) if legend: if title is '': From b928044f408944fad0d6d73ce054df2b1db289e5 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 08:23:08 +0000 Subject: [PATCH 04/51] sparse gp missing data --- GPy/core/sparse_gp.py | 25 +++++++++++-------- GPy/examples/dimensionality_reduction.py | 6 ++--- .../latent_function_inference/var_dtc.py | 6 ++--- GPy/models/bayesian_gplvm.py | 2 +- 4 files changed, 22 insertions(+), 17 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index bb3116ba..dfe2d38f 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -46,28 +46,33 @@ class SparseGP(GP): self.Z = Param('inducing inputs', Z) self.num_inducing = Z.shape[0] - self.q = NormalPosterior(X, X_variance) - - GP.__init__(self, self.q.mean, Y, kernel, likelihood, inference_method=inference_method, name=name) + if not (X_variance is None): + self.q = NormalPosterior(X, X_variance) + GP.__init__(self, self.q.mean, Y, kernel, likelihood, inference_method=inference_method, name=name) + else: + self.X = X + GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name) + self.add_parameter(self.Z, index=0) self.parameters_changed() def has_uncertain_inputs(self): - return self.q.has_uncertain_inputs() - + return hasattr(self, 'q') + def parameters_changed(self): if self.has_uncertain_inputs(): self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference_latent(self.kern, self.q, self.Z, self.likelihood, self.Y) - else: - self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, self.X_variance, self.Z, self.likelihood, self.Y) - self.likelihood.update_gradients(self.grad_dict.pop('partial_for_likelihood')) - if self.has_uncertain_inputs(): + # gradients + self.likelihood.update_gradients(self.grad_dict.pop('partial_for_likelihood')) self.kern.update_gradients_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict) self.Z.gradient = self.kern.gradients_Z_variational(posterior_variational=self.q, Z=self.Z, **self.grad_dict) else: + self.posterior, self._log_marginal_likelihood, self.grad_dict = self.inference_method.inference(self.kern, self.X, None, self.Z, self.likelihood, self.Y) + # gradients + self.likelihood.update_gradients(self.grad_dict.pop('partial_for_likelihood')) self.kern.update_gradients_sparse(X=self.X, Z=self.Z, **self.grad_dict) self.Z.gradient = self.kern.gradients_Z_sparse(X=self.X, Z=self.Z, **self.grad_dict) - + def _raw_predict(self, Xnew, X_variance_new=None, full_cov=False): """ Make a prediction for the latent function values diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index a2686d73..4345cb90 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -164,7 +164,7 @@ def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, _np.random.seed(0) data = GPy.util.datasets.oil() - kernel = GPy.kern.RBF(Q, 1., _np.random.uniform(0,1,(Q,)), ARD=True)# + GPy.kern.Bias(Q, _np.exp(-2)) + kernel = GPy.kern.RBF(Q, 1., 1./_np.random.uniform(0,1,(Q,)), ARD=True)# + GPy.kern.Bias(Q, _np.exp(-2)) Y = data['X'][:N] m = GPy.models.BayesianGPLVM(Y, Q, kernel=kernel, num_inducing=num_inducing, **k) m.data_labels = data['Y'][:N].argmax(axis=1) @@ -270,7 +270,7 @@ def bgplvm_simulation(optimize=True, verbose=1, from GPy import kern from GPy.models import BayesianGPLVM - D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 3, 10 + D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 5, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) @@ -293,7 +293,7 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, from GPy.models import BayesianGPLVM from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData - D1, D2, D3, N, num_inducing, Q = 15, 5, 8, 30, 5, 9 + D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 5, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 349cd72d..cc630bfa 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -45,7 +45,7 @@ class VarDTC(object): def inference(self, kern, X, X_variance, Z, likelihood, Y): """Inference for normal sparseGP""" uncertain_inputs = False - psi0, psi1, psi2 = _compute_psi(kern, X, X_variance, Z, uncertain_inputs) + psi0, psi1, psi2 = _compute_psi(kern, X, Z) return self._inference(kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs) def inference_latent(self, kern, posterior_variational, Z, likelihood, Y): @@ -205,7 +205,7 @@ class VarDTCMissingData(object): def inference(self, kern, X, X_variance, Z, likelihood, Y): """Inference for normal sparseGP""" uncertain_inputs = False - psi0, psi1, psi2 = _compute_psi(kern, X, X_variance, Z, uncertain_inputs) + psi0, psi1, psi2 = _compute_psi(kern, X, Z) return self._inference(kern, psi0, psi1, psi2, Z, likelihood, Y, uncertain_inputs) def inference_latent(self, kern, posterior_variational, Z, likelihood, Y): @@ -358,7 +358,7 @@ class VarDTCMissingData(object): return post, log_marginal, grad_dict -def _compute_psi(kern, X, X_variance, Z): +def _compute_psi(kern, X, Z): psi0 = kern.Kdiag(X) psi1 = kern.K(X, Z) psi2 = None diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 74b8abe0..8455c7a1 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -41,7 +41,7 @@ class BayesianGPLVM(SparseGP): if likelihood is None: likelihood = Gaussian() - self.q = NormalPosterior(X, X_variance) + #self.q = NormalPosterior(X, X_variance) self.variational_prior = NormalPrior() SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, X_variance, name, **kwargs) From e9957ed8961e7a62df8d8d6f022c388e64720760 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 08:23:46 +0000 Subject: [PATCH 05/51] docstrings in kern.py --- GPy/kern/_src/kern.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 2e412688..d4d54bbb 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -26,27 +26,41 @@ class Kern(Parameterized): raise NotImplementedError def Kdiag(self, Xa): raise NotImplementedError - def psi0(self,Z,variational_posterior): + def psi0(self, Z, variational_posterior): raise NotImplementedError - def psi1(self,Z,variational_posterior): + def psi1(self, Z, variational_posterior): raise NotImplementedError - def psi2(self,Z,variational_posterior): + def psi2(self, Z, variational_posterior): raise NotImplementedError def gradients_X(self, dL_dK, X, X2): raise NotImplementedError def gradients_X_diag(self, dL_dK, X): raise NotImplementedError + def update_gradients_full(self, dL_dK, X, X2): """Set the gradients of all parameters when doing full (N) inference.""" raise NotImplementedError + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): """ Set the gradients of all parameters when doing inference with uncertain inputs, using expectations of the kernel. + + The esential maths is + + dL_d{theta_i} = dL_dpsi0 * dpsi0_d{theta_i} + + dL_dpsi1 * dpsi1_d{theta_i} + + dL_dpsi2 * dpsi2_d{theta_i} """ raise NotImplementedError + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + """ + Returns the derivative of the objective wrt Z, using the chain rule + through the expectation variables. + """ raise NotImplementedError + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): """ Compute the gradients wrt the parameters of the variational @@ -106,7 +120,8 @@ class Kern(Parameterized): def prod(self, other, tensor=False): """ - Multiply two kernels (either on the same space, or on the tensor product of the input space). + Multiply two kernels (either on the same space, or on the tensor + product of the input space). :param other: the other kernel to be added :type other: GPy.kern From ec6e85b1c12cfaa42a91bb683d574f04d098d6c9 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 08:24:03 +0000 Subject: [PATCH 06/51] mucho changes to linear.py --- GPy/core/sparse_gp.py | 2 +- GPy/examples/dimensionality_reduction.py | 7 ++-- GPy/kern/_src/linear.py | 53 ++++++++++++------------ GPy/models/bayesian_gplvm.py | 2 +- 4 files changed, 32 insertions(+), 32 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 7677fea2..91d1f205 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -70,7 +70,7 @@ class SparseGP(GP): #gradients wrt Z self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) self.Z.gradient += self.kern.gradients_Z_expectations( - self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpis2'], Z=self.Z, variational_posterior=self.X) + self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) else: #gradients wrt kernel target = np.zeros(self.kern.size) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index a2686d73..eab9a965 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -274,6 +274,7 @@ def bgplvm_simulation(optimize=True, verbose=1, _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) + #k = kern.RBF(Q, ARD=True, lengthscale=10.) m = BayesianGPLVM(Y, Q, init="PCA", num_inducing=num_inducing, kernel=k) if optimize: @@ -281,7 +282,7 @@ def bgplvm_simulation(optimize=True, verbose=1, m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05) if plot: - m.q.plot("BGPLVM Latent Space 1D") + m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD('BGPLVM Simulation ARD Parameters') return m @@ -302,7 +303,7 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, m = BayesianGPLVM(Y.copy(), Q, init="random", num_inducing=num_inducing, kernel=k) m.inference_method = VarDTCMissingData() m.Y[inan] = _np.nan - m.q.variance *= .1 + m.X.variance *= .1 m.parameters_changed() m.Yreal = Y @@ -311,7 +312,7 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05) if plot: - m.q.plot("BGPLVM Latent Space 1D") + m.X.plot("BGPLVM Latent Space 1D") m.kern.plot_ARD('BGPLVM Simulation ARD Parameters') return m diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py index e503180a..44f71d85 100644 --- a/GPy/kern/_src/linear.py +++ b/GPy/kern/_src/linear.py @@ -22,22 +22,25 @@ class Linear(Kern): :param input_dim: the number of input dimensions :type input_dim: int :param variances: the vector of variances :math:`\sigma^2_i` - :type variances: array or list of the appropriate size (or float if there is only one variance parameter) - :param ARD: Auto Relevance Determination. If equal to "False", the kernel has only one variance parameter \sigma^2, otherwise there is one variance parameter per dimension. + :type variances: array or list of the appropriate size (or float if there + is only one variance parameter) + :param ARD: Auto Relevance Determination. If False, the kernel has only one + variance parameter \sigma^2, otherwise there is one variance + parameter per dimension. :type ARD: Boolean :rtype: kernel object + """ def __init__(self, input_dim, variances=None, ARD=False, name='linear'): super(Linear, self).__init__(input_dim, name) self.ARD = ARD - if ARD == False: + if not ARD: if variances is not None: variances = np.asarray(variances) assert variances.size == 1, "Only one variance needed for non-ARD kernel" else: variances = np.ones(1) - self._Xcache, self._X2cache = np.empty(shape=(2,)) else: if variances is not None: variances = np.asarray(variances) @@ -103,7 +106,6 @@ class Linear(Kern): #---------------------------------------# # PSI statistics # - # variational # #---------------------------------------# def psi0(self, Z, variational_posterior): @@ -118,18 +120,16 @@ class Linear(Kern): return np.dot(ZAinner, ZA.T) def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - mu, S = variational_posterior.mean, variational_posterior.variance + #psi1 + self.update_gradients_full(dL_dpsi1, variational_posterior.mean, Z) # psi0: tmp = dL_dpsi0[:, None] * self._mu2S(variational_posterior) - if self.ARD: grad = tmp.sum(0) - else: grad = np.atleast_1d(tmp.sum()) - #psi1 - self.update_gradients_full(dL_dpsi1, mu, Z) - grad += self.variances.gradient + if self.ARD: self.variances.gradient += tmp.sum(0) + else: self.variances.gradient += tmp.sum() #psi2 tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(variational_posterior, Z)[:, :, None, :] * (2. * Z)[None, None, :, :]) - if self.ARD: grad += tmp.sum(0).sum(0).sum(0) - else: grad += tmp.sum() + if self.ARD: self.variances.gradient += tmp.sum(0).sum(0).sum(0) + else: self.variances.gradient += tmp.sum() def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): #psi1 @@ -155,7 +155,7 @@ class Linear(Kern): #--------------------------------------------------# - def _weave_dpsi2_dmuS(self, dL_dpsi2, Z, pv, target_mu, target_S): + def _weave_dpsi2_dmuS(self, dL_dpsi2, Z, vp, target_mu, target_S): # Think N,num_inducing,num_inducing,input_dim ZA = Z * self.variances AZZA = ZA.T[:, None, :, None] * ZA[None, :, None, :] @@ -198,16 +198,15 @@ class Linear(Kern): weave_options = {'headers' : [''], 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], 'extra_link_args' : ['-lgomp']} - - mu = pv.mean + mu = vp.mean N,num_inducing,input_dim,mu = mu.shape[0],Z.shape[0],mu.shape[1],param_to_array(mu) weave.inline(code, support_code=support_code, libraries=['gomp'], arg_names=['N','num_inducing','input_dim','mu','AZZA','AZZA_2','target_mu','target_S','dL_dpsi2'], type_converters=weave.converters.blitz,**weave_options) - def _weave_dpsi2_dZ(self, dL_dpsi2, Z, pv, target): - AZA = self.variances*self._ZAinner(pv, Z) + def _weave_dpsi2_dZ(self, dL_dpsi2, Z, vp, target): + AZA = self.variances*self._ZAinner(vp, Z) code=""" int n,m,mm,q; #pragma omp parallel for private(n,mm,q) @@ -229,23 +228,23 @@ class Linear(Kern): 'extra_compile_args': ['-fopenmp -O3'], #-march=native'], 'extra_link_args' : ['-lgomp']} - N,num_inducing,input_dim = pv.mean.shape[0],Z.shape[0],pv.mean.shape[1] - mu = param_to_array(pv.mean) + N,num_inducing,input_dim = vp.mean.shape[0],Z.shape[0],vp.mean.shape[1] + mu = param_to_array(vp.mean) weave.inline(code, support_code=support_code, libraries=['gomp'], arg_names=['N','num_inducing','input_dim','AZA','target','dL_dpsi2'], type_converters=weave.converters.blitz,**weave_options) - def _mu2S(self, pv): - return np.square(pv.mean) + pv.variance + def _mu2S(self, vp): + return np.square(vp.mean) + vp.variance - def _ZAinner(self, pv, Z): + def _ZAinner(self, vp, Z): ZA = Z*self.variances - inner = (pv.mean[:, None, :] * pv.mean[:, :, None]) - diag_indices = np.diag_indices(pv.mean.shape[1], 2) - inner[:, diag_indices[0], diag_indices[1]] += pv.variance + inner = (vp.mean[:, None, :] * vp.mean[:, :, None]) + diag_indices = np.diag_indices(vp.mean.shape[1], 2) + inner[:, diag_indices[0], diag_indices[1]] += vp.variance - return np.dot(ZA, inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [num_inducing x N x input_dim]! + return np.dot(ZA, inner).swapaxes(0, 1) # NOTE: self.ZAinner \in [num_inducing x num_data x input_dim]! def input_sensitivity(self): if self.ARD: return self.variances diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 366995dc..5fed767b 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -66,7 +66,7 @@ class BayesianGPLVM(SparseGP): 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_q_variational(variational_posterior=self.X, Z=self.Z, **self.grad_dict) + self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **self.grad_dict) # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) From a9a4841e5f3ffea628821b27fadb988b961923fe Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 08:24:12 +0000 Subject: [PATCH 07/51] logexpNeg transformation --- GPy/core/parameterization/transformations.py | 32 ++++++++++++++++++-- 1 file changed, 30 insertions(+), 2 deletions(-) diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index 51194865..36291ca3 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -4,10 +4,10 @@ import numpy as np from domains import _POSITIVE,_NEGATIVE, _BOUNDED -import sys import weakref -_lim_val = -np.log(sys.float_info.epsilon) +_exp_lim_val = np.finfo(np.float64).max +_lim_val = np.log(_exp_lim_val)#-np.log(sys.float_info.epsilon) #=============================================================================== # Fixing constants @@ -34,6 +34,16 @@ class Transformation(object): def initialize(self, f): """ produce a sensible initial value for f(x)""" raise NotImplementedError + def plot(self, xlabel=r'transformed $\theta$', ylabel=r'$\theta$', axes=None, *args,**kw): + import sys + assert "matplotlib" in sys.modules, "matplotlib package has not been imported." + import matplotlib.pyplot as plt + from ...plotting.matplot_dep import base_plots + x = np.linspace(-8,8) + base_plots.meanplot(x, self.f(x),axes=axes*args,**kw) + axes = plt.gca() + axes.set_xlabel(xlabel) + axes.set_ylabel(ylabel) def __str__(self): raise NotImplementedError def __repr__(self): @@ -54,6 +64,24 @@ class Logexp(Transformation): return np.abs(f) def __str__(self): return '+ve' + + +class LogexpNeg(Transformation): + domain = _POSITIVE + def f(self, x): + return np.where(x>_lim_val, -x, -np.log(1. + np.exp(np.clip(x, -np.inf, _lim_val)))) + #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x))) + def finv(self, f): + return np.where(f>_lim_val, 0, np.log(np.exp(-f) - 1.)) + def gradfactor(self, f): + return np.where(f>_lim_val, -1, -1 + np.exp(-f)) + def initialize(self, f): + if np.any(f < 0.): + print "Warning: changing parameters to satisfy constraints" + return np.abs(f) + def __str__(self): + return '+ve' + class NegativeLogexp(Transformation): domain = _NEGATIVE From 102aabc6e50ac7ef04f80beeaa0df41414c6a42c Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 08:25:50 +0000 Subject: [PATCH 08/51] parameter inheritance structure --- GPy/core/parameterization/param.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index ccbc76d5..8eac69da 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -15,7 +15,7 @@ __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision __print_threshold__ = 5 ###### -class Param(Constrainable, ObservableArray, Gradcheckable, Indexable): +class Param(Constrainable, ObservableArray, Gradcheckable): """ Parameter object for GPy models. From 3c0f89bf53ea681cf8dcaf00df8d4c438fb32b63 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 08:33:45 +0000 Subject: [PATCH 09/51] vdtc_missing data corrections --- GPy/core/sparse_gp.py | 3 +-- GPy/inference/latent_function_inference/var_dtc.py | 12 ++++++------ 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index b1a3c12e..15a4e1f8 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -2,12 +2,11 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np -from ..util.linalg import mdot from gp import GP from parameterization.param import Param from ..inference.latent_function_inference import var_dtc from .. import likelihoods -from parameterization.variational import NormalPosterior +from parameterization.variational import VariationalPosterior class SparseGP(GP): """ diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 66ab3cbe..fec61204 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -204,14 +204,14 @@ class VarDTCMissingData(object): def inference(self, kern, X, Z, likelihood, Y): if isinstance(X, VariationalPosterior): uncertain_inputs = True - psi0 = kern.psi0(Z, X) - psi1 = kern.psi1(Z, X) - psi2 = kern.psi2(Z, X) + psi0_all = kern.psi0(Z, X) + psi1_all = kern.psi1(Z, X) + psi2_all = kern.psi2(Z, X) else: uncertain_inputs = False - psi0 = kern.Kdiag(X) - psi1 = kern.K(X, Z) - psi2 = None + psi0_all = kern.Kdiag(X) + psi1_all = kern.K(X, Z) + psi2_all = None Ys, traces = self._Y(Y) beta_all = 1./likelihood.variance From ed0c162597b343d037b24e4b26687642a766ee69 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 09:10:11 +0000 Subject: [PATCH 10/51] BayesianGPLVM init with paramschanged --- GPy/models/bayesian_gplvm.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 5fed767b..18a08e5d 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -49,6 +49,7 @@ class BayesianGPLVM(SparseGP): SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, name, **kwargs) self.add_parameter(self.X, index=0) + self.parameters_changed() def _getstate(self): """ From f6da5345d91f9ebba4da7e7aad7a3eff5eeb3b48 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 09:16:31 +0000 Subject: [PATCH 11/51] plot latent updated --- GPy/core/model.py | 209 ++---------------- GPy/kern/_src/rbf.py | 2 +- .../matplot_dep/dim_reduction_plots.py | 24 +- 3 files changed, 32 insertions(+), 203 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 21bcf0c7..6514d73a 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -20,14 +20,13 @@ class Model(Parameterized): self.optimization_runs = [] self.sampling_runs = [] self.preferred_optimizer = 'scg' - # self._set_params(self._get_params()) has been taken out as it should only be called on leaf nodes + def log_likelihood(self): raise NotImplementedError, "this needs to be implemented to use the model class" + def _log_likelihood_gradients(self): - # def dK_d(self, param, dL_dK, X, X2) g = np.zeros(self.size) try: - # [g.__setitem__(s, self.gradient_mapping[p]().flat) for p, s in itertools.izip(self._parameters_, self._param_slices_) if not p.is_fixed] [p._collect_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_) if not p.is_fixed] except ValueError: raise ValueError, 'Gradient for {} not defined, please specify gradients for parameters to optimize'.format(p.name) @@ -61,85 +60,6 @@ class Model(Parameterized): self.priors = state.pop() Parameterized._setstate(self, state) -# def set_prior(self, regexp, what): -# """ -# -# Sets priors on the model parameters. -# -# **Notes** -# -# Asserts that the prior is suitable for the constraint. If the -# wrong constraint is in place, an error is raised. If no -# constraint is in place, one is added (warning printed). -# -# For tied parameters, the prior will only be "counted" once, thus -# a prior object is only inserted on the first tied index -# -# :param regexp: regular expression of parameters on which priors need to be set. -# :type param: string, regexp, or integer array -# :param what: prior to set on parameter. -# :type what: GPy.core.Prior type -# -# """ -# if self.priors is None: -# self.priors = [None for i in range(self._get_params().size)] -# -# which = self.grep_param_names(regexp) -# -# # check tied situation -# tie_partial_matches = [tie for tie in self.tied_indices if (not set(tie).isdisjoint(set(which))) & (not set(tie) == set(which))] -# if len(tie_partial_matches): -# raise ValueError, "cannot place prior across partial ties" -# tie_matches = [tie for tie in self.tied_indices if set(which) == set(tie) ] -# if len(tie_matches) > 1: -# raise ValueError, "cannot place prior across multiple ties" -# elif len(tie_matches) == 1: -# which = which[:1] # just place a prior object on the first parameter -# -# -# # check constraints are okay -# -# if what.domain is _POSITIVE: -# constrained_positive_indices = [i for i, t in zip(self.constrained_indices, self.constraints) if t.domain is _POSITIVE] -# if len(constrained_positive_indices): -# constrained_positive_indices = np.hstack(constrained_positive_indices) -# else: -# constrained_positive_indices = np.zeros(shape=(0,)) -# bad_constraints = np.setdiff1d(self.all_constrained_indices(), constrained_positive_indices) -# assert not np.any(which[:, None] == bad_constraints), "constraint and prior incompatible" -# unconst = np.setdiff1d(which, constrained_positive_indices) -# if len(unconst): -# print "Warning: constraining parameters to be positive:" -# print '\n'.join([n for i, n in enumerate(self._get_param_names()) if i in unconst]) -# print '\n' -# self.constrain_positive(unconst) -# elif what.domain is _REAL: -# assert not np.any(which[:, None] == self.all_constrained_indices()), "constraint and prior incompatible" -# else: -# raise ValueError, "prior not recognised" -# -# # store the prior in a local list -# for w in which: -# self.priors[w] = what - - def get_gradient(self, name, return_names=False): - """ - Get model gradient(s) by name. The name is applied as a regular expression and all parameters that match that regular expression are returned. - - :param name: the name of parameters required (as a regular expression). - :type name: regular expression - :param return_names: whether or not to return the names matched (default False) - :type return_names: bool - """ - matches = self.grep_param_names(name) - if len(matches): - if return_names: - return self._log_likelihood_gradients()[matches], np.asarray(self._get_param_names())[matches].tolist() - else: - return self._log_likelihood_gradients()[matches] - else: - raise AttributeError, "no parameter matches %s" % name - def randomize(self): """ Randomize the model. @@ -183,7 +103,9 @@ class Model(Parameterized): :param messages: whether to display during optimisation :type messages: bool - .. note:: If num_processes is None, the number of workes in the multiprocessing pool is automatically set to the number of processors on the current machine. + .. note:: If num_processes is None, the number of workes in the + multiprocessing pool is automatically set to the number of processors + on the current machine. """ initial_parameters = self._get_params_transformed() @@ -227,45 +149,23 @@ class Model(Parameterized): self._set_params_transformed(initial_parameters) def ensure_default_constraints(self, warning=True): - """ + """ Ensure that any variables which should clearly be positive have been constrained somehow. The method performs a regular expression search on parameter names looking for the terms 'variance', 'lengthscale', 'precision' and 'kappa'. If any of these terms are present in the name the parameter is constrained positive. + + DEPRECATED. """ raise DeprecationWarning, 'parameters now have default constraints' - #positive_strings = ['variance', 'lengthscale', 'precision', 'kappa', 'sensitivity'] - # param_names = self._get_param_names() - -# for s in positive_strings: -# paramlist = self.grep_param_names(".*"+s) -# for param in paramlist: -# for p in param.flattened_parameters: -# rav_i = set(self._raveled_index_for(p)) -# for constraint in self.constraints.iter_properties(): -# rav_i -= set(self._constraint_indices(p, constraint)) -# rav_i -= set(np.nonzero(self._fixes_for(p)!=UNFIXED)[0]) -# ind = self._backtranslate_index(p, np.array(list(rav_i), dtype=int)) -# if ind.size != 0: -# p[np.unravel_index(ind, p.shape)].constrain_positive(warning=warning) -# if paramlist: -# self.__getitem__(None, paramlist).constrain_positive(warning=warning) -# 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 currently_constrained): -# to_make_positive.append(i) -# if len(to_make_positive): -# self.constrain_positive(np.asarray(to_make_positive), warning=warning) def objective_function(self, x): """ The objective function passed to the optimizer. It combines the likelihood and the priors. - + Failures are handled robustly. The algorithm will try several times to return the objective, and will raise the original exception if the objective cannot be computed. @@ -336,7 +236,9 @@ class Model(Parameterized): :messages: whether to display during optimisation :type messages: bool :param optimizer: which optimizer to use (defaults to self.preferred optimizer) - :type optimizer: string TODO: valid strings? + :type optimizer: string + + TODO: valid args """ if optimizer is None: optimizer = self.preferred_optimizer @@ -389,15 +291,15 @@ class Model(Parameterized): indices = np.r_[:self.size] which = (transformed_index[:,None]==indices[self._fixes_][None,:]).nonzero() transformed_index = (indices-(~self._fixes_).cumsum())[transformed_index[which[0]]] - + if transformed_index.size == 0: print "No free parameters to check" return - + # just check the global ratio dx = np.zeros_like(x) dx[transformed_index] = step * np.sign(np.random.uniform(-1, 1, transformed_index.size)) - + # evaulate around the point x f1 = self.objective_function(x + dx) f2 = self.objective_function(x - dx) @@ -405,7 +307,7 @@ class Model(Parameterized): dx = dx[transformed_index] gradient = gradient[transformed_index] - + numerical_gradient = (f1 - f2) / (2 * dx) global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gradient - numerical_gradient).mean() < tolerance) @@ -439,7 +341,7 @@ class Model(Parameterized): #print param_index, transformed_index else: transformed_index = param_index - + if param_index.size == 0: print "No free parameters to check" return @@ -472,82 +374,5 @@ class Model(Parameterized): print grad_string self._set_params_transformed(x) return ret - - def input_sensitivity(self): - """ - return an array describing the sesitivity of the model to each input - NB. Right now, we're basing this on the lengthscales (or - variances) of the kernel. TODO: proper sensitivity analysis - where we integrate across the model inputs and evaluate the - effect on the variance of the model output. """ - if not hasattr(self, 'kern'): - raise ValueError, "this model has no kernel" - - k = self.kern#[p for p in self.kern._parameters_ if hasattr(p, "ARD") and p.ARD] - from ..kern import RBF, Linear#, RBFInv - - if isinstance(k, RBF): - return 1. / k.lengthscale - #elif isinstance(k, RBFInv): - # return k.inv_lengthscale - elif isinstance(k, Linear): - return k.variances - else: - raise ValueError, "cannot determine sensitivity for this kernel" - - def pseudo_EM(self, stop_crit=.1, **kwargs): - """ - EM - like algorithm for Expectation Propagation and Laplace approximation - - :param stop_crit: convergence criterion - :type stop_crit: float - - .. Note: kwargs are passed to update_likelihood and optimize functions. - """ - assert isinstance(self.likelihood, likelihoods.EP) or isinstance(self.likelihood, likelihoods.EP_Mixed_Noise), "pseudo_EM is only available for EP likelihoods" - ll_change = stop_crit + 1. - iteration = 0 - last_ll = -np.inf - - convergence = False - alpha = 0 - stop = False - - # Handle **kwargs - ep_args = {} - for arg in kwargs.keys(): - if arg in ('epsilon', 'power_ep'): - ep_args[arg] = kwargs[arg] - del kwargs[arg] - - while not stop: - last_approximation = self.likelihood.copy() - last_params = self._get_params() - if len(ep_args) == 2: - self.update_likelihood_approximation(epsilon=ep_args['epsilon'], power_ep=ep_args['power_ep']) - elif len(ep_args) == 1: - if ep_args.keys()[0] == 'epsilon': - self.update_likelihood_approximation(epsilon=ep_args['epsilon']) - elif ep_args.keys()[0] == 'power_ep': - self.update_likelihood_approximation(power_ep=ep_args['power_ep']) - else: - self.update_likelihood_approximation() - new_ll = self.log_likelihood() - ll_change = new_ll - last_ll - - if ll_change < 0: - self.likelihood = last_approximation # restore previous likelihood approximation - self._set_params(last_params) # restore model parameters - print "Log-likelihood decrement: %s \nLast likelihood update discarded." % ll_change - stop = True - else: - self.optimize(**kwargs) - last_ll = self.log_likelihood() - if ll_change < stop_crit: - stop = True - iteration += 1 - if stop: - print "%s iterations." % iteration - self.update_likelihood_approximation() diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 7c43b18d..f98fa43e 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -225,7 +225,7 @@ class RBF(Stationary): if self.ARD: lengthscale2 = self.lengthscale **2 else: - lengthscale2 = np.ones(input_dim) * self.lengthscale2**2 + lengthscale2 = np.ones(input_dim) * self.lengthscale**2 code = """ double tmp; diff --git a/GPy/plotting/matplot_dep/dim_reduction_plots.py b/GPy/plotting/matplot_dep/dim_reduction_plots.py index 3f4ea9b0..10b352d3 100644 --- a/GPy/plotting/matplot_dep/dim_reduction_plots.py +++ b/GPy/plotting/matplot_dep/dim_reduction_plots.py @@ -2,6 +2,7 @@ import pylab as pb import numpy as np from latent_space_visualizations.controllers.imshow_controller import ImshowController,ImAnnotateController from ...util.misc import param_to_array +from ...core.parameterization.variational import VariationalPosterior from .base_plots import x_frame2D import itertools import Tango @@ -19,7 +20,7 @@ def most_significant_input_dimensions(model, which_indices): input_1, input_2 = 0, 1 else: try: - input_1, input_2 = np.argsort(model.input_sensitivity())[::-1][:2] + input_1, input_2 = np.argsort(model.kern.input_sensitivity())[::-1][:2] except: raise ValueError, "cannot automatically determine which dimensions to plot, please pass 'which_indices'" else: @@ -43,26 +44,29 @@ def plot_latent(model, labels=None, which_indices=None, labels = np.ones(model.num_data) input_1, input_2 = most_significant_input_dimensions(model, which_indices) - X = param_to_array(model.X) - # first, plot the output variance as a function of the latent space - Xtest, xx, yy, xmin, xmax = x_frame2D(X[:, [input_1, input_2]], resolution=resolution) - Xtest_full = np.zeros((Xtest.shape[0], model.X.shape[1])) + #fethch the data points X that we'd like to plot + X = model.X + if isinstance(X, VariationalPosterior): + X = param_to_array(X.mean) + else: + X = param_to_array(X) + + # create a function which computes the shading of latent space according to the output variance def plot_function(x): + Xtest_full = np.zeros((x.shape[0], model.X.shape[1])) Xtest_full[:, [input_1, input_2]] = x mu, var, low, up = model.predict(Xtest_full) var = var[:, :1] return np.log(var) + #Create an IMshow controller that can re-plot the latent space shading at a good resolution view = ImshowController(ax, plot_function, tuple(X[:, [input_1, input_2]].min(0)) + tuple(X[:, [input_1, input_2]].max(0)), resolution, aspect=aspect, interpolation='bilinear', cmap=pb.cm.binary) -# ax.imshow(var.reshape(resolution, resolution).T, -# extent=[xmin[0], xmax[0], xmin[1], xmax[1]], cmap=pb.cm.binary, interpolation='bilinear', origin='lower') - # make sure labels are in order of input: ulabels = [] for lab in labels: @@ -95,8 +99,8 @@ def plot_latent(model, labels=None, which_indices=None, if not np.all(labels == 1.) and legend: ax.legend(loc=0, numpoints=1) - ax.set_xlim(xmin[0], xmax[0]) - ax.set_ylim(xmin[1], xmax[1]) + #ax.set_xlim(xmin[0], xmax[0]) + #ax.set_ylim(xmin[1], xmax[1]) ax.grid(b=False) # remove the grid if present, it doesn't look good ax.set_aspect('auto') # set a nice aspect ratio From 60e8f88b9bbefc02c953e73d31e621d95dbcae82 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 10:38:09 +0000 Subject: [PATCH 12/51] lots of hacking on RBF --- GPy/examples/dimensionality_reduction.py | 8 +- GPy/kern/_src/rbf.py | 170 ++++++++++------------- 2 files changed, 80 insertions(+), 98 deletions(-) diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index fc0bc085..2044f08d 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -16,16 +16,16 @@ def bgplvm_test_model(optimize=False, verbose=1, plot=False, output_dim=200, nan output_dim = 1 input_dim = 3 else: - input_dim = 1 + input_dim = 2 output_dim = output_dim # generate GPLVM-like data X = _np.random.rand(num_inputs, input_dim) - #lengthscales = _np.random.rand(input_dim) - #k = (GPy.kern.RBF(input_dim, .5, lengthscales, ARD=True) + lengthscales = _np.random.rand(input_dim) + k = GPy.kern.RBF(input_dim, .5, lengthscales, ARD=True) ##+ GPy.kern.white(input_dim, 0.01) #) - k = GPy.kern.Linear(input_dim, ARD=1)# + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001) + #k = GPy.kern.Linear(input_dim, ARD=1)# + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim, 0.00001) K = k.K(X) Y = _np.random.multivariate_normal(_np.zeros(num_inputs), K, (output_dim,)).T diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index f98fa43e..3bf355be 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -35,87 +35,80 @@ class RBF(Stationary): # PSI statistics # #---------------------------------------# - def parameters_changed(self): - # reset cached results - self._Z, self._mu, self._S = np.empty(shape=(3, 1)) # cached versions of Z,mu,S - - def psi0(self, Z, variational_posterior): return self.Kdiag(variational_posterior.mean) def psi1(self, Z, variational_posterior): - mu = variational_posterior.mean - S = variational_posterior.variance - self._psi_computations(Z, mu, S) - return self._psi1 + _, _, _, psi1 = self._psi1computations(Z, variational_posterior) + return psi1 def psi2(self, Z, variational_posterior): - mu = variational_posterior.mean - S = variational_posterior.variance - self._psi_computations(Z, mu, S) - return self._psi2 + _, _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) + return psi2 def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - mu = variational_posterior.mean - S = variational_posterior.variance - self._psi_computations(Z, mu, S) l2 = self.lengthscale **2 #contributions from psi0: - self.variance.gradient += np.sum(dL_dpsi0) + self.variance.gradient = np.sum(dL_dpsi0) + self.lengthscale.gradient = 0. #from psi1 - self.variance.gradient += np.sum(dL_dpsi1 * self._psi1 / self.variance) - d_length = self._psi1[:,:,None] * ((self._psi1_dist_sq - 1.)/(self.lengthscale*self._psi1_denom) +1./self.lengthscale) + denom, _, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + d_length = psi1[:,:,None] * ((dist_sq - 1.)/(self.lengthscale*denom) +1./self.lengthscale) dpsi1_dlength = d_length * dL_dpsi1[:, :, None] if not self.ARD: self.lengthscale.gradient += dpsi1_dlength.sum() else: self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) + self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance #from psi2 - d_var = 2.*self._psi2 / self.variance - d_length = 2.*self._psi2[:, :, :, None] * (self._psi2_Zdist_sq * self._psi2_denom + self._psi2_mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * self._psi2_denom) + S = variational_posterior.variance + denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + d_length = 2.*psi2[:, :, :, None] * (Zdist_sq[None, :,:,:] * denom[:,None,None,:] + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom[:,None,None,:]) + #TODO: combine denom and l2 as denom_l2?? + #TODO: tidy the above! + #TODO: tensordot below? - self.variance.gradient += np.sum(dL_dpsi2 * d_var) dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] if not self.ARD: self.lengthscale.gradient += dpsi2_dlength.sum() else: self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) + self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - mu = variational_posterior.mean - S = variational_posterior.variance - self._psi_computations(Z, mu, S) l2 = self.lengthscale **2 #psi1 - denominator = (l2 * (self._psi1_denom)) - dpsi1_dZ = -self._psi1[:, :, None] * ((self._psi1_dist / denominator)) + denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + denominator = l2 * denom + dpsi1_dZ = -psi1[:, :, None] * (dist / denominator) grad = np.sum(dL_dpsi1[:, :, None] * dpsi1_dZ, 0) #psi2 - term1 = self._psi2_Zdist / l2 # num_inducing, num_inducing, input_dim - term2 = self._psi2_mudist / self._psi2_denom / l2 # N, num_inducing, num_inducing, input_dim - dZ = self._psi2[:, :, :, None] * (term1[None] + term2) + denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + term1 = Zdist / l2 # M, M, Q + term2 = mudist / denom[:,None,None,:] / l2 # N, M, M, Q + dZ = psi2[:, :, :, None] * (term1[None, :, :, :] + term2) #N,M,M,Q grad += 2*(dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0) return grad def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): - mu = variational_posterior.mean - S = variational_posterior.variance - self._psi_computations(Z, mu, S) l2 = self.lengthscale **2 #psi1 - tmp = self._psi1[:, :, None] / l2 / self._psi1_denom - grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * self._psi1_dist, 1) - grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (self._psi1_dist_sq - 1), 1) + denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + tmp = psi1[:, :, None] / l2 / denom + grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) + grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) #psi2 - tmp = self._psi2[:, :, :, None] / l2 / self._psi2_denom - grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * self._psi2_mudist).sum(1).sum(1) - grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*self._psi2_mudist_sq - 1)).sum(1).sum(1) + denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + tmp = psi2[:, :, :, None] / l2 / denom[:,None,None,:] + grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) + grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) return grad_mu, grad_S @@ -177,83 +170,72 @@ class RBF(Stationary): return target + #@cache_this TODO + def _psi1computations(self, Z, vp): + mu, S = vp.mean, vp.variance + l2 = self.lengthscale **2 + denom = S[:, None, :] / l2 + 1. # N,1,Q + dist = Z[None, :, :] - mu[:, None, :] # N,M,Q + dist_sq = np.square(dist) / l2 / denom # N,M,Q + exponent = -0.5 * np.sum(dist_sq + np.log(denom), -1)#N,M + psi1 = self.variance * np.exp(exponent) # N,M + return denom, dist, dist_sq, psi1 - def _psi_computations(self, Z, mu, S): - # here are the "statistics" for psi1 and psi2 - Z_changed = not fast_array_equal(Z, self._Z) - if Z_changed: - # Z has changed, compute Z specific stuff - self._psi2_Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q - self._psi2_Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q - self._psi2_Zdist_sq = np.square(self._psi2_Zdist / self.lengthscale) # M,M,Q - if Z_changed or not fast_array_equal(mu, self._mu) or not fast_array_equal(S, self._S): - # something's changed. recompute EVERYTHING - l2 = self.lengthscale **2 - # psi1 - self._psi1_denom = S[:, None, :] / l2 + 1. - self._psi1_dist = Z[None, :, :] - mu[:, None, :] - self._psi1_dist_sq = np.square(self._psi1_dist) / l2 / self._psi1_denom - self._psi1_exponent = -0.5 * np.sum(self._psi1_dist_sq + np.log(self._psi1_denom), -1) - self._psi1 = self.variance * np.exp(self._psi1_exponent) + #@cache_this TODO + def _psi2computations(self, Z, vp): + mu, S = vp.mean, vp.variance - # psi2 - self._psi2_denom = 2.*S[:, None, None, :] / l2 + 1. # N,M,M,Q - 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)/(l2*self._psi2_denom) - # self._psi2_exponent = np.sum(-self._psi2_Zdist_sq -self._psi2_mudist_sq -0.5*np.log(self._psi2_denom),-1) #N,M,M,Q - self._psi2 = np.square(self.variance) * np.exp(self._psi2_exponent) # N,M,M,Q + N, Q = mu.shape + M = Z.shape[0] - # store matrices for caching - self._Z, self._mu, self._S = Z, mu, S + #compute required distances + Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q + Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q + Zdist_sq = np.square(Zdist / self.lengthscale) # M,M,Q - def weave_psi2(self, mu, Zhat): - N, input_dim = mu.shape - num_inducing = Zhat.shape[0] + #allocate memory for the things we want to compute + mudist = np.empty((N, M, M, Q)) + mudist_sq = np.empty((N, M, M, Q)) + psi2 = np.empty((N, M, M)) - mudist = np.empty((N, num_inducing, num_inducing, input_dim)) - mudist_sq = np.empty((N, num_inducing, num_inducing, input_dim)) - psi2_exponent = np.zeros((N, num_inducing, num_inducing)) - psi2 = np.empty((N, num_inducing, num_inducing)) + l2 = self.lengthscale **2 + denom = 2.*S / l2 + 1. # N,Q + half_log_denom = 0.5 * np.log(denom) + denom_l2 = denom*l2 # TODO: Max and James: divide?? - psi2_Zdist_sq = self._psi2_Zdist_sq - _psi2_denom = self._psi2_denom.squeeze().reshape(N, self.input_dim) - half_log_psi2_denom = 0.5 * np.log(self._psi2_denom).squeeze().reshape(N, self.input_dim) - variance_sq = np.float64(np.square(self.variance)) - if self.ARD: - lengthscale2 = self.lengthscale **2 - else: - lengthscale2 = np.ones(input_dim) * self.lengthscale**2 + variance_sq = float(np.square(self.variance)) code = """ double tmp; + double exponent; - #pragma omp parallel for private(tmp) + #pragma omp parallel for private(tmp, exponentgg) for (int n=0; n Date: Wed, 26 Feb 2014 11:38:46 +0000 Subject: [PATCH 13/51] some work on ep, and some messing with wher ethe derivatives are computed (in the model, not the inference --- GPy/core/gp.py | 2 +- GPy/inference/latent_function_inference/ep.py | 15 +++++++++++ .../exact_gaussian_inference.py | 3 +-- GPy/util/multioutput.py | 25 +++++++++++++++++-- 4 files changed, 40 insertions(+), 5 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 81985773..50f9d342 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -70,7 +70,7 @@ class GP(Model): def parameters_changed(self): self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata) - self._dL_dK = grad_dict['dL_dK'] + self.kern.update_gradients_full(grad_dict['dL_dK']) def log_likelihood(self): return self._log_marginal_likelihood diff --git a/GPy/inference/latent_function_inference/ep.py b/GPy/inference/latent_function_inference/ep.py index 87c08221..1904d48c 100644 --- a/GPy/inference/latent_function_inference/ep.py +++ b/GPy/inference/latent_function_inference/ep.py @@ -23,11 +23,26 @@ class EP(object): self.old_mutilde, self.old_vtilde = None, None def inference(self, kern, X, likelihood, Y, Y_metadata=None): + num_data, output_dim = X.shape + assert output_dim ==1, "ep in 1D only (for now!)" K = kern.K(X) mu_tilde, tau_tilde = self.expectation_propagation() + Wi, LW, LWi, W_logdet = pdinv(K + np.diag(1./tau_tilde) + + alpha, _ = dpotrs(LW, mu_tilde, lower=1) + + log_marginal = 0.5*(-num_data * log_2_pi - W_logdet - np.sum(alpha * mu_tilde)) + + dL_dK = 0.5 * (tdot(alpha[:,None]) - Wi) + + #TODO: what abot derivatives of the likelihood parameters? + + return Posterior(woodbury_inv=Wi, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK} + + def expectation_propagation(self, K, Y, Y_metadata, likelihood) diff --git a/GPy/inference/latent_function_inference/exact_gaussian_inference.py b/GPy/inference/latent_function_inference/exact_gaussian_inference.py index 2c3bfa6a..922b52f4 100644 --- a/GPy/inference/latent_function_inference/exact_gaussian_inference.py +++ b/GPy/inference/latent_function_inference/exact_gaussian_inference.py @@ -49,8 +49,7 @@ class ExactGaussianInference(object): dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi) - kern.update_gradients_full(dL_dK, X) - + #TODO: does this really live here? likelihood.update_gradients(np.diag(dL_dK)) return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK} diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py index a57593a7..eb4d8d08 100644 --- a/GPy/util/multioutput.py +++ b/GPy/util/multioutput.py @@ -2,6 +2,27 @@ import numpy as np import warnings from .. import kern +def build_XY(input_list,output_list=None,index=None): + num_outputs = len(input_list) + _s = [0] + [ _x.shape[0] for _x in input_list ] + _s = np.cumsum(_s) + slices = [slice(a,b) for a,b in zip(_s[:-1],_s[1:])] + if output_list is not None: + assert num_outputs == len(output_list) + Y = np.vstack(output_list) + else: + Y = None + + if index is not None: + assert len(index) == num_outputs + I = np.vstack( [j*np.ones((_x.shape[0],1)) for _x,j in zip(input_list,index)] ) + else: + I = np.vstack( [j*np.ones((_x.shape[0],1)) for _x,j in zip(input_list,range(num_outputs))] ) + + X = np.vstack(input_list) + X = np.hstack([X,I]) + return X,Y,slices + def build_lcm(input_dim, num_outputs, CK = [], NC = [], W_columns=1,W=None,kappa=None): #TODO build_icm or build_lcm """ @@ -25,9 +46,9 @@ def build_lcm(input_dim, num_outputs, CK = [], NC = [], W_columns=1,W=None,kappa k.input_dim = input_dim + 1 warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.") - kernel = CK[0].prod(kern.coregionalize(num_outputs,W_columns,W,kappa),tensor=True) + kernel = CK[0].prod(kern.Coregionalize(num_outputs,W_columns,W,kappa),tensor=True) for k in CK[1:]: - k_coreg = kern.coregionalize(num_outputs,W_columns,W,kappa) + k_coreg = kern.Coregionalize(num_outputs,W_columns,W,kappa) kernel += k.prod(k_coreg,tensor=True) for k in NC: kernel += k From 95960c33c8d0957d734ab333aaf61ac0ebee6bf0 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 12:32:06 +0000 Subject: [PATCH 14/51] rbf with new parameter structure --- GPy/kern/_src/rbf.py | 77 +++++++++++++++++++++----------------------- 1 file changed, 36 insertions(+), 41 deletions(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 3bf355be..145bca5d 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -4,11 +4,7 @@ import numpy as np from scipy import weave -from kern import Kern -from ...util.linalg import tdot -from ...util.misc import fast_array_equal, param_to_array -from ...core.parameterization import Param -from ...core.parameterization.transformations import Logexp +from ...util.misc import param_to_array from stationary import Stationary class RBF(Stationary): @@ -66,7 +62,7 @@ class RBF(Stationary): #from psi2 S = variational_posterior.variance denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - d_length = 2.*psi2[:, :, :, None] * (Zdist_sq[None, :,:,:] * denom[:,None,None,:] + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom[:,None,None,:]) + d_length = 2.*psi2[:, :, :, None] * (Zdist_sq[None, :,:,:] * denom + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom) #TODO: combine denom and l2 as denom_l2?? #TODO: tidy the above! #TODO: tensordot below? @@ -91,7 +87,7 @@ class RBF(Stationary): #psi2 denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) term1 = Zdist / l2 # M, M, Q - term2 = mudist / denom[:,None,None,:] / l2 # N, M, M, Q + term2 = mudist / denom / l2 # N, M, M, Q dZ = psi2[:, :, :, None] * (term1[None, :, :, :] + term2) #N,M,M,Q grad += 2*(dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0) @@ -106,7 +102,7 @@ class RBF(Stationary): grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) #psi2 denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - tmp = psi2[:, :, :, None] / l2 / denom[:,None,None,:] + tmp = psi2[:, :, :, None] / l2 / denom grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) @@ -198,51 +194,50 @@ class RBF(Stationary): #allocate memory for the things we want to compute mudist = np.empty((N, M, M, Q)) mudist_sq = np.empty((N, M, M, Q)) + exponent = np.zeros((N,M,M)) psi2 = np.empty((N, M, M)) l2 = self.lengthscale **2 - denom = 2.*S / l2 + 1. # N,Q - half_log_denom = 0.5 * np.log(denom) - denom_l2 = denom*l2 # TODO: Max and James: divide?? - + denom = (2.*S[:,None,None,:] / l2) + 1. # N,Q + half_log_denom = 0.5 * np.log(denom[:,0,0,:]) + denom_l2 = denom[:,0,0,:]*l2 + variance_sq = float(np.square(self.variance)) code = """ - double tmp; - double exponent; + double tmp, exponent_tmp; - #pragma omp parallel for private(tmp, exponentgg) - for (int n=0; n #include From f6228abd4c25554357d515504311cc45687e37bc Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 13:18:29 +0000 Subject: [PATCH 15/51] gradients --- GPy/core/gp.py | 2 +- GPy/kern/_src/rbf.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 50f9d342..916aaa81 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -70,7 +70,7 @@ class GP(Model): def parameters_changed(self): self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata) - self.kern.update_gradients_full(grad_dict['dL_dK']) + self.kern.update_gradients_full(grad_dict['dL_dK'], self.X) def log_likelihood(self): return self._log_marginal_likelihood diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 145bca5d..666a79f4 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -61,8 +61,8 @@ class RBF(Stationary): #from psi2 S = variational_posterior.variance - denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - d_length = 2.*psi2[:, :, :, None] * (Zdist_sq[None, :,:,:] * denom + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom) + denom, _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * denom + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom) #TODO: combine denom and l2 as denom_l2?? #TODO: tidy the above! #TODO: tensordot below? From adce78ba5af6d9d2689fb2c0834270bcf031831d Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 13:18:37 +0000 Subject: [PATCH 16/51] maps import --- GPy/plotting/matplot_dep/maps.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/GPy/plotting/matplot_dep/maps.py b/GPy/plotting/matplot_dep/maps.py index 1c839f64..e941ab2d 100644 --- a/GPy/plotting/matplot_dep/maps.py +++ b/GPy/plotting/matplot_dep/maps.py @@ -4,7 +4,6 @@ import matplotlib.patches as patches from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection #from matplotlib import cm -import shapefile import re pb.ion() @@ -134,8 +133,10 @@ def plot_string_match(sf,regex,field,**kwargs): plot(shape_records,**kwargs) -def new_shape_string(sf,name,regex,field=2,type=shapefile.POINT): - +def new_shape_string(sf,name,regex,field=2,type=None): + import shapefile + if type is None: + type = shapefile.POINT newshp = shapefile.Writer(shapeType = sf.shapeType) newshp.autoBalance = 1 From 96d7bb3c8692469c131f4cd4150f8f6aaddeed5d Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 13:23:18 +0000 Subject: [PATCH 17/51] sparse GP no longer accepts X_variance --- GPy/core/sparse_gp.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 15a4e1f8..751a295d 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -31,7 +31,7 @@ class SparseGP(GP): """ - def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, X_variance=None, name='sparse gp'): + def __init__(self, X, Y, Z, kernel, likelihood, inference_method=None, name='sparse gp'): #pick a sensible inference method if inference_method is None: From 26aeb5e1dbb937367648e6c8ad339fe068bc5cff Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 13:37:15 +0000 Subject: [PATCH 18/51] lengthscale fixes --- GPy/kern/_src/stationary.py | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 2d0d284a..b3c648fc 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -73,16 +73,9 @@ class Stationary(Kern): """ if self.ARD: - if X2 is None: - Xl = X/self.lengthscale - Xsq = np.sum(np.square(Xl),1) - return np.sqrt(np.sqrt(-2.*tdot(Xl) +(Xsq[:,None] + Xsq[None,:]))) - else: - X1l = X/self.lengthscale - X2l = X2/self.lengthscale - X1sq = np.sum(np.square(X1l),1) - X2sq = np.sum(np.square(X2l),1) - return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) + if X2 is not None: + X2 = X2 / self.lengthscale + return self._unscaled_dist(X/self.lengthscale, X2) else: return self._unscaled_dist(X, X2)/self.lengthscale From 9867330861c6b515636c9c732771ed34c3ade677 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 14:30:28 +0000 Subject: [PATCH 19/51] moved plot functionality from add to kern --- GPy/kern/_src/add.py | 10 +--------- GPy/kern/_src/coregionalize.py | 11 ----------- GPy/kern/_src/kern.py | 8 ++++++++ GPy/kern/_src/rbf.py | 12 ++++++++---- 4 files changed, 17 insertions(+), 24 deletions(-) diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index 1466800c..77fe057d 100644 --- a/GPy/kern/_src/add.py +++ b/GPy/kern/_src/add.py @@ -184,20 +184,12 @@ class Add(Kern): target_S += b return target_mu, target_S - def plot(self, *args, **kwargs): - """ - See GPy.plotting.matplot_dep.plot - """ - assert "matplotlib" in sys.modules, "matplotlib package has not been imported." - from ..plotting.matplot_dep import kernel_plots - kernel_plots.plot(self,*args) - def input_sensitivity(self): in_sen = np.zeros((self.num_params, self.input_dim)) for i, [p, i_s] in enumerate(zip(self._parameters_, self.input_slices)): in_sen[i, i_s] = p.input_sensitivity() return in_sen - + def _getstate(self): """ Get the current state of the class, diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py index cafdd5ee..6679eba4 100644 --- a/GPy/kern/_src/coregionalize.py +++ b/GPy/kern/_src/coregionalize.py @@ -60,17 +60,6 @@ class Coregionalize(Kern): def K(self, X, X2=None): index = np.asarray(X, dtype=np.int) - #here's the old code (numpy) - #if index2 is None: - #index2 = index - #else: - #index2 = np.asarray(index2, dtype=np.int) - #false_target = target.copy() - #ii, jj = np.meshgrid(index, index2) - #ii, jj = ii.T, jj.T - #false_target += self.B[ii, jj] - - if X2 is None: target = np.empty((X.shape[0], X.shape[0]), dtype=np.float64) code=""" diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index d7d8f9ca..eb3291e0 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -68,6 +68,14 @@ class Kern(Parameterized): """ raise NotImplementedError + def plot(self, *args, **kwargs): + """ + See GPy.plotting.matplot_dep.plot + """ + assert "matplotlib" in sys.modules, "matplotlib package has not been imported." + from ..plotting.matplot_dep import kernel_plots + kernel_plots.plot(self,*args) + def plot_ARD(self, *args, **kw): """ See :class:`~GPy.plotting.matplot_dep.kernel_plots` diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 666a79f4..f5bafb48 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -178,6 +178,11 @@ class RBF(Stationary): return denom, dist, dist_sq, psi1 + #@cache_this(ignore_args=(1,)) + def _Z_distances(self, Z): + Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q + Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q + return Zhat, Zdist #@cache_this TODO def _psi2computations(self, Z, vp): @@ -187,8 +192,7 @@ class RBF(Stationary): M = Z.shape[0] #compute required distances - Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q - Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q + Zhat, Zdist = self._Z_distances(Z) Zdist_sq = np.square(Zdist / self.lengthscale) # M,M,Q #allocate memory for the things we want to compute @@ -201,7 +205,7 @@ class RBF(Stationary): denom = (2.*S[:,None,None,:] / l2) + 1. # N,Q half_log_denom = 0.5 * np.log(denom[:,0,0,:]) denom_l2 = denom[:,0,0,:]*l2 - + variance_sq = float(np.square(self.variance)) code = """ double tmp, exponent_tmp; @@ -237,7 +241,7 @@ class RBF(Stationary): } } """ - + support_code = """ #include #include From 65fd6dd24e0cd840106783ac9c82452298bb7c1a Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 15:46:14 +0000 Subject: [PATCH 20/51] observable pattern through and thorugh --- GPy/core/parameterization/array_core.py | 5 +- GPy/core/parameterization/param.py | 17 +++--- GPy/core/parameterization/parameter_core.py | 27 +++++---- GPy/core/parameterization/parameterized.py | 4 +- GPy/kern/_src/linear.py | 1 - GPy/kern/_src/rbf.py | 5 +- GPy/kern/_src/ssrbf.py | 5 +- GPy/kern/_src/stationary.py | 5 ++ GPy/models/gplvm.py | 2 +- GPy/models/sparse_gp_regression.py | 6 +- GPy/util/caching.py | 67 +++++++-------------- 11 files changed, 64 insertions(+), 80 deletions(-) diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index e8be0f77..e5b2ac07 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -59,11 +59,10 @@ class ObservableArray(np.ndarray, Observable): else: return s.size != 0 - def __setitem__(self, s, val, update=True): + def __setitem__(self, s, val): if self._s_not_empty(s): super(ObservableArray, self).__setitem__(s, val) - if update: - self._notify_observers() + self._notify_observers() def __getslice__(self, start, stop): return self.__getitem__(slice(start, stop)) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 8eac69da..53b2260e 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -137,8 +137,6 @@ class Param(Constrainable, ObservableArray, Gradcheckable): #=========================================================================== def _set_params(self, param, update=True): self.flat = param - #self._notify_tied_parameters() - self._notify_observers() def _get_params(self): return self.flat @@ -161,12 +159,10 @@ class Param(Constrainable, ObservableArray, Gradcheckable): try: new_arr._current_slice_ = s; new_arr._original_ = self.base is new_arr.base except AttributeError: pass # returning 0d array or float, double etc return new_arr - def __setitem__(self, s, val, update=True): - super(Param, self).__setitem__(s, val, update=update) - #self._notify_tied_parameters() - if update and self._s_not_empty(s): - self._notify_parameters_changed() - + def __setitem__(self, s, val): + super(Param, self).__setitem__(s, val) + #self._notify_observers() + #=========================================================================== # Index Operations: #=========================================================================== @@ -185,6 +181,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable): a = self._realshape_[i] + a internal_offset += a * extended_realshape[i] return internal_offset + def _raveled_index(self, slice_index=None): # return an index array on the raveled array, which is formed by the current_slice # of this object @@ -354,7 +351,7 @@ class ParamConcatenation(object): val = val._vals() ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; vals = self._vals(); vals[s] = val; del val - [numpy.place(p, ind[ps], vals[ps]) and update and p._notify_parameters_changed() + [numpy.place(p, ind[ps], vals[ps]) and update and p._notify_observers() for p, ps in zip(self.params, self._param_slices_)] def _vals(self): return numpy.hstack([p._get_params() for p in self.params]) @@ -363,7 +360,7 @@ class ParamConcatenation(object): #=========================================================================== def update_all_params(self): for p in self.params: - p._notify_parameters_changed() + p._notify_observers() def constrain(self, constraint, warning=True): [param.constrain(constraint, update=False) for param in self.params] diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 28d63b02..e1e5c21c 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -18,8 +18,13 @@ class Observable(object): def add_observer(self, observer, callble): self._observer_callables_[observer].append(callble) - def remove_observer(self, observer, callble): - del self._observer_callables_[observer][callble] + def remove_observer(self, observer, callble=None): + if callble is None: + del self._observer_callables_[observer] + else: + self._observer_callables_[observer].remove(callble) + if len(self._observer_callables_[observer]) == 0: + self.remove_observer(observer) def _notify_observers(self): [[callble(self) for callble in callables] @@ -72,9 +77,8 @@ class Parentable(object): return self._direct_parent_._highest_parent_ def _notify_parameters_changed(self): - if self.has_parent(): - self._direct_parent_._notify_parameters_changed() - + raise NotImplementedError, "shouldnt happen, abstract superclass" + class Nameable(Parentable): def __init__(self, name, *a, **kw): super(Nameable, self).__init__(*a, **kw) @@ -309,7 +313,7 @@ class Constrainable(Nameable, Indexable): print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) which.add(transform, self._raveled_index()) if update: - self._notify_parameters_changed() + self._notify_observers() def _remove_from_index_operations(self, which, transforms): if len(transforms) == 0: @@ -325,7 +329,7 @@ class Constrainable(Nameable, Indexable): return removed -class Parameterizable(Constrainable): +class Parameterizable(Constrainable, Observable): def __init__(self, *args, **kwargs): super(Parameterizable, self).__init__(*args, **kwargs) from GPy.core.parameterization.array_core import ParamList @@ -386,7 +390,7 @@ class Parameterizable(Constrainable): def _set_params(self, params, update=True): # don't overwrite this anymore! import itertools - [p._set_params(params[s], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + [p._set_params(params[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] self.parameters_changed() def copy(self): @@ -420,11 +424,10 @@ class Parameterizable(Constrainable): return s - def _notify_parameters_changed(self): + def _notify_parameters_changed(self, which): self.parameters_changed() - if self.has_parent(): - self._direct_parent_._notify_parameters_changed() - + self._notify_observers() + def parameters_changed(self): """ This method gets called when parameters have changed. diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index d463ed43..16ae43a7 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -11,7 +11,7 @@ from parameter_core import Constrainable, Pickleable, Parentable, Observable, Pa from transformations import __fixed__ from array_core import ParamList -class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable): +class Parameterized(Parameterizable, Pickleable, Gradcheckable): """ Parameterized class @@ -92,6 +92,7 @@ class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable): self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) + param.add_observer(self, self._notify_parameters_changed) self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" @@ -120,6 +121,7 @@ class Parameterized(Parameterizable, Pickleable, Observable, Gradcheckable): del self._parameters_[param._parent_index_] param._disconnect_parent() + param.remove_observer(self, self._notify_parameters_changed) self.constraints.shift_left(start, param.size) self._connect_fixes() self._connect_parameters() diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py index 44f71d85..2c82eb7a 100644 --- a/GPy/kern/_src/linear.py +++ b/GPy/kern/_src/linear.py @@ -9,7 +9,6 @@ from ...util.linalg import tdot from ...util.misc import fast_array_equal, param_to_array from ...core.parameterization import Param from ...core.parameterization.transformations import Logexp -from ...util.caching import cache_this class Linear(Kern): """ diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 666a79f4..5e973aee 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -6,6 +6,7 @@ import numpy as np from scipy import weave from ...util.misc import param_to_array from stationary import Stationary +from GPy.util.caching import Cache_this class RBF(Stationary): """ @@ -166,7 +167,7 @@ class RBF(Stationary): return target - #@cache_this TODO + @Cache_this(limit=1) def _psi1computations(self, Z, vp): mu, S = vp.mean, vp.variance l2 = self.lengthscale **2 @@ -179,7 +180,7 @@ class RBF(Stationary): - #@cache_this TODO + @Cache_this(limit=1) def _psi2computations(self, Z, vp): mu, S = vp.mean, vp.variance diff --git a/GPy/kern/_src/ssrbf.py b/GPy/kern/_src/ssrbf.py index 60df7225..cd921acb 100644 --- a/GPy/kern/_src/ssrbf.py +++ b/GPy/kern/_src/ssrbf.py @@ -6,7 +6,6 @@ from kern import Kern import numpy as np from ...util.linalg import tdot from ...util.config import * -from ...util.caching import cache_this from stationary import Stationary class SSRBF(Stationary): @@ -155,7 +154,7 @@ class SSRBF(Stationary): # Precomputations # #---------------------------------------# - @cache_this(1) + #@cache_this(1) def _K_computations(self, X, X2): """ K(X,X2) - X is NxQ @@ -175,7 +174,7 @@ class SSRBF(Stationary): self._K_dist2 = -2.*np.dot(X, X2.T) + (np.sum(np.square(X), axis=1)[:, None] + np.sum(np.square(X2), axis=1)[None, :]) self._K_dvar = np.exp(-0.5 * self._K_dist2) - @cache_this(1) + #@cache_this(1) def _psi_computations(self, Z, mu, S, gamma): """ Z - MxQ diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index b3c648fc..ebc6b880 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -9,6 +9,7 @@ from ...util.linalg import tdot from ... import util import numpy as np from scipy import integrate +from ...util.caching import Cache_this class Stationary(Kern): def __init__(self, input_dim, variance, lengthscale, ARD, name): @@ -39,15 +40,18 @@ class Stationary(Kern): def dK_dr(self, r): raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" + @Cache_this(limit=5, ignore_args=()) def K(self, X, X2=None): r = self._scaled_dist(X, X2) return self.K_of_r(r) + @Cache_this(limit=5, ignore_args=(0,)) def _dist(self, X, X2): if X2 is None: X2 = X return X[:, None, :] - X2[None, :, :] + @Cache_this(limit=5, ignore_args=(0,)) def _unscaled_dist(self, X, X2=None): """ Compute the square distance between each row of X and X2, or between @@ -61,6 +65,7 @@ class Stationary(Kern): X2sq = np.sum(np.square(X2),1) return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) + @Cache_this(limit=5, ignore_args=()) def _scaled_dist(self, X, X2=None): """ Efficiently compute the scaled distance, r. diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py index c2fd46fe..ba270dad 100644 --- a/GPy/models/gplvm.py +++ b/GPy/models/gplvm.py @@ -41,7 +41,7 @@ class GPLVM(GP): def parameters_changed(self): super(GPLVM, self).parameters_changed() - self.X.gradient = self.kern.gradients_X(self._dL_dK, self.X, None) + self.X.gradient = self.kern.gradients_X(self.dL_dK, self.X, None) def _getstate(self): return GP._getstate(self) diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py index 54c89a89..6a76df3f 100644 --- a/GPy/models/sparse_gp_regression.py +++ b/GPy/models/sparse_gp_regression.py @@ -8,6 +8,7 @@ from .. import likelihoods from .. import kern from ..inference.latent_function_inference import VarDTC from ..util.misc import param_to_array +from ..core.parameterization.variational import VariationalPosterior class SparseGPRegression(SparseGP): """ @@ -44,7 +45,10 @@ class SparseGPRegression(SparseGP): assert Z.shape[1] == input_dim likelihood = likelihoods.Gaussian() - + + if not (X_variance is None): + X = VariationalPosterior(X,X_variance) + SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC()) def _getstate(self): diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 55e546df..2c7be9bd 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -1,21 +1,23 @@ from ..core.parameterization.parameter_core import Observable class Cacher(object): - def __init__(self, operation, limit=5, reset_on_first=False): + def __init__(self, operation, limit=5, ignore_args=()): self.limit = int(limit) - self._reset_on_first = reset_on_first + self.ignore_args = ignore_args self.operation=operation self.cached_inputs = [] self.cached_outputs = [] self.inputs_changed = [] def __call__(self, *args): - if self._reset_on_first: - assert isinstance(args[0], Observable) - args[0].add_observer(self, self.reset) - cached_args = args + if len(self.ignore_args) != 0: + ca = [a for i,a in enumerate(args) if i not in self.ignore_args] + cached_args = [] + for a in ca: + if not any(a is ai for ai in cached_args): + cached_args.append(a) else: - cached_args = args[1:] + cached_args = args if not all([isinstance(arg, Observable) for arg in cached_args]): @@ -36,7 +38,7 @@ class Cacher(object): self.cached_inputs.append(cached_args) self.cached_outputs.append(self.operation(*args)) self.inputs_changed.append(False) - [a.add_observer(self, self.on_cache_changed) for a in args] + [a.add_observer(self, self.on_cache_changed) for a in cached_args] return self.cached_outputs[-1] def on_cache_changed(self, arg): @@ -48,42 +50,15 @@ class Cacher(object): self.cached_outputs = [] self.inputs_changed = [] - - - -def cache_this(limit=5, reset_on_self=False): - def limited_cache(f): - c = Cacher(f, limit, reset_on_first=reset_on_self) +class Cache_this(object): + def __init__(self, limit=5, ignore_args=()): + self.limit = limit + self.ignore_args = ignore_args + self.c = None + def __call__(self, f): def f_wrap(*args): - return c(*args) - f_wrap._cacher = c - return f_wrap - return limited_cache - - - - - - - - - - - - - #Xbase = X - #while Xbase is not None: - #try: - #i = self.cached_inputs.index(X) - #break - #except ValueError: - #Xbase = X.base - #continue - #self.inputs_changed[i] = True - - - - - - - + if self.c is None: + self.c = Cacher(f, self.limit, ignore_args=self.ignore_args) + return self.c(*args) + f_wrap._cacher = self + return f_wrap \ No newline at end of file From f1487c693576c69e45aa7ee744a7bb5898c0b826 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 26 Feb 2014 15:52:55 +0000 Subject: [PATCH 21/51] parameters changed notify added --- GPy/core/parameterization/parameter_core.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index e1e5c21c..b8dca23c 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -392,6 +392,7 @@ class Parameterizable(Constrainable, Observable): import itertools [p._set_params(params[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] self.parameters_changed() + self._notify_observers() def copy(self): """Returns a (deep) copy of the current model""" From 16e2fbd6cbe5c1f946755aff5cef2741b0939bea Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 15:54:43 +0000 Subject: [PATCH 22/51] pydot graphing half done --- GPy/core/parameterization/param.py | 9 +++++++++ GPy/core/parameterization/parameterized.py | 20 ++++++++++++++++++++ 2 files changed, 29 insertions(+) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 53b2260e..b48bec0d 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -60,6 +60,15 @@ class Param(Constrainable, ObservableArray, Gradcheckable): def __init__(self, name, input_array, default_constraint=None, *a, **kw): super(Param, self).__init__(name=name, default_constraint=default_constraint, *a, **kw) + def build_pydot(self,G): + import pydot + node = pydot.Node(id(self), shape='record', label=self.name) + G.add_node(node) + for o in self._observer_callables_.keys(): + print o, self.hirarchy_name() + + return node + def __array_finalize__(self, obj): # see InfoArray.__array_finalize__ for comments if obj is None: return diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 16ae43a7..83a62161 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -64,6 +64,26 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self._connect_parameters() del self._in_init_ + def build_pydot(self, G=None): + import pydot + iamroot = False + if G is None: + G = pydot.Dot(graph_type='digraph') + iamroot=True + node = pydot.Node(id(self), shape='record', label=self.name) + G.add_node(node) + for child in self._parameters_: + child_node = child.build_pydot(G) + G.add_edge(pydot.Edge(node, child_node)) + + for o in self._observer_callables_.keys(): + print id(o), self.hirarchy_name() + + if iamroot: + return G + return node + + def add_parameter(self, param, index=None): """ :param parameters: the parameters to add From 68f0af4deb139041fb68782df26e192b8248be38 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 26 Feb 2014 17:13:47 +0000 Subject: [PATCH 23/51] fixes in the plotting and in the dot graphing --- GPy/core/parameterization/param.py | 6 +++++- GPy/core/parameterization/parameterized.py | 6 +++++- GPy/plotting/matplot_dep/base_plots.py | 20 ++++++++++---------- GPy/plotting/matplot_dep/models_plots.py | 2 +- 4 files changed, 21 insertions(+), 13 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index b48bec0d..3b08eebd 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -65,7 +65,11 @@ class Param(Constrainable, ObservableArray, Gradcheckable): node = pydot.Node(id(self), shape='record', label=self.name) G.add_node(node) for o in self._observer_callables_.keys(): - print o, self.hirarchy_name() + label = o.name if hasattr(o, 'name') else str(o) + observed_node = pydot.Node(id(o), label=label) + G.add_node(observed_node) + edge = pydot.Edge(str(id(self)), str(id(o)), color='darkorange2', arrowhead='vee') + G.add_edge(edge) return node diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 83a62161..fb606278 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -77,7 +77,11 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): G.add_edge(pydot.Edge(node, child_node)) for o in self._observer_callables_.keys(): - print id(o), self.hirarchy_name() + label = o.name if hasattr(o, 'name') else str(o) + observed_node = pydot.Node(id(o), label=label) + G.add_node(observed_node) + edge = pydot.Edge(str(id(self)), str(id(o)), color='darkorange2', arrowhead='vee') + G.add_edge(edge) if iamroot: return G diff --git a/GPy/plotting/matplot_dep/base_plots.py b/GPy/plotting/matplot_dep/base_plots.py index a9d25223..e86ef6ca 100644 --- a/GPy/plotting/matplot_dep/base_plots.py +++ b/GPy/plotting/matplot_dep/base_plots.py @@ -18,17 +18,17 @@ def meanplot(x, mu, color=Tango.colorsHex['darkBlue'], ax=None, fignum=None, lin _, axes = ax_default(fignum, ax) #here's the mean return axes.plot(x,mu,color=color,linewidth=linewidth,**kw) - + def gpplot(x,mu,lower,upper,edgecol=Tango.colorsHex['darkBlue'],fillcol=Tango.colorsHex['lightBlue'],ax=None,fignum=None,xlabel='x',ylabel='y',**kwargs): - _, axes = ax_default(ax, fignum) - + _, axes = ax_default(fignum, ax) + mu = mu.flatten() x = x.flatten() lower = lower.flatten() upper = upper.flatten() plots = [] - + #here's the mean plots.append(meanplot(x, mu, edgecol, axes)) @@ -39,15 +39,15 @@ def gpplot(x,mu,lower,upper,edgecol=Tango.colorsHex['darkBlue'],fillcol=Tango.co plots.append(axes.fill(np.hstack((x,x[::-1])),np.hstack((upper,lower[::-1])),color=fillcol,**kwargs)) #this is the edge: - plots.append(meanplot(x, upper,color=edgecol,linewidth=0.2,axes=axes)) - plots.append(meanplot(x, lower,color=edgecol,linewidth=0.2,axes=axes)) - + plots.append(meanplot(x, upper,color=edgecol,linewidth=0.2,ax=axes)) + plots.append(meanplot(x, lower,color=edgecol,linewidth=0.2,ax=axes)) + axes.set_xlabel(xlabel) axes.set_ylabel(ylabel) - + return plots - - + + def removeRightTicks(ax=None): ax = ax or pb.gca() for i, line in enumerate(ax.get_yticklines()): diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 5123f514..d72d2a3e 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -86,7 +86,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', m, v, lower, upper = model.predict(Xgrid) Y = Y for d in which_data_ycols: - gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax, edgecol=linecol, fillcol=fillcol) + gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], ax=ax, edgecol=linecol, fillcol=fillcol) ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5) #optionally plot some samples From 023203dc3a449c6ef046774c8539a5aba5c2f4f2 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 27 Feb 2014 08:07:00 +0000 Subject: [PATCH 24/51] parent observer now static and always last --- GPy/core/parameterization/param.py | 2 ++ GPy/core/parameterization/parameter_core.py | 5 +++-- GPy/core/parameterization/parameterized.py | 1 - 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 53b2260e..96df98f4 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -161,6 +161,8 @@ class Param(Constrainable, ObservableArray, Gradcheckable): return new_arr def __setitem__(self, s, val): super(Param, self).__setitem__(s, val) + if self.has_parent(): + self._direct_parent_._notify_parameters_changed() #self._notify_observers() #=========================================================================== diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index b8dca23c..55fe660c 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -391,8 +391,7 @@ class Parameterizable(Constrainable, Observable): # don't overwrite this anymore! import itertools [p._set_params(params[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - self.parameters_changed() - self._notify_observers() + self._notify_parameters_changed() def copy(self): """Returns a (deep) copy of the current model""" @@ -428,6 +427,8 @@ class Parameterizable(Constrainable, Observable): def _notify_parameters_changed(self, which): self.parameters_changed() self._notify_observers() + if self.has_parent(): + self._direct_parent_._notify_parameters_changed() def parameters_changed(self): """ diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 16ae43a7..18cea61a 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -92,7 +92,6 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) - param.add_observer(self, self._notify_parameters_changed) self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" From 37d341cec4b243530c1367ce5fa71783f821e3e1 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 27 Feb 2014 08:10:54 +0000 Subject: [PATCH 25/51] parent observer now static and always last --- 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 55fe660c..94567761 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -424,7 +424,7 @@ class Parameterizable(Constrainable, Observable): return s - def _notify_parameters_changed(self, which): + def _notify_parameters_changed(self): self.parameters_changed() self._notify_observers() if self.has_parent(): From ebea658f5ccc9e329fc6d336dd503fd941266f22 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 27 Feb 2014 08:39:49 +0000 Subject: [PATCH 26/51] caching can handle None --- GPy/util/caching.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 2c7be9bd..47a44d58 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -12,16 +12,17 @@ class Cacher(object): def __call__(self, *args): if len(self.ignore_args) != 0: ca = [a for i,a in enumerate(args) if i not in self.ignore_args] - cached_args = [] - for a in ca: - if not any(a is ai for ai in cached_args): - cached_args.append(a) else: - cached_args = args - + ca = args + # this makes sure we only add an observer once, and that None can be in args + cached_args = [] + for a in ca: + if (not any(a is ai for ai in cached_args)) and a is not None: + cached_args.append(a) if not all([isinstance(arg, Observable) for arg in cached_args]): return self.operation(*args) + if cached_args in self.cached_inputs: i = self.cached_inputs.index(cached_args) if self.inputs_changed[i]: From a6d3fda23420b7669429ed6736266f182e66dfbf Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 27 Feb 2014 08:47:06 +0000 Subject: [PATCH 27/51] caching can handle None --- GPy/core/parameterization/array_core.py | 2 +- GPy/core/parameterization/parameter_core.py | 13 +++++++------ GPy/util/caching.py | 4 +++- 3 files changed, 11 insertions(+), 8 deletions(-) diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index e5b2ac07..a338ceed 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -4,7 +4,7 @@ __updated__ = '2013-12-16' import numpy as np -from parameter_core import Observable, Parameterizable +from parameter_core import Observable class ParamList(list): """ diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 94567761..6afa94cb 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -19,12 +19,13 @@ class Observable(object): self._observer_callables_[observer].append(callble) def remove_observer(self, observer, callble=None): - if callble is None: - del self._observer_callables_[observer] - else: - self._observer_callables_[observer].remove(callble) - if len(self._observer_callables_[observer]) == 0: - self.remove_observer(observer) + if observer in self._observer_callables_: + if callble is None: + del self._observer_callables_[observer] + elif callble in self._observer_callables_[observer]: + self._observer_callables_[observer].remove(callble) + if len(self._observer_callables_[observer]) == 0: + self.remove_observer(observer) def _notify_observers(self): [[callble(self) for callble in callables] diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 47a44d58..8e60cf26 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -19,8 +19,9 @@ class Cacher(object): for a in ca: if (not any(a is ai for ai in cached_args)) and a is not None: cached_args.append(a) - if not all([isinstance(arg, Observable) for arg in cached_args]): + print cached_args + import ipdb;ipdb.set_trace() return self.operation(*args) if cached_args in self.cached_inputs: @@ -46,6 +47,7 @@ class Cacher(object): self.inputs_changed = [any([a is arg for a in args]) or old_ic for args, old_ic in zip(self.cached_inputs, self.inputs_changed)] def reset(self, obj): + [[a.remove_observer(self, self.on_cache_changed) for a in args] for args in self.cached_inputs] [[a.remove_observer(self, self.reset) for a in args] for args in self.cached_inputs] self.cached_inputs = [] self.cached_outputs = [] From 283b0745aac0653ac18bc00821f4c8e7e59e6fa3 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 27 Feb 2014 09:16:36 +0000 Subject: [PATCH 28/51] added some caching --- GPy/kern/_src/linear.py | 20 ++++++++++---------- GPy/kern/_src/rbf.py | 1 + 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py index 2c82eb7a..ab5030eb 100644 --- a/GPy/kern/_src/linear.py +++ b/GPy/kern/_src/linear.py @@ -49,13 +49,8 @@ class Linear(Kern): self.variances = Param('variances', variances, Logexp()) self.add_parameter(self.variances) - self.variances.add_observer(self, self._on_changed) - def _on_changed(self, obj): - #TODO: move this to base class? isnt it jst for the caching? - self._notify_observers() - - #@cache_this(limit=3, reset_on_self=True) + @Cache_this(limit=2) def K(self, X, X2=None): if self.ARD: if X2 is None: @@ -66,7 +61,7 @@ class Linear(Kern): else: return self._dot_product(X, X2) * self.variances - #@cache_this(limit=3, reset_on_self=False) + @Cache_this(limit=1, ignore_args(0,)) def _dot_product(self, X, X2=None): if X2 is None: return tdot(X) @@ -113,6 +108,7 @@ class Linear(Kern): def psi1(self, Z, variational_posterior): return self.K(variational_posterior.mean, Z) #the variance, it does nothing + @Cache_this(limit=1) def psi2(self, Z, variational_posterior): ZA = Z * self.variances ZAinner = self._ZAinner(variational_posterior, Z) @@ -126,9 +122,11 @@ class Linear(Kern): if self.ARD: self.variances.gradient += tmp.sum(0) else: self.variances.gradient += tmp.sum() #psi2 - tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(variational_posterior, Z)[:, :, None, :] * (2. * Z)[None, None, :, :]) - if self.ARD: self.variances.gradient += tmp.sum(0).sum(0).sum(0) - else: self.variances.gradient += tmp.sum() + if self.ARD: + tmp = dL_dpsi2[:, :, :, None] * (self._ZAinner(variational_posterior, Z)[:, :, None, :] * Z[None, None, :, :]) + self.variances.gradient += 2.*tmp.sum(0).sum(0).sum(0) + else: + self.variances.gradient += 2.*np.sum(dL_dpsi2 * self.psi2(Z, variational_posterior))/self.variances def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): #psi1 @@ -234,9 +232,11 @@ class Linear(Kern): type_converters=weave.converters.blitz,**weave_options) + @Cache_this(limit=1, ignore_args=(0,)) def _mu2S(self, vp): return np.square(vp.mean) + vp.variance + @Cache_this(limit=1) def _ZAinner(self, vp, Z): ZA = Z*self.variances inner = (vp.mean[:, None, :] * vp.mean[:, :, None]) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 342b065e..cf5ea0c4 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -113,6 +113,7 @@ class RBF(Stationary): # Precomputations # #---------------------------------------# + #TODO: this function is unused, but it will be useful in the stationary class def _dL_dlengthscales_via_K(self, dL_dK, X, X2): """ A helper function for update_gradients_* methods From 2feb849bf75bb647a08b1884d07b1bec5cd222c0 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 27 Feb 2014 10:02:34 +0000 Subject: [PATCH 29/51] linear fix --- GPy/kern/_src/linear.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py index ab5030eb..1521026d 100644 --- a/GPy/kern/_src/linear.py +++ b/GPy/kern/_src/linear.py @@ -9,6 +9,7 @@ from ...util.linalg import tdot from ...util.misc import fast_array_equal, param_to_array from ...core.parameterization import Param from ...core.parameterization.transformations import Logexp +from ...util.caching import Cache_this class Linear(Kern): """ @@ -61,7 +62,7 @@ class Linear(Kern): else: return self._dot_product(X, X2) * self.variances - @Cache_this(limit=1, ignore_args(0,)) + @Cache_this(limit=1, ignore_args=(0,)) def _dot_product(self, X, X2=None): if X2 is None: return tdot(X) From 5cf792504a653b5ae6cc194b66f0eb7e04a3b273 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 27 Feb 2014 10:47:27 +0000 Subject: [PATCH 30/51] messing with caching --- GPy/kern/_src/stationary.py | 8 ++--- GPy/testing/kernel_tests.py | 19 +++++----- GPy/util/caching.py | 71 +++++++++++++++++++++++++++---------- 3 files changed, 64 insertions(+), 34 deletions(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index ebc6b880..8d8ae476 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -40,18 +40,18 @@ class Stationary(Kern): def dK_dr(self, r): raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" - @Cache_this(limit=5, ignore_args=()) + #@Cache_this(limit=5, ignore_args=()) def K(self, X, X2=None): r = self._scaled_dist(X, X2) return self.K_of_r(r) - @Cache_this(limit=5, ignore_args=(0,)) + #@Cache_this(limit=5, ignore_args=(0,)) def _dist(self, X, X2): if X2 is None: X2 = X return X[:, None, :] - X2[None, :, :] - @Cache_this(limit=5, ignore_args=(0,)) + #@Cache_this(limit=5, ignore_args=(0,)) def _unscaled_dist(self, X, X2=None): """ Compute the square distance between each row of X and X2, or between @@ -65,7 +65,7 @@ class Stationary(Kern): X2sq = np.sum(np.square(X2),1) return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) - @Cache_this(limit=5, ignore_args=()) + #@Cache_this(limit=5, ignore_args=()) def _scaled_dist(self, X, X2=None): """ Efficiently compute the scaled distance, r. diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index e5985145..d373a546 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -8,13 +8,6 @@ import sys verbose = True -try: - import sympy - SYMPY_AVAILABLE=True -except ImportError: - SYMPY_AVAILABLE=False - - class Kern_check_model(GPy.core.Model): """ This is a dummy model class used as a base class for checking that the @@ -70,14 +63,11 @@ class Kern_check_dKdiag_dtheta(Kern_check_model): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=None) self.add_parameter(self.kernel) - def parameters_changed(self): - self.kernel.update_gradients_diag(self.dL_dK, self.X) - def log_likelihood(self): return (np.diag(self.dL_dK)*self.kernel.Kdiag(self.X)).sum() def parameters_changed(self): - return self.kernel.update_gradients_diag(np.diag(self.dL_dK), self.X) + self.kernel.update_gradients_diag(np.diag(self.dL_dK), self.X) class Kern_check_dK_dX(Kern_check_model): """This class allows gradient checks for the gradient of a kernel with respect to X. """ @@ -99,6 +89,8 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX): def parameters_changed(self): self.X.gradient = self.kernel.gradients_X_diag(self.dL_dK, self.X) + + def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): """ This function runs on kernels to check the correctness of their @@ -217,11 +209,15 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): return pass_checks + class KernelTestsContinuous(unittest.TestCase): def setUp(self): self.X = np.random.randn(100,2) self.X2 = np.random.randn(110,2) + continuous_kerns = ['RBF', 'Linear'] + self.kernclasses = [getattr(GPy.kern, s) for s in continuous_kerns] + def test_Matern32(self): k = GPy.kern.Matern32(2) self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) @@ -234,6 +230,7 @@ class KernelTestsContinuous(unittest.TestCase): + if __name__ == "__main__": print "Running unit tests, please be (very) patient..." unittest.main() diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 8e60cf26..2899cb33 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -1,6 +1,13 @@ from ..core.parameterization.parameter_core import Observable class Cacher(object): + """ + + + + + """ + def __init__(self, operation, limit=5, ignore_args=()): self.limit = int(limit) self.ignore_args = ignore_args @@ -10,50 +17,75 @@ class Cacher(object): self.inputs_changed = [] def __call__(self, *args): + """ + A wrapper function for self.operation, + """ + + #ensure that specified arguments are ignored if len(self.ignore_args) != 0: - ca = [a for i,a in enumerate(args) if i not in self.ignore_args] + oa = [a for i,a in enumerate(args) if i not in self.ignore_args] else: - ca = args + oa = args + # this makes sure we only add an observer once, and that None can be in args - cached_args = [] - for a in ca: - if (not any(a is ai for ai in cached_args)) and a is not None: - cached_args.append(a) - if not all([isinstance(arg, Observable) for arg in cached_args]): - print cached_args - import ipdb;ipdb.set_trace() + observable_args = [] + for a in oa: + if (not any(a is ai for ai in observable_args)) and a is not None: + observable_args.append(a) + + #make sure that all the found argument really are observable: + #otherswise don't cache anything, pass args straight though + if not all([isinstance(arg, Observable) for arg in observable_args]): return self.operation(*args) - - if cached_args in self.cached_inputs: - i = self.cached_inputs.index(cached_args) + + #if the result is cached, return the cached computation + state = [all(a is b for a, b in zip(args, cached_i)) for cached_i in self.cached_inputs] + if any(state): + i = state.index(True) if self.inputs_changed[i]: + #(elements of) the args have changed since we last computed: update self.cached_outputs[i] = self.operation(*args) self.inputs_changed[i] = False return self.cached_outputs[i] else: + #first time we've seen these arguments: compute + + #first make sure the depth limit isn't exceeded if len(self.cached_inputs) == self.limit: args_ = self.cached_inputs.pop(0) - [a.remove_observer(self, self.on_cache_changed) for a in args_] + [a.remove_observer(self, self.on_cache_changed) for a in args_ if a is not None] self.inputs_changed.pop(0) self.cached_outputs.pop(0) - self.cached_inputs.append(cached_args) + #compute + self.cached_inputs.append(args) self.cached_outputs.append(self.operation(*args)) self.inputs_changed.append(False) - [a.add_observer(self, self.on_cache_changed) for a in cached_args] - return self.cached_outputs[-1] + [a.add_observer(self, self.on_cache_changed) for a in observable_args] + return self.cached_outputs[-1]#Max says return. def on_cache_changed(self, arg): + """ + A callback funtion, which sets local flags when the elements of some cached inputs change + + this function gets 'hooked up' to the inputs when we cache them, and upon their elements being changed we update here. + """ self.inputs_changed = [any([a is arg for a in args]) or old_ic for args, old_ic in zip(self.cached_inputs, self.inputs_changed)] def reset(self, obj): - [[a.remove_observer(self, self.on_cache_changed) for a in args] for args in self.cached_inputs] - [[a.remove_observer(self, self.reset) for a in args] for args in self.cached_inputs] + """ + Totally reset the cache + """ + [[a.remove_observer(self, self.on_cache_changed) for a in args if isinstance(a, Observable)] for args in self.cached_inputs] + [[a.remove_observer(self, self.reset) for a in args if isinstance(a, Observable)] for args in self.cached_inputs] self.cached_inputs = [] self.cached_outputs = [] self.inputs_changed = [] class Cache_this(object): + """ + A decorator which can be applied to bound methods in order to cache them + """ def __init__(self, limit=5, ignore_args=()): self.limit = limit self.ignore_args = ignore_args @@ -64,4 +96,5 @@ class Cache_this(object): self.c = Cacher(f, self.limit, ignore_args=self.ignore_args) return self.c(*args) f_wrap._cacher = self - return f_wrap \ No newline at end of file + f_wrap.__doc__ = "**cached**\n\n" + (f.__doc__ or "") + return f_wrap From d8c2f7813159ab94d2187d293ccf5994f1d02c8d Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Thu, 27 Feb 2014 16:28:42 +0000 Subject: [PATCH 31/51] [SSGPLVM] update SSGPLVM with new inferface and merge ssrbf into rbf --- GPy/kern/_src/rbf.py | 61 ++++++- GPy/kern/_src/rbf_psi_comp/__init__.py | 2 + GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py | 111 ++++++++++++ GPy/kern/_src/ssrbf.py | 178 ++++--------------- GPy/models/bayesian_gplvm.py | 31 ---- GPy/models/ss_gplvm.py | 2 +- 6 files changed, 206 insertions(+), 179 deletions(-) create mode 100644 GPy/kern/_src/rbf_psi_comp/__init__.py create mode 100644 GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index cf5ea0c4..7bf0adeb 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -7,6 +7,8 @@ from scipy import weave from ...util.misc import param_to_array from stationary import Stationary from GPy.util.caching import Cache_this +from ...core.parameterization import variational +from rbf_psi_comp import ssrbf_psi_comp class RBF(Stationary): """ @@ -36,14 +38,38 @@ class RBF(Stationary): return self.Kdiag(variational_posterior.mean) def psi1(self, Z, variational_posterior): - _, _, _, psi1 = self._psi1computations(Z, variational_posterior) + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + psi1, _, _, _, _, _, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + else: + _, _, _, psi1 = self._psi1computations(Z, variational_posterior) return psi1 def psi2(self, Z, variational_posterior): - _, _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + psi2, _, _, _, _, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + else: + _, _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) return psi2 def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + # Spike-and-Slab GPLVM + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + _, _dpsi1_dvariance, _, _, _, _, _dpsi1_dlengthscale = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + _, _dpsi2_dvariance, _, _, _, _, _dpsi2_dlengthscale = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + + #contributions from psi0: + self.variance.gradient = np.sum(dL_dpsi0) + + #from psi1 + self.variance.gradient += np.sum(dL_dpsi1 * _dpsi1_dvariance) + self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) + + + #from psi2 + self.variance.gradient += (dL_dpsi2 * _dpsi2_dvariance).sum() + self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) + return + l2 = self.lengthscale **2 #contributions from psi0: @@ -77,6 +103,19 @@ class RBF(Stationary): self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + # Spike-and-Slab GPLVM + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + _, _, _, _, _, _dpsi1_dZ, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + _, _, _, _, _, _dpsi2_dZ, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + + #psi1 + grad = (dL_dpsi1[:, :, None] * _dpsi1_dZ).sum(axis=0) + + #psi2 + grad += (dL_dpsi2[:, :, :, None] * _dpsi2_dZ).sum(axis=0).sum(axis=1) + + return grad + l2 = self.lengthscale **2 #psi1 @@ -95,6 +134,24 @@ class RBF(Stationary): return grad def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + # Spike-and-Slab GPLVM + if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): + ndata = variational_posterior.mean.shape[0] + + _, _, _dpsi1_dgamma, _dpsi1_dmu, _dpsi1_dS, _, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + _, _, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + + #psi1 + grad_mu = (dL_dpsi1[:, :, None] * _dpsi1_dmu).sum(axis=1) + grad_S = (dL_dpsi1[:, :, None] * _dpsi1_dS).sum(axis=1) + grad_gamma = (dL_dpsi1[:,:,None] * _dpsi1_dgamma).sum(axis=1) + #psi2 + grad_mu += (dL_dpsi2[:, :, :, None] * _dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1) + grad_S += (dL_dpsi2[:, :, :, None] * _dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1) + grad_gamma += (dL_dpsi2[:,:,:, None] * _dpsi2_dgamma).reshape(ndata,-1,self.input_dim).sum(axis=1) + + return grad_mu, grad_S, grad_gamma + l2 = self.lengthscale **2 #psi1 denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) diff --git a/GPy/kern/_src/rbf_psi_comp/__init__.py b/GPy/kern/_src/rbf_psi_comp/__init__.py new file mode 100644 index 00000000..4c0d373d --- /dev/null +++ b/GPy/kern/_src/rbf_psi_comp/__init__.py @@ -0,0 +1,2 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) diff --git a/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py b/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py new file mode 100644 index 00000000..f3d5ee6b --- /dev/null +++ b/GPy/kern/_src/rbf_psi_comp/ssrbf_psi_comp.py @@ -0,0 +1,111 @@ +# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Licensed under the BSD 3-clause license (see LICENSE.txt) + +""" +The package for the psi statistics computation +""" + +import numpy as np + +def _Z_distances(Z): + Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q + Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q + return Zhat, Zdist + +# def _psi1computations(self, Z, vp): +# mu, S = vp.mean, vp.variance +# l2 = lengthscale **2 +# denom = S[:, None, :] / l2 + 1. # N,1,Q +# dist = Z[None, :, :] - mu[:, None, :] # N,M,Q +# dist_sq = np.square(dist) / l2 / denom # N,M,Q +# exponent = -0.5 * np.sum(dist_sq + np.log(denom), -1)#N,M +# psi1 = self.variance * np.exp(exponent) # N,M +# return denom, dist, dist_sq, psi1 + +def _psi1computations(variance, lengthscale, Z, mu, S, gamma): + """ + Z - MxQ + mu - NxQ + S - NxQ + gamma - NxQ + """ + # here are the "statistics" for psi1 and psi2 + # Produced intermediate results: + # _psi1 NxM + # _dpsi1_dvariance NxM + # _dpsi1_dlengthscale NxMxQ + # _dpsi1_dZ NxMxQ + # _dpsi1_dgamma NxMxQ + # _dpsi1_dmu NxMxQ + # _dpsi1_dS NxMxQ + + lengthscale2 = np.square(lengthscale) + + # psi1 + _psi1_denom = S[:, None, :] / lengthscale2 + 1. # Nx1xQ + _psi1_denom_sqrt = np.sqrt(_psi1_denom) #Nx1xQ + _psi1_dist = Z[None, :, :] - mu[:, None, :] # NxMxQ + _psi1_dist_sq = np.square(_psi1_dist) / (lengthscale2 * _psi1_denom) # NxMxQ + _psi1_common = gamma[:,None,:] / (lengthscale2*_psi1_denom*_psi1_denom_sqrt) #Nx1xQ + _psi1_exponent1 = np.log(gamma[:,None,:]) -0.5 * (_psi1_dist_sq + np.log(_psi1_denom)) # NxMxQ + _psi1_exponent2 = np.log(1.-gamma[:,None,:]) -0.5 * (np.square(Z[None,:,:])/lengthscale2) # NxMxQ + _psi1_exponent = np.log(np.exp(_psi1_exponent1) + np.exp(_psi1_exponent2)) #NxMxQ + _psi1_exp_sum = _psi1_exponent.sum(axis=-1) #NxM + _psi1_exp_dist_sq = np.exp(-0.5*_psi1_dist_sq) # NxMxQ + _psi1_exp_Z = np.exp(-0.5*np.square(Z[None,:,:])/lengthscale2) # 1xMxQ + _psi1_q = variance * np.exp(_psi1_exp_sum[:,:,None] - _psi1_exponent) # NxMxQ + _psi1 = variance * np.exp(_psi1_exp_sum) # NxM + _dpsi1_dvariance = _psi1 / variance # NxM + _dpsi1_dgamma = _psi1_q * (_psi1_exp_dist_sq/_psi1_denom_sqrt-_psi1_exp_Z) # NxMxQ + _dpsi1_dmu = _psi1_q * (_psi1_exp_dist_sq * _psi1_dist * _psi1_common) # NxMxQ + _dpsi1_dS = _psi1_q * (_psi1_exp_dist_sq * _psi1_common * 0.5 * (_psi1_dist_sq - 1.)) # NxMxQ + _dpsi1_dZ = _psi1_q * (- _psi1_common * _psi1_dist * _psi1_exp_dist_sq - (1-gamma[:,None,:])/lengthscale2*Z[None,:,:]*_psi1_exp_Z) # NxMxQ + _dpsi1_dlengthscale = 2.*lengthscale*_psi1_q * (0.5*_psi1_common*(S[:,None,:]/lengthscale2+_psi1_dist_sq)*_psi1_exp_dist_sq + 0.5*(1-gamma[:,None,:])*np.square(Z[None,:,:]/lengthscale2)*_psi1_exp_Z) # NxMxQ + + return _psi1, _dpsi1_dvariance, _dpsi1_dgamma, _dpsi1_dmu, _dpsi1_dS, _dpsi1_dZ, _dpsi1_dlengthscale + +def _psi2computations(variance, lengthscale, Z, mu, S, gamma): + """ + Z - MxQ + mu - NxQ + S - NxQ + gamma - NxQ + """ + # here are the "statistics" for psi1 and psi2 + # Produced intermediate results: + # _psi2 NxMxM + # _psi2_dvariance NxMxM + # _psi2_dlengthscale NxMxMxQ + # _psi2_dZ NxMxMxQ + # _psi2_dgamma NxMxMxQ + # _psi2_dmu NxMxMxQ + # _psi2_dS NxMxMxQ + + lengthscale2 = np.square(lengthscale) + + _psi2_Zhat, _psi2_Zdist = _Z_distances(Z) + _psi2_Zdist_sq = np.square(_psi2_Zdist / lengthscale) # M,M,Q + _psi2_Z_sq_sum = (np.square(Z[:,None,:])+np.square(Z[None,:,:]))/lengthscale2 # MxMxQ + + # psi2 + _psi2_denom = 2.*S[:, None, None, :] / lengthscale2 + 1. # Nx1x1xQ + _psi2_denom_sqrt = np.sqrt(_psi2_denom) + _psi2_mudist = mu[:,None,None,:]-_psi2_Zhat #N,M,M,Q + _psi2_mudist_sq = np.square(_psi2_mudist)/(lengthscale2*_psi2_denom) + _psi2_common = gamma[:,None,None,:]/(lengthscale2 * _psi2_denom * _psi2_denom_sqrt) # Nx1x1xQ + _psi2_exponent1 = -_psi2_Zdist_sq -_psi2_mudist_sq -0.5*np.log(_psi2_denom)+np.log(gamma[:,None,None,:]) #N,M,M,Q + _psi2_exponent2 = np.log(1.-gamma[:,None,None,:]) - 0.5*(_psi2_Z_sq_sum) # NxMxMxQ + _psi2_exponent = np.log(np.exp(_psi2_exponent1) + np.exp(_psi2_exponent2)) + _psi2_exp_sum = _psi2_exponent.sum(axis=-1) #NxM + _psi2_q = np.square(variance) * np.exp(_psi2_exp_sum[:,:,:,None]-_psi2_exponent) # NxMxMxQ + _psi2_exp_dist_sq = np.exp(-_psi2_Zdist_sq -_psi2_mudist_sq) # NxMxMxQ + _psi2_exp_Z = np.exp(-0.5*_psi2_Z_sq_sum) # MxMxQ + _psi2 = np.square(variance) * np.exp(_psi2_exp_sum) # N,M,M + _dpsi2_dvariance = 2. * _psi2/variance # NxMxM + _dpsi2_dgamma = _psi2_q * (_psi2_exp_dist_sq/_psi2_denom_sqrt - _psi2_exp_Z) # NxMxMxQ + _dpsi2_dmu = _psi2_q * (-2.*_psi2_common*_psi2_mudist * _psi2_exp_dist_sq) # NxMxMxQ + _dpsi2_dS = _psi2_q * (_psi2_common * (2.*_psi2_mudist_sq - 1.) * _psi2_exp_dist_sq) # NxMxMxQ + _dpsi2_dZ = 2.*_psi2_q * (_psi2_common*(-_psi2_Zdist*_psi2_denom+_psi2_mudist)*_psi2_exp_dist_sq - (1-gamma[:,None,None,:])*Z[:,None,:]/lengthscale2*_psi2_exp_Z) # NxMxMxQ + _dpsi2_dlengthscale = 2.*lengthscale* _psi2_q * (_psi2_common*(S[:,None,None,:]/lengthscale2+_psi2_Zdist_sq*_psi2_denom+_psi2_mudist_sq)*_psi2_exp_dist_sq+(1-gamma[:,None,None,:])*_psi2_Z_sq_sum*0.5/lengthscale2*_psi2_exp_Z) # NxMxMxQ + + return _psi2, _dpsi2_dvariance, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _dpsi2_dZ, _dpsi2_dlengthscale diff --git a/GPy/kern/_src/ssrbf.py b/GPy/kern/_src/ssrbf.py index cd921acb..391ef1c7 100644 --- a/GPy/kern/_src/ssrbf.py +++ b/GPy/kern/_src/ssrbf.py @@ -7,6 +7,7 @@ import numpy as np from ...util.linalg import tdot from ...util.config import * from stationary import Stationary +from rbf_psi_comp import ssrbf_psi_comp class SSRBF(Stationary): """ @@ -54,101 +55,63 @@ class SSRBF(Stationary): # PSI statistics # #---------------------------------------# - def psi0(self, Z, posterior_variational): - ret = np.empty(posterior_variational.mean.shape[0]) + def psi0(self, Z, variational_posterior): + ret = np.empty(variational_posterior.mean.shape[0]) ret[:] = self.variance return ret - def psi1(self, Z, posterior_variational): - self._psi_computations(Z, posterior_variational.mean, posterior_variational.variance, posterior_variational.binary_prob) - return self._psi1 + def psi1(self, Z, variational_posterior): + _psi1, _, _, _, _, _, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + return _psi1 - def psi2(self, Z, posterior_variational): - self._psi_computations(Z, posterior_variational.mean, posterior_variational.variance, posterior_variational.binary_prob) - return self._psi2 + def psi2(self, Z, variational_posterior): + _psi2, _, _, _, _, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + return _psi2 - def dL_dpsi0_dmuSgamma(self, dL_dpsi0, Z, mu, S, gamma, target_mu, target_S, target_gamma): - pass - - - def dL_dpsi1_dmuSgamma(self, dL_dpsi1, Z, mu, S, gamma, target_mu, target_S, target_gamma): - self._psi_computations(Z, mu, S, gamma) - target_mu += (dL_dpsi1[:, :, None] * self._dpsi1_dmu).sum(axis=1) - target_S += (dL_dpsi1[:, :, None] * self._dpsi1_dS).sum(axis=1) - target_gamma += (dL_dpsi1[:,:,None] * self._dpsi1_dgamma).sum(axis=1) - - - def dL_dpsi2_dmuSgamma(self, dL_dpsi2, Z, mu, S, gamma, target_mu, target_S, target_gamma): - """Think N,num_inducing,num_inducing,input_dim """ - self._psi_computations(Z, mu, S, gamma) - target_mu += (dL_dpsi2[:, :, :, None] * self._dpsi2_dmu).reshape(mu.shape[0],-1,mu.shape[1]).sum(axis=1) - target_S += (dL_dpsi2[:, :, :, None] * self._dpsi2_dS).reshape(S.shape[0],-1,S.shape[1]).sum(axis=1) - target_gamma += (dL_dpsi2[:,:,:, None] *self._dpsi2_dgamma).reshape(gamma.shape[0],-1,gamma.shape[1]).sum(axis=1) - - def update_gradients_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, posterior_variational): - self._psi_computations(Z, posterior_variational.mean, posterior_variational.variance, posterior_variational.binary_prob) + def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + _, _dpsi1_dvariance, _, _, _, _, _dpsi1_dlengthscale = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + _, _dpsi2_dvariance, _, _, _, _, _dpsi2_dlengthscale = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) #contributions from psi0: self.variance.gradient = np.sum(dL_dpsi0) #from psi1 - self.variance.gradient += np.sum(dL_dpsi1 * self._dpsi1_dvariance) - self.lengthscale.gradient = (dL_dpsi1[:,:,None]*self._dpsi1_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) + self.variance.gradient += np.sum(dL_dpsi1 * _dpsi1_dvariance) + self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) #from psi2 - self.variance.gradient += (dL_dpsi2 * self._dpsi2_dvariance).sum() - self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * self._dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) - - #from Kmm - self._K_computations(Z, None) - dvardLdK = self._K_dvar * dL_dKmm - var_len3 = self.variance / (np.square(self.lengthscale)*self.lengthscale) - - self.variance.gradient += np.sum(dvardLdK) - self.lengthscale.gradient += (np.square(Z[:,None,:]-Z[None,:,:])*dvardLdK[:,:,None]).reshape(-1,self.input_dim).sum(axis=0)*var_len3 + self.variance.gradient += (dL_dpsi2 * _dpsi2_dvariance).sum() + self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) - - def gradients_Z_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, posterior_variational): - self._psi_computations(Z, posterior_variational.mean, posterior_variational.variance, posterior_variational.binary_prob) + def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + _, _, _, _, _, _dpsi1_dZ, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + _, _, _, _, _, _dpsi2_dZ, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) #psi1 - grad = (dL_dpsi1[:, :, None] * self._dpsi1_dZ).sum(axis=0) + grad = (dL_dpsi1[:, :, None] * _dpsi1_dZ).sum(axis=0) #psi2 - grad += (dL_dpsi2[:, :, :, None] * self._dpsi2_dZ).sum(axis=0).sum(axis=1) - - grad += self.gradients_X(dL_dKmm, Z, None) + grad += (dL_dpsi2[:, :, :, None] * _dpsi2_dZ).sum(axis=0).sum(axis=1) return grad - def gradients_q_variational(self, dL_dKmm, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, posterior_variational): - ndata = posterior_variational.mean.shape[0] - self._psi_computations(Z, posterior_variational.mean, posterior_variational.variance, posterior_variational.binary_prob) + def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): + ndata = variational_posterior.mean.shape[0] + + _, _, _dpsi1_dgamma, _dpsi1_dmu, _dpsi1_dS, _, _ = ssrbf_psi_comp._psi1computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + _, _, _dpsi2_dgamma, _dpsi2_dmu, _dpsi2_dS, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) + #psi1 - grad_mu = (dL_dpsi1[:, :, None] * self._dpsi1_dmu).sum(axis=1) - grad_S = (dL_dpsi1[:, :, None] * self._dpsi1_dS).sum(axis=1) - grad_gamma = (dL_dpsi1[:,:,None] * self._dpsi1_dgamma).sum(axis=1) + grad_mu = (dL_dpsi1[:, :, None] * _dpsi1_dmu).sum(axis=1) + grad_S = (dL_dpsi1[:, :, None] * _dpsi1_dS).sum(axis=1) + grad_gamma = (dL_dpsi1[:,:,None] * _dpsi1_dgamma).sum(axis=1) #psi2 - grad_mu += (dL_dpsi2[:, :, :, None] * self._dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1) - grad_S += (dL_dpsi2[:, :, :, None] * self._dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1) - grad_gamma += (dL_dpsi2[:,:,:, None] *self._dpsi2_dgamma).reshape(ndata,-1,self.input_dim).sum(axis=1) + grad_mu += (dL_dpsi2[:, :, :, None] * _dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1) + grad_S += (dL_dpsi2[:, :, :, None] * _dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1) + grad_gamma += (dL_dpsi2[:,:,:, None] * _dpsi2_dgamma).reshape(ndata,-1,self.input_dim).sum(axis=1) return grad_mu, grad_S, grad_gamma - - def gradients_X(self, dL_dK, X, X2=None): - #if self._X is None or X.base is not self._X.base or X2 is not None: - if X2==None: - _K_dist = X[:,None,:] - X[None,:,:] - _K_dist2 = np.square(_K_dist/self.lengthscale).sum(axis=-1) - dK_dX = self.variance*np.exp(-0.5 * self._K_dist2[:,:,None]) * (-2.*_K_dist/np.square(self.lengthscale)) - dL_dX = (dL_dK[:,:,None] * dK_dX).sum(axis=1) - else: - _K_dist = X[:,None,:] - X2[None,:,:] - _K_dist2 = np.square(_K_dist/self.lengthscale).sum(axis=-1) - dK_dX = self.variance*np.exp(-0.5 * self._K_dist2[:,:,None]) * (-_K_dist/np.square(self.lengthscale)) - dL_dX = (dL_dK[:,:,None] * dK_dX).sum(axis=1) - return dL_dX #---------------------------------------# # Precomputations # @@ -174,78 +137,3 @@ class SSRBF(Stationary): self._K_dist2 = -2.*np.dot(X, X2.T) + (np.sum(np.square(X), axis=1)[:, None] + np.sum(np.square(X2), axis=1)[None, :]) self._K_dvar = np.exp(-0.5 * self._K_dist2) - #@cache_this(1) - def _psi_computations(self, Z, mu, S, gamma): - """ - Z - MxQ - mu - NxQ - S - NxQ - gamma - NxQ - """ - # here are the "statistics" for psi1 and psi2 - # Produced intermediate results: - # _psi1 NxM - # _dpsi1_dvariance NxM - # _dpsi1_dlengthscale NxMxQ - # _dpsi1_dZ NxMxQ - # _dpsi1_dgamma NxMxQ - # _dpsi1_dmu NxMxQ - # _dpsi1_dS NxMxQ - # _psi2 NxMxM - # _psi2_dvariance NxMxM - # _psi2_dlengthscale NxMxMxQ - # _psi2_dZ NxMxMxQ - # _psi2_dgamma NxMxMxQ - # _psi2_dmu NxMxMxQ - # _psi2_dS NxMxMxQ - - lengthscale2 = np.square(self.lengthscale) - - _psi2_Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q - _psi2_Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q - _psi2_Zdist_sq = np.square(_psi2_Zdist / self.lengthscale) # M,M,Q - _psi2_Z_sq_sum = (np.square(Z[:,None,:])+np.square(Z[None,:,:]))/lengthscale2 # MxMxQ - - # psi1 - _psi1_denom = S[:, None, :] / lengthscale2 + 1. # Nx1xQ - _psi1_denom_sqrt = np.sqrt(_psi1_denom) #Nx1xQ - _psi1_dist = Z[None, :, :] - mu[:, None, :] # NxMxQ - _psi1_dist_sq = np.square(_psi1_dist) / (lengthscale2 * _psi1_denom) # NxMxQ - _psi1_common = gamma[:,None,:] / (lengthscale2*_psi1_denom*_psi1_denom_sqrt) #Nx1xQ - _psi1_exponent1 = np.log(gamma[:,None,:]) -0.5 * (_psi1_dist_sq + np.log(_psi1_denom)) # NxMxQ - _psi1_exponent2 = np.log(1.-gamma[:,None,:]) -0.5 * (np.square(Z[None,:,:])/lengthscale2) # NxMxQ - _psi1_exponent = np.log(np.exp(_psi1_exponent1) + np.exp(_psi1_exponent2)) #NxMxQ - _psi1_exp_sum = _psi1_exponent.sum(axis=-1) #NxM - _psi1_exp_dist_sq = np.exp(-0.5*_psi1_dist_sq) # NxMxQ - _psi1_exp_Z = np.exp(-0.5*np.square(Z[None,:,:])/lengthscale2) # 1xMxQ - _psi1_q = self.variance * np.exp(_psi1_exp_sum[:,:,None] - _psi1_exponent) # NxMxQ - self._psi1 = self.variance * np.exp(_psi1_exp_sum) # NxM - self._dpsi1_dvariance = self._psi1 / self.variance # NxM - self._dpsi1_dgamma = _psi1_q * (_psi1_exp_dist_sq/_psi1_denom_sqrt-_psi1_exp_Z) # NxMxQ - self._dpsi1_dmu = _psi1_q * (_psi1_exp_dist_sq * _psi1_dist * _psi1_common) # NxMxQ - self._dpsi1_dS = _psi1_q * (_psi1_exp_dist_sq * _psi1_common * 0.5 * (_psi1_dist_sq - 1.)) # NxMxQ - self._dpsi1_dZ = _psi1_q * (- _psi1_common * _psi1_dist * _psi1_exp_dist_sq - (1-gamma[:,None,:])/lengthscale2*Z[None,:,:]*_psi1_exp_Z) # NxMxQ - self._dpsi1_dlengthscale = 2.*self.lengthscale*_psi1_q * (0.5*_psi1_common*(S[:,None,:]/lengthscale2+_psi1_dist_sq)*_psi1_exp_dist_sq + 0.5*(1-gamma[:,None,:])*np.square(Z[None,:,:]/lengthscale2)*_psi1_exp_Z) # NxMxQ - - - # psi2 - _psi2_denom = 2.*S[:, None, None, :] / lengthscale2 + 1. # Nx1x1xQ - _psi2_denom_sqrt = np.sqrt(_psi2_denom) - _psi2_mudist = mu[:,None,None,:]-_psi2_Zhat #N,M,M,Q - _psi2_mudist_sq = np.square(_psi2_mudist)/(lengthscale2*_psi2_denom) - _psi2_common = gamma[:,None,None,:]/(lengthscale2 * _psi2_denom * _psi2_denom_sqrt) # Nx1x1xQ - _psi2_exponent1 = -_psi2_Zdist_sq -_psi2_mudist_sq -0.5*np.log(_psi2_denom)+np.log(gamma[:,None,None,:]) #N,M,M,Q - _psi2_exponent2 = np.log(1.-gamma[:,None,None,:]) - 0.5*(_psi2_Z_sq_sum) # NxMxMxQ - _psi2_exponent = np.log(np.exp(_psi2_exponent1) + np.exp(_psi2_exponent2)) - _psi2_exp_sum = _psi2_exponent.sum(axis=-1) #NxM - _psi2_q = np.square(self.variance) * np.exp(_psi2_exp_sum[:,:,:,None]-_psi2_exponent) # NxMxMxQ - _psi2_exp_dist_sq = np.exp(-_psi2_Zdist_sq -_psi2_mudist_sq) # NxMxMxQ - _psi2_exp_Z = np.exp(-0.5*_psi2_Z_sq_sum) # MxMxQ - self._psi2 = np.square(self.variance) * np.exp(_psi2_exp_sum) # N,M,M - self._dpsi2_dvariance = 2. * self._psi2/self.variance # NxMxM - self._dpsi2_dgamma = _psi2_q * (_psi2_exp_dist_sq/_psi2_denom_sqrt - _psi2_exp_Z) # NxMxMxQ - self._dpsi2_dmu = _psi2_q * (-2.*_psi2_common*_psi2_mudist * _psi2_exp_dist_sq) # NxMxMxQ - self._dpsi2_dS = _psi2_q * (_psi2_common * (2.*_psi2_mudist_sq - 1.) * _psi2_exp_dist_sq) # NxMxMxQ - self._dpsi2_dZ = 2.*_psi2_q * (_psi2_common*(-_psi2_Zdist*_psi2_denom+_psi2_mudist)*_psi2_exp_dist_sq - (1-gamma[:,None,None,:])*Z[:,None,:]/lengthscale2*_psi2_exp_Z) # NxMxMxQ - self._dpsi2_dlengthscale = 2.*self.lengthscale* _psi2_q * (_psi2_common*(S[:,None,None,:]/lengthscale2+_psi2_Zdist_sq*_psi2_denom+_psi2_mudist_sq)*_psi2_exp_dist_sq+(1-gamma[:,None,None,:])*_psi2_Z_sq_sum*0.5/lengthscale2*_psi2_exp_Z) # NxMxMxQ - \ No newline at end of file diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 18a08e5d..8763426a 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -150,37 +150,6 @@ class BayesianGPLVM(SparseGP): return dim_reduction_plots.plot_steepest_gradient_map(self,*args,**kwargs) -class BayesianGPLVMWithMissingData(BayesianGPLVM): - def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10, - Z=None, kernel=None, inference_method=None, likelihood=None, name='bayesian gplvm', **kwargs): - from ..util.subarray_and_sorting import common_subarrays - self.subarrays = common_subarrays(Y) - import ipdb;ipdb.set_trace() - BayesianGPLVM.__init__(self, Y, input_dim, X=X, X_variance=X_variance, init=init, num_inducing=num_inducing, Z=Z, kernel=kernel, inference_method=inference_method, likelihood=likelihood, name=name, **kwargs) - - - def parameters_changed(self): - super(BayesianGPLVM, self).parameters_changed() - self._log_marginal_likelihood -= self.KL_divergence() - - dL_dmu, dL_dS = self.dL_dmuS() - - # dL: - self.X.mean.gradient = dL_dmu - self.X.variance.gradient = dL_dS - - # dKL: - self.X.mean.gradient -= self.X.mean - self.X.variance.gradient -= (1. - (1. / (self.X.variance))) * 0.5 - -if __name__ == '__main__': - import numpy as np - X = np.random.randn(20,2) - W = np.linspace(0,1,10)[None,:] - Y = (X*W).sum(1) - missing = np.random.binomial(1,.1,size=Y.shape) - - pass def latent_cost_and_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2): """ diff --git a/GPy/models/ss_gplvm.py b/GPy/models/ss_gplvm.py index f21da605..94682c74 100644 --- a/GPy/models/ss_gplvm.py +++ b/GPy/models/ss_gplvm.py @@ -58,7 +58,7 @@ class SSGPLVM(SparseGP): super(SSGPLVM, self).parameters_changed() self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X) - self.X.mean.gradient, self.X.variance.gradient, self.X.binary_prob.gradient = self.kern.gradients_q_variational(posterior_variational=self.X, Z=self.Z, **self.grad_dict) + self.X.mean.gradient, self.X.variance.gradient, self.X.binary_prob.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **self.grad_dict) # update for the KL divergence self.variational_prior.update_gradients_KL(self.X) From 399da015e696d9ca05e5503fa03eb7bd2606bd03 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 08:56:52 +0000 Subject: [PATCH 32/51] rbf --- GPy/kern/_src/rbf.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index cf5ea0c4..0e5108ad 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -18,7 +18,7 @@ class RBF(Stationary): """ - def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='RBF'): + def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='rbf'): super(RBF, self).__init__(input_dim, variance, lengthscale, ARD, name) self.weave_options = {} @@ -200,7 +200,6 @@ class RBF(Stationary): #allocate memory for the things we want to compute mudist = np.empty((N, M, M, Q)) mudist_sq = np.empty((N, M, M, Q)) - exponent = np.zeros((N,M,M)) psi2 = np.empty((N, M, M)) l2 = self.lengthscale **2 From a35464b32f1a5a9b08bddee52ea9639d91ae5cdf Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 08:57:30 +0000 Subject: [PATCH 33/51] WARNING: switched caching off --- GPy/util/caching.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 2899cb33..76d030ca 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -38,6 +38,9 @@ class Cacher(object): if not all([isinstance(arg, Observable) for arg in observable_args]): return self.operation(*args) + # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING + return self.operation(*args) + #if the result is cached, return the cached computation state = [all(a is b for a, b in zip(args, cached_i)) for cached_i in self.cached_inputs] if any(state): From 999d2419dd6d6119dc17406053d24d0365512eb3 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 08:58:34 +0000 Subject: [PATCH 34/51] observer pattern now tested and fully operational. needed the good night rest : ) --- GPy/core/parameterization/array_core.py | 2 +- GPy/core/parameterization/lists_and_dicts.py | 18 ++++ GPy/core/parameterization/param.py | 1 + GPy/core/parameterization/parameter_core.py | 42 ++++---- GPy/core/parameterization/parameterized.py | 7 ++ GPy/testing/observable_tests.py | 108 +++++++++++++++++++ 6 files changed, 158 insertions(+), 20 deletions(-) create mode 100644 GPy/core/parameterization/lists_and_dicts.py create mode 100644 GPy/testing/observable_tests.py diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index a338ceed..208cd4fb 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -62,7 +62,7 @@ class ObservableArray(np.ndarray, Observable): def __setitem__(self, s, val): if self._s_not_empty(s): super(ObservableArray, self).__setitem__(s, val) - self._notify_observers() + self._notify_observers(self[s]) def __getslice__(self, start, stop): return self.__getitem__(slice(start, stop)) diff --git a/GPy/core/parameterization/lists_and_dicts.py b/GPy/core/parameterization/lists_and_dicts.py new file mode 100644 index 00000000..cdf9f5f6 --- /dev/null +++ b/GPy/core/parameterization/lists_and_dicts.py @@ -0,0 +1,18 @@ +''' +Created on 27 Feb 2014 + +@author: maxz +''' + +class ParamList(list): + """ + List to store ndarray-likes in. + It will look for 'is' instead of calling __eq__ on each element. + """ + def __contains__(self, other): + for el in self: + if el is other: + return True + return False + + pass diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 89d3a4e4..ca9905f7 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -172,6 +172,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable): try: new_arr._current_slice_ = s; new_arr._original_ = self.base is new_arr.base except AttributeError: pass # returning 0d array or float, double etc return new_arr + def __setitem__(self, s, val): super(Param, self).__setitem__(s, val) if self.has_parent(): diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 6afa94cb..58dd63d8 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -2,6 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED +import heapq __updated__ = '2013-12-16' @@ -11,25 +12,29 @@ def adjust_name_for_printing(name): return '' class Observable(object): + _updated = True def __init__(self, *args, **kwargs): - from collections import defaultdict - self._observer_callables_ = defaultdict(list) - - def add_observer(self, observer, callble): - self._observer_callables_[observer].append(callble) + self._observer_callables_ = [] + def add_observer(self, observer, callble, priority=0): + heapq.heappush(self._observer_callables_, (priority, observer, callble)) + def remove_observer(self, observer, callble=None): - if observer in self._observer_callables_: - if callble is None: - del self._observer_callables_[observer] - elif callble in self._observer_callables_[observer]: - self._observer_callables_[observer].remove(callble) - if len(self._observer_callables_[observer]) == 0: - self.remove_observer(observer) - - def _notify_observers(self): - [[callble(self) for callble in callables] - for callables in self._observer_callables_.itervalues()] + to_remove = [] + for p, obs, clble in self._observer_callables_: + if callble is not None: + if (obs == observer) and (callble == clble): + to_remove.append((p, obs, clble)) + else: + if obs is observer: + to_remove.append((p, obs, clble)) + for r in to_remove: + self._observer_callables_.remove(r) + + def _notify_observers(self, which=None): + if which is None: + which = self + [callble(which) for _, _, callble in heapq.nlargest(len(self._observer_callables_), self._observer_callables_)] class Pickleable(object): def _getstate(self): @@ -333,7 +338,7 @@ class Constrainable(Nameable, Indexable): class Parameterizable(Constrainable, Observable): def __init__(self, *args, **kwargs): super(Parameterizable, self).__init__(*args, **kwargs) - from GPy.core.parameterization.array_core import ParamList + from GPy.core.parameterization.lists_and_dicts import ParamList _parameters_ = ParamList() self._added_names_ = set() @@ -398,7 +403,7 @@ class Parameterizable(Constrainable, Observable): """Returns a (deep) copy of the current model""" import copy from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView - from .array_core import ParamList + from .lists_and_dicts import ParamList dc = dict() for k, v in self.__dict__.iteritems(): @@ -427,7 +432,6 @@ class Parameterizable(Constrainable, Observable): def _notify_parameters_changed(self): self.parameters_changed() - self._notify_observers() if self.has_parent(): self._direct_parent_._notify_parameters_changed() diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index f5fcc6ad..fe8c76e4 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -116,6 +116,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) + param.add_observer(self, self._pass_through_notify, -1) self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" @@ -169,6 +170,12 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self._param_slices_.append(slice(sizes[-2], sizes[-1])) self._add_parameter_name(p) + #=========================================================================== + # notification system + #=========================================================================== + def _pass_through_notify(self, which): + self._notify_observers(which) + #=========================================================================== # Pickling operations #=========================================================================== diff --git a/GPy/testing/observable_tests.py b/GPy/testing/observable_tests.py new file mode 100644 index 00000000..214a678f --- /dev/null +++ b/GPy/testing/observable_tests.py @@ -0,0 +1,108 @@ +''' +Created on 27 Feb 2014 + +@author: maxz +''' +import unittest +from GPy.core.parameterization.parameterized import Parameterized +from GPy.core.parameterization.param import Param +import numpy + + +class ParamTestParent(Parameterized): + parent_changed_count = 0 + def parameters_changed(self): + self.parent_changed_count += 1 + +class ParameterizedTest(Parameterized): + params_changed_count = 0 + def parameters_changed(self): + self.params_changed_count += 1 + +class Test(unittest.TestCase): + + def setUp(self): + self.parent = ParamTestParent('test parent') + self.par = ParameterizedTest('test model') + self.p = Param('test parameter', numpy.random.normal(1,2,(10,3))) + + self.par.add_parameter(self.p) + self.parent.add_parameter(self.par) + + self._observer_triggered = None + self._trigger_count = 0 + self._first = None + self._second = None + + def _trigger(self, which): + self._observer_triggered = float(which) + self._trigger_count += 1 + if self._first is not None: + self._second = self._trigger + else: + self._first = self._trigger + + def _trigger_priority(self, which): + if self._first is not None: + self._second = self._trigger_priority + else: + self._first = self._trigger_priority + + def test_observable(self): + self.par.add_observer(self, self._trigger, -1) + self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.p[0,1] = 3 # trigger observers + self.assertEqual(self._observer_triggered, 3, 'observer should have triggered') + self.assertEqual(self._trigger_count, 1, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 1, 'params changed once') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.par.remove_observer(self) + self.p[2,1] = 4 + self.assertEqual(self._observer_triggered, 3, 'observer should not have triggered') + self.assertEqual(self._trigger_count, 1, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 2, 'params changed second') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.par.add_observer(self, self._trigger, -1) + self.p[2,1] = 4 + self.assertEqual(self._observer_triggered, 4, 'observer should have triggered') + self.assertEqual(self._trigger_count, 2, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 3, 'params changed second') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.par.remove_observer(self, self._trigger) + self.p[0,1] = 3 + self.assertEqual(self._observer_triggered, 4, 'observer should not have triggered') + self.assertEqual(self._trigger_count, 2, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 4, 'params changed second') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + def test_set_params(self): + self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') + self.par._set_params(numpy.ones(self.par.size)) + self.assertEqual(self.par.params_changed_count, 1, 'now params changed') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + def test_priority(self): + self.par.add_observer(self, self._trigger, -1) + self.par.add_observer(self, self._trigger_priority, 0) + self.par._notify_observers(0) + self.assertEqual(self._first, self._trigger_priority, 'priority should be first') + self.assertEqual(self._second, self._trigger, 'priority should be first') + + self.par.remove_observer(self) + self._first = self._second = None + + self.par.add_observer(self, self._trigger, 1) + self.par.add_observer(self, self._trigger_priority, 0) + self.par._notify_observers(0) + self.assertEqual(self._first, self._trigger, 'priority should be second') + self.assertEqual(self._second, self._trigger_priority, 'priority should be second') + + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'Test.testName'] + unittest.main() \ No newline at end of file From 7ae9d03c4514396ecb3dd2c7e21be31684354c91 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 28 Feb 2014 11:01:54 +0000 Subject: [PATCH 35/51] efficiencies in stationary --- GPy/kern/_src/rbf.py | 4 --- GPy/kern/_src/stationary.py | 51 +++++++++++++++++++++++++++++-------- 2 files changed, 41 insertions(+), 14 deletions(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 7bf0adeb..d817b765 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -311,7 +311,3 @@ class RBF(Stationary): type_converters=weave.converters.blitz, **self.weave_options) return denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 - - def input_sensitivity(self): - if self.ARD: return 1./self.lengthscale - else: return (1./self.lengthscale).repeat(self.input_dim) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 8d8ae476..ae4cd879 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -12,6 +12,35 @@ from scipy import integrate from ...util.caching import Cache_this class Stationary(Kern): + """ + Stationary kernels (covaraince functions). + + Stationary covariance fucntion depend only on r, where r is defined as + + r = \sqrt{ \sum_{q=1}^Q (x_q - x'_q)^2 } + + The covaraince function k(x, x' can then be written k(r). + + In this implementation, r is scaled by the lengthscales parameter(s): + + r = \sqrt{ \sum_{q=1}^Q \frac{(x_q - x'_q)^2}{\ell_q^2} }. + + By default, there's only one lengthscale: seaprate lengthscales for each + dimension can be enables by setting ARD=True. + + To implement a stationary covaraince function using this class, one need + only define the covariance function k(r), and it derivative. + + ... + def K_of_r(self, r): + return foo + def dK_dr(self, r): + return bar + + The lengthscale(s) and variance parameters are added to the structure automatically. + + """ + def __init__(self, input_dim, variance, lengthscale, ARD, name): super(Stationary, self).__init__(input_dim, name) self.ARD = ARD @@ -20,11 +49,11 @@ class Stationary(Kern): lengthscale = np.ones(1) else: lengthscale = np.asarray(lengthscale) - assert lengthscale.size == 1, "Only lengthscale needed for non-ARD kernel" + assert lengthscale.size == 1, "Only 1 lengthscale needed for non-ARD kernel" else: if lengthscale is not None: lengthscale = np.asarray(lengthscale) - assert lengthscale.size in [1, input_dim], "Bad lengthscales" + assert lengthscale.size in [1, input_dim], "Bad number of lengthscales" if lengthscale.size != input_dim: lengthscale = np.ones(input_dim)*lengthscale else: @@ -35,10 +64,10 @@ class Stationary(Kern): self.add_parameters(self.variance, self.lengthscale) def K_of_r(self, r): - raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" + raise NotImplementedError, "implement the covariance function as a fn of r to use this class" def dK_dr(self, r): - raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" + raise NotImplementedError, "implement derivative of the covariance function wrt r to use this class" #@Cache_this(limit=5, ignore_args=()) def K(self, X, X2=None): @@ -84,7 +113,6 @@ class Stationary(Kern): else: return self._unscaled_dist(X, X2)/self.lengthscale - def Kdiag(self, X): ret = np.empty(X.shape[0]) ret[:] = self.variance @@ -98,15 +126,18 @@ class Stationary(Kern): r = self._scaled_dist(X, X2) K = self.K_of_r(r) - rinv = self._inv_dist(X, X2) dL_dr = self.dK_dr(r) * dL_dK if self.ARD: - x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 - self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0) + #rinv = self._inv_dist(X, X2) + #x_xl3 = np.square(self._dist(X, X2)) # TODO: this is rather high memory? Should we loop instead? + #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 + self.lengthscale.gradient = np.zeros(self.input_dim) + tmp = dL_dr*self._inv_dist(X, X2) + if X2 is None: X2 = X + [np.copyto(self.lengthscale.gradient[q:q+1], -np.sum(tmp * np.square(X[:,q:q+1] - X2[:,q:q+1].T))/self.lengthscale[q]**3) for q in xrange(self.input_dim)] else: - x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 - self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum() + self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale self.variance.gradient = np.sum(K * dL_dK)/self.variance From 2771e3f71f53a421a153d8aab81c57d21c755aee Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 11:33:31 +0000 Subject: [PATCH 36/51] observer pattern has a handle to trigger only > min_priority observers --- GPy/core/model.py | 14 -- GPy/core/parameterization/array_core.py | 1 + GPy/core/parameterization/param.py | 16 ++- GPy/core/parameterization/parameter_core.py | 136 ++++++++++++++----- GPy/core/parameterization/parameterized.py | 39 ++---- GPy/core/parameterization/priors.py | 30 ++++ GPy/core/parameterization/transformations.py | 8 +- GPy/testing/observable_tests.py | 27 +++- 8 files changed, 181 insertions(+), 90 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 6514d73a..0925a199 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -60,20 +60,6 @@ class Model(Parameterized): self.priors = state.pop() Parameterized._setstate(self, state) - def randomize(self): - """ - Randomize the model. - Make this draw from the prior if one exists, else draw from N(0,1) - """ - # first take care of all parameters (from N(0,1)) - # x = self._get_params_transformed() - x = np.random.randn(self.size_transformed) - x = self._untransform_params(x) - # now draw from prior where possible - [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None] - self._set_params(x) - # self._set_params_transformed(self._get_params_transformed()) # makes sure all of the tied parameters get the same init (since there's only one prior object...) - def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs): """ Perform random restarts of the model, and set the model to the best diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 208cd4fb..e9e5ca8c 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -66,6 +66,7 @@ class ObservableArray(np.ndarray, Observable): def __getslice__(self, start, stop): return self.__getitem__(slice(start, stop)) + def __setslice__(self, start, stop, val): return self.__setitem__(slice(start, stop), val) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index ca9905f7..14cba600 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -3,8 +3,8 @@ import itertools import numpy -from parameter_core import Constrainable, Gradcheckable, Indexable, Parentable, adjust_name_for_printing -from array_core import ObservableArray, ParamList +from parameter_core import OptimizationHandlable, Gradcheckable, adjust_name_for_printing +from array_core import ObservableArray ###### printing __constraints_name__ = "Constraint" @@ -15,7 +15,7 @@ __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision __print_threshold__ = 5 ###### -class Param(Constrainable, ObservableArray, Gradcheckable): +class Param(OptimizationHandlable, ObservableArray, Gradcheckable): """ Parameter object for GPy models. @@ -148,8 +148,11 @@ class Param(Constrainable, ObservableArray, Gradcheckable): #=========================================================================== # get/set parameters #=========================================================================== - def _set_params(self, param, update=True): + def _set_params(self, param, trigger_parent=True): self.flat = param + if trigger_parent: min_priority = None + else: min_priority = -numpy.inf + self._notify_observers(None, min_priority) def _get_params(self): return self.flat @@ -175,9 +178,6 @@ class Param(Constrainable, ObservableArray, Gradcheckable): def __setitem__(self, s, val): super(Param, self).__setitem__(s, val) - if self.has_parent(): - self._direct_parent_._notify_parameters_changed() - #self._notify_observers() #=========================================================================== # Index Operations: @@ -205,6 +205,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable): ind = self._indices(slice_index) if ind.ndim < 2: ind = ind[:, None] return numpy.asarray(numpy.apply_along_axis(lambda x: numpy.sum(extended_realshape * x), 1, ind), dtype=int) + def _expand_index(self, slice_index=None): # this calculates the full indexing arrays from the slicing objects given by get_item for _real..._ attributes # it basically translates slices to their respective index arrays and turns negative indices around @@ -346,6 +347,7 @@ class ParamConcatenation(object): See :py:class:`GPy.core.parameter.Param` for more details on constraining. """ # self.params = params + from lists_and_dicts import ParamList self.params = ParamList([]) for p in params: for p in p.flattened_parameters: diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 58dd63d8..b4483b4d 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -31,10 +31,24 @@ class Observable(object): for r in to_remove: self._observer_callables_.remove(r) - def _notify_observers(self, which=None): + def _notify_observers(self, which=None, min_priority=None): + """ + Notifies all observers. Which is the element, which kicked off this + notification loop. + + NOTE: notifies only observers with priority p > min_priority! + ^^^^^^^^^^^^^^^^ + + :param which: object, which started this notification loop + :param min_priority: only notify observers with priority > min_priority + if min_priority is None, notify all observers in order + """ if which is None: which = self - [callble(which) for _, _, callble in heapq.nlargest(len(self._observer_callables_), self._observer_callables_)] + if min_priority is None: + [callble(which) for _, _, callble in heapq.nlargest(len(self._observer_callables_), self._observer_callables_)] + else: + [callble(which) for p, _, callble in heapq.nlargest(len(self._observer_callables_), self._observer_callables_) if p > min_priority] class Pickleable(object): def _getstate(self): @@ -210,9 +224,9 @@ class Constrainable(Nameable, Indexable): #=========================================================================== # Prior Operations #=========================================================================== - def set_prior(self, prior, warning=True, update=True): + def set_prior(self, prior, warning=True): repriorized = self.unset_priors() - self._add_to_index_operations(self.priors, repriorized, prior, warning, update) + self._add_to_index_operations(self.priors, repriorized, prior, warning) def unset_priors(self, *priors): return self._remove_from_index_operations(self.priors, priors) @@ -238,7 +252,7 @@ class Constrainable(Nameable, Indexable): # Constrain operations -> done #=========================================================================== - def constrain(self, transform, warning=True, update=True): + def constrain(self, transform, warning=True): """ :param transform: the :py:class:`GPy.core.transformations.Transformation` to constrain the this parameter to. @@ -248,9 +262,9 @@ class Constrainable(Nameable, Indexable): :py:class:`GPy.core.transformations.Transformation`. """ if isinstance(transform, Transformation): - self._set_params(transform.initialize(self._get_params()), update=False) + self._set_params(transform.initialize(self._get_params()), trigger_parent=True) reconstrained = self.unconstrain() - self._add_to_index_operations(self.constraints, reconstrained, transform, warning, update) + self._add_to_index_operations(self.constraints, reconstrained, transform, warning) def unconstrain(self, *transforms): """ @@ -261,30 +275,30 @@ class Constrainable(Nameable, Indexable): """ return self._remove_from_index_operations(self.constraints, transforms) - def constrain_positive(self, warning=True, update=True): + def constrain_positive(self, warning=True): """ :param warning: print a warning if re-constraining parameters. Constrain this parameter to the default positive constraint. """ - self.constrain(Logexp(), warning=warning, update=update) + self.constrain(Logexp(), warning=warning) - def constrain_negative(self, warning=True, update=True): + def constrain_negative(self, warning=True): """ :param warning: print a warning if re-constraining parameters. Constrain this parameter to the default negative constraint. """ - self.constrain(NegativeLogexp(), warning=warning, update=update) + self.constrain(NegativeLogexp(), warning=warning) - def constrain_bounded(self, lower, upper, warning=True, update=True): + def constrain_bounded(self, lower, upper, warning=True): """ :param lower, upper: the limits to bound this parameter to :param warning: print a warning if re-constraining parameters. Constrain this parameter to lie within the given range. """ - self.constrain(Logistic(lower, upper), warning=warning, update=update) + self.constrain(Logistic(lower, upper), warning=warning) def unconstrain_positive(self): """ @@ -314,12 +328,10 @@ class Constrainable(Nameable, Indexable): for p in self._parameters_: p._parent_changed(parent) - def _add_to_index_operations(self, which, reconstrained, transform, warning, update): + def _add_to_index_operations(self, which, reconstrained, transform, warning): if warning and reconstrained.size > 0: print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) which.add(transform, self._raveled_index()) - if update: - self._notify_observers() def _remove_from_index_operations(self, which, transforms): if len(transforms) == 0: @@ -334,8 +346,69 @@ class Constrainable(Nameable, Indexable): return removed +class OptimizationHandlable(Constrainable, Observable): + def _get_params_transformed(self): + # transformed parameters (apply transformation rules) + p = self._get_params() + [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + if self._has_fixes(): + return p[self._fixes_] + return p + + def _set_params_transformed(self, p): + # inverse apply transformations for parameters and set the resulting parameters + self._set_params(self._untransform_params(p)) + + def _untransform_params(self, p): + p = p.copy() + if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp + [np.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + return p + + def _get_params(self): + # don't overwrite this anymore! + if not self.size: + return np.empty(shape=(0,), dtype=np.float64) + return np.hstack([x._get_params() for x in self._parameters_ if x.size > 0]) -class Parameterizable(Constrainable, Observable): + def _set_params(self, params, trigger_parent=True): + # don't overwrite this anymore! + raise NotImplementedError, "This needs to be implemented seperately" + + #=========================================================================== + # Optimization handles: + #=========================================================================== + def _get_param_names(self): + n = np.array([p.hirarchy_name() + '[' + str(i) + ']' for p in self.flattened_parameters for i in p._indices()]) + return n + def _get_param_names_transformed(self): + n = self._get_param_names() + if self._has_fixes(): + return n[self._fixes_] + return n + + #=========================================================================== + # Randomizeable + #=========================================================================== + def randomize(self): + """ + Randomize the model. + Make this draw from the prior if one exists, else draw from N(0,1) + """ + import numpy as np + # first take care of all parameters (from N(0,1)) + # x = self._get_params_transformed() + x = np.random.randn(self.size_transformed) + x = self._untransform_params(x) + # now draw from prior where possible + [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None] + self._set_params(x) + # self._set_params_transformed(self._get_params_transformed()) # makes sure all of the tied parameters get the same init (since there's only one prior object...) + + +import numpy as np + +class Parameterizable(OptimizationHandlable): def __init__(self, *args, **kwargs): super(Parameterizable, self).__init__(*args, **kwargs) from GPy.core.parameterization.lists_and_dicts import ParamList @@ -382,23 +455,21 @@ class Parameterizable(Constrainable, Observable): import itertools [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + def _set_params(self, params, trigger_parent=True): + import itertools + [p._set_params(params[s], trigger_parent=False) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + if trigger_parent: min_priority = None + else: min_priority = -np.inf + self._notify_observers(None, min_priority) + def _set_gradient(self, g): import itertools [p._set_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - def _get_params(self): - import numpy as np - # don't overwrite this anymore! - if not self.size: - return np.empty(shape=(0,), dtype=np.float64) - return np.hstack([x._get_params() for x in self._parameters_ if x.size > 0]) - - def _set_params(self, params, update=True): - # don't overwrite this anymore! - import itertools - [p._set_params(params[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - self._notify_parameters_changed() - + + #=========================================================================== + # TODO: not working yet + #=========================================================================== def copy(self): """Returns a (deep) copy of the current model""" import copy @@ -429,11 +500,6 @@ class Parameterizable(Constrainable, Observable): s.add_parameter(p) return s - - def _notify_parameters_changed(self): - self.parameters_changed() - if self.has_parent(): - self._direct_parent_._notify_parameters_changed() def parameters_changed(self): """ diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index fe8c76e4..0093c6f3 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -58,6 +58,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self._in_init_ = True self._parameters_ = ParamList() self.size = sum(p.size for p in self._parameters_) + self.add_observer(self, self._parameters_changed_notification, -100) if not self._has_fixes(): self._fixes_ = None self._param_slices_ = [] @@ -65,7 +66,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): del self._in_init_ def build_pydot(self, G=None): - import pydot + import pydot # @UnresolvedImport iamroot = False if G is None: G = pydot.Dot(graph_type='digraph') @@ -116,7 +117,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) - param.add_observer(self, self._pass_through_notify, -1) + param.add_observer(self, self._pass_through_notify_observers, -np.inf) self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" @@ -173,9 +174,10 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): #=========================================================================== # notification system #=========================================================================== - def _pass_through_notify(self, which): + def _parameters_changed_notification(self, which): + self.parameters_changed() + def _pass_through_notify_observers(self, which): self._notify_observers(which) - #=========================================================================== # Pickling operations #=========================================================================== @@ -244,32 +246,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): g[self._offset_for(p) + numpy.array(list(i))] += g[self._raveled_index_for(t)] if self._has_fixes(): return g[self._fixes_] return g - #=========================================================================== - # Optimization handles: - #=========================================================================== - def _get_param_names(self): - n = numpy.array([p.hirarchy_name() + '[' + str(i) + ']' for p in self.flattened_parameters for i in p._indices()]) - return n - def _get_param_names_transformed(self): - n = self._get_param_names() - if self._has_fixes(): - return n[self._fixes_] - return n - def _get_params_transformed(self): - # transformed parameters (apply transformation rules) - p = self._get_params() - [numpy.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - if self._has_fixes(): - return p[self._fixes_] - return p - def _set_params_transformed(self, p): - # inverse apply transformations for parameters and set the resulting parameters - self._set_params(self._untransform_params(p)) - def _untransform_params(self, p): - p = p.copy() - if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp - [numpy.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - return p + #=========================================================================== # Indexable Handling #=========================================================================== @@ -304,6 +281,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): this is not in the global view of things! """ return numpy.r_[:self.size] + #=========================================================================== # Fixing parameters: #=========================================================================== @@ -311,6 +289,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): if self._has_fixes(): return self._fixes_[self._raveled_index_for(param)] return numpy.ones(self.size, dtype=bool)[self._raveled_index_for(param)] + #=========================================================================== # Convenience for fixed, tied checking of param: #=========================================================================== diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py index 906fe003..29adc923 100644 --- a/GPy/core/parameterization/priors.py +++ b/GPy/core/parameterization/priors.py @@ -64,6 +64,36 @@ class Gaussian(Prior): return np.random.randn(n) * self.sigma + self.mu +class Uniform(Prior): + domain = _REAL + _instances = [] + def __new__(cls, lower, upper): # Singleton: + if cls._instances: + cls._instances[:] = [instance for instance in cls._instances if instance()] + for instance in cls._instances: + if instance().lower == lower and instance().upper == upper: + return instance() + o = super(Prior, cls).__new__(cls, lower, upper) + cls._instances.append(weakref.ref(o)) + return cls._instances[-1]() + + def __init__(self, lower, upper): + self.lower = float(lower) + self.upper = float(upper) + + def __str__(self): + return "[" + str(np.round(self.lower)) + ', ' + str(np.round(self.upper)) + ']' + + def lnpdf(self, x): + region = (x>=self.lower) * (x<=self.upper) + return region + + def lnpdf_grad(self, x): + return np.zeros(x.shape) + + def rvs(self, n): + return np.random.uniform(self.lower, self.upper, size=n) + class LogGaussian(Prior): """ Implementation of the univariate *log*-Gaussian probability function, coupled with random variables. diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index 36291ca3..60fcc469 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -6,8 +6,11 @@ import numpy as np from domains import _POSITIVE,_NEGATIVE, _BOUNDED import weakref +import sys +#_lim_val = -np.log(sys.float_info.epsilon) + _exp_lim_val = np.finfo(np.float64).max -_lim_val = np.log(_exp_lim_val)#-np.log(sys.float_info.epsilon) +_lim_val = np.log(_exp_lim_val)# #=============================================================================== # Fixing constants @@ -35,7 +38,6 @@ class Transformation(object): """ produce a sensible initial value for f(x)""" raise NotImplementedError def plot(self, xlabel=r'transformed $\theta$', ylabel=r'$\theta$', axes=None, *args,**kw): - import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." import matplotlib.pyplot as plt from ...plotting.matplot_dep import base_plots @@ -52,7 +54,7 @@ class Transformation(object): class Logexp(Transformation): domain = _POSITIVE def f(self, x): - return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -np.inf, _lim_val)))) + return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -_lim_val, _lim_val)))) #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x))) def finv(self, f): return np.where(f>_lim_val, f, np.log(np.exp(f) - 1.)) diff --git a/GPy/testing/observable_tests.py b/GPy/testing/observable_tests.py index 214a678f..6b4f1a87 100644 --- a/GPy/testing/observable_tests.py +++ b/GPy/testing/observable_tests.py @@ -18,16 +18,26 @@ class ParameterizedTest(Parameterized): params_changed_count = 0 def parameters_changed(self): self.params_changed_count += 1 + def _set_params(self, params, trigger_parent=True): + Parameterized._set_params(self, params, trigger_parent=trigger_parent) class Test(unittest.TestCase): def setUp(self): self.parent = ParamTestParent('test parent') self.par = ParameterizedTest('test model') + self.par2 = ParameterizedTest('test model 2') self.p = Param('test parameter', numpy.random.normal(1,2,(10,3))) self.par.add_parameter(self.p) + self.par.add_parameter(Param('test1', numpy.random.normal(0,1,(1,)))) + self.par.add_parameter(Param('test2', numpy.random.normal(0,1,(1,)))) + + self.par2.add_parameter(Param('par2 test1', numpy.random.normal(0,1,(1,)))) + self.par2.add_parameter(Param('par2 test2', numpy.random.normal(0,1,(1,)))) + self.parent.add_parameter(self.par) + self.parent.add_parameter(self.par2) self._observer_triggered = None self._trigger_count = 0 @@ -84,7 +94,22 @@ class Test(unittest.TestCase): self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') self.par._set_params(numpy.ones(self.par.size)) self.assertEqual(self.par.params_changed_count, 1, 'now params changed') - self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) + + self.parent._set_params(numpy.ones(self.parent.size) * 2) + self.assertEqual(self.par.params_changed_count, 2, 'now params changed') + self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) + + + def test_priority_notify(self): + self.assertEqual(self.par.params_changed_count, 0) + self.par._notify_observers(0, None) + self.assertEqual(self.par.params_changed_count, 1) + self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) + + self.par._notify_observers(0, -numpy.inf) + self.assertEqual(self.par.params_changed_count, 2) + self.assertEqual(self.parent.parent_changed_count, 1) def test_priority(self): self.par.add_observer(self, self._trigger, -1) From 1d1123fcae9bd99ff60bf158287f5f138fe61559 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 11:39:05 +0000 Subject: [PATCH 37/51] plotting with uncertain inputs --- GPy/plotting/matplot_dep/models_plots.py | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index d72d2a3e..4ca4441e 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -56,10 +56,13 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - - X, Y = param_to_array(model.X, model.Y) - if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): X_variance = model.X_variance - + + if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): + X = model.X.mean + X_variance = param_to_array(model.X.variance) + else: + X = model.X + X, Y = param_to_array(X, model.Y) if hasattr(model, 'Z'): Z = param_to_array(model.Z) #work out what the inputs are for plotting (1D or 2D) @@ -98,10 +101,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #add error bars for uncertain (if input uncertainty is being modelled) - #if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): - # ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), - # xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()), - # ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): + ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), + xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()), + ecolor='k', fmt=None, elinewidth=.5, alpha=.5) #set the limits of the plot to some sensible values From 83f5b9377ac95d97f1b26441567542cfdb2ac8be Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 11:46:09 +0000 Subject: [PATCH 38/51] caching in place again and working : ) --- GPy/inference/latent_function_inference/var_dtc.py | 4 ++-- GPy/kern/_src/rbf.py | 4 ++-- GPy/kern/_src/stationary.py | 10 +++++----- GPy/util/caching.py | 2 +- 4 files changed, 10 insertions(+), 10 deletions(-) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index fec61204..2b7ca7ad 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -60,7 +60,7 @@ class VarDTC(object): _, output_dim = Y.shape #see whether we've got a different noise variance for each datum - beta = 1./np.squeeze(likelihood.variance) + beta = 1./max(1e-6, np.squeeze(likelihood.variance)) # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! #self.YYTfactor = self.get_YYTfactor(Y) @@ -214,7 +214,7 @@ class VarDTCMissingData(object): psi2_all = None Ys, traces = self._Y(Y) - beta_all = 1./likelihood.variance + beta_all = 1./max(1e-6, likelihood.variance) het_noise = beta_all.size != 1 import itertools diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index d4a60077..38022bd4 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -159,7 +159,7 @@ class RBF(Stationary): grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) #psi2 - denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + denom, _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) tmp = psi2[:, :, :, None] / l2 / denom grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) @@ -237,7 +237,7 @@ class RBF(Stationary): return denom, dist, dist_sq, psi1 - #@cache_this(ignore_args=(1,)) + @Cache_this(limit=1, ignore_args=(0,)) def _Z_distances(self, Z): Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index ae4cd879..bc51d850 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -69,18 +69,18 @@ class Stationary(Kern): def dK_dr(self, r): raise NotImplementedError, "implement derivative of the covariance function wrt r to use this class" - #@Cache_this(limit=5, ignore_args=()) + @Cache_this(limit=5, ignore_args=()) def K(self, X, X2=None): r = self._scaled_dist(X, X2) return self.K_of_r(r) - #@Cache_this(limit=5, ignore_args=(0,)) + @Cache_this(limit=5, ignore_args=(0,)) def _dist(self, X, X2): if X2 is None: X2 = X return X[:, None, :] - X2[None, :, :] - #@Cache_this(limit=5, ignore_args=(0,)) + @Cache_this(limit=5, ignore_args=(0,)) def _unscaled_dist(self, X, X2=None): """ Compute the square distance between each row of X and X2, or between @@ -94,7 +94,7 @@ class Stationary(Kern): X2sq = np.sum(np.square(X2),1) return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) - #@Cache_this(limit=5, ignore_args=()) + @Cache_this(limit=5, ignore_args=()) def _scaled_dist(self, X, X2=None): """ Efficiently compute the scaled distance, r. @@ -147,7 +147,7 @@ class Stationary(Kern): diagonal, where we return zero (the distance on the diagonal is zero). This term appears in derviatives. """ - dist = self._scaled_dist(X, X2) + dist = self._scaled_dist(X, X2).copy() if X2 is None: nondiag = util.diag.offdiag_view(dist) nondiag[:] = 1./nondiag diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 76d030ca..a2017407 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -39,7 +39,7 @@ class Cacher(object): return self.operation(*args) # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING - return self.operation(*args) + # return self.operation(*args) #if the result is cached, return the cached computation state = [all(a is b for a, b in zip(args, cached_i)) for cached_i in self.cached_inputs] From 7aad39e70e0f786ae8a1b9b0690b06747c593e95 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 28 Feb 2014 12:06:28 +0000 Subject: [PATCH 39/51] non essential tidying in stationary --- GPy/kern/_src/stationary.py | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index ae4cd879..b2868772 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -74,16 +74,10 @@ class Stationary(Kern): r = self._scaled_dist(X, X2) return self.K_of_r(r) - #@Cache_this(limit=5, ignore_args=(0,)) - def _dist(self, X, X2): - if X2 is None: - X2 = X - return X[:, None, :] - X2[None, :, :] - #@Cache_this(limit=5, ignore_args=(0,)) def _unscaled_dist(self, X, X2=None): """ - Compute the square distance between each row of X and X2, or between + Compute the Euclidean distance between each row of X and X2, or between each pair of rows of X if X2 is None. """ if X2 is None: @@ -99,7 +93,7 @@ class Stationary(Kern): """ Efficiently compute the scaled distance, r. - r = \sum_{q=1}^Q (x_q - x'q)^2/l_q^2 + r = \sqrt( \sum_{q=1}^Q (x_q - x'q)^2/l_q^2 ) Note that if thre is only one lengthscale, l comes outside the sum. In this case we compute the unscaled distance first (in a separate @@ -129,10 +123,10 @@ class Stationary(Kern): dL_dr = self.dK_dr(r) * dL_dK if self.ARD: - #rinv = self._inv_dist(X, X2) - #x_xl3 = np.square(self._dist(X, X2)) # TODO: this is rather high memory? Should we loop instead? + #rinv = self._inv_dis# this is rather high memory? Should we loop instead?t(X, X2) + #d = X[:, None, :] - X2[None, :, :] + #x_xl3 = np.square(d) #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 - self.lengthscale.gradient = np.zeros(self.input_dim) tmp = dL_dr*self._inv_dist(X, X2) if X2 is None: X2 = X [np.copyto(self.lengthscale.gradient[q:q+1], -np.sum(tmp * np.square(X[:,q:q+1] - X2[:,q:q+1].T))/self.lengthscale[q]**3) for q in xrange(self.input_dim)] @@ -162,7 +156,9 @@ class Stationary(Kern): r = self._scaled_dist(X, X2) invdist = self._inv_dist(X, X2) dL_dr = self.dK_dr(r) * dL_dK - #The high-memory numpy way: ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2 + #The high-memory numpy way: + #d = X[:, None, :] - X2[None, :, :] + #ret = np.sum((invdist*dL_dr)[:,:,None]*d,1)/self.lengthscale**2 #if X2 is None: #ret *= 2. @@ -245,7 +241,7 @@ class Matern52(Stationary): .. math:: - k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } + k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) """ def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat52'): super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, name) @@ -256,7 +252,7 @@ class Matern52(Stationary): def dK_dr(self, r): return self.variance*(10./3*r -5.*r -5.*np.sqrt(5.)/3*r**2)*np.exp(-np.sqrt(5.)*r) - def Gram_matrix(self,F,F1,F2,F3,lower,upper): + def Gram_matrix(self, F, F1, F2, F3, lower, upper): """ Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1. From 82a25d691b8f614b3cbccee1d2e0d95b0a5afd50 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 12:08:16 +0000 Subject: [PATCH 40/51] fixed caching bug with args having Nones --- GPy/util/caching.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index a2017407..250efe11 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -1,4 +1,5 @@ from ..core.parameterization.parameter_core import Observable +import itertools class Cacher(object): """ @@ -40,9 +41,9 @@ class Cacher(object): # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING # return self.operation(*args) - + #if the result is cached, return the cached computation - state = [all(a is b for a, b in zip(args, cached_i)) for cached_i in self.cached_inputs] + state = [all(a is b for a, b in itertools.izip_longest(args, cached_i)) for cached_i in self.cached_inputs] if any(state): i = state.index(True) if self.inputs_changed[i]: From af50fa3e57e969e6d2571e78aca8743a631aefb3 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 12:09:03 +0000 Subject: [PATCH 41/51] prediction code need updating, started with woodbury vector, but how to predict variance in sparse gp with uncertain inputs? --- GPy/core/sparse_gp.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 751a295d..b06ffbc7 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -85,11 +85,11 @@ class SparseGP(GP): self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) - def _raw_predict(self, Xnew, X_variance_new=None, full_cov=False): + def _raw_predict(self, Xnew, full_cov=False): """ Make a prediction for the latent function values """ - if X_variance_new is None: + if not isinstance(Xnew, VariationalPosterior): Kx = self.kern.K(self.Z, Xnew) mu = np.dot(Kx.T, self.posterior.woodbury_vector) if full_cov: @@ -100,13 +100,13 @@ class SparseGP(GP): Kxx = self.kern.Kdiag(Xnew) var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T else: - Kx = self.kern.psi1(self.Z, Xnew, X_variance_new) - mu = np.dot(Kx, self.Cpsi1V) + Kx = self.kern.psi1(self.Z, Xnew) + mu = np.dot(Kx, self.posterior.woodbury_vector) if full_cov: raise NotImplementedError, "TODO" else: - Kxx = self.kern.psi0(self.Z, Xnew, X_variance_new) - psi2 = self.kern.psi2(self.Z, Xnew, X_variance_new) + Kxx = self.kern.psi0(self.Z, Xnew) + psi2 = self.kern.psi2(self.Z, Xnew) var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1) return mu, var From ab7dff9a3db8cc7e047de251535791f8e4a755ca Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 28 Feb 2014 12:59:49 +0000 Subject: [PATCH 42/51] no longer caching denom in psi2_rbf --- GPy/kern/_src/rbf.py | 83 ++++++++------------------------------------ 1 file changed, 14 insertions(+), 69 deletions(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 38022bd4..88f88761 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -48,7 +48,7 @@ class RBF(Stationary): if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): psi2, _, _, _, _, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) else: - _, _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) + _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) return psi2 def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): @@ -80,19 +80,16 @@ class RBF(Stationary): denom, _, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) d_length = psi1[:,:,None] * ((dist_sq - 1.)/(self.lengthscale*denom) +1./self.lengthscale) dpsi1_dlength = d_length * dL_dpsi1[:, :, None] - if not self.ARD: - self.lengthscale.gradient += dpsi1_dlength.sum() - else: + if self.ARD: self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) + else: + self.lengthscale.gradient += dpsi1_dlength.sum() self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance #from psi2 S = variational_posterior.variance - denom, _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * denom + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom) - #TODO: combine denom and l2 as denom_l2?? - #TODO: tidy the above! - #TODO: tensordot below? + _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2) / (2.*S[:,None,None,:] + l2)*self.lengthscale dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] if not self.ARD: @@ -125,9 +122,11 @@ class RBF(Stationary): grad = np.sum(dL_dpsi1[:, :, None] * dpsi1_dZ, 0) #psi2 - denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) term1 = Zdist / l2 # M, M, Q - term2 = mudist / denom / l2 # N, M, M, Q + S = variational_posterior.variance + term2 = mudist / (2.*S[:,None,None,:] + l2) # N, M, M, Q + dZ = psi2[:, :, :, None] * (term1[None, :, :, :] + term2) #N,M,M,Q grad += 2*(dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0) @@ -159,8 +158,9 @@ class RBF(Stationary): grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) #psi2 - denom, _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - tmp = psi2[:, :, :, None] / l2 / denom + _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + S = variational_posterior.variance + tmp = psi2[:, :, :, None] / (2.*S[:,None,None,:] + l2) grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) @@ -170,61 +170,6 @@ class RBF(Stationary): # Precomputations # #---------------------------------------# - #TODO: this function is unused, but it will be useful in the stationary class - def _dL_dlengthscales_via_K(self, dL_dK, X, X2): - """ - A helper function for update_gradients_* methods - - Computes the derivative of the objective L wrt the lengthscales via - - dL_dl = sum_{i,j}(dL_dK_{ij} dK_dl) - - assumes self._K_computations has just been called. - - This is only valid if self.ARD=True - """ - target = np.zeros(self.input_dim) - dvardLdK = self._K_dvar * dL_dK - var_len3 = self.variance / np.power(self.lengthscale, 3) - if X2 is None: - # save computation for the symmetrical case - dvardLdK = dvardLdK + dvardLdK.T - code = """ - int q,i,j; - double tmp; - for(q=0; q Date: Fri, 28 Feb 2014 13:37:26 +0000 Subject: [PATCH 43/51] einsumming in stationary --- GPy/kern/_src/stationary.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index a88904da..5f88c3b0 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -13,13 +13,13 @@ from ...util.caching import Cache_this class Stationary(Kern): """ - Stationary kernels (covaraince functions). + Stationary kernels (covariance functions). Stationary covariance fucntion depend only on r, where r is defined as r = \sqrt{ \sum_{q=1}^Q (x_q - x'_q)^2 } - The covaraince function k(x, x' can then be written k(r). + The covariance function k(x, x' can then be written k(r). In this implementation, r is scaled by the lengthscales parameter(s): @@ -28,7 +28,7 @@ class Stationary(Kern): By default, there's only one lengthscale: seaprate lengthscales for each dimension can be enables by setting ARD=True. - To implement a stationary covaraince function using this class, one need + To implement a stationary covariance function using this class, one need only define the covariance function k(r), and it derivative. ... @@ -74,6 +74,11 @@ class Stationary(Kern): r = self._scaled_dist(X, X2) return self.K_of_r(r) + @Cache_this(limit=3, ignore_args=()) + def dK_dr_via_X(self, X, X2): + #a convenience function, so we can cache dK_dr + return self.dK_dr(self._scaled_dist(X, X2)) + @Cache_this(limit=5, ignore_args=(0,)) def _unscaled_dist(self, X, X2=None): """ @@ -117,11 +122,11 @@ class Stationary(Kern): self.lengthscale.gradient = 0. def update_gradients_full(self, dL_dK, X, X2=None): - r = self._scaled_dist(X, X2) - K = self.K_of_r(r) - dL_dr = self.dK_dr(r) * dL_dK + self.variance.gradient = np.einsum('ij,ij,i', self.K(X, X2), dL_dK, 1./self.variance) + #now the lengthscale gradient(s) + dL_dr = self.dK_dr_via_X(X, X2) * dL_dK if self.ARD: #rinv = self._inv_dis# this is rather high memory? Should we loop instead?t(X, X2) #d = X[:, None, :] - X2[None, :, :] @@ -129,11 +134,11 @@ 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 - [np.copyto(self.lengthscale.gradient[q:q+1], -np.sum(tmp * np.square(X[:,q:q+1] - X2[:,q:q+1].T))/self.lengthscale[q]**3) for q in xrange(self.input_dim)] + self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)]) else: + r = self._scaled_dist(X, X2) self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale - self.variance.gradient = np.sum(K * dL_dK)/self.variance def _inv_dist(self, X, X2=None): """ @@ -153,9 +158,8 @@ class Stationary(Kern): """ Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X """ - r = self._scaled_dist(X, X2) invdist = self._inv_dist(X, X2) - dL_dr = self.dK_dr(r) * dL_dK + dL_dr = self.dK_dr_via_X(X, X2) * dL_dK #The high-memory numpy way: #d = X[:, None, :] - X2[None, :, :] #ret = np.sum((invdist*dL_dr)[:,:,None]*d,1)/self.lengthscale**2 @@ -168,7 +172,7 @@ class Stationary(Kern): tmp *= 2. X2 = X ret = np.empty(X.shape, dtype=np.float64) - [np.copyto(ret[:,q], np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), 1)) for q in xrange(self.input_dim)] + [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] ret /= self.lengthscale**2 return ret From c87bda9e49623808a8ad236740ac2d34ae75bab0 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 28 Feb 2014 14:20:17 +0000 Subject: [PATCH 44/51] einsumming in rbf for speed --- GPy/kern/_src/rbf.py | 115 ++++++++++++++++++++++++------------------- 1 file changed, 63 insertions(+), 52 deletions(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 88f88761..baa5b932 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -70,34 +70,39 @@ class RBF(Stationary): self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) return - l2 = self.lengthscale **2 + elif isinstance(variational_posterior, variational.NormalPosterior): + + l2 = self.lengthscale **2 - #contributions from psi0: - self.variance.gradient = np.sum(dL_dpsi0) - self.lengthscale.gradient = 0. + #contributions from psi0: + self.variance.gradient = np.sum(dL_dpsi0) + self.lengthscale.gradient = 0. + + #from psi1 + denom, _, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + d_length = psi1[:,:,None] * ((dist_sq - 1.)/(self.lengthscale*denom) +1./self.lengthscale) + dpsi1_dlength = d_length * dL_dpsi1[:, :, None] + if self.ARD: + self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) + else: + self.lengthscale.gradient += dpsi1_dlength.sum() + self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance + + #from psi2 + S = variational_posterior.variance + _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2) / (2.*S[:,None,None,:] + l2)*self.lengthscale + + dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] + if not self.ARD: + self.lengthscale.gradient += dpsi2_dlength.sum() + else: + self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) + + self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance - #from psi1 - denom, _, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) - d_length = psi1[:,:,None] * ((dist_sq - 1.)/(self.lengthscale*denom) +1./self.lengthscale) - dpsi1_dlength = d_length * dL_dpsi1[:, :, None] - if self.ARD: - self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) else: - self.lengthscale.gradient += dpsi1_dlength.sum() - self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance - - #from psi2 - S = variational_posterior.variance - _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2) / (2.*S[:,None,None,:] + l2)*self.lengthscale - - dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] - if not self.ARD: - self.lengthscale.gradient += dpsi2_dlength.sum() - else: - self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) - - self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance + raise ValueError, "unknown distriubtion received for psi-statistics" def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): # Spike-and-Slab GPLVM @@ -112,25 +117,26 @@ class RBF(Stationary): grad += (dL_dpsi2[:, :, :, None] * _dpsi2_dZ).sum(axis=0).sum(axis=1) return grad - - l2 = self.lengthscale **2 - #psi1 - denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) - denominator = l2 * denom - dpsi1_dZ = -psi1[:, :, None] * (dist / denominator) - grad = np.sum(dL_dpsi1[:, :, None] * dpsi1_dZ, 0) + elif isinstance(variational_posterior, variational.NormalPosterior): + + l2 = self.lengthscale **2 - #psi2 - Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - term1 = Zdist / l2 # M, M, Q - S = variational_posterior.variance - term2 = mudist / (2.*S[:,None,None,:] + l2) # N, M, M, Q + #psi1 + denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + grad = np.einsum('ij,ij,ijk,ijk->jk', dL_dpsi1, psi1, dist, -1./(denom*l2)) - dZ = psi2[:, :, :, None] * (term1[None, :, :, :] + term2) #N,M,M,Q - grad += 2*(dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0) + #psi2 + Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + term1 = Zdist / l2 # M, M, Q + S = variational_posterior.variance + term2 = mudist / (2.*S[:,None,None,:] + l2) # N, M, M, Q - return grad + grad += 2.*np.einsum('ijk,ijk,ijkl->kl', dL_dpsi2, psi2, term1[None,:,:,:] + term2) + + return grad + else: + raise ValueError, "unknown distriubtion received for psi-statistics" def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): # Spike-and-Slab GPLVM @@ -150,19 +156,24 @@ class RBF(Stationary): grad_gamma += (dL_dpsi2[:,:,:, None] * _dpsi2_dgamma).reshape(ndata,-1,self.input_dim).sum(axis=1) return grad_mu, grad_S, grad_gamma + + elif isinstance(variational_posterior, variational.NormalPosterior): - l2 = self.lengthscale **2 - #psi1 - denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) - tmp = psi1[:, :, None] / l2 / denom - grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) - grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) - #psi2 - _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - S = variational_posterior.variance - tmp = psi2[:, :, :, None] / (2.*S[:,None,None,:] + l2) - grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) - grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) + l2 = self.lengthscale **2 + #psi1 + denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + tmp = psi1[:, :, None] / l2 / denom + grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) + grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) + #psi2 + _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + S = variational_posterior.variance + tmp = psi2[:, :, :, None] / (2.*S[:,None,None,:] + l2) + grad_mu += -2.*np.einsum('ijk,ijkl,ijkl->il', dL_dpsi2, tmp , mudist) + grad_S += np.einsum('ijk,ijkl,ijkl->il', dL_dpsi2 , tmp , (2.*mudist_sq - 1)) + + else: + raise ValueError, "unknown distriubtion received for psi-statistics" return grad_mu, grad_S From 47e4026141b0712777eda3713b150f43d2756c11 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 16:18:47 +0000 Subject: [PATCH 45/51] hierarchy edits. adding removing parameters withing hierarchy --- GPy/core/model.py | 4 +- GPy/core/parameterization/index_operations.py | 12 +++--- GPy/core/parameterization/param.py | 43 +++++++++++++------ GPy/core/parameterization/parameter_core.py | 42 ++++++++++-------- GPy/core/parameterization/parameterized.py | 35 ++++++++++----- GPy/examples/dimensionality_reduction.py | 10 ++--- .../latent_function_inference/var_dtc.py | 5 +-- GPy/kern/_src/kern.py | 10 +++-- GPy/models/sparse_gp_regression.py | 4 +- GPy/testing/parameterized_tests.py | 3 +- GPy/util/caching.py | 2 +- 11 files changed, 106 insertions(+), 64 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 0925a199..6fd80d76 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -265,7 +265,7 @@ class Model(Parameterized): and numerical gradients is within of unity. """ - x = self._get_params_transformed().copy() + x = self._get_params_transformed() if not verbose: # make sure only to test the selected parameters @@ -283,7 +283,7 @@ class Model(Parameterized): return # just check the global ratio - dx = np.zeros_like(x) + dx = np.zeros(x.shape()) dx[transformed_index] = step * np.sign(np.random.uniform(-1, 1, transformed_index.size)) # evaulate around the point x diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index b5399741..6450c41c 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -194,9 +194,13 @@ class ParameterIndexOperationsView(object): def shift_right(self, start, size): - raise NotImplementedError, 'Shifting only supported in original ParamIndexOperations' - + self._param_index_ops.shift_right(start+self._offset, size) + def shift_left(self, start, size): + self._param_index_ops.shift_left(start+self._offset, size) + self._offset -= size + self._size -= size + def clear(self): for i, ind in self.items(): self._param_index_ops.remove(i, ind+self._offset) @@ -232,9 +236,7 @@ class ParameterIndexOperationsView(object): def __getitem__(self, prop): ind = self._filter_index(self._param_index_ops[prop]) - if ind.size > 0: - return ind - raise KeyError, prop + return ind def __str__(self, *args, **kwargs): import pprint diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 14cba600..d52442d1 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -269,7 +269,7 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): return [t._short() for t in self._tied_to_] or [''] def __repr__(self, *args, **kwargs): name = "\033[1m{x:s}\033[0;0m:\n".format( - x=self.hirarchy_name()) + x=self.hierarchy_name()) return name + super(Param, self).__repr__(*args, **kwargs) def _ties_for(self, rav_index): # size = sum(p.size for p in self._tied_to_) @@ -303,12 +303,12 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): gen = map(lambda x: " ".join(map(str, x)), gen) return reduce(lambda a, b:max(a, len(b)), gen, len(header)) def _max_len_values(self): - return reduce(lambda a, b:max(a, len("{x:=.{0}g}".format(__precision__, x=b))), self.flat, len(self.hirarchy_name())) + return reduce(lambda a, b:max(a, len("{x:=.{0}g}".format(__precision__, x=b))), self.flat, len(self.hierarchy_name())) def _max_len_index(self, ind): return reduce(lambda a, b:max(a, len(str(b))), ind, len(__index_name__)) def _short(self): # short string to print - name = self.hirarchy_name() + name = self.hierarchy_name() if self._realsize_ < 2: return name ind = self._indices() @@ -331,8 +331,8 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): if lp is None: lp = self._max_len_names(prirs, __tie_name__) sep = '-' header_format = " {i:{5}^{2}s} | \033[1m{x:{5}^{1}s}\033[0;0m | {c:{5}^{0}s} | {p:{5}^{4}s} | {t:{5}^{3}s}" - if only_name: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hirarchy_name(), c=sep*lc, i=sep*li, t=sep*lt, p=sep*lp) # nice header for printing - else: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hirarchy_name(), c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # nice header for printing + if only_name: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hierarchy_name(), c=sep*lc, i=sep*li, t=sep*lt, p=sep*lp) # nice header for printing + else: header = header_format.format(lc, lx, li, lt, lp, ' ', 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] + [" {i!s:^{3}s} | {x: >{1}.{2}g} | {c:^{0}s} | {p:^{5}s} | {t:^{4}s} ".format(lc, lx, __precision__, li, lt, lp, 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)]) # return all the constraints with right indices # except: return super(Param, self).__str__() @@ -356,6 +356,21 @@ class ParamConcatenation(object): self._param_sizes = [p.size for p in self.params] startstops = numpy.cumsum([0] + self._param_sizes) self._param_slices_ = [slice(start, stop) for start,stop in zip(startstops, startstops[1:])] + + parents = dict() + for p in self.params: + if p.has_parent(): + parent = p._direct_parent_ + level = 0 + while parent is not None: + if parent in parents: + parents[parent] = max(level, parents[parent]) + else: + parents[parent] = level + level += 1 + parent = parent._direct_parent_ + import operator + self.parents = map(lambda x: x[0], sorted(parents.iteritems(), key=operator.itemgetter(1))) #=========================================================================== # Get/set items, enable broadcasting #=========================================================================== @@ -369,24 +384,26 @@ class ParamConcatenation(object): val = val._vals() ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; vals = self._vals(); vals[s] = val; del val - [numpy.place(p, ind[ps], vals[ps]) and update and p._notify_observers() + [numpy.place(p, ind[ps], vals[ps]) for p, ps in zip(self.params, self._param_slices_)] + if update: + self.update_all_params() def _vals(self): return numpy.hstack([p._get_params() for p in self.params]) #=========================================================================== # parameter operations: #=========================================================================== def update_all_params(self): - for p in self.params: - p._notify_observers() - + for par in self.parents: + par._notify_observers(-numpy.inf) + def constrain(self, constraint, warning=True): - [param.constrain(constraint, update=False) for param in self.params] + [param.constrain(constraint, trigger_parent=False) for param in self.params] self.update_all_params() constrain.__doc__ = Param.constrain.__doc__ def constrain_positive(self, warning=True): - [param.constrain_positive(warning, update=False) for param in self.params] + [param.constrain_positive(warning, trigger_parent=False) for param in self.params] self.update_all_params() constrain_positive.__doc__ = Param.constrain_positive.__doc__ @@ -396,12 +413,12 @@ class ParamConcatenation(object): fix = constrain_fixed def constrain_negative(self, warning=True): - [param.constrain_negative(warning, update=False) for param in self.params] + [param.constrain_negative(warning, trigger_parent=False) for param in self.params] self.update_all_params() constrain_negative.__doc__ = Param.constrain_negative.__doc__ def constrain_bounded(self, lower, upper, warning=True): - [param.constrain_bounded(lower, upper, warning, update=False) for param in self.params] + [param.constrain_bounded(lower, upper, warning, trigger_parent=False) for param in self.params] self.update_all_params() constrain_bounded.__doc__ = Param.constrain_bounded.__doc__ diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index b4483b4d..4b1b16e0 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -6,6 +6,11 @@ import heapq __updated__ = '2013-12-16' +class HierarchyError(Exception): + """ + Gets thrown when something is wrong with the parameter hierarchy + """ + def adjust_name_for_printing(name): if name is not None: return name.replace(" ", "_").replace(".", "_").replace("-", "").replace("+", "").replace("!", "").replace("*", "").replace("/", "") @@ -114,11 +119,11 @@ class Nameable(Parentable): self._name = name if self.has_parent(): self._direct_parent_._name_changed(self, from_name) - def hirarchy_name(self, adjust_for_printing=True): + def hierarchy_name(self, adjust_for_printing=True): if adjust_for_printing: adjust = lambda x: adjust_name_for_printing(x) else: adjust = lambda x: x if self.has_parent(): - return self._direct_parent_.hirarchy_name() + "." + adjust(self.name) + return self._direct_parent_.hierarchy_name() + "." + adjust(self.name) return adjust(self.name) @@ -175,7 +180,7 @@ class Constrainable(Nameable, Indexable): #=========================================================================== # Fixing Parameters: #=========================================================================== - def constrain_fixed(self, value=None, warning=True): + def constrain_fixed(self, value=None, warning=True, trigger_parent=True): """ Constrain this paramter to be fixed to the current value it carries. @@ -183,7 +188,7 @@ class Constrainable(Nameable, Indexable): """ if value is not None: self[:] = value - self.constrain(__fixed__, warning=warning) + self.constrain(__fixed__, warning=warning, trigger_parent=trigger_parent) rav_i = self._highest_parent_._raveled_index_for(self) self._highest_parent_._set_fixed(rav_i) fix = constrain_fixed @@ -224,7 +229,7 @@ class Constrainable(Nameable, Indexable): #=========================================================================== # Prior Operations #=========================================================================== - def set_prior(self, prior, warning=True): + def set_prior(self, prior, warning=True, trigger_parent=True): repriorized = self.unset_priors() self._add_to_index_operations(self.priors, repriorized, prior, warning) @@ -252,7 +257,7 @@ class Constrainable(Nameable, Indexable): # Constrain operations -> done #=========================================================================== - def constrain(self, transform, warning=True): + def constrain(self, transform, warning=True, trigger_parent=True): """ :param transform: the :py:class:`GPy.core.transformations.Transformation` to constrain the this parameter to. @@ -262,7 +267,7 @@ class Constrainable(Nameable, Indexable): :py:class:`GPy.core.transformations.Transformation`. """ if isinstance(transform, Transformation): - self._set_params(transform.initialize(self._get_params()), trigger_parent=True) + self._set_params(transform.initialize(self._get_params()), trigger_parent=trigger_parent) reconstrained = self.unconstrain() self._add_to_index_operations(self.constraints, reconstrained, transform, warning) @@ -275,30 +280,30 @@ class Constrainable(Nameable, Indexable): """ return self._remove_from_index_operations(self.constraints, transforms) - def constrain_positive(self, warning=True): + def constrain_positive(self, warning=True, trigger_parent=True): """ :param warning: print a warning if re-constraining parameters. Constrain this parameter to the default positive constraint. """ - self.constrain(Logexp(), warning=warning) + self.constrain(Logexp(), warning=warning, trigger_parent=trigger_parent) - def constrain_negative(self, warning=True): + def constrain_negative(self, warning=True, trigger_parent=True): """ :param warning: print a warning if re-constraining parameters. Constrain this parameter to the default negative constraint. """ - self.constrain(NegativeLogexp(), warning=warning) + self.constrain(NegativeLogexp(), warning=warning, trigger_parent=trigger_parent) - def constrain_bounded(self, lower, upper, warning=True): + def constrain_bounded(self, lower, upper, warning=True, trigger_parent=True): """ :param lower, upper: the limits to bound this parameter to :param warning: print a warning if re-constraining parameters. Constrain this parameter to lie within the given range. """ - self.constrain(Logistic(lower, upper), warning=warning) + self.constrain(Logistic(lower, upper), warning=warning, trigger_parent=trigger_parent) def unconstrain_positive(self): """ @@ -359,6 +364,9 @@ class OptimizationHandlable(Constrainable, Observable): # inverse apply transformations for parameters and set the resulting parameters self._set_params(self._untransform_params(p)) + def _size_transformed(self): + return self.size - self.constraints[__fixed__].size + def _untransform_params(self, p): p = p.copy() if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp @@ -373,13 +381,13 @@ class OptimizationHandlable(Constrainable, Observable): def _set_params(self, params, trigger_parent=True): # don't overwrite this anymore! - raise NotImplementedError, "This needs to be implemented seperately" + raise NotImplementedError, "This needs to be implemented in Param and Parametrizable" #=========================================================================== # Optimization handles: #=========================================================================== def _get_param_names(self): - n = np.array([p.hirarchy_name() + '[' + str(i) + ']' for p in self.flattened_parameters for i in p._indices()]) + n = np.array([p.hierarchy_name() + '[' + str(i) + ']' for p in self.flattened_parameters for i in p._indices()]) return n def _get_param_names_transformed(self): n = self._get_param_names() @@ -398,7 +406,7 @@ class OptimizationHandlable(Constrainable, Observable): import numpy as np # first take care of all parameters (from N(0,1)) # x = self._get_params_transformed() - x = np.random.randn(self.size_transformed) + x = np.random.randn(self._size_transformed()) x = self._untransform_params(x) # now draw from prior where possible [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None] @@ -435,7 +443,7 @@ class Parameterizable(OptimizationHandlable): if pname in self._added_names_: del self.__dict__[pname] self._add_parameter_name(param) - else: + elif pname not in dir(self): self.__dict__[pname] = param self._added_names_.add(pname) diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 0093c6f3..56d785c3 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -7,7 +7,7 @@ import cPickle import itertools from re import compile, _pattern_type from param import ParamConcatenation -from parameter_core import Constrainable, Pickleable, Parentable, Observable, Parameterizable, adjust_name_for_printing, Gradcheckable +from parameter_core import Pickleable, Parameterizable, adjust_name_for_printing, Gradcheckable from transformations import __fixed__ from array_core import ParamList @@ -105,6 +105,14 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.remove_parameter(param) self.add_parameter(param, index) elif param not in self._parameters_: + if param.has_parent(): + parent = param._direct_parent_ + while parent is not None: + if parent is self: + from GPy.core.parameterization.parameter_core import HierarchyError + raise HierarchyError, "You cannot add a parameter twice into the hirarchy" + parent = parent._direct_parent_ + param._direct_parent_.remove_parameter(param) # make sure the size is set if index is None: self.constraints.update(param.constraints, self.size) @@ -117,13 +125,16 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) + param.add_observer(self, self._pass_through_notify_observers, -np.inf) + self.size += param.size + + self._connect_parameters() + self._notify_parent_change() + self._connect_fixes() else: raise RuntimeError, """Parameter exists already added and no copy made""" - self._connect_parameters() - self._notify_parent_change() - self._connect_fixes() def add_parameters(self, *parameters): @@ -146,12 +157,19 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): del self._parameters_[param._parent_index_] param._disconnect_parent() - param.remove_observer(self, self._notify_parameters_changed) + param.remove_observer(self, self._pass_through_notify_observers) self.constraints.shift_left(start, param.size) + self._connect_fixes() self._connect_parameters() self._notify_parent_change() + parent = self._direct_parent_ + while parent is not None: + parent._connect_fixes() + parent._connect_parameters() + parent._notify_parent_change() + parent = parent._direct_parent_ def _connect_parameters(self): # connect parameterlist to this parameterized object @@ -351,7 +369,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): # Printing: #=========================================================================== def _short(self): - return self.hirarchy_name() + return self.hierarchy_name() @property def flattened_parameters(self): return [xi for x in self._parameters_ for xi in x.flattened_parameters] @@ -359,11 +377,6 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): def _parameter_sizes_(self): return [x.size for x in self._parameters_] @property - def size_transformed(self): - if self._has_fixes(): - return sum(self._fixes_) - return self.size - @property def parameter_shapes(self): return [xi for x in self._parameters_ for xi in x.parameter_shapes] @property diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 2044f08d..9ebb54a2 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -187,10 +187,10 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False): _np.random.seed(1234) x = _np.linspace(0, 4 * _np.pi, N)[:, None] - s1 = _np.vectorize(lambda x: -_np.sin(_np.exp(x))) + s1 = _np.vectorize(lambda x: _np.sin(x)) s2 = _np.vectorize(lambda x: _np.cos(x)**2) s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x))) - sS = _np.vectorize(lambda x: x*_np.sin(x)) + sS = _np.vectorize(lambda x: _np.cos(x)) s1 = s1(x) s2 = s2(x) @@ -202,7 +202,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False): s3 -= s3.mean(); s3 /= s3.std(0) sS -= sS.mean(); sS /= sS.std(0) - S1 = _np.hstack([s1, s2, sS]) + S1 = _np.hstack([s1, sS]) S2 = _np.hstack([s2, s3, sS]) S3 = _np.hstack([s3, sS]) @@ -270,7 +270,7 @@ def bgplvm_simulation(optimize=True, verbose=1, from GPy import kern from GPy.models import BayesianGPLVM - D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 5, 9 + D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) @@ -294,7 +294,7 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, from GPy.models import BayesianGPLVM from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData - D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 5, 9 + D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 7, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 2b7ca7ad..64707298 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -60,8 +60,7 @@ class VarDTC(object): _, output_dim = Y.shape #see whether we've got a different noise variance for each datum - beta = 1./max(1e-6, np.squeeze(likelihood.variance)) - + beta = 1./np.fmax(likelihood.variance, 1e-6) # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! #self.YYTfactor = self.get_YYTfactor(Y) #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) @@ -214,7 +213,7 @@ class VarDTCMissingData(object): psi2_all = None Ys, traces = self._Y(Y) - beta_all = 1./max(1e-6, likelihood.variance) + beta_all = 1./np.fmax(likelihood.variance, 1e-6) het_noise = beta_all.size != 1 import itertools diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index eb3291e0..14e6ae49 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -112,10 +112,12 @@ class Kern(Parameterized): """ assert isinstance(other, Kern), "only kernels can be added to kernels..." from add import Add - return Add([self, other], tensor) - - def __call__(self, X, X2=None): - return self.K(X, X2) + kernels = [] + if not tensor and isinstance(self, Add): kernels.extend(self._parameters_) + else: kernels.append(self) + if not tensor and isinstance(other, Add): kernels.extend(other._parameters_) + else: kernels.append(other) + return Add(kernels, tensor) def __mul__(self, other): """ Here we overload the '*' operator. See self.prod for more information""" diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py index 6a76df3f..4980d26a 100644 --- a/GPy/models/sparse_gp_regression.py +++ b/GPy/models/sparse_gp_regression.py @@ -8,7 +8,7 @@ from .. import likelihoods from .. import kern from ..inference.latent_function_inference import VarDTC from ..util.misc import param_to_array -from ..core.parameterization.variational import VariationalPosterior +from ..core.parameterization.variational import NormalPosterior class SparseGPRegression(SparseGP): """ @@ -47,7 +47,7 @@ class SparseGPRegression(SparseGP): likelihood = likelihoods.Gaussian() if not (X_variance is None): - X = VariationalPosterior(X,X_variance) + X = NormalPosterior(X,X_variance) SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC()) diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index 6f13d294..0da3f3ae 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -6,6 +6,7 @@ Created on Feb 13, 2014 import unittest import GPy import numpy as np +from GPy.core.parameterization.parameter_core import HierarchyError class Test(unittest.TestCase): @@ -65,7 +66,7 @@ class Test(unittest.TestCase): self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1]) def test_add_parameter_already_in_hirarchy(self): - self.test1.add_parameter(self.white._parameters_[0]) + self.assertRaises(HierarchyError, self.test1.add_parameter, self.white._parameters_[0]) def test_default_constraints(self): self.assertIs(self.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 250efe11..a2434c80 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -40,7 +40,7 @@ class Cacher(object): return self.operation(*args) # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING - # return self.operation(*args) + return self.operation(*args) #if the result is cached, return the cached computation state = [all(a is b for a, b in itertools.izip_longest(args, cached_i)) for cached_i in self.cached_inputs] From eae3c28dc04d5b1d18046dfda437b16f20049152 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 16:20:52 +0000 Subject: [PATCH 46/51] checkgrad --- 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 6fd80d76..0e3913c8 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -283,7 +283,7 @@ class Model(Parameterized): return # just check the global ratio - dx = np.zeros(x.shape()) + dx = np.zeros(x.shape) dx[transformed_index] = step * np.sign(np.random.uniform(-1, 1, transformed_index.size)) # evaulate around the point x From 20e02e63a976a7c2cf50a1f80ca3a45983f9ec3c Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 16:58:24 +0000 Subject: [PATCH 47/51] global gradient test done and some parameterized fixes --- GPy/core/model.py | 8 +++- GPy/core/parameterization/array_core.py | 13 ------ GPy/core/parameterization/index_operations.py | 42 +------------------ GPy/core/parameterization/lists_and_dicts.py | 19 ++++++++- GPy/core/parameterization/param.py | 9 ++-- GPy/core/parameterization/parameter_core.py | 9 ++-- GPy/core/parameterization/parameterized.py | 32 ++------------ 7 files changed, 39 insertions(+), 93 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 0e3913c8..d27cbc69 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -226,6 +226,11 @@ class Model(Parameterized): TODO: valid args """ + if self.is_fixed: + raise RuntimeError, "Cannot optimize, when everything is fixed" + if self.size == 0: + raise RuntimeError, "Model without parameters cannot be minimized" + if optimizer is None: optimizer = self.preferred_optimizer @@ -294,9 +299,8 @@ class Model(Parameterized): dx = dx[transformed_index] gradient = gradient[transformed_index] - numerical_gradient = (f1 - f2) / (2 * dx) global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) - return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gradient - numerical_gradient).mean() < tolerance) + return (np.abs(1. - global_ratio) < tolerance) else: # check the gradient of each parameter individually, and do some pretty printing try: diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index e9e5ca8c..cf353ead 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -6,19 +6,6 @@ __updated__ = '2013-12-16' import numpy as np from parameter_core import Observable -class ParamList(list): - """ - List to store ndarray-likes in. - It will look for 'is' instead of calling __eq__ on each element. - """ - def __contains__(self, other): - for el in self: - if el is other: - return True - return False - - pass - class ObservableArray(np.ndarray, Observable): """ An ndarray which reports changes to its observers. diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index 6450c41c..a9f3768e 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -5,47 +5,7 @@ Created on Oct 2, 2013 ''' import numpy from numpy.lib.function_base import vectorize -from param import Param -from collections import defaultdict - -class ParamDict(defaultdict): - def __init__(self): - """ - Default will be self._default, if not set otherwise - """ - defaultdict.__init__(self, self.default_factory) - - def __getitem__(self, key): - try: - return defaultdict.__getitem__(self, key) - except KeyError: - for a in self.iterkeys(): - if numpy.all(a==key) and a._parent_index_==key._parent_index_: - return defaultdict.__getitem__(self, a) - raise - - def __contains__(self, key): - if defaultdict.__contains__(self, key): - return True - for a in self.iterkeys(): - if numpy.all(a==key) and a._parent_index_==key._parent_index_: - return True - return False - - def __setitem__(self, key, value): - if isinstance(key, Param): - for a in self.iterkeys(): - if numpy.all(a==key) and a._parent_index_==key._parent_index_: - return super(ParamDict, self).__setitem__(a, value) - defaultdict.__setitem__(self, key, value) - -class SetDict(ParamDict): - def default_factory(self): - return set() - -class IntArrayDict(ParamDict): - def default_factory(self): - return numpy.int_([]) +from lists_and_dicts import IntArrayDict class ParameterIndexOperations(object): ''' diff --git a/GPy/core/parameterization/lists_and_dicts.py b/GPy/core/parameterization/lists_and_dicts.py index cdf9f5f6..5b13b3b5 100644 --- a/GPy/core/parameterization/lists_and_dicts.py +++ b/GPy/core/parameterization/lists_and_dicts.py @@ -4,7 +4,24 @@ Created on 27 Feb 2014 @author: maxz ''' -class ParamList(list): +from collections import defaultdict +class DefaultArrayDict(defaultdict): + def __init__(self): + """ + Default will be self._default, if not set otherwise + """ + defaultdict.__init__(self, self.default_factory) + +class SetDict(DefaultArrayDict): + def default_factory(self): + return set() + +class IntArrayDict(DefaultArrayDict): + def default_factory(self): + import numpy as np + return np.int_([]) + +class ArrayList(list): """ List to store ndarray-likes in. It will look for 'is' instead of calling __eq__ on each element. diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index d52442d1..22610a70 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -50,7 +50,7 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): obj._realsize_ = obj.size obj._realndim_ = obj.ndim obj._updated_ = False - from index_operations import SetDict + from lists_and_dicts import SetDict obj._tied_to_me_ = SetDict() obj._tied_to_ = [] obj._original_ = True @@ -232,7 +232,8 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): #=========================================================================== @property def is_fixed(self): - return self._highest_parent_._is_fixed(self) + from transformations import __fixed__ + return self.constraints[__fixed__].size == self.size #def round(self, decimals=0, out=None): # view = super(Param, self).round(decimals, out).view(Param) # view.__array_finalize__(self) @@ -347,8 +348,8 @@ class ParamConcatenation(object): See :py:class:`GPy.core.parameter.Param` for more details on constraining. """ # self.params = params - from lists_and_dicts import ParamList - self.params = ParamList([]) + from lists_and_dicts import ArrayList + self.params = ArrayList([]) for p in params: for p in p.flattened_parameters: if p not in self.params: diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 4b1b16e0..e7344fa5 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -335,6 +335,7 @@ class Constrainable(Nameable, Indexable): def _add_to_index_operations(self, which, reconstrained, transform, warning): if warning and reconstrained.size > 0: + # TODO: figure out which parameters have changed and only print those print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) which.add(transform, self._raveled_index()) @@ -419,8 +420,8 @@ import numpy as np class Parameterizable(OptimizationHandlable): def __init__(self, *args, **kwargs): super(Parameterizable, self).__init__(*args, **kwargs) - from GPy.core.parameterization.lists_and_dicts import ParamList - _parameters_ = ParamList() + from GPy.core.parameterization.lists_and_dicts import ArrayList + _parameters_ = ArrayList() self._added_names_ = set() def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): @@ -482,7 +483,7 @@ class Parameterizable(OptimizationHandlable): """Returns a (deep) copy of the current model""" import copy from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView - from .lists_and_dicts import ParamList + from .lists_and_dicts import ArrayList dc = dict() for k, v in self.__dict__.iteritems(): @@ -496,7 +497,7 @@ class Parameterizable(OptimizationHandlable): dc['_direct_parent_'] = None dc['_parent_index_'] = None - dc['_parameters_'] = ParamList() + dc['_parameters_'] = ArrayList() dc['constraints'].clear() dc['priors'].clear() dc['size'] = 0 diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 56d785c3..6fd60442 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -9,7 +9,7 @@ from re import compile, _pattern_type from param import ParamConcatenation from parameter_core import Pickleable, Parameterizable, adjust_name_for_printing, Gradcheckable from transformations import __fixed__ -from array_core import ParamList +from lists_and_dicts import ArrayList class Parameterized(Parameterizable, Pickleable, Gradcheckable): """ @@ -56,7 +56,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): def __init__(self, name=None, *a, **kw): super(Parameterized, self).__init__(name=name, parent=None, parent_index=None, *a, **kw) self._in_init_ = True - self._parameters_ = ParamList() + self._parameters_ = ArrayList() self.size = sum(p.size for p in self._parameters_) self.add_observer(self, self._parameters_changed_notification, -100) if not self._has_fixes(): @@ -265,16 +265,6 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): if self._has_fixes(): return g[self._fixes_] return g - #=========================================================================== - # Indexable Handling - #=========================================================================== - def _backtranslate_index(self, param, ind): - # translate an index in parameterized indexing into the index of param - ind = ind - self._offset_for(param) - ind = ind[ind >= 0] - internal_offset = param._internal_offset() - ind = ind[ind < param.size + internal_offset] - return ind def _offset_for(self, param): # get the offset in the parameterized index array for param if param.has_parent(): @@ -300,35 +290,21 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): """ return numpy.r_[:self.size] - #=========================================================================== - # Fixing parameters: - #=========================================================================== - def _fixes_for(self, param): - if self._has_fixes(): - return self._fixes_[self._raveled_index_for(param)] - return numpy.ones(self.size, dtype=bool)[self._raveled_index_for(param)] - #=========================================================================== # Convenience for fixed, tied checking of param: #=========================================================================== - def fixed_indices(self): - return np.array([x.is_fixed for x in self._parameters_]) - def _is_fixed(self, param): - # returns if the whole param is fixed - if not self._has_fixes(): - return False - return not self._fixes_[self._raveled_index_for(param)].any() - # return not self._fixes_[self._offset_for(param): self._offset_for(param)+param._realsize_].any() @property def is_fixed(self): for p in self._parameters_: if not p.is_fixed: return False return True + def _get_original(self, param): # if advanced indexing is activated it happens that the array is a copy # you can retrieve the original param through this method, by passing # the copy here return self._parameters_[param._parent_index_] + #=========================================================================== # Get/set parameters: #=========================================================================== From 024b92996e5b92dbe1ad69a55a0cc2546e224f27 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 28 Feb 2014 17:18:34 +0000 Subject: [PATCH 48/51] caching switched on --- GPy/util/caching.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index a2434c80..250efe11 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -40,7 +40,7 @@ class Cacher(object): return self.operation(*args) # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING - return self.operation(*args) + # return self.operation(*args) #if the result is cached, return the cached computation state = [all(a is b for a, b in itertools.izip_longest(args, cached_i)) for cached_i in self.cached_inputs] From 8b197b79a079583dbbc4c3ff8f7ddf2ff086e9c5 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 28 Feb 2014 17:35:00 +0000 Subject: [PATCH 49/51] stability in stationary) --- GPy/kern/_src/stationary.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/GPy/kern/_src/stationary.py b/GPy/kern/_src/stationary.py index 5f88c3b0..44e17d8a 100644 --- a/GPy/kern/_src/stationary.py +++ b/GPy/kern/_src/stationary.py @@ -87,7 +87,9 @@ class Stationary(Kern): """ if X2 is None: Xsq = np.sum(np.square(X),1) - return np.sqrt(-2.*tdot(X) + (Xsq[:,None] + Xsq[None,:])) + r2 = -2.*tdot(X) + (Xsq[:,None] + Xsq[None,:]) + util.diag.view(r2)[:,]= 0. # force diagnoal to be zero: sometime numerically a little negative + return np.sqrt(r2) else: X1sq = np.sum(np.square(X),1) X2sq = np.sum(np.square(X2),1) From 32c6237672c8aef7616ec2b90b737799febbc441 Mon Sep 17 00:00:00 2001 From: Zhenwen Dai Date: Fri, 28 Feb 2014 18:11:00 +0000 Subject: [PATCH 50/51] changes on rbf --- GPy/kern/_src/rbf.py | 1 - 1 file changed, 1 deletion(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 7bf0adeb..24b70671 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -257,7 +257,6 @@ class RBF(Stationary): #allocate memory for the things we want to compute mudist = np.empty((N, M, M, Q)) mudist_sq = np.empty((N, M, M, Q)) - exponent = np.zeros((N,M,M)) psi2 = np.empty((N, M, M)) l2 = self.lengthscale **2 From 496624af784c18c1ff65f5b0ea2c4dab34c2eb3b Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 28 Feb 2014 21:23:47 +0000 Subject: [PATCH 51/51] weaving a faster rbf --- GPy/kern/_src/rbf.py | 52 +++++++++++++++++++++++++++++++++++++++----- 1 file changed, 47 insertions(+), 5 deletions(-) diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index baa5b932..007bac77 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -91,13 +91,11 @@ class RBF(Stationary): #from psi2 S = variational_posterior.variance _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2) / (2.*S[:,None,None,:] + l2)*self.lengthscale - dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] if not self.ARD: - self.lengthscale.gradient += dpsi2_dlength.sum() + self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2).sum() else: - self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) + self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2) self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance @@ -224,7 +222,7 @@ class RBF(Stationary): code = """ double tmp, exponent_tmp; - //#pragma omp parallel for private(tmp, exponent_tmp) + #pragma omp parallel for private(tmp, exponent_tmp) for (int n=0; nl', dL_dpsi2, psi2, Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2, 1./(2.*S + l2))*self.lengthscale + + result = np.zeros(self.input_dim) + code = """ + double tmp; + for(int q=0; q + #include + """ + N,Q = S.shape + M = psi2.shape[-1] + + S = param_to_array(S) + weave.inline(code, support_code=support_code, libraries=['gomp'], + arg_names=['psi2', 'dL_dpsi2', 'N', 'M', 'Q', 'mudist_sq', 'l2', 'Zdist_sq', 'S', 'result'], + type_converters=weave.converters.blitz, **self.weave_options) + + return 2.*result*self.lengthscale