From d3321251ef0fcaf1d995dea69ed7374cd77db9a0 Mon Sep 17 00:00:00 2001 From: Mike Croucher Date: Mon, 7 Sep 2015 16:34:53 +0100 Subject: [PATCH 1/3] Switched to scipy.special.log1p@ --- GPy/likelihoods/link_functions.py | 2 +- GPy/testing/link_function_tests.py | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/GPy/likelihoods/link_functions.py b/GPy/likelihoods/link_functions.py index 30ad32ad..4947fdb8 100644 --- a/GPy/likelihoods/link_functions.py +++ b/GPy/likelihoods/link_functions.py @@ -141,7 +141,7 @@ class Log_ex_1(GPTransformation): """ def transf(self,f): - return scipy.log1p(safe_exp(f)) + return scipy.special.log1p(safe_exp(f)) def dtransf_df(self,f): ef = safe_exp(f) diff --git a/GPy/testing/link_function_tests.py b/GPy/testing/link_function_tests.py index a4b631f8..2f8fc5a8 100644 --- a/GPy/testing/link_function_tests.py +++ b/GPy/testing/link_function_tests.py @@ -1,5 +1,5 @@ import numpy as np -import scipy as sp +import scipy from scipy.special import cbrt from GPy.models import GradientChecker _lim_val = np.finfo(np.float64).max @@ -92,8 +92,8 @@ class LinkFunctionTests(np.testing.TestCase): link = Log_ex_1() lim_of_inf = _lim_val_exp - np.testing.assert_almost_equal(np.log1p(np.exp(self.mid_f)), link.transf(self.mid_f)) - assert np.isinf(np.log1p(np.exp(np.log(self.f_upper_lim)))) + np.testing.assert_almost_equal(scipy.special.log1p(np.exp(self.mid_f)), link.transf(self.mid_f)) + assert np.isinf(scipy.special.log1p(np.exp(np.log(self.f_upper_lim)))) #Check the clipping works np.testing.assert_almost_equal(link.transf(self.f_lower_lim), 0, decimal=5) #Need to look at most significant figures here rather than the decimals From 415d99d62d820211bc0afdcfe969f912f8fbcc91 Mon Sep 17 00:00:00 2001 From: Mike Croucher Date: Mon, 7 Sep 2015 16:41:14 +0100 Subject: [PATCH 2/3] Used scipy.log1p since it gives more consistent results cross-platform --- GPy/testing/link_function_tests.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/testing/link_function_tests.py b/GPy/testing/link_function_tests.py index 2f8fc5a8..9f41f736 100644 --- a/GPy/testing/link_function_tests.py +++ b/GPy/testing/link_function_tests.py @@ -97,13 +97,13 @@ class LinkFunctionTests(np.testing.TestCase): #Check the clipping works np.testing.assert_almost_equal(link.transf(self.f_lower_lim), 0, decimal=5) #Need to look at most significant figures here rather than the decimals - np.testing.assert_approx_equal(link.transf(self.f_upper_lim), np.log1p(_lim_val), significant=5) + np.testing.assert_approx_equal(link.transf(self.f_upper_lim), scipy.special.log1p(_lim_val), significant=5) self.check_overflow(link, lim_of_inf) #Check that it would otherwise fail beyond_lim_of_inf = lim_of_inf + 10.0 old_err_state = np.seterr(over='ignore') - self.assertTrue(np.isinf(np.log1p(np.exp(beyond_lim_of_inf)))) + self.assertTrue(np.isinf(scipy.special.log1p(np.exp(beyond_lim_of_inf)))) np.seterr(**old_err_state) From 929cf0a4890e418ecec0b000ed7fefa2372bc082 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 7 Sep 2015 16:52:59 +0100 Subject: [PATCH 3/3] [more coverage] and predictive var fixes --- GPy/core/gp.py | 40 +++++++++++++------ GPy/models/gplvm.py | 2 + .../matplot_dep/dim_reduction_plots.py | 2 +- GPy/testing/run_coverage.sh | 2 +- 4 files changed, 32 insertions(+), 14 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index ad082b3c..9a199faa 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -108,9 +108,15 @@ class GP(Model): # The predictive variable to be used to predict using the posterior object's # woodbury_vector and woodbury_inv is defined as predictive_variable + # as long as the posterior has the right woodbury entries. + # It is the input variable used for the covariance between + # X_star and the posterior of the GP. # This is usually just a link to self.X (full GP) or self.Z (sparse GP). # Make sure to name this variable and the predict functions will "just work" - # as long as the posterior has the right woodbury entries. + # In maths the predictive variable is: + # K_{xx} - K_{xp}W_{pp}^{-1}K_{px} + # W_{pp} := \texttt{Woodbury inv} + # p := _predictive_variable self._predictive_variable = self.X @@ -213,7 +219,7 @@ class GP(Model): Kxx = kern.K(Xnew) if self.posterior.woodbury_inv.ndim == 2: var = Kxx - np.dot(Kx.T, np.dot(self.posterior.woodbury_inv, Kx)) - elif self.posterior.woodbury_inv.ndim == 3: + elif self.posterior.woodbury_inv.ndim == 3: # Missing data var = np.empty((Kxx.shape[0],Kxx.shape[1],self.posterior.woodbury_inv.shape[2])) from ..util.linalg import mdot for i in range(var.shape[2]): @@ -223,7 +229,7 @@ class GP(Model): Kxx = kern.Kdiag(Xnew) if self.posterior.woodbury_inv.ndim == 2: var = (Kxx - np.sum(np.dot(self.posterior.woodbury_inv.T, Kx) * Kx, 0))[:,None] - elif self.posterior.woodbury_inv.ndim == 3: + elif self.posterior.woodbury_inv.ndim == 3: # Missing data var = np.empty((Kxx.shape[0],self.posterior.woodbury_inv.shape[2])) for i in range(var.shape[1]): var[:, i] = (Kxx - (np.sum(np.dot(self.posterior.woodbury_inv[:, :, i].T, Kx) * Kx, 0))) @@ -364,11 +370,15 @@ class GP(Model): var_jac = dK2_dXdX - np.einsum('qim,miq->iq', dK_dXnew_full.T.dot(wi), dK_dXnew_full) return var_jac - if self.posterior.woodbury_inv.ndim == 3: - var_jac = [] - for d in range(self.posterior.woodbury_inv.shape[2]): - var_jac.append(compute_cov_inner(self.posterior.woodbury_inv[:, :, d])) - var_jac = np.concatenate(var_jac) + if self.posterior.woodbury_inv.ndim == 3: # Missing data: + if full_cov: + var_jac = np.empty((Xnew.shape[0],Xnew.shape[0],Xnew.shape[1],self.output_dim)) + for d in range(self.posterior.woodbury_inv.shape[2]): + var_jac[:, :, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d]) + else: + var_jac = np.empty((Xnew.shape[0],Xnew.shape[1],self.output_dim)) + for d in range(self.posterior.woodbury_inv.shape[2]): + var_jac[:, :, d] = compute_cov_inner(self.posterior.woodbury_inv[:, :, d]) else: var_jac = compute_cov_inner(self.posterior.woodbury_inv) return mean_jac, var_jac @@ -391,10 +401,11 @@ class GP(Model): mu_jac, var_jac = self.predict_jacobian(Xnew, kern, full_cov=False) mumuT = np.einsum('iqd,ipd->iqp', mu_jac, mu_jac) + Sigma = np.zeros(mumuT.shape) if var_jac.ndim == 3: - Sigma = np.einsum('iqd,ipd->iqp', var_jac, var_jac) + Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = var_jac.sum(-1) else: - Sigma = self.output_dim*np.einsum('iq,ip->iqp', var_jac, var_jac) + Sigma[(slice(None), )+np.diag_indices(Xnew.shape[1], 2)] = self.output_dim*var_jac G = 0. if mean: G += mumuT @@ -412,8 +423,13 @@ class GP(Model): """ G = self.predict_wishard_embedding(Xnew, kern, mean, covariance) from ..util.linalg import jitchol - return np.array([np.sqrt(np.exp(2*np.sum(np.log(np.diag(jitchol(G[n, :, :])))))) for n in range(Xnew.shape[0])]) - #return np.array([np.sqrt(np.linalg.det(G[n, :, :])) for n in range(Xnew.shape[0])]) + mag = np.empty(Xnew.shape[0]) + for n in range(Xnew.shape[0]): + try: + mag[n] = np.sqrt(np.exp(2*np.sum(np.log(np.diag(jitchol(G[n, :, :])))))) + except: + mag[n] = np.sqrt(np.linalg.det(G[n, :, :])) + return mag def posterior_samples_f(self,X,size=10, full_cov=True): """ diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py index d4f4f564..17d42e5a 100644 --- a/GPy/models/gplvm.py +++ b/GPy/models/gplvm.py @@ -36,8 +36,10 @@ class GPLVM(GP): likelihood = Gaussian() super(GPLVM, self).__init__(X, Y, kernel, likelihood, name='GPLVM') + self.X = Param('latent_mean', X) self.link_parameter(self.X, index=0) + self._predictive_variable = self.X def parameters_changed(self): super(GPLVM, self).parameters_changed() diff --git a/GPy/plotting/matplot_dep/dim_reduction_plots.py b/GPy/plotting/matplot_dep/dim_reduction_plots.py index a36f168d..ff0887b8 100644 --- a/GPy/plotting/matplot_dep/dim_reduction_plots.py +++ b/GPy/plotting/matplot_dep/dim_reduction_plots.py @@ -119,7 +119,7 @@ def plot_latent(model, labels=None, which_indices=None, Xtest_full[:, [input_1, input_2]] = x _, var = model.predict(Xtest_full, **predict_kwargs) var = var[:, :1] - return np.log(var) + return 2*np.sqrt(var) #Create an IMshow controller that can re-plot the latent space shading at a good resolution if plot_limits is None: diff --git a/GPy/testing/run_coverage.sh b/GPy/testing/run_coverage.sh index f2e52230..f9b33015 100755 --- a/GPy/testing/run_coverage.sh +++ b/GPy/testing/run_coverage.sh @@ -1 +1 @@ -nosetests . --with-coverage --logging-level=INFO --cover-html --cover-html-dir=coverage --cover-package=GPy --cover-erase +nosetests . --with-coverage --logging-level=INFO --cover-html --cover-html-dir=coverage --cover-package=GPy --cover-erase --cover-omit=GPy.examples