From c76f1a4d6dcc70115e0eff4520f8257f395292a9 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Tue, 11 Feb 2014 14:06:42 +0000 Subject: [PATCH 01/40] Minor reorganising --- GPy/core/parameterization/param.py | 8 ++++---- GPy/inference/latent_function_inference/laplace.py | 5 ++--- GPy/testing/likelihood_tests.py | 4 ++-- 3 files changed, 8 insertions(+), 9 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 583c6425..772cbd42 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -157,14 +157,14 @@ class Param(ObservableArray, Constrainable, Gradcheckable): #=========================================================================== def tie_to(self, param): """ - :param param: the parameter object to tie this parameter to. + :param param: the parameter object to tie this parameter to. Can be ParamConcatenation (retrieved by regexp search) - + Tie this parameter to the given parameter. Broadcasting is not allowed, but you can tie a whole dimension to one parameter: self[:,0].tie_to(other), where other is a one-value parameter. - + Note: For now only one parameter can have ties, so all of a parameter will be removed, when re-tieing! """ @@ -534,7 +534,7 @@ class ParamConcatenation(object): def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance) #checkgrad.__doc__ = Gradcheckable.checkgrad.__doc__ - + __lt__ = lambda self, val: self._vals() < val __le__ = lambda self, val: self._vals() <= val __eq__ = lambda self, val: self._vals() == val diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index bc81a86a..4edb9a1d 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -92,12 +92,11 @@ class LaplaceInference(object): iteration = 0 while difference > self._mode_finding_tolerance and iteration < self._mode_finding_max_iter: W = -likelihood.d2logpdf_df2(f, Y, extra_data=Y_metadata) - - W_f = W*f grad = likelihood.dlogpdf_df(f, Y, extra_data=Y_metadata) + W_f = W*f + b = W_f + grad # R+W p46 line 6. - #W12BiW12Kb, B_logdet = self._compute_B_statistics(K, W.copy(), np.dot(K, b), likelihood.log_concave) W12BiW12, _, _ = self._compute_B_statistics(K, W, likelihood.log_concave) W12BiW12Kb = np.dot(W12BiW12, np.dot(K, b)) diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index 7f48ac95..9920d648 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -1,10 +1,10 @@ import numpy as np import unittest import GPy -from GPy.models import GradientChecker +from ..models import GradientChecker import functools import inspect -from GPy.likelihoods import link_functions +from ..likelihoods import link_functions from ..core.parameterization import Param from functools import partial #np.random.seed(300) From 610dbed38a313303946696c4c54b6770b100127b Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Tue, 11 Feb 2014 15:23:00 +0000 Subject: [PATCH 02/40] Fixed copy bug of observable array --- GPy/core/parameterization/array_core.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 4b5b7700..fe2a951e 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -25,7 +25,7 @@ class ParamList(list): if el is other: return True return False - + pass class C(np.ndarray): __array_priority__ = 1. @@ -80,8 +80,8 @@ class ObservableArray(ListArray, 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) + return self.__setitem__(slice(start, stop), val) def __copy__(self, *args): - return ObservableArray(self.base.base.copy(*args)) + return ObservableArray(self.view(np.ndarray).copy()) def copy(self, *args): return self.__copy__(*args) From 48c68640ef218c93e0e19b5f7ae05f52ab0d2a40 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 11 Feb 2014 16:30:12 +0000 Subject: [PATCH 03/40] renaming dtc --- GPy/inference/latent_function_inference/__init__.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py index 5184f0b4..139c562f 100644 --- a/GPy/inference/latent_function_inference/__init__.py +++ b/GPy/inference/latent_function_inference/__init__.py @@ -26,3 +26,4 @@ etc. from exact_gaussian_inference import ExactGaussianInference from laplace import LaplaceInference expectation_propagation = 'foo' # TODO +from dtc import DTC From 1b29c0939e13699775b403da4b3610052a099f0a Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 11 Feb 2014 16:44:59 +0000 Subject: [PATCH 04/40] renaming dtc again --- GPy/inference/latent_function_inference/{DTC.py => dtc.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename GPy/inference/latent_function_inference/{DTC.py => dtc.py} (100%) diff --git a/GPy/inference/latent_function_inference/DTC.py b/GPy/inference/latent_function_inference/dtc.py similarity index 100% rename from GPy/inference/latent_function_inference/DTC.py rename to GPy/inference/latent_function_inference/dtc.py From 54b7581c39688e2fefaae3da77216fd7f184f77a Mon Sep 17 00:00:00 2001 From: James Hensman Date: Tue, 11 Feb 2014 20:05:36 +0000 Subject: [PATCH 05/40] changes to DTC --- GPy/core/sparse_gp.py | 4 ++-- GPy/inference/latent_function_inference/dtc.py | 12 ++++++------ 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 1879145a..73cc8755 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -38,9 +38,9 @@ class SparseGP(GP): if inference_method is None: if isinstance(likelihood, likelihoods.Gaussian): inference_method = varDTC.VarDTC() - else: + else: #inference_method = ?? - raise NotImplementedError, "what to do what to do?" + raise NotImplementedError, "what to do what to do?" print "defaulting to ", inference_method, "for latent function inference" self.Z = Param('inducing inputs', Z) diff --git a/GPy/inference/latent_function_inference/dtc.py b/GPy/inference/latent_function_inference/dtc.py index cc6e8606..c6f6844c 100644 --- a/GPy/inference/latent_function_inference/dtc.py +++ b/GPy/inference/latent_function_inference/dtc.py @@ -52,20 +52,20 @@ class DTC(object): b, _ = dtrtrs(LA, tmp*beta, lower=1) tmp, _ = dtrtrs(LA, b, lower=1, trans=1) v, _ = dtrtrs(L, tmp, lower=1, trans=1) - tmp = tdrtrs(LA, Li, lower=1, trans=0) + tmp, _ = dtrtrs(LA, Li, lower=1, trans=0) P = tdot(tmp.T) #compute log marginal log_marginal = -0.5*num_data*output_dim*np.log(2*np.pi) + \ -np.sum(np.log(np.diag(LA)))*output_dim + \ 0.5*num_data*output_dim*np.log(beta) + \ - -0.5*beta*np.sum(np.square(Y)) + + -0.5*beta*np.sum(np.square(Y)) + \ 0.5*np.sum(np.square(b)) # Compute dL_dKmm tmp, _ = dtrtrs(L, A_I, lower=1, trans=1) dL_dK, _ = dtrtrs(L, tmp.T, lower=1, trans=0) - tmp, _ = dtrtrs(LA, tmp.T. lower=1, trans=1) + tmp, _ = dtrtrs(LA, tmp.T, lower=1, trans=1) dL_dK -= tdot(tmp.T) dL_dK *= output_dim dL_dK -= tdot(v) @@ -79,17 +79,17 @@ class DTC(object): #compute dL_dR Uv = np.dot(U, v) - dL_dR = 0.5*(np.sum(U*np.dot(P, U.T), 1) - beta * np.sum(np.square(Y, 1)) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) + dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - beta * np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) )*beta**2 - grad_dict = {'dL_dKmm': dL_dKmm, 'dL_dKdiag':np.zeros_like(Knn), 'dL_dKnm':dL_dU} + grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn), 'dL_dKnm':dL_dU.T} #update gradients kern.update_gradients_sparse(X=X, Z=Z, **grad_dict) likelihood.update_gradients(dL_dR) #construct a posterior object - post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=Lm) + post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L) return post, log_marginal, grad_dict From 34dec0ade3d55f358a34f3c8587c9499d0c955f6 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 12 Feb 2014 11:13:24 +0000 Subject: [PATCH 06/40] derivatives working in DTC --- .../latent_function_inference/dtc.py | 21 +++++++------------ 1 file changed, 7 insertions(+), 14 deletions(-) diff --git a/GPy/inference/latent_function_inference/dtc.py b/GPy/inference/latent_function_inference/dtc.py index c6f6844c..85792d66 100644 --- a/GPy/inference/latent_function_inference/dtc.py +++ b/GPy/inference/latent_function_inference/dtc.py @@ -40,9 +40,9 @@ class DTC(object): Kmmi, L, Li, _ = pdinv(Kmm) # Compute A - LiUT, _ = dtrtrs(L, U.T*np.sqrt(beta), lower=1) - A_I = tdot(LiUT) - A = A_I + np.eye(num_inducing) + #LiUT, _ = dtrtrs(L, U.T*np.sqrt(beta), lower=1) + LiUT = np.dot(Li, U.T)*np.sqrt(beta) + A = tdot(LiUT) + np.eye(num_inducing) # factor A LA = jitchol(A) @@ -63,24 +63,17 @@ class DTC(object): 0.5*np.sum(np.square(b)) # Compute dL_dKmm - tmp, _ = dtrtrs(L, A_I, lower=1, trans=1) - dL_dK, _ = dtrtrs(L, tmp.T, lower=1, trans=0) - tmp, _ = dtrtrs(LA, tmp.T, lower=1, trans=1) - dL_dK -= tdot(tmp.T) - dL_dK *= output_dim - dL_dK -= tdot(v) - dL_dK /=2. + vvT_P = tdot(v.reshape(-1,1)) + P + dL_dK = 0.5*(Kmmi - vvT_P) # Compute dL_dU - vvT_P = tdot(v.reshape(-1,1)) + P vY = np.dot(v.reshape(-1,1),Y.T) - dL_dU = vY + np.dot(vvT_P, U.T) + dL_dU = vY - np.dot(vvT_P, U.T) dL_dU *= beta #compute dL_dR Uv = np.dot(U, v) - dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - beta * np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1) - )*beta**2 + dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./beta + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1))*beta**2 grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':np.zeros_like(Knn), 'dL_dKnm':dL_dU.T} From a29cf87e53b28864381edaca71355c10dee64512 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Feb 2014 11:31:12 +0000 Subject: [PATCH 07/40] deleted listarray --- GPy/core/parameterization/array_core.py | 15 +-------------- GPy/core/parameterization/param.py | 4 +--- 2 files changed, 2 insertions(+), 17 deletions(-) diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 1d300d80..dba813a5 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 ListArray(np.ndarray): - """ - ndarray which can be stored in lists and checked if it is in. - WARNING: This overrides the functionality of x==y!!! - Use numpy.equal(x,y) for element-wise equality testing. - """ - - def __new__(cls, input_array): - obj = np.asanyarray(input_array).view(cls) - return obj - #def __eq__(self, other): - # return other is self - class ParamList(list): def __contains__(self, other): @@ -29,7 +16,7 @@ class ParamList(list): pass -class ObservableArray(ListArray, Observable): +class ObservableArray(np.ndarray, Observable): """ An ndarray which reports changes to its observers. The observers can add themselves with a callable, which diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index b76b3858..5d987ae6 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -46,6 +46,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable): """ __array_priority__ = -1 # Never give back Param _fixes_ = None + _parameters_ = [] def __new__(cls, name, input_array, default_constraint=None): obj = numpy.atleast_1d(super(Param, cls).__new__(cls, input_array=input_array)) obj._current_slice_ = (slice(obj.shape[0]),) @@ -144,9 +145,6 @@ class Param(ObservableArray, Constrainable, Gradcheckable): # from_name = self.name # self.name = new_name # self._direct_parent_._name_changed(self, from_name) - @property - def _parameters_(self): - return [] def _collect_gradient(self, target): target[:] = self.gradient.flat #=========================================================================== From 87acda0034ee82aee600bc08028332b5dabfedf0 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 12 Feb 2014 12:29:52 +0000 Subject: [PATCH 08/40] here's fitc --- .../latent_function_inference/__init__.py | 1 + .../latent_function_inference/fitc.py | 264 +++++------------- 2 files changed, 66 insertions(+), 199 deletions(-) diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py index 139c562f..c89f771b 100644 --- a/GPy/inference/latent_function_inference/__init__.py +++ b/GPy/inference/latent_function_inference/__init__.py @@ -27,3 +27,4 @@ from exact_gaussian_inference import ExactGaussianInference from laplace import LaplaceInference expectation_propagation = 'foo' # TODO from dtc import DTC +from fitc import FITC diff --git a/GPy/inference/latent_function_inference/fitc.py b/GPy/inference/latent_function_inference/fitc.py index 82d0973f..476715e7 100644 --- a/GPy/inference/latent_function_inference/fitc.py +++ b/GPy/inference/latent_function_inference/fitc.py @@ -1,225 +1,91 @@ -# Copyright (c) 2012, GPy authors (see AUTHORS.txt). +# Copyright (c) 2012, James Hensman # Licensed under the BSD 3-clause license (see LICENSE.txt) -class VarDTC(object): +from posterior import Posterior +from ...util.linalg import jitchol, tdot, dtrtrs, dpotri, pdinv +import numpy as np +log_2_pi = np.log(2*np.pi) + +class FITC(object): + """ + An object for inference when the likelihood is Gaussian, but we want to do sparse inference. + + The function self.inference returns a Posterior object, which summarizes + the posterior. + + """ def __init__(self): self.const_jitter = 1e-6 def inference(self, kern, X, X_variance, Z, likelihood, Y): + assert X_variance is None, "cannot use X_variance with FITC. Try varDTC." num_inducing, _ = Z.shape num_data, output_dim = Y.shape - #make sure we're not using variational uncertain inputs - assert = X_variance is None, "variational inducing inputs only for use with varDTC inference" - #Note: we can;t do the variational thing after making the GITC conditional approximation because K~ appears in the log determinant. + #make sure the noise is not hetero + sigma_n = np.squeeze(likelihood.variance) + if sigma_n.size <1: + raise NotImplementedError, "no hetero noise with this implementation of FITC" - #see whether we've got a different noise variance for each datum - sigma2 = np.squeeze(likelihood.variance) - het_noise = False - if sigma2.size <1: - het_noise = True - - # kernel computations, using BGPLVM notation Kmm = kern.K(Z) - Kdiag = kern.Kdiag(X) - Kmn = kern.K(X, Z) + Knn = kern.Kdiag(X) + Knm = kern.K(X, Z) + U = Knm - #factor Kmm - Lm = jitchol(Kmm) - V = dtrtrs(Lm, Kmn, lower=1) + #factor Kmm + Kmmi, L, Li, _ = pdinv(Kmm) - #compute effective noise - g_sn2 = sigma2 + Kdiag - np.sum(V*V, 0) + #compute beta_star, the effective noise precision + LiUT = np.dot(Li, U.T) + sigma_star = Knn + sigma_n - np.sum(np.square(LiUT),0) + beta_star = 1./sigma_star - #compute and factor B - tmp = Kmn / np.sqrt(g_snd) - tmp, _ = dtrtrs(Lm, tmp, lower=1) - A = tdot(tmp) - B = np.eye(num_inducing) + A - Bi, Lb, LBi, log_det_B = pdinv(B) + # Compute and factor A + A = tdot(LiUT*np.sqrt(beta_star)) + np.eye(num_inducing) + LA = jitchol(A) - #compute posterior parameters - tmp, _ = dtrtrs() - woodbury_vec = + # back substutue to get b, P, v + URiy = np.dot(U.T*beta_star,Y) + tmp, _ = dtrtrs(L, URiy, lower=1) + b, _ = dtrtrs(LA, tmp, lower=1) + tmp, _ = dtrtrs(LA, b, lower=1, trans=1) + v, _ = dtrtrs(L, tmp, lower=1, trans=1) + tmp, _ = dtrtrs(LA, Li, lower=1, trans=0) + P = tdot(tmp.T) - + #compute log marginal + log_marginal = -0.5*num_data*output_dim*np.log(2*np.pi) + \ + -np.sum(np.log(np.diag(LA)))*output_dim + \ + 0.5*output_dim*np.sum(np.log(beta_star)) + \ + -0.5*np.sum(np.square(Y.T*np.sqrt(beta_star))) + \ + 0.5*np.sum(np.square(b)) + #compute dL_dR + Uv = np.dot(U, v) + dL_dR = 0.5*(np.sum(U*np.dot(U,P), 1) - 1./beta_star + np.sum(np.square(Y), 1) - 2.*np.sum(Uv*Y, 1) + np.sum(np.square(Uv), 1))*beta_star**2 - psi1V = np.dot(self.psi1, self.V_star) - Lmi_psi1V, _ = dtrtrs(Lm, self.psi1V, lower=1, trans=0) - LBi_Lmi_psi1V, _ = dtrtrs(LB, Lmi_psi1V, lower=1, trans=0) + # Compute dL_dKmm + vvT_P = tdot(v.reshape(-1,1)) + P + dL_dK = 0.5*(Kmmi - vvT_P) + KiU = np.dot(Kmmi, U.T) + dL_dK += np.dot(KiU*dL_dR, KiU.T) - Kmmipsi1 = np.dot(self.Lmi.T,Lmipsi1) - b_psi1_Ki = self.beta_star * Kmmipsi1.T - Ki_pbp_Ki = np.dot(Kmmipsi1,b_psi1_Ki) - Kmmi = np.dot(self.Lmi.T,self.Lmi) - LBiLmi = np.dot(self.LBi,self.Lmi) - LBL_inv = np.dot(LBiLmi.T,LBiLmi) - VVT = np.outer(self.V_star,self.V_star) - VV_p_Ki = np.dot(VVT,Kmmipsi1.T) - Ki_pVVp_Ki = np.dot(Kmmipsi1,VV_p_Ki) - psi1beta = self.psi1*self.beta_star.T - H = self.Kmm + mdot(self.psi1,psi1beta.T) - LH = jitchol(H) - LHi = chol_inv(LH) - Hi = np.dot(LHi.T,LHi) + # Compute dL_dU + vY = np.dot(v.reshape(-1,1),Y.T) + dL_dU = vY - np.dot(vvT_P, U.T) + dL_dU *= beta_star + dL_dU -= 2.*KiU*dL_dR - betapsi1TLmiLBi = np.dot(psi1beta.T,LBiLmi.T) - alpha = np.array([np.dot(a.T,a) for a in betapsi1TLmiLBi])[:,None] - gamma_1 = mdot(VVT,self.psi1.T,Hi) - pHip = mdot(self.psi1.T,Hi,self.psi1) - gamma_2 = mdot(self.beta_star*pHip,self.V_star) - gamma_3 = self.V_star * gamma_2 + grad_dict = {'dL_dKmm': dL_dK, 'dL_dKdiag':dL_dR, 'dL_dKnm':dL_dU.T} - self._dL_dpsi0 = -0.5 * self.beta_star#dA_dpsi0: logdet(self.beta_star) - self._dL_dpsi0 += .5 * self.V_star**2 #dA_psi0: yT*beta_star*y - self._dL_dpsi0 += .5 *alpha #dC_dpsi0 - self._dL_dpsi0 += 0.5*mdot(self.beta_star*pHip,self.V_star)**2 - self.V_star * mdot(self.V_star.T,pHip*self.beta_star).T #dD_dpsi0 + #update gradients + kern.update_gradients_sparse(X=X, Z=Z, **grad_dict) + likelihood.update_gradients(dL_dR) - self._dL_dpsi1 = b_psi1_Ki.copy() #dA_dpsi1: logdet(self.beta_star) - self._dL_dpsi1 += -np.dot(psi1beta.T,LBL_inv) #dC_dpsi1 - self._dL_dpsi1 += gamma_1 - mdot(psi1beta.T,Hi,self.psi1,gamma_1) #dD_dpsi1 + #construct a posterior object + post = Posterior(woodbury_inv=Kmmi-P, woodbury_vector=v, K=Kmm, mean=None, cov=None, K_chol=L) - self._dL_dKmm = -0.5 * np.dot(Kmmipsi1,b_psi1_Ki) #dA_dKmm: logdet(self.beta_star) - self._dL_dKmm += .5*(LBL_inv - Kmmi) + mdot(LBL_inv,psi1beta,Kmmipsi1.T) #dC_dKmm - self._dL_dKmm += -.5 * mdot(Hi,self.psi1,gamma_1) #dD_dKmm + return post, log_marginal, grad_dict - self._dpsi1_dtheta = 0 - self._dpsi1_dX = 0 - self._dKmm_dtheta = 0 - self._dKmm_dX = 0 - self._dpsi1_dX_jkj = 0 - self._dpsi1_dtheta_jkj = 0 - - for i,V_n,alpha_n,gamma_n,gamma_k in zip(range(self.num_data),self.V_star,alpha,gamma_2,gamma_3): - K_pp_K = np.dot(Kmmipsi1[:,i:(i+1)],Kmmipsi1[:,i:(i+1)].T) - _dpsi1 = (-V_n**2 - alpha_n + 2.*gamma_k - gamma_n**2) * Kmmipsi1.T[i:(i+1),:] - _dKmm = .5*(V_n**2 + alpha_n + gamma_n**2 - 2.*gamma_k) * K_pp_K #Diag_dD_dKmm - self._dpsi1_dtheta += self.kern._param_grad_helper(_dpsi1,self.X[i:i+1,:],self.Z) - self._dKmm_dtheta += self.kern._param_grad_helper(_dKmm,self.Z) - self._dKmm_dX += self.kern.gradients_X(_dKmm ,self.Z) - self._dpsi1_dX += self.kern.gradients_X(_dpsi1.T,self.Z,self.X[i:i+1,:]) - - # the partial derivative vector for the likelihood - if self.likelihood.num_params == 0: - # save computation here. - self.partial_for_likelihood = None - elif self.likelihood.is_heteroscedastic: - raise NotImplementedError, "heteroscedatic derivates not implemented." - else: - # likelihood is not heterscedatic - dbstar_dnoise = self.likelihood.precision * (self.beta_star**2 * self.Diag0[:,None] - self.beta_star) - Lmi_psi1 = mdot(self.Lmi,self.psi1) - LBiLmipsi1 = np.dot(self.LBi,Lmi_psi1) - aux_0 = np.dot(self._LBi_Lmi_psi1V.T,LBiLmipsi1) - aux_1 = self.likelihood.Y.T * np.dot(self._LBi_Lmi_psi1V.T,LBiLmipsi1) - aux_2 = np.dot(LBiLmipsi1.T,self._LBi_Lmi_psi1V) - - dA_dnoise = 0.5 * self.input_dim * (dbstar_dnoise/self.beta_star).sum() - 0.5 * self.input_dim * np.sum(self.likelihood.Y**2 * dbstar_dnoise) - dC_dnoise = -0.5 * np.sum(mdot(self.LBi.T,self.LBi,Lmi_psi1) * Lmi_psi1 * dbstar_dnoise.T) - - dD_dnoise_1 = mdot(self.V_star*LBiLmipsi1.T,LBiLmipsi1*dbstar_dnoise.T*self.likelihood.Y.T) - alpha = mdot(LBiLmipsi1,self.V_star) - alpha_ = mdot(LBiLmipsi1.T,alpha) - dD_dnoise_2 = -0.5 * self.input_dim * np.sum(alpha_**2 * dbstar_dnoise ) - - dD_dnoise_1 = mdot(self.V_star.T,self.psi1.T,self.Lmi.T,self.LBi.T,self.LBi,self.Lmi,self.psi1,dbstar_dnoise*self.likelihood.Y) - dD_dnoise_2 = 0.5*mdot(self.V_star.T,self.psi1.T,Hi,self.psi1,dbstar_dnoise*self.psi1.T,Hi,self.psi1,self.V_star) - dD_dnoise = dD_dnoise_1 + dD_dnoise_2 - - self.partial_for_likelihood = dA_dnoise + dC_dnoise + dD_dnoise - - def log_likelihood(self): - """ Compute the (lower bound on the) log marginal likelihood """ - A = -0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) + 0.5 * np.sum(np.log(self.beta_star)) - 0.5 * np.sum(self.V_star * self.likelihood.Y) - C = -self.output_dim * (np.sum(np.log(np.diag(self.LB)))) - D = 0.5 * np.sum(np.square(self._LBi_Lmi_psi1V)) - return A + C + D - - def _log_likelihood_gradients(self): - pass - return np.hstack((self.dL_dZ().flatten(), self.dL_dtheta(), self.likelihood._gradients(partial=self.partial_for_likelihood))) - - def dL_dtheta(self): - dL_dtheta = self.kern.dKdiag_dtheta(self._dL_dpsi0,self.X) - dL_dtheta += self.kern._param_grad_helper(self._dL_dpsi1,self.X,self.Z) - dL_dtheta += self.kern._param_grad_helper(self._dL_dKmm,X=self.Z) - dL_dtheta += self._dKmm_dtheta - dL_dtheta += self._dpsi1_dtheta - return dL_dtheta - - def dL_dZ(self): - dL_dZ = self.kern.gradients_X(self._dL_dpsi1.T,self.Z,self.X) - dL_dZ += self.kern.gradients_X(self._dL_dKmm,X=self.Z) - dL_dZ += self._dpsi1_dX - dL_dZ += self._dKmm_dX - return dL_dZ - - def _raw_predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): - assert X_variance_new is None, "FITC model is not defined for handling uncertain inputs." - - if self.likelihood.is_heteroscedastic: - Iplus_Dprod_i = 1./(1.+ self.Diag0 * self.likelihood.precision.flatten()) - self.Diag = self.Diag0 * Iplus_Dprod_i - self.P = Iplus_Dprod_i[:,None] * self.psi1.T - self.RPT0 = np.dot(self.Lmi,self.psi1) - self.L = np.linalg.cholesky(np.eye(self.num_inducing) + np.dot(self.RPT0,((1. - Iplus_Dprod_i)/self.Diag0)[:,None]*self.RPT0.T)) - self.R,info = dtrtrs(self.L,self.Lmi,lower=1) - self.RPT = np.dot(self.R,self.P.T) - self.Sigma = np.diag(self.Diag) + np.dot(self.RPT.T,self.RPT) - self.w = self.Diag * self.likelihood.v_tilde - self.Gamma = np.dot(self.R.T, np.dot(self.RPT,self.likelihood.v_tilde)) - self.mu = self.w + np.dot(self.P,self.Gamma) - - """ - Make a prediction for the generalized FITC model - - Arguments - --------- - X : Input prediction data - Nx1 numpy array (floats) - """ - # q(u|f) = N(u| R0i*mu_u*f, R0i*C*R0i.T) - - # Ci = I + (RPT0)Di(RPT0).T - # C = I - [RPT0] * (input_dim+[RPT0].T*[RPT0])^-1*[RPT0].T - # = I - [RPT0] * (input_dim + self.Qnn)^-1 * [RPT0].T - # = I - [RPT0] * (U*U.T)^-1 * [RPT0].T - # = I - V.T * V - U = np.linalg.cholesky(np.diag(self.Diag0) + self.Qnn) - V,info = dtrtrs(U,self.RPT0.T,lower=1) - C = np.eye(self.num_inducing) - np.dot(V.T,V) - mu_u = np.dot(C,self.RPT0)*(1./self.Diag0[None,:]) - #self.C = C - #self.RPT0 = np.dot(self.R0,self.Knm.T) P0.T - #self.mu_u = mu_u - #self.U = U - # q(u|y) = N(u| R0i*mu_H,R0i*Sigma_H*R0i.T) - mu_H = np.dot(mu_u,self.mu) - self.mu_H = mu_H - Sigma_H = C + np.dot(mu_u,np.dot(self.Sigma,mu_u.T)) - # q(f_star|y) = N(f_star|mu_star,sigma2_star) - Kx = self.kern.K(self.Z, Xnew, which_parts=which_parts) - KR0T = np.dot(Kx.T,self.Lmi.T) - mu_star = np.dot(KR0T,mu_H) - if full_cov: - Kxx = self.kern.K(Xnew,which_parts=which_parts) - var = Kxx + np.dot(KR0T,np.dot(Sigma_H - np.eye(self.num_inducing),KR0T.T)) - else: - Kxx = self.kern.Kdiag(Xnew,which_parts=which_parts) - var = (Kxx + np.sum(KR0T.T*np.dot(Sigma_H - np.eye(self.num_inducing),KR0T.T),0))[:,None] - return mu_star[:,None],var - else: - raise NotImplementedError, "Heteroscedastic case not implemented." - """ - Kx = self.kern.K(self.Z, Xnew) - mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V) - if full_cov: - Kxx = self.kern.K(Xnew) - var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx) #NOTE this won't work for plotting - else: - Kxx = self.kern.Kdiag(Xnew) - var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0) - return mu,var[:,None] - """ From cb0fb9d5934975fea90ef098a47771860abddc90 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Feb 2014 13:08:13 +0000 Subject: [PATCH 09/40] deleted list array --- GPy/core/parameterization/array_core.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 6ece51e9..a65b8c24 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -4,10 +4,13 @@ __updated__ = '2013-12-16' import numpy as np -from parameter_core import Observable +from parameter_core import Observable, Constrainable, Gradcheckable 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: @@ -26,7 +29,7 @@ class ObservableArray(np.ndarray, Observable): __array_priority__ = -1 # Never give back ObservableArray def __new__(cls, input_array): cls.__name__ = "ObservableArray\n " - obj = super(ObservableArray, cls).__new__(cls, input_array).view(cls) + obj = np.atleast_1d(input_array).view(cls) obj._observers_ = {} return obj def __array_finalize__(self, obj): From e4fcde3299a4e9af2d14584d3c8c3040fffea2b8 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Wed, 12 Feb 2014 15:27:47 +0000 Subject: [PATCH 10/40] fixed plotting bug --- GPy/plotting/matplot_dep/inference_plots.py | 1 + GPy/plotting/matplot_dep/models_plots.py | 18 ++++++++++-------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/GPy/plotting/matplot_dep/inference_plots.py b/GPy/plotting/matplot_dep/inference_plots.py index f9bb464a..6a3a8a93 100644 --- a/GPy/plotting/matplot_dep/inference_plots.py +++ b/GPy/plotting/matplot_dep/inference_plots.py @@ -2,6 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import pylab as pb +import sys #import numpy as np #import Tango #from base_plots import gpplot, x_frame1D, x_frame2D diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 0aa1d4a4..177b9a95 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -94,14 +94,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. - #add inducing inputs (if a sparse model is used) - if hasattr(model,"Z"): - #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] - Zu = param_to_array(model.Z[:,free_dims]) - z_height = ax.get_ylim()[0] - ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) - - + #add error bars for uncertain (if input uncertainty is being modelled) if hasattr(model,"has_uncertain_inputs"): ax.errorbar(model.X[which_data, free_dims], model.likelihood.data[which_data, 0], @@ -115,6 +108,15 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', ax.set_xlim(xmin, xmax) ax.set_ylim(ymin, ymax) + #add inducing inputs (if a sparse model is used) + if hasattr(model,"Z"): + #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] + Zu = param_to_array(model.Z[:,free_dims]) + z_height = ax.get_ylim()[0] + ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) + + + #2D plotting elif len(free_dims) == 2: From c788c463d8bb8553f7077ee0d296ad68b1ef587e Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Wed, 12 Feb 2014 15:53:38 +0000 Subject: [PATCH 11/40] Added metadata --- GPy/core/gp.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index ab725897..214c2324 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -26,7 +26,7 @@ class GP(Model): """ - def __init__(self, X, Y, kernel, likelihood, inference_method=None, name='gp'): + def __init__(self, X, Y, kernel, likelihood, inference_method=None, Y_metadata=None, name='gp'): super(GP, self).__init__(name) assert X.ndim == 2 @@ -38,6 +38,12 @@ class GP(Model): assert Y.shape[0] == self.num_data _, self.output_dim = self.Y.shape + if Y_metadata is not None: + assert Y_metadata.shape == self.Y.shape + self.Y_metadata = ObservableArray(Y_metadata) + else: + self.Y_metadata = None + assert isinstance(kernel, kern.kern) self.kern = kernel @@ -58,7 +64,7 @@ class GP(Model): self.parameters_changed() def parameters_changed(self): - self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y) + 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'] def log_likelihood(self): From 46ce76dee83580ceb3c59293a71aa891e1d8e026 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Wed, 12 Feb 2014 16:48:57 +0000 Subject: [PATCH 12/40] Fixed bernoulli likelihood divide by 0 and log of 0 --- GPy/examples/classification.py | 14 ++++++++----- .../latent_function_inference/__init__.py | 8 ++++---- .../latent_function_inference/laplace.py | 3 ++- GPy/likelihoods/bernoulli.py | 20 +++++++++++++++---- GPy/testing/likelihood_tests.py | 8 ++------ 5 files changed, 33 insertions(+), 20 deletions(-) diff --git a/GPy/examples/classification.py b/GPy/examples/classification.py index f9aaddd1..8637cc35 100644 --- a/GPy/examples/classification.py +++ b/GPy/examples/classification.py @@ -87,18 +87,22 @@ def toy_linear_1d_classification_laplace(seed=default_seed, optimize=True, plot= Y = data['Y'][:, 0:1] Y[Y.flatten() == -1] = 0 - bern_noise_model = GPy.likelihoods.bernoulli() - laplace_likelihood = GPy.likelihoods.Laplace(Y.copy(), bern_noise_model) + likelihood = GPy.likelihoods.Bernoulli() + laplace_inf = GPy.inference.latent_function_inference.Laplace() + kernel = GPy.kern.rbf(1) # Model definition - m = GPy.models.GPClassification(data['X'], Y, likelihood=laplace_likelihood) - print m + m = GPy.core.GP(data['X'], Y, kernel=kernel, likelihood=likelihood, inference_method=laplace_inf) # Optimize if optimize: #m.update_likelihood_approximation() # Parameters optimization: - m.optimize('bfgs', messages=1) + try: + m.optimize('scg', messages=1) + except Exception as e: + return m + #m.pseudo_EM() # Plot diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py index c89f771b..b114b8b9 100644 --- a/GPy/inference/latent_function_inference/__init__.py +++ b/GPy/inference/latent_function_inference/__init__.py @@ -9,12 +9,12 @@ prior over a finite set of points f. This prior is where K is the kernel matrix. We also have a likelihood (see GPy.likelihoods) which defines how the data are -related to the latent function: p(y | f). If the likelihood is also a Gaussian, -the inference over f is tractable (see exact_gaussian_inference.py). +related to the latent function: p(y | f). If the likelihood is also a Gaussian, +the inference over f is tractable (see exact_gaussian_inference.py). If the likelihood object is something other than Gaussian, then exact inference is not tractable. We then resort to a Laplace approximation (laplace.py) or -expectation propagation (ep.py). +expectation propagation (ep.py). The inference methods return a "Posterior" instance, which is a simple structure which contains a summary of the posterior. The model classes can then @@ -24,7 +24,7 @@ etc. """ from exact_gaussian_inference import ExactGaussianInference -from laplace import LaplaceInference +from laplace import Laplace expectation_propagation = 'foo' # TODO from dtc import DTC from fitc import FITC diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index 4edb9a1d..1b7bdad8 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -17,7 +17,7 @@ from posterior import Posterior import warnings from scipy import optimize -class LaplaceInference(object): +class Laplace(object): def __init__(self): """ @@ -52,6 +52,7 @@ class LaplaceInference(object): f_hat, Ki_fhat = self.rasm_mode(K, Y, likelihood, Ki_f_init, Y_metadata=Y_metadata) + self.f_hat = f_hat #Compute hessian and other variables at mode log_marginal, woodbury_vector, woodbury_inv, dL_dK, dL_dthetaL = self.mode_computations(f_hat, Ki_fhat, K, Y, likelihood, kern, Y_metadata) diff --git a/GPy/likelihoods/bernoulli.py b/GPy/likelihoods/bernoulli.py index 00626cd3..388ce173 100644 --- a/GPy/likelihoods/bernoulli.py +++ b/GPy/likelihoods/bernoulli.py @@ -116,7 +116,8 @@ class Bernoulli(Likelihood): Each y_i must be in {0, 1} """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape - objective = (link_f**y) * ((1.-link_f)**(1.-y)) + #objective = (link_f**y) * ((1.-link_f)**(1.-y)) + objective = np.where(y, link_f, 1.-link_f) return np.exp(np.sum(np.log(objective))) def logpdf_link(self, link_f, y, extra_data=None): @@ -136,7 +137,9 @@ class Bernoulli(Likelihood): """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape #objective = y*np.log(link_f) + (1.-y)*np.log(link_f) + state = np.seterr(divide='ignore') objective = np.where(y==1, np.log(link_f), np.log(1-link_f)) + np.seterr(**state) return np.sum(objective) def dlogpdf_dlink(self, link_f, y, extra_data=None): @@ -155,7 +158,10 @@ class Bernoulli(Likelihood): :rtype: Nx1 array """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape - grad = (y/link_f) - (1.-y)/(1-link_f) + #grad = (y/link_f) - (1.-y)/(1-link_f) + state = np.seterr(divide='ignore') + grad = np.where(y, 1./link_f, -1./(1-link_f)) + np.seterr(**state) return grad def d2logpdf_dlink2(self, link_f, y, extra_data=None): @@ -180,7 +186,10 @@ class Bernoulli(Likelihood): (the distribution for y_i depends only on link(f_i) not on link(f_(j!=i)) """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape - d2logpdf_dlink2 = -y/(link_f**2) - (1-y)/((1-link_f)**2) + #d2logpdf_dlink2 = -y/(link_f**2) - (1-y)/((1-link_f)**2) + state = np.seterr(divide='ignore') + d2logpdf_dlink2 = np.where(y, -1./np.square(link_f), -1./np.square(1.-link_f)) + np.seterr(**state) return d2logpdf_dlink2 def d3logpdf_dlink3(self, link_f, y, extra_data=None): @@ -199,7 +208,10 @@ class Bernoulli(Likelihood): :rtype: Nx1 array """ assert np.atleast_1d(link_f).shape == np.atleast_1d(y).shape - d3logpdf_dlink3 = 2*(y/(link_f**3) - (1-y)/((1-link_f)**3)) + #d3logpdf_dlink3 = 2*(y/(link_f**3) - (1-y)/((1-link_f)**3)) + state = np.seterr(divide='ignore') + d3logpdf_dlink3 = np.where(y, 2./(link_f**3), -2./((1.-link_f)**3)) + np.seterr(**state) return d3logpdf_dlink3 def samples(self, gp): diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index 9920d648..458831a0 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -516,9 +516,8 @@ class TestNoiseModels(object): Y = Y/Y.max() white_var = 1e-6 kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) - laplace_likelihood = GPy.inference.latent_function_inference.LaplaceInference() + laplace_likelihood = GPy.inference.latent_function_inference.Laplace() m = GPy.core.GP(X.copy(), Y.copy(), kernel, likelihood=model, inference_method=laplace_likelihood) - m.ensure_default_constraints() m['white'].constrain_fixed(white_var) #Set constraints @@ -555,7 +554,6 @@ class TestNoiseModels(object): kernel = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) ep_inf = GPy.inference.latent_function_inference.EP() m = GPy.core.GP(X.copy(), Y.copy(), kernel=kernel, likelihood=model, inference_method=ep_inf) - m.ensure_default_constraints() m['white'].constrain_fixed(white_var) for param_num in range(len(param_names)): @@ -644,13 +642,11 @@ class LaplaceTests(unittest.TestCase): m1['variance'] = initial_var_guess m1['variance'].constrain_bounded(1e-4, 10) m1['rbf'].constrain_bounded(1e-4, 10) - m1.ensure_default_constraints() m1.randomize() gauss_distr2 = GPy.likelihoods.Gaussian(variance=initial_var_guess) - laplace_inf = GPy.inference.latent_function_inference.LaplaceInference() + laplace_inf = GPy.inference.latent_function_inference.Laplace() m2 = GPy.core.GP(X, Y.copy(), kernel=kernel2, likelihood=gauss_distr2, inference_method=laplace_inf) - m2.ensure_default_constraints() m2['white'].constrain_fixed(1e-6) m2['rbf'].constrain_bounded(1e-4, 10) m2['variance'].constrain_bounded(1e-4, 10) From a264cdaa980c95af378c071e0b733c0a4f19c7fd Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Feb 2014 17:11:55 +0000 Subject: [PATCH 13/40] redid constraints --- GPy/core/model.py | 46 ++-- GPy/core/parameterization/array_core.py | 2 +- GPy/core/parameterization/index_operations.py | 105 ++++++++- GPy/core/parameterization/param.py | 14 +- GPy/core/parameterization/parameter_core.py | 93 ++++++-- GPy/core/parameterization/parameterized.py | 209 +++++++----------- GPy/core/parameterization/transformations.py | 5 +- GPy/models/bayesian_gplvm.py | 6 +- 8 files changed, 290 insertions(+), 190 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index aacf20be..748a5e08 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -17,10 +17,10 @@ import itertools # import numdifftools as ndt class Model(Parameterized): - _fail_count = 0 # Count of failed optimization steps (see objective) - _allowed_failures = 10 # number of allowed failures + _fail_count = 0 # Count of failed optimization steps (see objective) + _allowed_failures = 10 # number of allowed failures def __init__(self, name): - super(Model, self).__init__(name)#Parameterized.__init__(self) + super(Model, self).__init__(name) # Parameterized.__init__(self) self.priors = [] self._priors = ParameterIndexOperations() self.optimization_runs = [] @@ -30,10 +30,10 @@ class Model(Parameterized): 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) + # 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] + # [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) @@ -100,7 +100,7 @@ class Model(Parameterized): 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 + which = which[:1] # just place a prior object on the first parameter # check constraints are okay @@ -168,14 +168,14 @@ class Model(Parameterized): 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 = self._get_params_transformed() x = np.random.randn(self.size_transformed) x = self._untransform_params(x) # now draw from prior where possible if self.priors is not None and len(self.priors): [np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) 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...) + # 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): """ @@ -220,8 +220,8 @@ class Model(Parameterized): job = pool.apply_async(opt_wrapper, args=(self,), kwds=kwargs) jobs.append(job) - pool.close() # signal that no more data coming in - pool.join() # wait for all the tasks to complete + pool.close() # signal that no more data coming in + pool.join() # wait for all the tasks to complete except KeyboardInterrupt: print "Ctrl+c received, terminating and joining pool." pool.terminate() @@ -378,7 +378,7 @@ class Model(Parameterized): def optimize_SGD(self, momentum=0.1, learning_rate=0.01, iterations=20, **kwargs): # assert self.Y.shape[1] > 1, "SGD only works with D > 1" - sgd = SGD.StochasticGD(self, iterations, learning_rate, momentum, **kwargs) # @UndefinedVariable + sgd = SGD.StochasticGD(self, iterations, learning_rate, momentum, **kwargs) # @UndefinedVariable sgd.run() self.optimization_runs.append(sgd) @@ -412,7 +412,7 @@ class Model(Parameterized): gradient = self.objective_function_gradients(x) numerical_gradient = (f1 - f2) / (2 * dx) - global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient==0, 1e-32, gradient))) + 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) else: @@ -444,18 +444,18 @@ class Model(Parameterized): return gradient = self.objective_function_gradients(x) - np.where(gradient==0, 1e-312, gradient) + np.where(gradient == 0, 1e-312, gradient) ret = True for i, ind in enumerate(param_list): xx = x.copy() - xx[ind] += step + xx[i] += step f1 = self.objective_function(xx) - xx[ind] -= 2.*step + xx[i] -= 2.*step f2 = self.objective_function(xx) numerical_gradient = (f1 - f2) / (2 * step) - ratio = (f1 - f2) / (2 * step * gradient[ind]) - difference = np.abs((f1 - f2) / 2 / step - gradient[ind]) + ratio = (f1 - f2) / (2 * step * gradient[i]) + difference = np.abs((f1 - f2) / 2 / step - gradient[i]) if (np.abs(1. - ratio) < tolerance) or np.abs(difference) < tolerance: formatted_name = "\033[92m {0} \033[0m".format(names[ind]) @@ -466,7 +466,7 @@ class Model(Parameterized): r = '%.6f' % float(ratio) d = '%.6f' % float(difference) - g = '%.6f' % gradient[ind] + g = '%.6f' % gradient[i] ng = '%.6f' % float(numerical_gradient) grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4]) print grad_string @@ -517,10 +517,10 @@ class Model(Parameterized): alpha = 0 stop = False - #Handle **kwargs + # Handle **kwargs ep_args = {} for arg in kwargs.keys(): - if arg in ('epsilon','power_ep'): + if arg in ('epsilon', 'power_ep'): ep_args[arg] = kwargs[arg] del kwargs[arg] @@ -528,7 +528,7 @@ class Model(Parameterized): 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']) + 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']) @@ -540,8 +540,8 @@ class Model(Parameterized): 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 + 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: diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index a65b8c24..975b78f4 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, Constrainable, Gradcheckable +from parameter_core import Observable, Parameterizable class ParamList(list): """ diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index 8abb31e9..da7f1d9b 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -58,7 +58,7 @@ class ParameterIndexOperations(object): index array, for multi-param handling. ''' def __init__(self): - self._properties = ParamDict() + self._properties = IntArrayDict() #self._reverse = collections.defaultdict(list) def __getstate__(self): @@ -71,16 +71,19 @@ class ParameterIndexOperations(object): def iteritems(self): return self._properties.iteritems() + def items(self): + return self._properties.items() + def properties(self): return self._properties.keys() - def iter_properties(self): + def iterproperties(self): return self._properties.iterkeys() def shift(self, start, size): for ind in self.iterindices(): toshift = ind>=start - if len(toshift) > 0: + if toshift.size > 0: ind[toshift] += size def clear(self): @@ -96,7 +99,7 @@ class ParameterIndexOperations(object): return self._properties.values() def properties_for(self, index): - return vectorize(lambda i: [prop for prop in self.iter_properties() if i in self._properties[prop]], otypes=[list])(index) + return vectorize(lambda i: [prop for prop in self.iterproperties() if i in self[prop]], otypes=[list])(index) def add(self, prop, indices): try: @@ -114,8 +117,13 @@ class ParameterIndexOperations(object): del self._properties[prop] return removed.astype(int) return numpy.array([]).astype(int) + def __getitem__(self, prop): return self._properties[prop] + + def __str__(self, *args, **kwargs): + import pprint + return pprint.pformat(dict(self._properties)) def combine_indices(arr1, arr2): return numpy.union1d(arr1, arr2) @@ -126,5 +134,94 @@ def remove_indices(arr, to_remove): def index_empty(index): return numpy.size(index) == 0 +class ParameterIndexOperationsView(object): + def __init__(self, param_index_operations, offset, size): + self._param_index_ops = param_index_operations + self._offset = offset + self._size = size + + def __getstate__(self): + return [self._param_index_ops, self._offset, self._size] + def __setstate__(self, state): + self._param_index_ops = state[0] + self._offset = state[1] + self._size = state[2] + + + def _filter_index(self, ind): + return ind[(ind >= self._offset) * (ind < (self._offset + self._size))] - self._offset + + + def iteritems(self): + for i, ind in self._param_index_ops.iteritems(): + ind2 = self._filter_index(ind) + if ind2.size > 0: + yield i, ind2 + + def items(self): + return [[i,v] for i,v in self.iteritems()] + + def properties(self): + return [i for i in self.iterproperties()] + + + def iterproperties(self): + for i, _ in self.iteritems(): + yield i + + + def shift(self, start, size): + raise NotImplementedError, 'Shifting only supported in original ParamIndexOperations' + + + def clear(self): + for i, ind in self.items(): + self._param_index_ops.remove(i, ind+self._offset) + + + def size(self): + return reduce(lambda a,b: a+b.size, self.iterindices(), 0) + + + def iterindices(self): + for _, ind in self.iteritems(): + yield ind + + + def indices(self): + [ind for ind in self.iterindices()] + + + def properties_for(self, index): + return vectorize(lambda i: [prop for prop in self.iterproperties() if i in self[prop]], otypes=[list])(index) + + + def add(self, prop, indices): + self._param_index_ops.add(prop, indices+self._offset) + + + def remove(self, prop, indices): + removed = self._param_index_ops.remove(prop, indices+self._offset) + if removed.size > 0: + return removed - self._size + return removed + + + def __getitem__(self, prop): + ind = self._filter_index(self._param_index_ops[prop]) + if ind.size > 0: + return ind + raise KeyError, prop + + def __str__(self, *args, **kwargs): + import pprint + return pprint.pformat(dict(self.iteritems())) + + def update(self, parameter_index_view): + for i, v in parameter_index_view.iteritems(): + self.add(i, v) + + pass + diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 0f0ef667..56b8aab6 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -3,7 +3,7 @@ import itertools import numpy -from parameter_core import Constrainable, Gradcheckable, adjust_name_for_printing +from parameter_core import Constrainable, Gradcheckable, Indexable, Parameterizable, adjust_name_for_printing from array_core import ObservableArray, ParamList ###### printing @@ -14,13 +14,7 @@ __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision __print_threshold__ = 5 ###### -class Float(numpy.float64, Constrainable): - def __init__(self, f, base): - super(Float,self).__init__(f) - self._base = base - - -class Param(ObservableArray, Constrainable, Gradcheckable): +class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameterizable): """ Parameter object for GPy models. @@ -364,7 +358,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable): return [self.shape] @property def _constraints_str(self): - return [' '.join(map(lambda c: str(c[0]) if c[1].size == self._realsize_ else "{" + str(c[0]) + "}", self._highest_parent_._constraints_iter_items(self)))] + return [' '.join(map(lambda c: str(c[0]) if c[1].size == self._realsize_ else "{" + str(c[0]) + "}", self.constraints.iteritems()))] @property def _ties_str(self): return [t._short() for t in self._tied_to_] or [''] @@ -390,7 +384,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable): else: ties[i, matches[0]] = numpy.take(tt_rav_index, matches[1], mode='wrap') return map(lambda a: sum(a, []), zip(*[[[tie.flatten()] if tx != None else [] for tx in t] for t, tie in zip(ties, self._tied_to_)])) def _constraints_for(self, rav_index): - return self._highest_parent_._constraints_for(self, rav_index) + return self.constraints.properties_for(rav_index) def _indices(self, slice_index=None): # get a int-array containing all indices in the first axis. if slice_index is None: diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 1075d808..1c507688 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -1,7 +1,7 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -from transformations import Logexp, NegativeLogexp, Logistic +from transformations import Transformation, Logexp, NegativeLogexp, Logistic __updated__ = '2013-12-16' @@ -10,6 +10,11 @@ def adjust_name_for_printing(name): return name.replace(" ", "_").replace(".", "_").replace("-","").replace("+","").replace("!","").replace("*","").replace("/","") return '' +#=============================================================================== +# Printing: +__fixed__ = "fixed" +#=============================================================================== + class Observable(object): _observers_ = {} def add_observer(self, observer, callble): @@ -20,7 +25,23 @@ class Observable(object): def _notify_observers(self): [callble(self) for callble in self._observers_.itervalues()] - +class Parameterizable(object): + def __init__(self, *args, **kwargs): + from GPy.core.parameterization.array_core import ParamList + _parameters_ = ParamList() + + def parameter_names(self): + return [p.name for p in self._parameters_] + + def parameters_changed(self): + """ + This method gets called when parameters have changed. + Another way of listening to param changes is to + add self as a listener to the param, such that + updates get passed through. See :py:function:``GPy.core.param.Observable.add_observer`` + """ + pass + class Pickleable(object): def _getstate(self): """ @@ -52,7 +73,7 @@ class Parentable(object): super(Parentable,self).__init__() self._direct_parent_ = direct_parent self._parent_index_ = parent_index - + def has_parent(self): return self._direct_parent_ is not None @@ -89,11 +110,22 @@ class Gradcheckable(Parentable): def _checkgrad(self, param): raise NotImplementedError, "Need log likelihood to check gradient against" +class Indexable(object): + def _raveled_index(self): + raise NotImplementedError, "Need to be able to get the raveled Index" + + def _internal_offset(self): + return 0 + + def _offset_for(self, param): + raise NotImplementedError, "shouldnt happen, offset required from non parameterization object?" -class Constrainable(Nameable): +class Constrainable(Nameable, Indexable, Parameterizable): def __init__(self, name, default_constraint=None): super(Constrainable,self).__init__(name) self._default_constraint_ = default_constraint + from index_operations import ParameterIndexOperations + self.constraints = ParameterIndexOperations() #=========================================================================== # Fixing Parameters: #=========================================================================== @@ -105,17 +137,28 @@ class Constrainable(Nameable): """ if value is not None: self[:] = value - self._highest_parent_._fix(self,warning) + self.constrain(__fixed__, warning=warning) + self._highest_parent_._set_fixed(self._raveled_index()) fix = constrain_fixed def unconstrain_fixed(self): """ This parameter will no longer be fixed. """ - self._highest_parent_._unfix(self) + unconstrained = self.unconstrain(__fixed__) + self._highest_parent_._set_unfixed(unconstrained) unfix = unconstrain_fixed #=========================================================================== # Constrain operations -> done #=========================================================================== + def _parent_changed(self, parent): + c = self.constraints + from index_operations import ParameterIndexOperationsView + self.constraints = ParameterIndexOperationsView(parent.constraints, parent._offset_for(self), self.size) + self.constraints.update(c) + del c + for p in self._parameters_: + p._parent_changed(parent) + def constrain(self, transform, warning=True, update=True): """ :param transform: the :py:class:`GPy.core.transformations.Transformation` @@ -125,15 +168,21 @@ class Constrainable(Nameable): Constrain the parameter to the given :py:class:`GPy.core.transformations.Transformation`. """ - if self.has_parent(): - self._highest_parent_._add_constrain(self, transform, warning) - if update: - self._highest_parent_.parameters_changed() - else: - for p in self._parameters_: - self._add_constrain(p, transform, warning) - if update: - self.parameters_changed() + if isinstance(transform, Transformation): + self._set_params(transform.initialize(self._get_params()), update=False) + reconstrained = self.unconstrain() + self.constraints.add(transform, self._raveled_index()) + if reconstrained.size > 0: + print "WARNING: reconstraining parameters {}".format(self.parameter_names) + if update: + self._highest_parent_.parameters_changed() + # if self.has_parent(): + # self._highest_parent_._add_constrain(self, transform, warning) + # else: + # for p in self._parameters_: + # self._add_constrain(p, transform, warning) + # if update: + # self.parameters_changed() def constrain_positive(self, warning=True, update=True): """ @@ -167,12 +216,14 @@ class Constrainable(Nameable): remove all :py:class:`GPy.core.transformations.Transformation` transformats of this parameter object. """ - if self.has_parent(): - self._highest_parent_._remove_constrain(self, *transforms) - else: - for p in self._parameters_: - self._remove_constrain(p, *transforms) - + if len(transforms) == 0: + transforms = self.constraints.properties() + import numpy as np + removed = np.empty((0,),dtype=int) + for t in transforms: + removed = np.intersect1d(removed, self.constraints.remove(t, self._raveled_index())) + return removed + def unconstrain_positive(self): """ Remove positive constraint of this parameter. diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 85a1e179..a30a1e3b 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -8,16 +8,9 @@ import cPickle import itertools from re import compile, _pattern_type from param import ParamConcatenation, Param -from parameter_core import Constrainable, Pickleable, Observable, adjust_name_for_printing, Gradcheckable -from index_operations import ParameterIndexOperations,\ - index_empty +from parameter_core import Constrainable, Pickleable, Observable, adjust_name_for_printing, Gradcheckable, __fixed__ from array_core import ParamList -#=============================================================================== -# Printing: -__fixed__ = "fixed" -#=============================================================================== - #=============================================================================== # constants FIXED = False @@ -69,7 +62,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): def __init__(self, name=None): super(Parameterized, self).__init__(name=name) self._in_init_ = True - self._constraints_ = None#ParameterIndexOperations() self._parameters_ = ParamList() self.size = sum(p.size for p in self._parameters_) if not self._has_fixes(): @@ -79,11 +71,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_ = set() del self._in_init_ - @property - def constraints(self): - if self._constraints_ is None: - self._constraints_ = ParameterIndexOperations() - return self._constraints_ #=========================================================================== # Parameter connection for model creation: #=========================================================================== @@ -128,12 +115,14 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): Add all parameters to this param class, you can insert parameters at any given index using the :func:`list.insert` syntax """ + # if param.has_parent(): + # raise AttributeError, "parameter {} already in another model, create new object (or copy) for adding".format(param._short()) if param in self._parameters_ and index is not None: # make sure fixes and constraints are indexed right if self._has_fixes(): - param_slice = slice(self._offset_for(param),self._offset_for(param)+param.size) + param_slice = slice(self._offset_for(param), self._offset_for(param) + param.size) dest_index = sum((p.size for p in self._parameters_[:index])) - dest_slice = slice(dest_index,dest_index+param.size) + dest_slice = slice(dest_index, dest_index + param.size) fixes_param = self._fixes_[param_slice].copy() self._fixes_[param_slice] = self._fixes_[dest_slice] self._fixes_[dest_slice] = fixes_param @@ -164,22 +153,18 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): ins = sum((p.size for p in self._parameters_[:index])) if self._has_fixes(): self._fixes_ = np.r_[self._fixes_[:ins], fixes_param, self._fixes[ins:]] elif not np.all(fixes_param): - self._fixes_ = np.ones(self.size+param.size, dtype=bool) - self._fixes_[ins:ins+param.size] = fixes_param + self._fixes_ = np.ones(self.size + param.size, dtype=bool) + self._fixes_[ins:ins + param.size] = fixes_param self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" self._connect_parameters() - # make sure the constraints are pulled over: - if hasattr(param, "_constraints_") and param._constraints_ is not None: - for t, ind in param._constraints_.iteritems(): - - self.constraints.add(t, ind+self._offset_for(param)) - param._constraints_.clear() + for p in self._parameters_: + p._parent_changed(self) if param._default_constraint_ is not None: - self._add_constrain(param, param._default_constraint_, False) - if self._has_fixes() and np.all(self._fixes_): # ==UNFIXED - self._fixes_= None + param.constrain(param._default_constraint_, False) + if self._has_fixes() and np.all(self._fixes_): # ==UNFIXED + self._fixes_ = None def add_parameters(self, *parameters): """ @@ -202,30 +187,22 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): or p in names_params_indices)]) self._connect_parameters() - def parameters_changed(self): - """ - This method gets called when parameters have changed. - Another way of listening to param changes is to - add self as a listener to the param, such that - updates get passed through. See :py:function:``GPy.core.param.Observable.add_observer`` - """ - # will be called as soon as parameters have changed - pass - def _connect_parameters(self): # connect parameterlist to this parameterized object # This just sets up the right connection for the params objects # to be used as parameters + # it also sets the constraints for each parameter to the constraints + # of their respective parents if not hasattr(self, "_parameters_") or len(self._parameters_) < 1: # no parameters for this class return sizes = [0] self._param_slices_ = [] - for i,p in enumerate(self._parameters_): + for i, p in enumerate(self._parameters_): p._direct_parent_ = self p._parent_index_ = i not_unique = [] - sizes.append(p.size+sizes[-1]) + sizes.append(p.size + sizes[-1]) self._param_slices_.append(slice(sizes[-2], sizes[-1])) pname = adjust_name_for_printing(p.name) # and makes sure to not delete programmatically added parameters @@ -237,7 +214,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): elif not (pname in not_unique): self.__dict__[pname] = p self._added_names_.add(pname) - #=========================================================================== # Pickling operations #=========================================================================== @@ -255,16 +231,16 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): cPickle.dump(self, f, protocol) def copy(self): """Returns a (deep) copy of the current model """ - #dc = dict() - #for k, v in self.__dict__.iteritems(): - #if k not in ['_highest_parent_', '_direct_parent_']: - #dc[k] = copy.deepcopy(v) + # dc = dict() + # for k, v in self.__dict__.iteritems(): + # if k not in ['_highest_parent_', '_direct_parent_']: + # dc[k] = copy.deepcopy(v) - #dc = copy.deepcopy(self.__dict__) - #dc['_highest_parent_'] = None - #dc['_direct_parent_'] = None - #s = self.__class__.new() - #s.__dict__ = dc + # dc = copy.deepcopy(self.__dict__) + # dc['_highest_parent_'] = None + # dc['_direct_parent_'] = None + # s = self.__class__.new() + # s.__dict__ = dc return copy.deepcopy(self) def __getstate__(self): if self._has_get_set_state(): @@ -272,8 +248,8 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self.__dict__ def __setstate__(self, state): if self._has_get_set_state(): - self._setstate(state) # set state - #self._set_params(self._get_params()) # restore all values + self._setstate(state) # set state + # self._set_params(self._get_params()) # restore all values return self.__dict__ = state def _has_get_set_state(self): @@ -289,7 +265,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): """ return [ self._fixes_, - self._constraints_, + self.constraints, self._parameters_, self._name, self._added_names_, @@ -300,7 +276,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._name = state.pop() self._parameters_ = state.pop() self._connect_parameters() - self._constraints_ = state.pop() + self.constraints = state.pop() self._fixes_ = state.pop() self.parameters_changed() #=========================================================================== @@ -310,9 +286,9 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): if self.has_parent(): return g x = self._get_params() - [numpy.put(g, i, g[i]*c.gradfactor(x[i])) for c,i in self.constraints.iteritems() if c != __fixed__] + [numpy.put(g, i, g[i] * c.gradfactor(x[i])) for c, i in self.constraints.iteritems() if c != __fixed__] for p in self.flattened_parameters: - for t,i in p._tied_to_me_.iteritems(): + for t, i in p._tied_to_me_.iteritems(): 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 @@ -320,7 +296,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): # Optimization handles: #=========================================================================== def _get_param_names(self): - n = numpy.array([p.name_hirarchical+'['+str(i)+']' for p in self.flattened_parameters for i in p._indices()]) + n = numpy.array([p.name_hirarchical + '[' + 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() @@ -331,16 +307,16 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): # don't overwrite this anymore! if not self.size: return np.empty(shape=(0,), dtype=np.float64) - return numpy.hstack([x._get_params() for x in self._parameters_ if x.size>0]) + return numpy.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! - [p._set_params(params[s], update=update) for p,s in itertools.izip(self._parameters_,self._param_slices_)] + [p._set_params(params[s], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] self.parameters_changed() 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__] + [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 @@ -350,7 +326,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): 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__] + [numpy.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] return p def _name_changed(self, param, old_name): if hasattr(self, old_name) and old_name in self._added_names_: @@ -365,7 +341,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): #=========================================================================== 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 - self._offset_for(param) ind = ind[ind >= 0] internal_offset = param._internal_offset() ind = ind[ind < param.size + internal_offset] @@ -377,7 +353,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self._param_slices_[param._direct_parent_._get_original(param)._parent_index_].start return self._offset_for(param._direct_parent_) + param._direct_parent_._offset_for(param) return 0 - + def _raveled_index_for(self, param): """ get the raveled index for a param @@ -387,7 +363,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): if isinstance(param, ParamConcatenation): return numpy.hstack((self._raveled_index_for(p) for p in param.params)) return param._raveled_index() + self._offset_for(param) - + def _raveled_index(self): """ get the raveled index for this object, @@ -404,7 +380,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): except AttributeError: pass self._fixes_[param_or_index] = FIXED - if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED + if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED def _set_unfixed(self, param_or_index): if not self._has_fixes(): self._fixes_ = numpy.ones(self.size, dtype=bool) try: @@ -415,18 +391,18 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): for constr, ind in self.constraints.iteritems(): if constr is __fixed__: self._fixes_[ind] = FIXED - if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED + if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED 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)] - def _fix(self, param, warning=True): - f = self._add_constrain(param, __fixed__, warning) - self._set_fixed(f) - def _unfix(self, param): - if self._has_fixes(): - f = self._remove_constrain(param, __fixed__) - self._set_unfixed(f) + # def _fix(self, param, warning=True): + # f = self._add_constrain(param, __fixed__, warning) + # self._set_fixed(f) + # def _unfix(self, param): + # if self._has_fixes(): + # f = self._remove_constrain(param, __fixed__) + # self._set_unfixed(f) #=========================================================================== # Convenience for fixed, tied checking of param: #=========================================================================== @@ -437,7 +413,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): 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() + # 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_: @@ -455,54 +431,33 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): #=========================================================================== # Constraint Handling: #=========================================================================== - def _add_constrain(self, param, transform, warning=True): - rav_i = self._raveled_index_for(param) - reconstrained = self._remove_constrain(param, index=rav_i) # remove constraints before - # if removing constraints before adding new is not wanted, just delete the above line! - self.constraints.add(transform, rav_i) - param = self._get_original(param) - if not (transform == __fixed__): - param._set_params(transform.initialize(param._get_params()), update=False) - if warning and any(reconstrained): - # if you want to print the whole params object, which was reconstrained use: - # m = str(param[self._backtranslate_index(param, reconstrained)]) - print "Warning: re-constraining parameters:\n{}".format(param._short()) - return rav_i - def _remove_constrain(self, param, *transforms, **kwargs): - if not transforms: - transforms = self.constraints.properties() - removed_indices = numpy.array([]).astype(int) - if "index" in kwargs: index = kwargs['index'] - else: index = self._raveled_index_for(param) - for constr in transforms: - removed = self.constraints.remove(constr, index) - if constr is __fixed__: - self._set_unfixed(removed) - removed_indices = numpy.union1d(removed_indices, removed) - return removed_indices - # convienience for iterating over items - def _constraints_iter_items(self, param): - for constr, ind in self.constraints.iteritems(): - ind = self._backtranslate_index(param, ind) - if not index_empty(ind): - yield constr, ind - def _constraints_iter(self, param): - for constr, _ in self._constraints_iter_items(param): - yield constr - def _contraints_iter_indices(self, param): - # iterate through all constraints belonging to param - for _, ind in self._constraints_iter_items(param): - yield ind - def _constraint_indices(self, param, constraint): - # indices in model range for param and constraint - return self._backtranslate_index(param, self.constraints[constraint]) + self._offset_for(param) - def _constraints_for(self, param, rav_index): - # constraint for param given its internal rav_index - return self.constraints.properties_for(rav_index+self._offset_for(param)) - def _constraints_for_collect(self, param, rav_index): - # constraint for param given its internal rav_index - cs = self._constraints_for(param, rav_index) - return set(itertools.chain(*cs)) + #=========================================================================== + # def _add_constrain(self, param, transform, warning=True): + # rav_i = self._raveled_index_for(param) + # reconstrained = self._remove_constrain(param, index=rav_i) # remove constraints before + # # if removing constraints before adding new is not wanted, just delete the above line! + # self.constraints.add(transform, rav_i) + # param = self._get_original(param) + # if not (transform == __fixed__): + # param._set_params(transform.initialize(param._get_params()), update=False) + # if warning and any(reconstrained): + # # if you want to print the whole params object, which was reconstrained use: + # # m = str(param[self._backtranslate_index(param, reconstrained)]) + # print "Warning: re-constraining parameters:\n{}".format(param._short()) + # return rav_i + # def _remove_constrain(self, param, *transforms, **kwargs): + # if not transforms: + # transforms = self.constraints.properties() + # removed_indices = numpy.array([]).astype(int) + # if "index" in kwargs: index = kwargs['index'] + # else: index = self._raveled_index_for(param) + # for constr in transforms: + # removed = self.constraints.remove(constr, index) + # if constr is __fixed__: + # self._set_unfixed(removed) + # removed_indices = numpy.union1d(removed_indices, removed) + # return removed_indices + #=========================================================================== #=========================================================================== # Get/set parameters: #=========================================================================== @@ -539,7 +494,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): # def __getattribute__(self, name): # #try: # return object.__getattribute__(self, name) - #except AttributeError: + # except AttributeError: # _, a, tb = sys.exc_info() # try: # return self.__getitem__(name) @@ -592,22 +547,22 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return [','.join(x._ties_str) for x in self.flattened_parameters] def __str__(self, header=True): - name = adjust_name_for_printing(self.name) + "." + name = adjust_name_for_printing(self.name) + "." constrs = self._constraints_str; ts = self._ties_str desc = self._description_str; names = self.parameter_names nl = max([len(str(x)) for x in names + [name]]) sl = max([len(str(x)) for x in desc + ["Value"]]) - cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]]) + cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]]) tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]]) format_spec = " \033[1m{{name:<{0}s}}\033[0;0m | {{desc:^{1}s}} | {{const:^{2}s}} | {{t:^{3}s}}".format(nl, sl, cl, tl) to_print = [] for n, d, c, t in itertools.izip(names, desc, constrs, ts): to_print.append(format_spec.format(name=n, desc=d, const=c, t=t)) - #to_print = [format_spec.format(p=p, const=c, t=t) if isinstance(p, Param) else p.__str__(header=False) for p, c, t in itertools.izip(self._parameters_, constrs, ts)] - sep = '-'*(nl+sl+cl+tl+8*2+3) + # to_print = [format_spec.format(p=p, const=c, t=t) if isinstance(p, Param) else p.__str__(header=False) for p, c, t in itertools.izip(self._parameters_, constrs, ts)] + sep = '-' * (nl + sl + cl + tl + 8 * 2 + 3) if header: header = " {{0:<{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}}".format(nl, sl, cl, tl).format(name, "Value", "Constraint", "Tied to") - #header += '\n' + sep + # header += '\n' + sep to_print.insert(0, header) return '\n'.format(sep).join(to_print) pass diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index c4cab1e9..e17da404 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -27,6 +27,8 @@ class Transformation(object): raise NotImplementedError def __str__(self): raise NotImplementedError + def __repr__(self): + return self.__class__.__name__ class Logexp(Transformation): domain = _POSITIVE @@ -56,7 +58,7 @@ class NegativeLogexp(Transformation): return -self.logexp.initialize(f) # np.abs(f) def __str__(self): return '-ve' - + class LogexpClipped(Logexp): max_bound = 1e100 min_bound = 1e-10 @@ -94,7 +96,6 @@ class LogexpClipped(Logexp): def __str__(self): return '+ve_c' - class Exponent(Transformation): # TODO: can't allow this to go to zero, need to set a lower bound. Similar with negative Exponent below. See old MATLAB code. domain = _POSITIVE diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 62d9a5a9..acb04350 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -23,7 +23,7 @@ class BayesianGPLVM(SparseGP, GPLVM): """ def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10, - Z=None, kernel=None, inference_method=None, likelihood=Gaussian(), name='bayesian gplvm', **kwargs): + Z=None, kernel=None, inference_method=None, likelihood=None, name='bayesian gplvm', **kwargs): if X == None: X = self.initialise_latent(init, input_dim, Y) self.init = init @@ -37,7 +37,9 @@ class BayesianGPLVM(SparseGP, GPLVM): if kernel is None: kernel = kern.rbf(input_dim) # + kern.white(input_dim) - + + if likelihood is None: + likelihood = Gaussian() self.q = Normal(X, X_variance) SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, X_variance, name, **kwargs) self.add_parameter(self.q, index=0) From 646400f49e3587e9baace963ebcfc4f028c47a6b Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 12 Feb 2014 17:30:20 +0000 Subject: [PATCH 14/40] fixed gradchecker and fixes for paramterized --- GPy/core/model.py | 2 +- GPy/core/parameterization/index_operations.py | 2 +- GPy/core/parameterization/parameter_core.py | 4 +++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 748a5e08..f2632c93 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -437,7 +437,7 @@ class Model(Parameterized): else: param_list = self._raveled_index_for(target_param) if self._has_fixes(): - param_list = np.intersect1d(param_list, np.r_[:self.size][self._fixes_], True) + param_list = np.intersect1d(np.r_[:self.size][self._fixes_], param_list, True) if param_list.size == 0: print "No free parameters to check" diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index da7f1d9b..d2549dad 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -205,7 +205,7 @@ class ParameterIndexOperationsView(object): def remove(self, prop, indices): removed = self._param_index_ops.remove(prop, indices+self._offset) if removed.size > 0: - return removed - self._size + return removed - self._size + 1 return removed diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 1c507688..6afea470 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -145,7 +145,9 @@ class Constrainable(Nameable, Indexable, Parameterizable): This parameter will no longer be fixed. """ unconstrained = self.unconstrain(__fixed__) + import ipdb;ipdb.set_trace() self._highest_parent_._set_unfixed(unconstrained) + unfix = unconstrain_fixed #=========================================================================== # Constrain operations -> done @@ -221,7 +223,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): import numpy as np removed = np.empty((0,),dtype=int) for t in transforms: - removed = np.intersect1d(removed, self.constraints.remove(t, self._raveled_index())) + removed = np.union1d(removed, self.constraints.remove(t, self._raveled_index())) return removed def unconstrain_positive(self): From d2e8807a88a1a6830c6fcf3a6ecb251579eaea78 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Wed, 12 Feb 2014 18:02:24 +0000 Subject: [PATCH 15/40] Fixed some examples and tests, and stated that Y metadata doesnt need to be the same size as Y --- GPy/core/gp.py | 1 - GPy/examples/non_gaussian.py | 40 ++++++++++++++------------------- GPy/testing/likelihood_tests.py | 4 ++++ 3 files changed, 21 insertions(+), 24 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 214c2324..ded4a9f5 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -39,7 +39,6 @@ class GP(Model): _, self.output_dim = self.Y.shape if Y_metadata is not None: - assert Y_metadata.shape == self.Y.shape self.Y_metadata = ObservableArray(Y_metadata) else: self.Y_metadata = None diff --git a/GPy/examples/non_gaussian.py b/GPy/examples/non_gaussian.py index 23122691..2a5e0c42 100644 --- a/GPy/examples/non_gaussian.py +++ b/GPy/examples/non_gaussian.py @@ -42,38 +42,35 @@ def student_t_approx(optimize=True, plot=True): kernel4 = GPy.kern.rbf(X.shape[1]) + GPy.kern.white(X.shape[1]) #Gaussian GP model on clean data - m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1) - # optimize - m1.ensure_default_constraints() - m1['white'] = 1e-5 - m1['white'].constrain_fixed('white') - m1.randomize() + #m1 = GPy.models.GPRegression(X, Y.copy(), kernel=kernel1) + ## optimize + #m1['white'].constrain_fixed(1e-5) + #m1.randomize() - #Gaussian GP model on corrupt data - m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2) - m2.ensure_default_constraints() - m1['white'] = 1e-5 - m1['white'].constrain_fixed('white') - m2.randomize() + ##Gaussian GP model on corrupt data + #m2 = GPy.models.GPRegression(X, Yc.copy(), kernel=kernel2) + #m1['white'].constrain_fixed(1e-5) + #m2.randomize() #Student t GP model on clean data t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd) - laplace_inf = GPy.inference.latent_function_inference.LaplaceInference() + laplace_inf = GPy.inference.latent_function_inference.Laplace() m3 = GPy.core.GP(X, Y.copy(), kernel3, likelihood=t_distribution, inference_method=laplace_inf) - m3.ensure_default_constraints() m3['t_noise'].constrain_bounded(1e-6, 10.) - m3['white'] = 1e-5 - m3['white'].constrain_fixed() + m3['white'].constrain_fixed(1e-5) m3.randomize() + debug = True + print m3 + if debug: + m3.optimize(messages=1) + return m3 #Student t GP model on corrupt data t_distribution = GPy.likelihoods.StudentT(deg_free=deg_free, sigma2=edited_real_sd) - laplace_inf = GPy.inference.latent_function_inference.LaplaceInference() + laplace_inf = GPy.inference.latent_function_inference.Laplace() m4 = GPy.core.GP(X, Yc.copy(), kernel4, likelihood=t_distribution, inference_method=laplace_inf) - m4.ensure_default_constraints() m4['t_noise'].constrain_bounded(1e-6, 10.) - m4['white'] = 1e-5 - m4['white'].constrain_fixed() + m4['white'].constrain_fixed(1e-5) m4.randomize() if optimize: @@ -156,7 +153,6 @@ def boston_example(optimize=True, plot=True): #Gaussian GP print "Gauss GP" mgp = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelgp.copy()) - mgp.ensure_default_constraints() mgp.constrain_fixed('white', 1e-5) mgp['rbf_len'] = rbf_len mgp['noise'] = noise @@ -174,7 +170,6 @@ def boston_example(optimize=True, plot=True): g_distribution = GPy.likelihoods.noise_model_constructors.gaussian(variance=noise, N=N, D=D) g_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), g_distribution) mg = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=g_likelihood) - mg.ensure_default_constraints() mg.constrain_positive('noise_variance') mg.constrain_fixed('white', 1e-5) mg['rbf_len'] = rbf_len @@ -194,7 +189,6 @@ def boston_example(optimize=True, plot=True): t_distribution = GPy.likelihoods.noise_model_constructors.student_t(deg_free=df, sigma2=noise) stu_t_likelihood = GPy.likelihoods.Laplace(Y_train.copy(), t_distribution) mstu_t = GPy.models.GPRegression(X_train.copy(), Y_train.copy(), kernel=kernelstu.copy(), likelihood=stu_t_likelihood) - mstu_t.ensure_default_constraints() mstu_t.constrain_fixed('white', 1e-5) mstu_t.constrain_bounded('t_noise', 0.0001, 1000) mstu_t['rbf_len'] = rbf_len diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index 458831a0..a70073e4 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -540,6 +540,10 @@ class TestNoiseModels(object): #import ipdb; ipdb.set_trace() #NOTE this test appears to be stochastic for some likelihoods (student t?) # appears to all be working in test mode right now... + + if not m.checkgrad(): + import ipdb; ipdb.set_trace() # XXX BREAKPOINT + assert m.checkgrad(step=step) ########### From 3c68a9688d8beb5b85eb2484aa630ffd198a7ca1 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 13 Feb 2014 08:53:14 +0000 Subject: [PATCH 16/40] general bugfixing --- GPy/core/gp.py | 5 +++-- GPy/core/sparse_gp.py | 20 ++++++++++--------- .../latent_function_inference/__init__.py | 1 + GPy/kern/parts/rbf.py | 2 +- GPy/likelihoods/gaussian.py | 11 +++++++--- GPy/models/sparse_gp_regression.py | 4 ++-- 6 files changed, 26 insertions(+), 17 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index ab725897..b0b25319 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -55,7 +55,8 @@ class GP(Model): self.add_parameter(self.kern) self.add_parameter(self.likelihood) - self.parameters_changed() + if self.__class__ is GP: + self.parameters_changed() def parameters_changed(self): self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y) @@ -111,7 +112,7 @@ class GP(Model): This is to allow for different normalizations of the output dimensions. """ - # normalize X values + #predict the latent function values mu, var = self._raw_predict(Xnew, full_cov=full_cov, which_parts=which_parts) # now push through likelihood diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index bf28d062..bbbea2fe 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -39,7 +39,7 @@ class SparseGP(GP): if isinstance(likelihood, likelihoods.Gaussian): inference_method = varDTC.VarDTC() else: - #inference_method = ?? + #inference_method = ?? raise NotImplementedError, "what to do what to do?" print "defaulting to ", inference_method, "for latent function inference" @@ -51,19 +51,21 @@ class SparseGP(GP): self.X_variance = X_variance GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name) - self.add_parameter(self.Z, index=0) + self.parameters_changed() + def parameters_changed(self): 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) - #The derivative of the bound wrt the inducing inputs Z - self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) - if self.X_variance is None: - self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) - else: - self.Z.gradient += self.kern.dpsi1_dZ(self.grad_dict['dL_dpsi1'], self.Z, self.X, self.X_variance) - self.Z.gradient += self.kern.dpsi2_dZ(self.grad_dict['dL_dpsi2'], self.Z, self.X, self.X_variance) + #The derivative of the bound wrt the inducing inputs Z ( unless they're all fixed) + if not self.Z.is_fixed: + self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) + if self.X_variance is None: + self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) + else: + self.Z.gradient += self.kern.dpsi1_dZ(self.grad_dict['dL_dpsi1'], self.Z, self.X, self.X_variance) + self.Z.gradient += self.kern.dpsi2_dZ(self.grad_dict['dL_dpsi2'], self.Z, self.X, self.X_variance) def _raw_predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): """ diff --git a/GPy/inference/latent_function_inference/__init__.py b/GPy/inference/latent_function_inference/__init__.py index c89f771b..dd5d93e0 100644 --- a/GPy/inference/latent_function_inference/__init__.py +++ b/GPy/inference/latent_function_inference/__init__.py @@ -26,5 +26,6 @@ etc. from exact_gaussian_inference import ExactGaussianInference from laplace import LaplaceInference expectation_propagation = 'foo' # TODO +from varDTC import VarDTC from dtc import DTC from fitc import FITC diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index 4b38bd0f..4a8ebd4d 100644 --- a/GPy/kern/parts/rbf.py +++ b/GPy/kern/parts/rbf.py @@ -53,7 +53,7 @@ class RBF(Kernpart): self.variance = Param('variance', variance, Logexp()) - self.lengthscale = Param('lengthscale', lengthscale) + self.lengthscale = Param('lengthscale', lengthscale, Logexp()) self.lengthscale.add_observer(self, self.update_lengthscale) self.update_lengthscale(self.lengthscale) diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py index 1a596eff..b82750ac 100644 --- a/GPy/likelihoods/gaussian.py +++ b/GPy/likelihoods/gaussian.py @@ -88,12 +88,9 @@ class Gaussian(Likelihood): return mu, var, low, up def predictive_mean(self, mu, sigma): - #new_sigma2 = self.predictive_variance(mu, sigma) - #return new_sigma2*(mu/sigma**2 + self.gp_link.transf(mu)/self.variance) return mu def predictive_variance(self, mu, sigma, predictive_mean=None): - #what on earth was the sum of precisions doing here? JH return self.variance + sigma**2 def pdf_link(self, link_f, y, extra_data=None): @@ -305,3 +302,11 @@ class Gaussian(Likelihood): gp = gp.flatten() Ysim = np.array([np.random.normal(self.gp_link.transf(gpj), scale=np.sqrt(self.variance), size=1) for gpj in gp]) return Ysim.reshape(orig_shape) + + def log_predictive_density(self, y_test, mu_star, var_star): + """ + assumes independence + """ + v = var_star + self.variance + return -0.5*np.log(2*np.pi) -0.5*np.log(v) - 0.5*np.square(y_test - mu_star)/v + diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py index 8740a1f5..61defb7d 100644 --- a/GPy/models/sparse_gp_regression.py +++ b/GPy/models/sparse_gp_regression.py @@ -6,6 +6,7 @@ import numpy as np from ..core import SparseGP from .. import likelihoods from .. import kern +from ..inference.latent_function_inference import VarDTC class SparseGPRegression(SparseGP): """ @@ -43,8 +44,7 @@ class SparseGPRegression(SparseGP): likelihood = likelihoods.Gaussian() - SparseGP.__init__(self, X, Y, Z, kernel, likelihood, X_variance=X_variance) - #self.ensure_default_constraints() + SparseGP.__init__(self, X, Y, Z, kernel, likelihood, X_variance=X_variance, inference_method=VarDTC()) def _getstate(self): return SparseGP._getstate(self) From 83d338964419e8c178ebf5feaabb1f2e2f78bdba Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 13 Feb 2014 11:34:23 +0000 Subject: [PATCH 17/40] an ugly hack to work around the 'stickiness' of ObservableArray. TODO: remove this hack --- GPy/inference/latent_function_inference/dtc.py | 4 ++++ GPy/inference/latent_function_inference/fitc.py | 4 ++++ 2 files changed, 8 insertions(+) diff --git a/GPy/inference/latent_function_inference/dtc.py b/GPy/inference/latent_function_inference/dtc.py index 85792d66..bcd0aab8 100644 --- a/GPy/inference/latent_function_inference/dtc.py +++ b/GPy/inference/latent_function_inference/dtc.py @@ -22,6 +22,10 @@ class DTC(object): def inference(self, kern, X, X_variance, Z, likelihood, Y): assert X_variance is None, "cannot use X_variance with DTC. Try varDTC." + #TODO: MAX! fix this! + from ...util.misc import param_to_array + Y = param_to_array(Y) + num_inducing, _ = Z.shape num_data, output_dim = Y.shape diff --git a/GPy/inference/latent_function_inference/fitc.py b/GPy/inference/latent_function_inference/fitc.py index 476715e7..3ad51155 100644 --- a/GPy/inference/latent_function_inference/fitc.py +++ b/GPy/inference/latent_function_inference/fitc.py @@ -19,6 +19,10 @@ class FITC(object): def inference(self, kern, X, X_variance, Z, likelihood, Y): assert X_variance is None, "cannot use X_variance with FITC. Try varDTC." + + #TODO: MAX! fix this! + from ...util.misc import param_to_array + Y = param_to_array(Y) num_inducing, _ = Z.shape num_data, output_dim = Y.shape From 80a734e1534976e6ef2b671101bc387f3156a176 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Thu, 13 Feb 2014 13:53:49 +0000 Subject: [PATCH 18/40] variouschanges --- GPy/core/parameterization/transformations.py | 3 ++- .../latent_function_inference/dtc.py | 5 ++--- .../latent_function_inference/posterior.py | 20 ++++++++++++++----- .../latent_function_inference/varDTC.py | 2 +- 4 files changed, 20 insertions(+), 10 deletions(-) diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index e17da404..06d3643a 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -33,7 +33,8 @@ class Transformation(object): class Logexp(Transformation): domain = _POSITIVE def f(self, x): - return np.where(x>_lim_val, x, np.log(1. + np.exp(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, f, np.log(np.exp(f) - 1.)) def gradfactor(self, f): diff --git a/GPy/inference/latent_function_inference/dtc.py b/GPy/inference/latent_function_inference/dtc.py index bcd0aab8..dbbff6d0 100644 --- a/GPy/inference/latent_function_inference/dtc.py +++ b/GPy/inference/latent_function_inference/dtc.py @@ -44,9 +44,8 @@ class DTC(object): Kmmi, L, Li, _ = pdinv(Kmm) # Compute A - #LiUT, _ = dtrtrs(L, U.T*np.sqrt(beta), lower=1) - LiUT = np.dot(Li, U.T)*np.sqrt(beta) - A = tdot(LiUT) + np.eye(num_inducing) + LiUTbeta = np.dot(Li, U.T)*np.sqrt(beta) + A = tdot(LiUTbeta) + np.eye(num_inducing) # factor A LA = jitchol(A) diff --git a/GPy/inference/latent_function_inference/posterior.py b/GPy/inference/latent_function_inference/posterior.py index c3aa9b36..09ac96e8 100644 --- a/GPy/inference/latent_function_inference/posterior.py +++ b/GPy/inference/latent_function_inference/posterior.py @@ -2,7 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np -from ...util.linalg import pdinv, dpotrs, tdot, dtrtrs, dpotri, symmetrify +from ...util.linalg import pdinv, dpotrs, tdot, dtrtrs, dpotri, symmetrify, jitchol, dtrtri class Posterior(object): """ @@ -80,8 +80,8 @@ class Posterior(object): @property def covariance(self): if self._covariance is None: - LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1) - self._covariance = self._K - tdot(LiK.T) + #LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1) + self._covariance = self._K - self._K.dot(self.woodbury_inv).dot(self._K) return self._covariance @property @@ -93,20 +93,30 @@ class Posterior(object): @property def woodbury_chol(self): if self._woodbury_chol is None: - #try computing woodbury chol from cov + #compute woodbury chol from if self._woodbury_inv is not None: _, _, self._woodbury_chol, _ = pdinv(self._woodbury_inv) + #Li = jitchol(self._woodbury_inv) + #self._woodbury_chol, _ = dtrtri(Li) + #W, _, _, _, = pdinv(self._woodbury_inv) + #symmetrify(W) + #self._woodbury_chol = jitchol(W) + #try computing woodbury chol from cov elif self._covariance is not None: + raise NotImplementedError, "TODO: check code here" B = self._K - self._covariance tmp, _ = dpotrs(self.K_chol, B) self._woodbury_inv, _ = dpotrs(self.K_chol, tmp.T) _, _, self._woodbury_chol, _ = pdinv(self._woodbury_inv) + else: + raise ValueError, "insufficient information to compute posterior" return self._woodbury_chol @property def woodbury_inv(self): if self._woodbury_inv is None: - self._woodbury_inv, _ = dpotri(self.woodbury_chol) + self._woodbury_inv, _ = dpotri(self.woodbury_chol, lower=0) + #self._woodbury_inv, _ = dpotrs(self.woodbury_chol, np.eye(self.woodbury_chol.shape[0]), lower=1) symmetrify(self._woodbury_inv) return self._woodbury_inv diff --git a/GPy/inference/latent_function_inference/varDTC.py b/GPy/inference/latent_function_inference/varDTC.py index 08329b5a..d7f770c8 100644 --- a/GPy/inference/latent_function_inference/varDTC.py +++ b/GPy/inference/latent_function_inference/varDTC.py @@ -34,7 +34,7 @@ class VarDTC(object): Note that L may have fewer columns than Y. """ N, D = Y.shape - if (N>D): + if (N>=D): return param_to_array(Y) else: return jitchol(tdot(Y)) From f71af385056aa1948f62a04a32c7791d6c20a1ab Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 13 Feb 2014 15:06:05 +0000 Subject: [PATCH 19/40] fixing now works, removing parameters needs fixing --- GPy/__init__.py | 1 + GPy/core/parameterization/index_operations.py | 15 ++- GPy/core/parameterization/parameter_core.py | 91 ++++++++----- GPy/core/parameterization/parameterized.py | 127 +++--------------- GPy/core/parameterization/transformations.py | 9 ++ GPy/kern/parts/rbf.py | 2 +- GPy/plotting/matplot_dep/visualize.py | 2 +- 7 files changed, 100 insertions(+), 147 deletions(-) diff --git a/GPy/__init__.py b/GPy/__init__.py index 41369a50..94128dfa 100644 --- a/GPy/__init__.py +++ b/GPy/__init__.py @@ -4,6 +4,7 @@ import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) import core +from core.parameterization import transformations import models import mappings import inference diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index d2549dad..091b6372 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -57,9 +57,11 @@ class ParameterIndexOperations(object): You can give an offset to set an offset for the given indices in the index array, for multi-param handling. ''' - def __init__(self): + def __init__(self, constraints=None): self._properties = IntArrayDict() - #self._reverse = collections.defaultdict(list) + if constraints is not None: + for t, i in constraints.iteritems(): + self.add(t, i) def __getstate__(self): return self._properties#, self._reverse @@ -191,7 +193,7 @@ class ParameterIndexOperationsView(object): def indices(self): - [ind for ind in self.iterindices()] + return [ind for ind in self.iterindices()] def properties_for(self, index): @@ -206,6 +208,8 @@ class ParameterIndexOperationsView(object): removed = self._param_index_ops.remove(prop, indices+self._offset) if removed.size > 0: return removed - self._size + 1 + if self[prop].size == 0: + del self[prop] return removed @@ -222,6 +226,9 @@ class ParameterIndexOperationsView(object): def update(self, parameter_index_view): for i, v in parameter_index_view.iteritems(): self.add(i, v) - + + + def copy(self): + return ParameterIndexOperations(dict(self.iteritems())) pass diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 6afea470..cfee60bd 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -1,7 +1,7 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -from transformations import Transformation, Logexp, NegativeLogexp, Logistic +from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED __updated__ = '2013-12-16' @@ -10,11 +10,6 @@ def adjust_name_for_printing(name): return name.replace(" ", "_").replace(".", "_").replace("-","").replace("+","").replace("!","").replace("*","").replace("/","") return '' -#=============================================================================== -# Printing: -__fixed__ = "fixed" -#=============================================================================== - class Observable(object): _observers_ = {} def add_observer(self, observer, callble): @@ -119,6 +114,14 @@ class Indexable(object): def _offset_for(self, param): raise NotImplementedError, "shouldnt happen, offset required from non parameterization object?" + + def _raveled_index_for(self, param): + """ + get the raveled index for a param + that is an int array, containing the indexes for the flattened + param inside this parameterized logic. + """ + raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization object?" class Constrainable(Nameable, Indexable, Parameterizable): def __init__(self, name, default_constraint=None): @@ -126,6 +129,9 @@ class Constrainable(Nameable, Indexable, Parameterizable): self._default_constraint_ = default_constraint from index_operations import ParameterIndexOperations self.constraints = ParameterIndexOperations() + if self._default_constraint_ is not None: + self.constrain(self._default_constraint_) + #=========================================================================== # Fixing Parameters: #=========================================================================== @@ -138,17 +144,40 @@ class Constrainable(Nameable, Indexable, Parameterizable): if value is not None: self[:] = value self.constrain(__fixed__, warning=warning) - self._highest_parent_._set_fixed(self._raveled_index()) + rav_i = self._highest_parent_._raveled_index_for(self) + self._highest_parent_._set_fixed(rav_i) fix = constrain_fixed + def unconstrain_fixed(self): """ This parameter will no longer be fixed. """ unconstrained = self.unconstrain(__fixed__) - import ipdb;ipdb.set_trace() - self._highest_parent_._set_unfixed(unconstrained) - + self._highest_parent_._set_unfixed(unconstrained) unfix = unconstrain_fixed + + def _set_fixed(self, index): + import numpy as np + if not self._has_fixes(): self._fixes_ = np.ones(self.size, dtype=bool) + self._fixes_[index] = FIXED + if np.all(self._fixes_): self._fixes_ = None # ==UNFIXED + + def _set_unfixed(self, index): + import numpy as np + if not self._has_fixes(): self._fixes_ = np.ones(self.size, dtype=bool) + #rav_i = self._raveled_index_for(param)[index] + self._fixes_[index] = UNFIXED + if np.all(self._fixes_): self._fixes_ = None # ==UNFIXED + + def _connect_fixes(self): + import numpy as np + fixed_indices = self.constraints[__fixed__] + if fixed_indices.size > 0: + self._fixes_ = np.ones(self.size, dtype=bool) * UNFIXED + self._fixes_[fixed_indices] = FIXED + else: + self._fixes_ = None + #=========================================================================== # Constrain operations -> done #=========================================================================== @@ -174,18 +203,29 @@ class Constrainable(Nameable, Indexable, Parameterizable): self._set_params(transform.initialize(self._get_params()), update=False) reconstrained = self.unconstrain() self.constraints.add(transform, self._raveled_index()) - if reconstrained.size > 0: + if warning and reconstrained.size > 0: print "WARNING: reconstraining parameters {}".format(self.parameter_names) if update: self._highest_parent_.parameters_changed() - # if self.has_parent(): - # self._highest_parent_._add_constrain(self, transform, warning) - # else: - # for p in self._parameters_: - # self._add_constrain(p, transform, warning) - # if update: - # self.parameters_changed() + def unconstrain(self, *transforms): + """ + :param transforms: The transformations to unconstrain from. + + remove all :py:class:`GPy.core.transformations.Transformation` + transformats of this parameter object. + """ + if len(transforms) == 0: + transforms = self.constraints.properties() + import numpy as np + removed = np.empty((0,),dtype=int) + for t in transforms: + unconstrained = self.constraints.remove(t, self._raveled_index()) + removed = np.union1d(removed, unconstrained) + if t is __fixed__: + self._highest_parent_._set_unfixed(unconstrained) + return removed + def constrain_positive(self, warning=True, update=True): """ :param warning: print a warning if re-constraining parameters. @@ -211,21 +251,6 @@ class Constrainable(Nameable, Indexable, Parameterizable): """ self.constrain(Logistic(lower, upper), warning=warning, update=update) - def unconstrain(self, *transforms): - """ - :param transforms: The transformations to unconstrain from. - - remove all :py:class:`GPy.core.transformations.Transformation` - transformats of this parameter object. - """ - if len(transforms) == 0: - transforms = self.constraints.properties() - import numpy as np - removed = np.empty((0,),dtype=int) - for t in transforms: - removed = np.union1d(removed, self.constraints.remove(t, self._raveled_index())) - return removed - def unconstrain_positive(self): """ Remove positive constraint of this parameter. diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index a30a1e3b..80bf2959 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -8,15 +8,10 @@ import cPickle import itertools from re import compile, _pattern_type from param import ParamConcatenation, Param -from parameter_core import Constrainable, Pickleable, Observable, adjust_name_for_printing, Gradcheckable, __fixed__ +from parameter_core import Constrainable, Pickleable, Observable, adjust_name_for_printing, Gradcheckable +from transformations import __fixed__, FIXED, UNFIXED from array_core import ParamList -#=============================================================================== -# constants -FIXED = False -UNFIXED = True -#=============================================================================== - class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): """ Parameterized class @@ -71,37 +66,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_ = set() del self._in_init_ - #=========================================================================== - # Parameter connection for model creation: - #=========================================================================== -# def set_as_parameter(self, name, array, gradient, index=None, gradient_parent=None): -# """ -# :param name: name of the param (in print and plots), can be callable without parameters -# :type name: str, callable -# :param array: array which the param consists of -# :type array: array-like -# :param gradient: gradient method of the param -# :type gradient: callable -# :param index: (optional) index of the param when printing -# -# (:param gradient_parent: connect these parameters to this class, but tell -# updates to highest_parent, this is needed when parameterized classes -# contain parameterized classes, but want to access the parameters -# of their children) -# -# -# Set array (e.g. self.X) as param with name and gradient. -# I.e: self.set_as_parameter('curvature', self.lengthscale, self.dK_dlengthscale) -# -# Note: the order in which parameters are added can be adjusted by -# giving an index, of where to put this param in printing -# """ -# if index is None: -# self._parameters_.append(Param(name, array, gradient)) -# else: -# self._parameters_.insert(index, Param(name, array, gradient)) -# self._connect_parameters(gradient_parent=gradient_parent) - def _has_fixes(self): return hasattr(self, "_fixes_") and self._fixes_ is not None @@ -118,53 +82,22 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): # if param.has_parent(): # raise AttributeError, "parameter {} already in another model, create new object (or copy) for adding".format(param._short()) if param in self._parameters_ and index is not None: - # make sure fixes and constraints are indexed right - if self._has_fixes(): - param_slice = slice(self._offset_for(param), self._offset_for(param) + param.size) - dest_index = sum((p.size for p in self._parameters_[:index])) - dest_slice = slice(dest_index, dest_index + param.size) - fixes_param = self._fixes_[param_slice].copy() - self._fixes_[param_slice] = self._fixes_[dest_slice] - self._fixes_[dest_slice] = fixes_param - - del self._parameters_[param._parent_index_] - self._parameters_.insert(index, param) + self.remove_parameter(param) + self.add_parameter(param, index) elif param not in self._parameters_: # make sure the size is set - if not hasattr(self, 'size'): - self.size = sum(p.size for p in self._parameters_) if index is None: self._parameters_.append(param) - - # make sure fixes and constraints are indexed right - if param._has_fixes(): fixes_param = param._fixes_.copy() - else: fixes_param = numpy.ones(param.size, dtype=bool) - if self._has_fixes(): self._fixes_ = np.r_[self._fixes_, fixes_param] - elif param._has_fixes(): self._fixes_ = np.r_[np.ones(self.size, dtype=bool), fixes_param] - else: start = sum(p.size for p in self._parameters_[:index]) self.constraints.shift(start, param.size) self._parameters_.insert(index, param) - - # make sure fixes and constraints are indexed right - if param._has_fixes(): fixes_param = param._fixes_.copy() - else: fixes_param = numpy.ones(param.size, dtype=bool) - ins = sum((p.size for p in self._parameters_[:index])) - if self._has_fixes(): self._fixes_ = np.r_[self._fixes_[:ins], fixes_param, self._fixes[ins:]] - elif not np.all(fixes_param): - self._fixes_ = np.ones(self.size + param.size, dtype=bool) - self._fixes_[ins:ins + param.size] = fixes_param self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" self._connect_parameters() for p in self._parameters_: p._parent_changed(self) - if param._default_constraint_ is not None: - param.constrain(param._default_constraint_, False) - if self._has_fixes() and np.all(self._fixes_): # ==UNFIXED - self._fixes_ = None def add_parameters(self, *parameters): """ @@ -173,18 +106,20 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): """ [self.add_parameter(p) for p in parameters] - def remove_parameter(self, *names_params_indices): + def remove_parameter(self, param): """ - :param names_params_indices: mix of parameter_names, param objects, or indices - to remove from being a param of this parameterized object. - - note: if it is a string object it will not (!) be regexp-matched - automatically. + :param param: param object to remove from being a parameter of this parameterized object. """ - self._parameters_ = ParamList([p for p in self._parameters_ - if not (p._parent_index_ in names_params_indices - or p.name in names_params_indices - or p in names_params_indices)]) + if not param in self._parameters_: + raise RuntimeError, "Parameter {} does not belong to this object, remove parameters directly from their respective parents".format(param._short()) + del self._parameters_[param._parent_index_] + self.size -= param.size + constr = param.constraints.copy() + param.constraints.clear() + param.constraints = constr + param._direct_parent_ = None + param._parent_index_ = None + param._connect_fixes() self._connect_parameters() def _connect_parameters(self): @@ -214,6 +149,8 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): elif not (pname in not_unique): self.__dict__[pname] = p self._added_names_.add(pname) + self._connect_fixes() + #=========================================================================== # Pickling operations #=========================================================================== @@ -337,7 +274,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_.add(pname) self.__dict__[pname] = param #=========================================================================== - # Index Handling + # Indexable Handling #=========================================================================== def _backtranslate_index(self, param, ind): # translate an index in parameterized indexing into the index of param @@ -373,36 +310,10 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): #=========================================================================== # Fixing parameters: #=========================================================================== - def _set_fixed(self, param_or_index): - if not self._has_fixes(): self._fixes_ = numpy.ones(self.size, dtype=bool) - try: - param_or_index = self._raveled_index_for(param_or_index) - except AttributeError: - pass - self._fixes_[param_or_index] = FIXED - if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED - def _set_unfixed(self, param_or_index): - if not self._has_fixes(): self._fixes_ = numpy.ones(self.size, dtype=bool) - try: - param_or_index = self._raveled_index_for(param_or_index) - except AttributeError: - pass - self._fixes_[param_or_index] = UNFIXED - for constr, ind in self.constraints.iteritems(): - if constr is __fixed__: - self._fixes_[ind] = FIXED - if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED 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)] - # def _fix(self, param, warning=True): - # f = self._add_constrain(param, __fixed__, warning) - # self._set_fixed(f) - # def _unfix(self, param): - # if self._has_fixes(): - # f = self._remove_constrain(param, __fixed__) - # self._set_unfixed(f) #=========================================================================== # Convenience for fixed, tied checking of param: #=========================================================================== diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index e17da404..9582ad4a 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -6,8 +6,17 @@ import numpy as np from domains import _POSITIVE,_NEGATIVE, _BOUNDED import sys import weakref + _lim_val = -np.log(sys.float_info.epsilon) +#=============================================================================== +# Fixing constants +__fixed__ = "fixed" +FIXED = False +UNFIXED = True +#=============================================================================== + + class Transformation(object): domain = None _instance = None diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index 4b38bd0f..4a8ebd4d 100644 --- a/GPy/kern/parts/rbf.py +++ b/GPy/kern/parts/rbf.py @@ -53,7 +53,7 @@ class RBF(Kernpart): self.variance = Param('variance', variance, Logexp()) - self.lengthscale = Param('lengthscale', lengthscale) + self.lengthscale = Param('lengthscale', lengthscale, Logexp()) self.lengthscale.add_observer(self, self.update_lengthscale) self.update_lengthscale(self.lengthscale) diff --git a/GPy/plotting/matplot_dep/visualize.py b/GPy/plotting/matplot_dep/visualize.py index 99e8a0da..fb085f39 100644 --- a/GPy/plotting/matplot_dep/visualize.py +++ b/GPy/plotting/matplot_dep/visualize.py @@ -4,7 +4,6 @@ import GPy import numpy as np import matplotlib as mpl import time -import Image try: import visual visual_available = True @@ -324,6 +323,7 @@ class image_show(matplotlib_show): else: self.vals = 255*(self.vals - self.vals.min())/(self.vals.max() - self.vals.min()) if not self.palette == []: # applying using an image palette (e.g. if the image has been quantized) + from PIL import Image self.vals = Image.fromarray(self.vals.astype('uint8')) self.vals.putpalette(self.palette) # palette is a list, must be loaded before calling this function From 13212abd0b0b8fb40aebd367c16ccbc334aabeaa Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 13 Feb 2014 16:44:45 +0000 Subject: [PATCH 20/40] parameter adding and removing now fully functional according to tests, including fixes --- GPy/core/parameterization/index_operations.py | 13 ++- GPy/core/parameterization/parameter_core.py | 11 ++- GPy/core/parameterization/parameterized.py | 14 ++- GPy/testing/parameterized_tests.py | 93 +++++++++++++++++++ 4 files changed, 123 insertions(+), 8 deletions(-) create mode 100644 GPy/testing/parameterized_tests.py diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index 091b6372..b816e05f 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -57,6 +57,7 @@ class ParameterIndexOperations(object): You can give an offset to set an offset for the given indices in the index array, for multi-param handling. ''' + _offset = 0 def __init__(self, constraints=None): self._properties = IntArrayDict() if constraints is not None: @@ -120,6 +121,14 @@ class ParameterIndexOperations(object): return removed.astype(int) return numpy.array([]).astype(int) + def update(self, parameter_index_view, offset=0): + for i, v in parameter_index_view.iteritems(): + self.add(i, v+offset) + + + def copy(self): + return ParameterIndexOperations(dict(self.iteritems())) + def __getitem__(self, prop): return self._properties[prop] @@ -223,9 +232,9 @@ class ParameterIndexOperationsView(object): import pprint return pprint.pformat(dict(self.iteritems())) - def update(self, parameter_index_view): + def update(self, parameter_index_view, offset=0): for i, v in parameter_index_view.iteritems(): - self.add(i, v) + self.add(i, v+offset) def copy(self): diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index cfee60bd..65504652 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -72,6 +72,13 @@ class Parentable(object): def has_parent(self): return self._direct_parent_ is not None + def _notify_parent_change(self): + for p in self._parameters_: + p._parent_changed(self) + + def _parent_changed(self): + raise NotImplementedError, "shouldnt happen, Parentable objects need to be able to change their parent" + @property def _highest_parent_(self): if self._direct_parent_ is None: @@ -182,11 +189,9 @@ class Constrainable(Nameable, Indexable, Parameterizable): # Constrain operations -> done #=========================================================================== def _parent_changed(self, parent): - c = self.constraints from index_operations import ParameterIndexOperationsView self.constraints = ParameterIndexOperationsView(parent.constraints, parent._offset_for(self), self.size) - self.constraints.update(c) - del c + self._fixes_ = None for p in self._parameters_: p._parent_changed(parent) diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 80bf2959..a976eb93 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -87,17 +87,20 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): elif param not in self._parameters_: # make sure the size is set if index is None: + self.constraints.update(param.constraints, self.size) self._parameters_.append(param) else: start = sum(p.size for p in self._parameters_[:index]) self.constraints.shift(start, param.size) + self.constraints.update(param.constraints, start) self._parameters_.insert(index, param) self.size += param.size else: raise RuntimeError, """Parameter exists already added and no copy made""" self._connect_parameters() - for p in self._parameters_: - p._parent_changed(self) + self._notify_parent_change() + self._connect_fixes() + def add_parameters(self, *parameters): """ @@ -120,7 +123,13 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): param._direct_parent_ = None param._parent_index_ = None param._connect_fixes() + param._notify_parent_change() + pname = adjust_name_for_printing(param.name) + if pname in self._added_names_: + del self.__dict__[pname] self._connect_parameters() + #self._notify_parent_change() + self._connect_fixes() def _connect_parameters(self): # connect parameterlist to this parameterized object @@ -149,7 +158,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): elif not (pname in not_unique): self.__dict__[pname] = p self._added_names_.add(pname) - self._connect_fixes() #=========================================================================== # Pickling operations diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py new file mode 100644 index 00000000..ff57606a --- /dev/null +++ b/GPy/testing/parameterized_tests.py @@ -0,0 +1,93 @@ +''' +Created on Feb 13, 2014 + +@author: maxzwiessele +''' +import unittest +import GPy +import numpy as np + +class Test(unittest.TestCase): + + def setUp(self): + self.rbf = GPy.kern.rbf(1) + self.white = GPy.kern.white(1) + from GPy.core.parameterization import Param + from GPy.core.parameterization.transformations import Logistic + self.param = Param('param', np.random.rand(25,2), Logistic(0, 1)) + + self.test1 = GPy.core.Parameterized("test model") + self.test1.add_parameter(self.white) + self.test1.add_parameter(self.rbf, 0) + self.test1.add_parameter(self.param) + + def test_add_parameter(self): + self.assertEquals(self.rbf._parent_index_, 0) + self.assertEquals(self.white._parent_index_, 1) + pass + + def test_fixes(self): + self.white.fix(warning=False) + self.test1.remove_parameter(self.test1.param) + self.assertTrue(self.test1._has_fixes()) + + from GPy.core.parameterization.transformations import FIXED, UNFIXED + self.assertListEqual(self.test1._fixes_.tolist(),[UNFIXED,UNFIXED,FIXED]) + + self.test1.add_parameter(self.white, 0) + self.assertListEqual(self.test1._fixes_.tolist(),[FIXED,UNFIXED,UNFIXED]) + + + def test_remove_parameter(self): + from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__ + self.white.fix() + self.test1.remove_parameter(self.white) + self.assertIs(self.test1._fixes_,None) + + self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) + self.assertIs(self.white.constraints,self.white.white.constraints._param_index_ops) + self.assertEquals(self.white.white.constraints._offset, 0) + self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) + + self.test1.add_parameter(self.white, 0) + self.assertIs(self.test1.constraints, self.white.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) + self.assertListEqual(self.test1.constraints[__fixed__].tolist(), [0]) + self.assertIs(self.white._fixes_,None) + self.assertListEqual(self.test1._fixes_.tolist(),[FIXED] + [UNFIXED] * 52) + self.test1.remove_parameter(self.white) + self.assertIs(self.test1._fixes_,None) + self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) + self.assertIs(self.white.constraints,self.white.white.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) + + def test_add_parameter_already_in_hirarchy(self): + self.test1.add_parameter(self.white._parameters_[0]) + + def test_default_constraints(self): + self.assertIs(self.rbf.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops) + self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) + self.assertListEqual(self.rbf.constraints.indices()[0].tolist(), range(2)) + from GPy.core.parameterization.transformations import Logexp + kern = self.rbf+self.white + self.assertListEqual(kern.constraints[Logexp()].tolist(), range(3)) + + def test_constraints(self): + self.rbf.constrain(GPy.transformations.Square(), False) + self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), range(2)) + self.assertListEqual(self.test1.constraints[GPy.transformations.Logexp()].tolist(), [2]) + + self.test1.remove_parameter(self.rbf) + self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), []) + + def test_constraints_views(self): + self.assertEqual(self.white.constraints._offset, 2) + self.assertEqual(self.rbf.constraints._offset, 0) + self.assertEqual(self.param.constraints._offset, 3) + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'Test.test_add_parameter'] + unittest.main() \ No newline at end of file From c84372127f8165a3b8c28682d6ec060864d42055 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 13 Feb 2014 16:49:17 +0000 Subject: [PATCH 21/40] checkgrad was changing parameters --- GPy/core/model.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/core/model.py b/GPy/core/model.py index f2632c93..2e2e63ed 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -470,6 +470,7 @@ class Model(Parameterized): ng = '%.6f' % float(numerical_gradient) grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4]) print grad_string + self._set_params_transformed(x) return ret def input_sensitivity(self): From c78ddde6debb0c6cb2b591aa46313596423fc058 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 13 Feb 2014 20:24:45 +0000 Subject: [PATCH 22/40] gradcheck now fully functional --- GPy/core/model.py | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 2e2e63ed..1252198a 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -410,7 +410,13 @@ class Model(Parameterized): f1 = self.objective_function(x + dx) f2 = self.objective_function(x - dx) gradient = self.objective_function_gradients(x) - + + param_list = self._raveled_index_for(target_param) + if self._has_fixes(): + param_list = np.intersect1d(np.r_[:self.size][self._fixes_], param_list, True) + + + numerical_gradient = (f1 - f2) / (2 * dx) global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) @@ -448,14 +454,14 @@ class Model(Parameterized): ret = True for i, ind in enumerate(param_list): xx = x.copy() - xx[i] += step + xx[ind] += step f1 = self.objective_function(xx) - xx[i] -= 2.*step + xx[ind] -= 2.*step f2 = self.objective_function(xx) - + print ind numerical_gradient = (f1 - f2) / (2 * step) - ratio = (f1 - f2) / (2 * step * gradient[i]) - difference = np.abs((f1 - f2) / 2 / step - gradient[i]) + ratio = (f1 - f2) / (2 * step * gradient[ind]) + difference = np.abs((f1 - f2) / 2 / step - gradient[ind]) if (np.abs(1. - ratio) < tolerance) or np.abs(difference) < tolerance: formatted_name = "\033[92m {0} \033[0m".format(names[ind]) @@ -466,7 +472,7 @@ class Model(Parameterized): r = '%.6f' % float(ratio) d = '%.6f' % float(difference) - g = '%.6f' % gradient[i] + g = '%.6f' % gradient[ind] ng = '%.6f' % float(numerical_gradient) grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4]) print grad_string From 9af4c34f90ebdefeb402603496f34e8584e7f318 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 13 Feb 2014 21:14:08 +0000 Subject: [PATCH 23/40] gradcheck fixes are not easy --- GPy/core/model.py | 53 +++++++++++++-------- GPy/core/parameterization/param.py | 2 +- GPy/core/parameterization/parameter_core.py | 8 ++-- GPy/core/parameterization/parameterized.py | 8 +--- 4 files changed, 42 insertions(+), 29 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 1252198a..e7c6c7ec 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -411,12 +411,21 @@ class Model(Parameterized): f2 = self.objective_function(x - dx) gradient = self.objective_function_gradients(x) - param_list = self._raveled_index_for(target_param) - if self._has_fixes(): - param_list = np.intersect1d(np.r_[:self.size][self._fixes_], param_list, True) - - + if target_param is None: + transformed_index = range(len(x)) + else: + transformed_index = self._raveled_index_for(target_param) + if self._has_fixes(): + indices = np.r_[:self.size] + which = (transformed_index[:,None]==indices[self._fixes_]).nonzero()[0] + indices -= (~self._fixes_).cumsum() + transformed_index = indices[which] + + if transformed_index.size == 0: + print "No free parameters to check" + return + 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))) @@ -439,40 +448,46 @@ class Model(Parameterized): separator = '-' * len(header_string[0]) print '\n'.join([header_string[0], separator]) if target_param is None: - param_list = range(len(x)) + param_index = range(len(x)) + transformed_index = param_index else: - param_list = self._raveled_index_for(target_param) + param_index = self._raveled_index_for(target_param) if self._has_fixes(): - param_list = np.intersect1d(np.r_[:self.size][self._fixes_], param_list, True) - - if param_list.size == 0: + indices = np.r_[:self.size] + which = (param_index[:,None]==indices[self._fixes_][None,:]).nonzero() + transformed_index = (indices-(~self._fixes_).cumsum())[which[1]] + param_index = indices[which[0]] + print param_index, transformed_index + else: + transformed_index = param_index + + if param_index.size == 0: print "No free parameters to check" return gradient = self.objective_function_gradients(x) np.where(gradient == 0, 1e-312, gradient) ret = True - for i, ind in enumerate(param_list): + for nind, xind in itertools.izip(param_index, transformed_index): xx = x.copy() - xx[ind] += step + xx[xind] += step f1 = self.objective_function(xx) - xx[ind] -= 2.*step + xx[xind] -= 2.*step f2 = self.objective_function(xx) - print ind numerical_gradient = (f1 - f2) / (2 * step) - ratio = (f1 - f2) / (2 * step * gradient[ind]) - difference = np.abs((f1 - f2) / 2 / step - gradient[ind]) + ratio = (f1 - f2) / (2 * step * gradient[xind]) + difference = np.abs((f1 - f2) / 2 / step - gradient[xind]) if (np.abs(1. - ratio) < tolerance) or np.abs(difference) < tolerance: - formatted_name = "\033[92m {0} \033[0m".format(names[ind]) + formatted_name = "\033[92m {0} \033[0m".format(names[nind]) ret &= True else: - formatted_name = "\033[91m {0} \033[0m".format(names[ind]) + formatted_name = "\033[91m {0} \033[0m".format(names[nind]) ret &= False r = '%.6f' % float(ratio) d = '%.6f' % float(difference) - g = '%.6f' % gradient[ind] + g = '%.6f' % gradient[xind] ng = '%.6f' % float(numerical_gradient) grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4]) print grad_string diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 56b8aab6..940cd1a0 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -348,7 +348,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri def _description_str(self): if self.size <= 1: return ["%f" % self] else: return [str(self.shape)] - def _parameter_names(self, add_name): + def parameter_names(self, add_name=False): return [self.name] @property def flattened_parameters(self): diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 65504652..7ad228f7 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -25,8 +25,10 @@ class Parameterizable(object): from GPy.core.parameterization.array_core import ParamList _parameters_ = ParamList() - def parameter_names(self): - return [p.name for p in self._parameters_] + def parameter_names(self, add_name=False): + if add_name: + return [adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] + return [xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] def parameters_changed(self): """ @@ -209,7 +211,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): reconstrained = self.unconstrain() self.constraints.add(transform, self._raveled_index()) if warning and reconstrained.size > 0: - print "WARNING: reconstraining parameters {}".format(self.parameter_names) + print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) if update: self._highest_parent_.parameters_changed() diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index a976eb93..593c37e9 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -434,11 +434,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self._direct_parent_.hirarchy_name() + adjust_name_for_printing(self.name) else: return adjust_name_for_printing(self.name) - def _parameter_names(self, add_name=False): - if add_name: - return [adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x._parameter_names(add_name=True)] - return [xi for x in self._parameters_ for xi in x._parameter_names(add_name=True)] - parameter_names = property(_parameter_names, doc="Names for all parameters handled by this parameterization object -- will add hirarchy name entries for printing") + #parameter_names = property(parameter_names, doc="Names for all parameters handled by this parameterization object -- will add hirarchy name entries for printing") def _collect_gradient(self, target): [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] @property @@ -468,7 +464,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): name = adjust_name_for_printing(self.name) + "." constrs = self._constraints_str; ts = self._ties_str - desc = self._description_str; names = self.parameter_names + desc = self._description_str; names = self.parameter_names() nl = max([len(str(x)) for x in names + [name]]) sl = max([len(str(x)) for x in desc + ["Value"]]) cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]]) From ffd09c7820f82d6bc8f3b46ab4b59ee34952b7e6 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Thu, 13 Feb 2014 21:59:08 +0000 Subject: [PATCH 24/40] =?UTF-8?q?checkgrad=20=20=20=20(=E2=95=AF=C2=B0?= =?UTF-8?q?=E2=96=A1=C2=B0=EF=BC=89=E2=95=AF=EF=B8=B5=20=E2=94=BB=E2=94=81?= =?UTF-8?q?=E2=94=BB?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- GPy/core/model.py | 4 ++-- GPy/core/parameterization/array_core.py | 2 +- GPy/core/parameterization/param.py | 7 ++++++- GPy/examples/regression.py | 1 - 4 files changed, 9 insertions(+), 5 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index e7c6c7ec..6bc3edb9 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -455,8 +455,8 @@ class Model(Parameterized): if self._has_fixes(): indices = np.r_[:self.size] which = (param_index[:,None]==indices[self._fixes_][None,:]).nonzero() - transformed_index = (indices-(~self._fixes_).cumsum())[which[1]] - param_index = indices[which[0]] + param_index = param_index[which[0]] + transformed_index = (indices-(~self._fixes_).cumsum())[param_index] print param_index, transformed_index else: transformed_index = param_index diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 975b78f4..78b9e8f2 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -28,8 +28,8 @@ class ObservableArray(np.ndarray, Observable): """ __array_priority__ = -1 # Never give back ObservableArray def __new__(cls, input_array): - cls.__name__ = "ObservableArray\n " obj = np.atleast_1d(input_array).view(cls) + cls.__name__ = "ObservableArray\n " obj._observers_ = {} return obj def __array_finalize__(self, obj): diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 940cd1a0..97fcfd32 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -43,6 +43,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri _parameters_ = [] def __new__(cls, name, input_array, default_constraint=None): obj = numpy.atleast_1d(super(Param, cls).__new__(cls, input_array=input_array)) + cls.__name__ = "Param" obj._current_slice_ = (slice(obj.shape[0]),) obj._realshape_ = obj.shape obj._realsize_ = obj.size @@ -57,7 +58,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri def __init__(self, name, input_array, default_constraint=None): super(Param, self).__init__(name=name, default_constraint=default_constraint) - + def __array_finalize__(self, obj): # see InfoArray.__array_finalize__ for comments if obj is None: return @@ -75,6 +76,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri self._original_ = getattr(obj, '_original_', None) self._name = getattr(obj, 'name', None) self.gradient = getattr(obj, 'gradient', None) + self.constraints = getattr(obj, 'constraints', None) def __array_wrap__(self, out_arr, context=None): return out_arr.view(numpy.ndarray) @@ -391,6 +393,9 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri slice_index = self._current_slice_ if isinstance(slice_index, (tuple, list)): clean_curr_slice = [s for s in slice_index if numpy.any(s != Ellipsis)] + for i in range(self._realndim_-len(clean_curr_slice)): + i+=len(clean_curr_slice) + clean_curr_slice += range(self._realshape_[i]) if (all(isinstance(n, (numpy.ndarray, list, tuple)) for n in clean_curr_slice) and len(set(map(len, clean_curr_slice))) <= 1): return numpy.fromiter(itertools.izip(*clean_curr_slice), diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index f8d3d5a9..c073369f 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -473,7 +473,6 @@ def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True): Z = np.random.uniform(-3., 3., (7, 1)) k = GPy.kern.rbf(1) - import ipdb;ipdb.set_trace() # create simple GP Model - no input uncertainty on this one m = GPy.models.SparseGPRegression(X, Y, kernel=GPy.kern.rbf(1), Z=Z) From 8d7273b8ff0bc2f91ebebc7ca886251226d10c3f Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 09:57:18 +0000 Subject: [PATCH 25/40] non verbose checkgrad adjusted to new system --- GPy/core/model.py | 28 +++++++++++++++------------- 1 file changed, 15 insertions(+), 13 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 6bc3edb9..c944f943 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -403,32 +403,34 @@ class Model(Parameterized): x = self._get_params_transformed().copy() if not verbose: - # just check the global ratio - dx = step * np.sign(np.random.uniform(-1, 1, x.size)) - - # evaulate around the point x - f1 = self.objective_function(x + dx) - f2 = self.objective_function(x - dx) - gradient = self.objective_function_gradients(x) - + # make sure only to test the selected parameters if target_param is None: transformed_index = range(len(x)) else: transformed_index = self._raveled_index_for(target_param) if self._has_fixes(): indices = np.r_[:self.size] - which = (transformed_index[:,None]==indices[self._fixes_]).nonzero()[0] - indices -= (~self._fixes_).cumsum() - transformed_index = indices[which] + 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) + gradient = self.objective_function_gradients(x) + + 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) else: # check the gradient of each parameter individually, and do some pretty printing @@ -457,7 +459,7 @@ class Model(Parameterized): which = (param_index[:,None]==indices[self._fixes_][None,:]).nonzero() param_index = param_index[which[0]] transformed_index = (indices-(~self._fixes_).cumsum())[param_index] - print param_index, transformed_index + #print param_index, transformed_index else: transformed_index = param_index From 3375606039782cbea9edabedf6d04d1e857aad5e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Fri, 14 Feb 2014 11:38:34 +0000 Subject: [PATCH 26/40] lots of F-ordering nonsense. Seems to work though --- GPy/core/sparse_gp.py | 2 +- .../latent_function_inference/posterior.py | 2 +- .../latent_function_inference/varDTC.py | 22 ++- GPy/util/linalg.py | 171 ++++++++++-------- 4 files changed, 109 insertions(+), 88 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index bbbea2fe..09e84c12 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -2,7 +2,7 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np -from ..util.linalg import mdot, tdot, symmetrify, backsub_both_sides, chol_inv, dtrtrs, dpotrs, dpotri +from ..util.linalg import mdot, tdot, symmetrify, backsub_both_sides, dtrtrs, dpotrs, dpotri from gp import GP from parameterization.param import Param from ..inference.latent_function_inference import varDTC diff --git a/GPy/inference/latent_function_inference/posterior.py b/GPy/inference/latent_function_inference/posterior.py index 09ac96e8..f28bf9d1 100644 --- a/GPy/inference/latent_function_inference/posterior.py +++ b/GPy/inference/latent_function_inference/posterior.py @@ -115,7 +115,7 @@ class Posterior(object): @property def woodbury_inv(self): if self._woodbury_inv is None: - self._woodbury_inv, _ = dpotri(self.woodbury_chol, lower=0) + self._woodbury_inv, _ = dpotri(self.woodbury_chol, lower=1) #self._woodbury_inv, _ = dpotrs(self.woodbury_chol, np.eye(self.woodbury_chol.shape[0]), lower=1) symmetrify(self._woodbury_inv) return self._woodbury_inv diff --git a/GPy/inference/latent_function_inference/varDTC.py b/GPy/inference/latent_function_inference/varDTC.py index d7f770c8..90ca7a6a 100644 --- a/GPy/inference/latent_function_inference/varDTC.py +++ b/GPy/inference/latent_function_inference/varDTC.py @@ -91,7 +91,7 @@ class VarDTC(object): tmp = psi1 * (np.sqrt(beta.reshape(num_data, 1))) else: tmp = psi1 * (np.sqrt(beta)) - tmp, _ = dtrtrs(Lm, np.asfortranarray(tmp.T), lower=1) + tmp, _ = dtrtrs(Lm, tmp.T, lower=1) A = tdot(tmp) # factor B @@ -100,14 +100,16 @@ class VarDTC(object): LB = jitchol(B) # 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) + #self.YYTfactor = self.get_YYTfactor(Y) + #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) + VVT_factor = beta*param_to_array(Y) trYYT = self.get_trYYT(Y) + psi1Vf = np.dot(psi1.T, VVT_factor) # back substutue C into psi1Vf - tmp, info1 = dtrtrs(Lm, np.asfortranarray(psi1Vf), lower=1, trans=0) - _LBi_Lmi_psi1Vf, _ = dtrtrs(LB, np.asfortranarray(tmp), lower=1, trans=0) + tmp, info1 = dtrtrs(Lm, psi1Vf, lower=1, trans=0) + _LBi_Lmi_psi1Vf, _ = dtrtrs(LB, tmp, lower=1, trans=0) tmp, info2 = dtrtrs(LB, _LBi_Lmi_psi1Vf, lower=1, trans=1) Cpsi1Vf, info3 = dtrtrs(Lm, tmp, lower=1, trans=1) @@ -152,8 +154,8 @@ class VarDTC(object): raise NotImplementedError, "heteroscedatic derivates with uncertain inputs not implemented" else: LBi = chol_inv(LB) - Lmi_psi1, nil = dtrtrs(Lm, np.asfortranarray(psi1.T), lower=1, trans=0) - _LBi_Lmi_psi1, _ = dtrtrs(LB, np.asfortranarray(Lmi_psi1), lower=1, trans=0) + Lmi_psi1, nil = dtrtrs(Lm, psi1.T, lower=1, trans=0) + _LBi_Lmi_psi1, _ = dtrtrs(LB, Lmi_psi1, lower=1, trans=0) partial_for_likelihood = -0.5 * beta + 0.5 * likelihood.V**2 partial_for_likelihood += 0.5 * output_dim * (psi0 - np.sum(Lmi_psi1**2,0))[:,None] * beta**2 @@ -195,12 +197,14 @@ class VarDTC(object): if VVT_factor.shape[1] == Y.shape[1]: woodbury_vector = Cpsi1Vf # == Cpsi1V else: + print 'foobar' psi1V = np.dot(Y.T*beta, psi1).T - tmp, _ = dtrtrs(Lm, np.asfortranarray(psi1V), lower=1, trans=0) + tmp, _ = dtrtrs(Lm, psi1V, lower=1, trans=0) tmp, _ = dpotrs(LB, tmp, lower=1) woodbury_vector, _ = dtrtrs(Lm, tmp, lower=1, trans=1) - Bi, _ = dpotri(LB, lower=0) + Bi, _ = dpotri(LB, lower=1) symmetrify(Bi) + Bi = dpotri(LB, lower=1)[0] woodbury_inv = backsub_both_sides(Lm, np.eye(num_inducing) - Bi) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index 44f3700d..6bcddc3e 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -17,7 +17,8 @@ import os from config import * if np.all(np.float64((scipy.__version__).split('.')[:2]) >= np.array([0, 12])): - import scipy.linalg.lapack as lapack + #import scipy.linalg.lapack.clapack as lapack + from scipy.linalg import lapack else: from scipy.linalg.lapack import flapack as lapack @@ -41,42 +42,103 @@ else: _blas_available = False warnings.warn("warning: caught this exception:" + str(e)) -def dtrtri(L, lower=0): +def force_F_ordered_symmetric(A): """ - Wrapper for lapack dtrtrs function + return a F ordered version of A, assuming A is symmetric + """ + if A.flags['F_CONTIGUOUS']: + return A + if A.flags['C_CONTIGUOUS']: + return A.T + else: + return np.asfortranarray(A) + +def force_F_ordered(A): + """ + return a F ordered version of A, assuming A is triangular + """ + if A.flags['F_CONTIGUOUS']: + return A + print "why are your arrays not F order?" + return np.asfortranarray(A) + +def jitchol(A, maxtries=5): + A = force_F_ordered_symmetric(A) + L, info = lapack.dpotrf(A, lower=1) + if info == 0: + return L + else: + if maxtries==0: + raise linalg.LinAlgError, "not positive definite, even with jitter." + diagA = np.diag(A) + if np.any(diagA <= 0.): + raise linalg.LinAlgError, "not pd: non-positive diagonal elements" + jitter = diagA.mean() * 1e-6 + + return jitchol(A+np.eye(A.shape[0])*jitter, maxtries-1) + +#def jitchol(A, maxtries=5): +# A = np.ascontiguousarray(A) +# L, info = lapack.dpotrf(A, lower=1) +# if info == 0: +# return L +# else: +# diagA = np.diag(A) +# if np.any(diagA <= 0.): +# raise linalg.LinAlgError, "not pd: non-positive diagonal elements" +# jitter = diagA.mean() * 1e-6 +# while maxtries > 0 and np.isfinite(jitter): +# print 'Warning: adding jitter of {:.10e}'.format(jitter) +# try: +# return linalg.cholesky(A + np.eye(A.shape[0]).T * jitter, lower=True) +# except: +# jitter *= 10 +# finally: +# maxtries -= 1 +# raise linalg.LinAlgError, "not positive definite, even with jitter." +# + + + + +def dtrtri(L, lower=1): + """ + Wrapper for lapack dtrtri function Inverse of L :param L: Triangular Matrix L :param lower: is matrix lower (true) or upper (false) :returns: Li, info """ + L = force_F_ordered(L) return lapack.dtrtri(L, lower=lower) -def dtrtrs(A, B, lower=0, trans=0, unitdiag=0): +def dtrtrs(A, B, lower=1, trans=0, unitdiag=0): """ Wrapper for lapack dtrtrs function - :param A: Matrix A + :param A: Matrix A(triangular) :param B: Matrix B :param lower: is matrix lower (true) or upper (false) :returns: """ + A = force_F_ordered(A) + #Note: B does not seem to need to be F ordered! return lapack.dtrtrs(A, B, lower=lower, trans=trans, unitdiag=unitdiag) -def dpotrs(A, B, lower=0): +def dpotrs(A, B, lower=1): """ Wrapper for lapack dpotrs function - :param A: Matrix A :param B: Matrix B :param lower: is matrix lower (true) or upper (false) :returns: - """ + A = force_F_ordered(A) return lapack.dpotrs(A, B, lower=lower) -def dpotri(A, lower=0): +def dpotri(A, lower=1): """ Wrapper for lapack dpotri function @@ -85,7 +147,12 @@ def dpotri(A, lower=0): :returns: A inverse """ - return lapack.dpotri(A, lower=lower) + assert lower==1, "scipy linalg behaviour is very weird. please use lower, fortran ordered arrays" + + A = force_F_ordered(A) + R, info = lapack.dpotri(A, lower=0) + symmetrify(R) + return R, info def pddet(A): """ @@ -133,60 +200,8 @@ def _mdot_r(a, b): b = b[0] return np.dot(a, b) -def jitchol(A, maxtries=5): - A = np.asfortranarray(A) - L, info = lapack.dpotrf(A, lower=1) - if info == 0: - return L - else: - diagA = np.diag(A) - if np.any(diagA <= 0.): - raise linalg.LinAlgError, "not pd: non-positive diagonal elements" - jitter = diagA.mean() * 1e-6 - while maxtries > 0 and np.isfinite(jitter): - print 'Warning: adding jitter of {:.10e}'.format(jitter) - try: - return linalg.cholesky(A + np.eye(A.shape[0]).T * jitter, lower=True) - except: - jitter *= 10 - finally: - maxtries -= 1 - raise linalg.LinAlgError, "not positive definite, even with jitter." - - - -def jitchol_old(A, maxtries=5): - """ - :param A: An almost pd square matrix - - :rval L: the Cholesky decomposition of A - - .. note: - - Adds jitter to K, to enforce positive-definiteness - if stuff breaks, please check: - np.allclose(sp.linalg.cholesky(XXT, lower = True), np.triu(sp.linalg.cho_factor(XXT)[0]).T) - - """ - try: - return linalg.cholesky(A, lower=True) - except linalg.LinAlgError: - diagA = np.diag(A) - if np.any(diagA < 0.): - raise linalg.LinAlgError, "not pd: negative diagonal elements" - jitter = diagA.mean() * 1e-6 - for i in range(1, maxtries + 1): - print '\rWarning: adding jitter of {:.10e} '.format(jitter), - try: - return linalg.cholesky(A + np.eye(A.shape[0]).T * jitter, lower=True) - except: - jitter *= 10 - - raise linalg.LinAlgError, "not positive definite, even with jitter." - def pdinv(A, *args): """ - :param A: A DxD pd numpy array :rval Ai: the inverse of A @@ -201,7 +216,7 @@ def pdinv(A, *args): """ L = jitchol(A, *args) logdet = 2.*np.sum(np.log(np.diag(L))) - Li = chol_inv(L) + Li = dtrtri(L) Ai, _ = lapack.dpotri(L) # Ai = np.tril(Ai) + np.tril(Ai,-1).T symmetrify(Ai) @@ -209,7 +224,7 @@ def pdinv(A, *args): return Ai, L, Li, logdet -def chol_inv(L): +def dtrtri(L): """ Inverts a Cholesky lower triangular matrix @@ -218,7 +233,8 @@ def chol_inv(L): """ - return lapack.dtrtri(L, lower=True)[0] + L = force_F_ordered(L) + return lapack.dtrtri(L, lower=1)[0] def multiple_pdinv(A): @@ -234,7 +250,7 @@ def multiple_pdinv(A): N = A.shape[-1] chols = [jitchol(A[:, :, i]) for i in range(N)] halflogdets = [np.sum(np.log(np.diag(L[0]))) for L in chols] - invs = [lapack.dpotri(L[0], True)[0] for L in chols] + invs = [dpotri(L[0], True)[0] for L in chols] invs = [np.triu(I) + np.triu(I, 1).T for I in invs] return np.dstack(invs), np.array(halflogdets) @@ -425,7 +441,8 @@ def tdot_blas(mat, out=None): symmetrify(out, upper=True) - return out + + return np.ascontiguousarray(out) def tdot(*args, **kwargs): if _blas_available: @@ -557,24 +574,24 @@ def cholupdate(L, x): def backsub_both_sides(L, X, transpose='left'): """ Return L^-T * X * L^-1, assumuing X is symmetrical and L is lower cholesky""" if transpose == 'left': - tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=1) - return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=1)[0].T + tmp, _ = dtrtrs(L, X, lower=1, trans=1) + return dtrtrs(L, tmp.T, lower=1, trans=1)[0].T else: - tmp, _ = lapack.dtrtrs(L, np.asfortranarray(X), lower=1, trans=0) - return lapack.dtrtrs(L, np.asfortranarray(tmp.T), lower=1, trans=0)[0].T + tmp, _ = dtrtrs(L, X, lower=1, trans=0) + return dtrtrs(L, tmp.T, lower=1, trans=0)[0].T def PCA(Y, input_dim): """ -Principal component analysis: maximum likelihood solution by SVD + Principal component analysis: maximum likelihood solution by SVD -:param Y: NxD np.array of data -:param input_dim: int, dimension of projection + :param Y: NxD np.array of data + :param input_dim: int, dimension of projection -:rval X: - Nxinput_dim np.array of dimensionality reduced data -:rval W: - input_dimxD mapping from X to Y + :rval X: - Nxinput_dim np.array of dimensionality reduced data + :rval W: - input_dimxD mapping from X to Y -""" + """ if not np.allclose(Y.mean(axis=0), 0.0): print "Y is not zero mean, centering it locally (GPy.util.linalg.PCA)" From b1c98c2c3d901dd6fe630c7ee42c4249db9e2da2 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 13:18:48 +0000 Subject: [PATCH 27/40] Priors added --- GPy/__init__.py | 3 +- GPy/core/__init__.py | 1 - GPy/core/model.py | 138 +++++------ GPy/core/parameterization/index_operations.py | 6 +- GPy/core/parameterization/param.py | 224 +++--------------- GPy/core/parameterization/parameter_core.py | 113 +++++++-- GPy/core/parameterization/parameterized.py | 71 ++---- GPy/core/parameterization/priors.py | 6 +- 8 files changed, 206 insertions(+), 356 deletions(-) diff --git a/GPy/__init__.py b/GPy/__init__.py index 94128dfa..b6bf81b7 100644 --- a/GPy/__init__.py +++ b/GPy/__init__.py @@ -4,7 +4,7 @@ import warnings warnings.filterwarnings("ignore", category=DeprecationWarning) import core -from core.parameterization import transformations +from core.parameterization import transformations, priors import models import mappings import inference @@ -15,7 +15,6 @@ import testing from numpy.testing import Tester from nose.tools import nottest import kern -from core import priors import plotting @nottest diff --git a/GPy/core/__init__.py b/GPy/core/__init__.py index 39887284..839529d6 100644 --- a/GPy/core/__init__.py +++ b/GPy/core/__init__.py @@ -2,7 +2,6 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) from model import * -from parameterization import priors from parameterization.parameterized import * from gp import GP from sparse_gp import SparseGP diff --git a/GPy/core/model.py b/GPy/core/model.py index c944f943..58008a99 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -21,8 +21,6 @@ class Model(Parameterized): _allowed_failures = 10 # number of allowed failures def __init__(self, name): super(Model, self).__init__(name) # Parameterized.__init__(self) - self.priors = [] - self._priors = ParameterIndexOperations() self.optimization_runs = [] self.sampling_runs = [] self.preferred_optimizer = 'scg' @@ -67,66 +65,66 @@ 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 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): """ @@ -146,22 +144,6 @@ class Model(Parameterized): else: raise AttributeError, "no parameter matches %s" % name - def log_prior(self): - """evaluate the prior""" - if self.priors is not None: - return np.sum([p.lnpdf(x) for p, x in zip(self.priors, self._get_params()) if p is not None]) - else: - return 0. - - def _log_prior_gradients(self): - """evaluate the gradients of the priors""" - if self.priors is None: - return 0. - x = self._get_params() - ret = np.zeros(x.size) - [np.put(ret, i, p.lnpdf_grad(xx)) for i, (p, xx) in enumerate(zip(self.priors, x)) if not p is None] - return ret - def randomize(self): """ Randomize the model. diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index b816e05f..d5061be7 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -107,7 +107,7 @@ class ParameterIndexOperations(object): def add(self, prop, indices): try: self._properties[prop] = combine_indices(self._properties[prop], indices) - except KeyError: + except KeyError: self._properties[prop] = indices def remove(self, prop, indices): @@ -215,10 +215,10 @@ class ParameterIndexOperationsView(object): def remove(self, prop, indices): removed = self._param_index_ops.remove(prop, indices+self._offset) - if removed.size > 0: - return removed - self._size + 1 if self[prop].size == 0: del self[prop] + if removed.size > 0: + return removed - self._size + 1 return removed diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 97fcfd32..a52f2e18 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -10,6 +10,7 @@ from array_core import ObservableArray, ParamList __constraints_name__ = "Constraint" __index_name__ = "Index" __tie_name__ = "Tied to" +__priors_name__ = "Prior" __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision used, sublassing numpy ndarray after all __print_threshold__ = 5 ###### @@ -77,9 +78,10 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri self._name = getattr(obj, 'name', None) self.gradient = getattr(obj, 'gradient', None) self.constraints = getattr(obj, 'constraints', None) + self.priors = getattr(obj, 'priors', None) - def __array_wrap__(self, out_arr, context=None): - return out_arr.view(numpy.ndarray) + #def __array_wrap__(self, out_arr, context=None): + # return out_arr.view(numpy.ndarray) #=========================================================================== # Pickling operations #=========================================================================== @@ -118,154 +120,17 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri #=========================================================================== # get/set parameters #=========================================================================== - def _set_params(self, param, update=True): self.flat = param - self._notify_tied_parameters() + #self._notify_tied_parameters() self._notify_observers() def _get_params(self): return self.flat -# @property -# def name(self): -# """ -# Name of this parameter. -# This can be a callable without parameters. The callable will be called -# every time the name property is accessed. -# """ -# if callable(self.name): -# return self.name() -# return self.name -# @name.setter -# def name(self, new_name): -# from_name = self.name -# self.name = new_name -# self._direct_parent_._name_changed(self, from_name) + def _collect_gradient(self, target): target[:] = self.gradient.flat - #=========================================================================== - # Tying operations -> bugged, TODO - #=========================================================================== - def tie_to(self, param): - """ - :param param: the parameter object to tie this parameter to. - Can be ParamConcatenation (retrieved by regexp search) - Tie this parameter to the given parameter. - Broadcasting is not allowed, but you can tie a whole dimension to - one parameter: self[:,0].tie_to(other), where other is a one-value - parameter. - - Note: For now only one parameter can have ties, so all of a parameter - will be removed, when re-tieing! - """ - #Note: this method will tie to the parameter which is the last in - # the chain of ties. Thus, if you tie to a tied parameter, - # this tie will be created to the parameter the param is tied - # to. - - assert isinstance(param, Param), "Argument {1} not of type {0}".format(Param, param.__class__) - param = numpy.atleast_1d(param) - if param.size != 1: - raise NotImplementedError, "Broadcast tying is not implemented yet" - try: - if self._original_: - self[:] = param - else: # this happens when indexing created a copy of the array - self._direct_parent_._get_original(self)[self._current_slice_] = param - except ValueError: - raise ValueError("Trying to tie {} with shape {} to {} with shape {}".format(self.name, self.shape, param.name, param.shape)) - if param is self: - raise RuntimeError, 'Cyclic tieing is not allowed' -# if len(param._tied_to_) > 0: -# if (self._direct_parent_._get_original(self) is param._direct_parent_._get_original(param) -# and len(set(self._raveled_index())&set(param._tied_to_[0]._raveled_index()))!=0): -# raise RuntimeError, 'Cyclic tieing is not allowed' -# self.tie_to(param._tied_to_[0]) -# return - if not param in self._direct_parent_._get_original(self)._tied_to_: - self._direct_parent_._get_original(self)._tied_to_ += [param] - param._add_tie_listener(self) - self._highest_parent_._set_fixed(self) - cs = self._highest_parent_._constraints_for(param, param._raveled_index()) - for cs in self._highest_parent_._constraints_for(param, param._raveled_index()): - [self.constrain(c, warning=False) for c in cs] -# for t in self._tied_to_me_.keys(): -# if t is not self: -# t.untie(self) -# t.tie_to(param) - - def untie(self, *ties): - """ - remove all ties. - """ - [t._direct_parent_._get_original(t)._remove_tie_listener(self) for t in self._tied_to_] - new_ties = [] - for t in self._direct_parent_._get_original(self)._tied_to_: - for tied in t._tied_to_me_.keys(): - if t._parent_index_ is tied._parent_index_: - new_ties.append(tied) - self._direct_parent_._get_original(self)._tied_to_ = new_ties - self._direct_parent_._get_original(self)._highest_parent_._set_unfixed(self) -# self._direct_parent_._remove_tie(self, *params) - def _notify_tied_parameters(self): - for tied, ind in self._tied_to_me_.iteritems(): - tied._on_tied_parameter_changed(self.base, list(ind)) - def _add_tie_listener(self, tied_to_me): - for t in self._tied_to_me_.keys(): - if tied_to_me._parent_index_ is t._parent_index_: - t_rav_i = t._raveled_index() - tr_rav_i = tied_to_me._raveled_index() - new_index = list(set(t_rav_i) | set(tr_rav_i)) - tmp = t._direct_parent_._get_original(t)[numpy.unravel_index(new_index, t._realshape_)] - self._tied_to_me_[tmp] = self._tied_to_me_[t] | set(self._raveled_index()) - del self._tied_to_me_[t] - return - self._tied_to_me_[tied_to_me] = set(self._raveled_index()) - def _remove_tie_listener(self, to_remove): - for t in self._tied_to_me_.keys(): - if t._parent_index_ == to_remove._parent_index_: - t_rav_i = t._raveled_index() - tr_rav_i = to_remove._raveled_index() - import ipdb;ipdb.set_trace() - new_index = list(set(t_rav_i) - set(tr_rav_i)) - if new_index: - tmp = t._direct_parent_._get_original(t)[numpy.unravel_index(new_index, t._realshape_)] - self._tied_to_me_[tmp] = self._tied_to_me_[t] - del self._tied_to_me_[t] - if len(self._tied_to_me_[tmp]) == 0: - del self._tied_to_me_[tmp] - else: - del self._tied_to_me_[t] - def _on_tied_parameter_changed(self, val, ind): - if not self._updated_: # not fast_array_equal(self, val[ind]): - val = numpy.atleast_1d(val) - self._updated_ = True - if self._original_: - self.__setitem__(slice(None), val[ind], update=False) - else: # this happens when indexing created a copy of the array - self._direct_parent_._get_original(self).__setitem__(self._current_slice_, val[ind], update=False) - self._notify_tied_parameters() - self._updated_ = False - #=========================================================================== - # Prior Operations - #=========================================================================== - def set_prior(self, prior): - """ - :param prior: prior to be set for this parameter - - Set prior for this parameter. - """ - if not hasattr(self._highest_parent_, '_set_prior'): - raise AttributeError("Parent of type {} does not support priors".format(self._highest_parent_.__class__)) - self._highest_parent_._set_prior(self, prior) - def unset_prior(self, *priors): - """ - :param priors: priors to remove from this parameter - - Remove all priors from this parameter - """ - self._highest_parent_._remove_prior(self, *priors) #=========================================================================== # Array operations -> done #=========================================================================== @@ -283,6 +148,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri self._notify_tied_parameters() if update: self._highest_parent_.parameters_changed() + #=========================================================================== # Index Operations: #=========================================================================== @@ -328,8 +194,9 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri return numpy.r_[a] return numpy.r_[:b] return itertools.imap(f, itertools.izip_longest(slice_index[:self._realndim_], self._realshape_, fillvalue=slice(self.size))) + #=========================================================================== - # Convienience + # Convenience #=========================================================================== @property def is_fixed(self): @@ -338,11 +205,10 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri view = super(Param, self).round(decimals, out).view(Param) view.__array_finalize__(self) return view - def _has_fixes(self): - return False round.__doc__ = numpy.round.__doc__ def _get_original(self, param): return self + #=========================================================================== # Printing -> done #=========================================================================== @@ -362,6 +228,9 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri def _constraints_str(self): return [' '.join(map(lambda c: str(c[0]) if c[1].size == self._realsize_ else "{" + str(c[0]) + "}", self.constraints.iteritems()))] @property + def _priors_str(self): + return [' '.join(map(lambda c: str(c[0]) if c[1].size == self._realsize_ else "{" + str(c[0]) + "}", self.priors.iteritems()))] + @property def _ties_str(self): return [t._short() for t in self._tied_to_] or [''] @property @@ -385,8 +254,6 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri if len(ind) != 1: ties[i, matches[0][ind_rav_matches]] = numpy.take(tt_rav_index, matches[1], mode='wrap')[ind_rav_matches] else: ties[i, matches[0]] = numpy.take(tt_rav_index, matches[1], mode='wrap') return map(lambda a: sum(a, []), zip(*[[[tie.flatten()] if tx != None else [] for tx in t] for t, tie in zip(ties, self._tied_to_)])) - def _constraints_for(self, rav_index): - return self.constraints.properties_for(rav_index) def _indices(self, slice_index=None): # get a int-array containing all indices in the first axis. if slice_index is None: @@ -404,6 +271,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri return numpy.fromiter(itertools.product(*expanded_index), dtype=[('', int)] * self._realndim_, count=reduce(lambda a, b: a * b.size, expanded_index, 1)).view((int, self._realndim_)) def _max_len_names(self, gen, header): + 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.name_hirarchical)) @@ -418,21 +286,26 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri if ind.size > 4: indstr = ','.join(map(str, ind[:2])) + "..." + ','.join(map(str, ind[-2:])) else: indstr = ','.join(map(str, ind)) return name + '[' + indstr + ']' - def __str__(self, constr_matrix=None, indices=None, ties=None, lc=None, lx=None, li=None, lt=None): + def __str__(self, constr_matrix=None, indices=None, prirs=None, ties=None, lc=None, lx=None, li=None, lp=None, lt=None, only_name=False): filter_ = self._current_slice_ vals = self.flat if indices is None: indices = self._indices(filter_) ravi = self._raveled_index(filter_) - if constr_matrix is None: constr_matrix = self._constraints_for(ravi) + if constr_matrix is None: constr_matrix = self.constraints.properties_for(ravi) + if prirs is None: prirs = self.priors.properties_for(ravi) if ties is None: ties = self._ties_for(ravi) ties = [' '.join(map(lambda x: x._short(), t)) for t in ties] if lc is None: lc = self._max_len_names(constr_matrix, __constraints_name__) if lx is None: lx = self._max_len_values() if li is None: li = self._max_len_index(indices) if lt is None: lt = self._max_len_names(ties, __tie_name__) - header = " {i:^{2}s} | \033[1m{x:^{1}s}\033[0;0m | {c:^{0}s} | {t:^{3}s}".format(lc, lx, li, lt, x=self.name_hirarchical, c=__constraints_name__, i=__index_name__, t=__tie_name__) # nice header for printing + 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.name_hirarchical, 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.name_hirarchical, 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} | {t:^{4}s} ".format(lc, lx, __precision__, li, lt, x=x, c=" ".join(map(str, c)), t=(t or ''), i=i) for i, x, c, t in itertools.izip(indices, vals, constr_matrix, ties)]) # return all the constraints with right indices + 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__() class ParamConcatenation(object): @@ -538,53 +411,20 @@ class ParamConcatenation(object): def __str__(self, *args, **kwargs): def f(p): ind = p._raveled_index() - return p._constraints_for(ind), p._ties_for(ind) + return p.constraints.properties_for(ind), p._ties_for(ind), p.priors.properties_for(ind) params = self.params - constr_matrices, ties_matrices = zip(*map(f, params)) + constr_matrices, ties_matrices, prior_matrices = zip(*map(f, params)) indices = [p._indices() for p in params] lc = max([p._max_len_names(cm, __constraints_name__) for p, cm in itertools.izip(params, constr_matrices)]) lx = max([p._max_len_values() for p in params]) li = max([p._max_len_index(i) for p, i in itertools.izip(params, indices)]) lt = max([p._max_len_names(tm, __tie_name__) for p, tm in itertools.izip(params, ties_matrices)]) - strings = [p.__str__(cm, i, tm, lc, lx, li, lt) for p, cm, i, tm in itertools.izip(params,constr_matrices,indices,ties_matrices)] + lp = max([p._max_len_names(pm, __constraints_name__) for p, pm in itertools.izip(params, prior_matrices)]) + strings = [] + start = True + for p, cm, i, tm, pm in itertools.izip(params,constr_matrices,indices,ties_matrices,prior_matrices): + strings.append(p.__str__(constr_matrix=cm, indices=i, prirs=pm, ties=tm, lc=lc, lx=lx, li=li, lp=lp, lt=lt, only_name=(1-start))) + start = False return "\n".join(strings) - return "\n{}\n".format(" -"+"- | -".join(['-'*l for l in [li,lx,lc,lt]])).join(strings) def __repr__(self): - return "\n".join(map(repr,self.params)) - -if __name__ == '__main__': - - - from GPy.core.parameterized import Parameterized - from GPy.core.parameter import Param - - #X = numpy.random.randn(2,3,1,5,2,4,3) - X = numpy.random.randn(3,2) - print "random done" - p = Param("q_mean", X) - p1 = Param("q_variance", numpy.random.rand(*p.shape)) - p2 = Param("Y", numpy.random.randn(p.shape[0], 1)) - - p3 = Param("variance", numpy.random.rand()) - p4 = Param("lengthscale", numpy.random.rand(2)) - - m = Parameterized() - rbf = Parameterized(name='rbf') - - rbf.add_parameter(p3,p4) - m.add_parameter(p,p1,rbf) - - print "setting params" - #print m.q_v[3:5,[1,4,5]] - print "constraining variance" - #m[".*variance"].constrain_positive() - #print "constraining rbf" - #m.rbf_l.constrain_positive() - #m.q_variance[1,[0,5,11,19,2]].tie_to(m.rbf_v) - #m.rbf_v.tie_to(m.rbf_l[0]) - #m.rbf_l[0].tie_to(m.rbf_l[1]) - #m.q_v.tie_to(m.rbf_v) -# m.rbf_l.tie_to(m.rbf_va) - # pt = numpy.array(params._get_params_transformed()) - # ptr = numpy.random.randn(*pt.shape) -# params.X.tie_to(params.rbf_v) + return "\n".join(map(repr,self.params)) \ No newline at end of file diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 7ad228f7..7e5fea73 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -29,7 +29,25 @@ class Parameterizable(object): if add_name: return [adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] return [xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] - + + def _collect_gradient(self, target): + import itertools + [p._collect_gradient(target[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], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + self.parameters_changed() + + def parameters_changed(self): """ This method gets called when parameters have changed. @@ -130,17 +148,18 @@ class Indexable(object): that is an int array, containing the indexes for the flattened param inside this parameterized logic. """ - raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization object?" - + raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization object?" + class Constrainable(Nameable, Indexable, Parameterizable): def __init__(self, name, default_constraint=None): super(Constrainable,self).__init__(name) self._default_constraint_ = default_constraint from index_operations import ParameterIndexOperations self.constraints = ParameterIndexOperations() + self.priors = ParameterIndexOperations() if self._default_constraint_ is not None: self.constrain(self._default_constraint_) - + #=========================================================================== # Fixing Parameters: #=========================================================================== @@ -186,16 +205,41 @@ class Constrainable(Nameable, Indexable, Parameterizable): self._fixes_[fixed_indices] = FIXED else: self._fixes_ = None + + def _has_fixes(self): + return hasattr(self, "_fixes_") and self._fixes_ is not None + #=========================================================================== + # Prior Operations + #=========================================================================== + def set_prior(self, prior, warning=True, update=True): + repriorized = self.unset_priors() + self._add_to_index_operations(self.priors, repriorized, prior, warning, update) + + def unset_priors(self, *priors): + return self._remove_from_index_operations(self.priors, priors) + + def log_prior(self): + """evaluate the prior""" + import numpy as np + if self.priors.size > 0: + x = self._get_params() + return np.sum([p.lnpdf(x[ind]) for p, ind in self.priors.iteritems()]) + return 0. + + def _log_prior_gradients(self): + """evaluate the gradients of the priors""" + import numpy as np + if self.priors.size > 0: + x = self._get_params() + ret = np.zeros(x.size) + [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.iteritems()] + return ret + return 0. + #=========================================================================== # Constrain operations -> done #=========================================================================== - def _parent_changed(self, parent): - from index_operations import ParameterIndexOperationsView - self.constraints = ParameterIndexOperationsView(parent.constraints, parent._offset_for(self), self.size) - self._fixes_ = None - for p in self._parameters_: - p._parent_changed(parent) def constrain(self, transform, warning=True, update=True): """ @@ -209,11 +253,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): if isinstance(transform, Transformation): self._set_params(transform.initialize(self._get_params()), update=False) reconstrained = self.unconstrain() - self.constraints.add(transform, self._raveled_index()) - if warning and reconstrained.size > 0: - print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) - if update: - self._highest_parent_.parameters_changed() + self._add_to_index_operations(self.constraints, reconstrained, transform, warning, update) def unconstrain(self, *transforms): """ @@ -222,16 +262,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): remove all :py:class:`GPy.core.transformations.Transformation` transformats of this parameter object. """ - if len(transforms) == 0: - transforms = self.constraints.properties() - import numpy as np - removed = np.empty((0,),dtype=int) - for t in transforms: - unconstrained = self.constraints.remove(t, self._raveled_index()) - removed = np.union1d(removed, unconstrained) - if t is __fixed__: - self._highest_parent_._set_unfixed(unconstrained) - return removed + return self._remove_from_index_operations(self.constraints, transforms) def constrain_positive(self, warning=True, update=True): """ @@ -277,3 +308,35 @@ class Constrainable(Nameable, Indexable, Parameterizable): Remove (lower, upper) bounded constrain from this parameter/ """ self.unconstrain(Logistic(lower, upper)) + + def _parent_changed(self, parent): + from index_operations import ParameterIndexOperationsView + self.constraints = ParameterIndexOperationsView(parent.constraints, parent._offset_for(self), self.size) + self.priors = ParameterIndexOperationsView(parent.priors, parent._offset_for(self), self.size) + self._fixes_ = None + for p in self._parameters_: + p._parent_changed(parent) + + def _add_to_index_operations(self, which, reconstrained, transform, warning, update): + 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._highest_parent_.parameters_changed() + + def _remove_from_index_operations(self, which, transforms): + if len(transforms) == 0: + transforms = which.properties() + import numpy as np + removed = np.empty((0, ), dtype=int) + for t in transforms: + unconstrained = which.remove(t, self._raveled_index()) + removed = np.union1d(removed, unconstrained) + if t is __fixed__: + self._highest_parent_._set_unfixed(unconstrained) + + return removed + + + + diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 593c37e9..c8a841c0 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -66,9 +66,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_ = set() del self._in_init_ - def _has_fixes(self): - return hasattr(self, "_fixes_") and self._fixes_ is not None - def add_parameter(self, param, index=None): """ :param parameters: the parameters to add @@ -88,11 +85,14 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): # make sure the size is set if index is None: self.constraints.update(param.constraints, self.size) + self.priors.update(param.priors, self.size) self._parameters_.append(param) else: start = sum(p.size for p in self._parameters_[:index]) self.constraints.shift(start, param.size) + self.priors.shift(start, param.size) self.constraints.update(param.constraints, start) + self.priors.update(param.priors, start) self._parameters_.insert(index, param) self.size += param.size else: @@ -210,6 +210,7 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): """ return [ self._fixes_, + self.priors, self.constraints, self._parameters_, self._name, @@ -220,9 +221,10 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): self._added_names_ = state.pop() self._name = state.pop() self._parameters_ = state.pop() - self._connect_parameters() self.constraints = state.pop() + self.priors = state.pop() self._fixes_ = state.pop() + self._connect_parameters() self.parameters_changed() #=========================================================================== # Gradient control @@ -248,16 +250,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): if self._has_fixes(): return n[self._fixes_] return n - def _get_params(self): - # don't overwrite this anymore! - if not self.size: - return np.empty(shape=(0,), dtype=np.float64) - return numpy.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! - [p._set_params(params[s], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - self.parameters_changed() def _get_params_transformed(self): # transformed parameters (apply transformation rules) p = self._get_params() @@ -348,36 +340,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self._direct_parent_.hirarchy_name() + adjust_name_for_printing(self.name) + "." return '' #=========================================================================== - # Constraint Handling: - #=========================================================================== - #=========================================================================== - # def _add_constrain(self, param, transform, warning=True): - # rav_i = self._raveled_index_for(param) - # reconstrained = self._remove_constrain(param, index=rav_i) # remove constraints before - # # if removing constraints before adding new is not wanted, just delete the above line! - # self.constraints.add(transform, rav_i) - # param = self._get_original(param) - # if not (transform == __fixed__): - # param._set_params(transform.initialize(param._get_params()), update=False) - # if warning and any(reconstrained): - # # if you want to print the whole params object, which was reconstrained use: - # # m = str(param[self._backtranslate_index(param, reconstrained)]) - # print "Warning: re-constraining parameters:\n{}".format(param._short()) - # return rav_i - # def _remove_constrain(self, param, *transforms, **kwargs): - # if not transforms: - # transforms = self.constraints.properties() - # removed_indices = numpy.array([]).astype(int) - # if "index" in kwargs: index = kwargs['index'] - # else: index = self._raveled_index_for(param) - # for constr in transforms: - # removed = self.constraints.remove(constr, index) - # if constr is __fixed__: - # self._set_unfixed(removed) - # removed_indices = numpy.union1d(removed_indices, removed) - # return removed_indices - #=========================================================================== - #=========================================================================== # Get/set parameters: #=========================================================================== def grep_param_names(self, regexp): @@ -434,9 +396,6 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): return self._direct_parent_.hirarchy_name() + adjust_name_for_printing(self.name) else: return adjust_name_for_printing(self.name) - #parameter_names = property(parameter_names, doc="Names for all parameters handled by this parameterization object -- will add hirarchy name entries for printing") - def _collect_gradient(self, target): - [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] @property def flattened_parameters(self): return [xi for x in self._parameters_ for xi in x.flattened_parameters] @@ -455,6 +414,9 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): def _constraints_str(self): return [cs for p in self._parameters_ for cs in p._constraints_str] @property + def _priors_str(self): + return [cs for p in self._parameters_ for cs in p._priors_str] + @property def _description_str(self): return [xi for x in self._parameters_ for xi in x._description_str] @property @@ -463,20 +425,23 @@ class Parameterized(Constrainable, Pickleable, Observable, Gradcheckable): def __str__(self, header=True): name = adjust_name_for_printing(self.name) + "." - constrs = self._constraints_str; ts = self._ties_str + constrs = self._constraints_str; + ts = self._ties_str + prirs = self._priors_str desc = self._description_str; names = self.parameter_names() nl = max([len(str(x)) for x in names + [name]]) sl = max([len(str(x)) for x in desc + ["Value"]]) cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]]) tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]]) - format_spec = " \033[1m{{name:<{0}s}}\033[0;0m | {{desc:^{1}s}} | {{const:^{2}s}} | {{t:^{3}s}}".format(nl, sl, cl, tl) + pl = max([len(str(x)) if x else 0 for x in prirs + ["Prior"]]) + format_spec = " \033[1m{{name:<{0}s}}\033[0;0m | {{desc:^{1}s}} | {{const:^{2}s}} | {{pri:^{3}s}} | {{t:^{4}s}}".format(nl, sl, cl, pl, tl) to_print = [] - for n, d, c, t in itertools.izip(names, desc, constrs, ts): - to_print.append(format_spec.format(name=n, desc=d, const=c, t=t)) + for n, d, c, t, p in itertools.izip(names, desc, constrs, ts, prirs): + to_print.append(format_spec.format(name=n, desc=d, const=c, t=t, pri=p)) # to_print = [format_spec.format(p=p, const=c, t=t) if isinstance(p, Param) else p.__str__(header=False) for p, c, t in itertools.izip(self._parameters_, constrs, ts)] - sep = '-' * (nl + sl + cl + tl + 8 * 2 + 3) + sep = '-' * (nl + sl + cl + + pl + tl + 8 * 2 + 3) if header: - header = " {{0:<{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}}".format(nl, sl, cl, tl).format(name, "Value", "Constraint", "Tied to") + header = " {{0:<{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}} | {{4:^{4}s}}".format(nl, sl, cl, pl, tl).format(name, "Value", "Constraint", "Prior", "Tied to") # header += '\n' + sep to_print.insert(0, header) return '\n'.format(sep).join(to_print) diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py index f1208f18..3155be64 100644 --- a/GPy/core/parameterization/priors.py +++ b/GPy/core/parameterization/priors.py @@ -9,15 +9,16 @@ from domains import _REAL, _POSITIVE import warnings import weakref -class Prior: +class Prior(object): domain = None def pdf(self, x): return np.exp(self.lnpdf(x)) def plot(self): + import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." - from ..plotting.matplot_dep import priors_plots + from ...plotting.matplot_dep import priors_plots priors_plots.univariate_plot(self) @@ -150,6 +151,7 @@ class MultivariateGaussian: return np.random.multivariate_normal(self.mu, self.var, n) def plot(self): + import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ..plotting.matplot_dep import priors_plots priors_plots.multivariate_plot(self) From 0919544eb7b26501cd1304a6bcbfe0784defa9e7 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 13:21:44 +0000 Subject: [PATCH 28/40] oops index operations had an assignment error --- GPy/core/parameterization/index_operations.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index d5061be7..bfd0bf21 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -215,8 +215,6 @@ class ParameterIndexOperationsView(object): def remove(self, prop, indices): removed = self._param_index_ops.remove(prop, indices+self._offset) - if self[prop].size == 0: - del self[prop] if removed.size > 0: return removed - self._size + 1 return removed From f21e8e2394140cced9768fa74ff7c981e604a157 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 13:43:43 +0000 Subject: [PATCH 29/40] roundtrip error fixed for likelihood tests --- GPy/core/model.py | 3 +-- GPy/core/parameterization/param.py | 16 ++++++++-------- .../latent_function_inference/laplace.py | 2 +- 3 files changed, 10 insertions(+), 11 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 58008a99..55083aaf 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -154,8 +154,7 @@ class Model(Parameterized): x = np.random.randn(self.size_transformed) x = self._untransform_params(x) # now draw from prior where possible - if self.priors is not None and len(self.priors): - [np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) if not p is None] + [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...) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index a52f2e18..2751bb14 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -145,7 +145,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri return new_arr def __setitem__(self, s, val, update=True): super(Param, self).__setitem__(s, val, update=update) - self._notify_tied_parameters() + #self._notify_tied_parameters() if update: self._highest_parent_.parameters_changed() @@ -201,11 +201,11 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri @property def is_fixed(self): return self._highest_parent_._is_fixed(self) - def round(self, decimals=0, out=None): - view = super(Param, self).round(decimals, out).view(Param) - view.__array_finalize__(self) - return view - round.__doc__ = numpy.round.__doc__ + #def round(self, decimals=0, out=None): + # view = super(Param, self).round(decimals, out).view(Param) + # view.__array_finalize__(self) + # return view + #round.__doc__ = numpy.round.__doc__ def _get_original(self, param): return self @@ -337,8 +337,8 @@ class ParamConcatenation(object): def __setitem__(self, s, val, update=True): 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 p._notify_tied_parameters() - for p, ps in zip(self.params, self._param_slices_)] + [numpy.place(p, ind[ps], vals[ps])# and p._notify_tied_parameters() + for p, ps in zip(self.params, self._param_slices_)] if update: self.params[0]._highest_parent_.parameters_changed() def _vals(self): diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index 1b7bdad8..1c524d71 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -216,7 +216,7 @@ class Laplace(object): """ if not log_concave: #print "Under 1e-10: {}".format(np.sum(W < 1e-6)) - W[W < 1e-6] = 1e-6 # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur + W.__setitem__(W < 1e-6, 1e-6, update=False) # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur # If the likelihood is non-log-concave. We wan't to say that there is a negative variance # To cause the posterior to become less certain than the prior and likelihood, # This is a property only held by non-log-concave likelihoods From 6419401d60116bab030692b80c0bb525fd291c08 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 14:40:32 +0000 Subject: [PATCH 30/40] parameters changed more structured now, parameters changed goes from bottom to top, when calling _notify_parameters_changed() --- GPy/core/parameterization/param.py | 6 +- GPy/core/parameterization/parameter_core.py | 83 +++++++++++---------- 2 files changed, 48 insertions(+), 41 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 2751bb14..46fa3377 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -147,7 +147,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri super(Param, self).__setitem__(s, val, update=update) #self._notify_tied_parameters() if update: - self._highest_parent_.parameters_changed() + self._notify_parameters_changed() #=========================================================================== # Index Operations: @@ -340,14 +340,14 @@ class ParamConcatenation(object): [numpy.place(p, ind[ps], vals[ps])# and p._notify_tied_parameters() for p, ps in zip(self.params, self._param_slices_)] if update: - self.params[0]._highest_parent_.parameters_changed() + self.params[0]._notify_parameters_changed() def _vals(self): return numpy.hstack([p._get_params() for p in self.params]) #=========================================================================== # parameter operations: #=========================================================================== def update_all_params(self): - self.params[0]._highest_parent_.parameters_changed() + self.params[0]._notify_parameters_changed() 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 7e5fea73..9741b4ef 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -19,43 +19,6 @@ class Observable(object): del self._observers_[observer] def _notify_observers(self): [callble(self) for callble in self._observers_.itervalues()] - -class Parameterizable(object): - def __init__(self, *args, **kwargs): - from GPy.core.parameterization.array_core import ParamList - _parameters_ = ParamList() - - def parameter_names(self, add_name=False): - if add_name: - return [adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] - return [xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] - - def _collect_gradient(self, target): - import itertools - [p._collect_gradient(target[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], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - self.parameters_changed() - - - def parameters_changed(self): - """ - This method gets called when parameters have changed. - Another way of listening to param changes is to - add self as a listener to the param, such that - updates get passed through. See :py:function:``GPy.core.param.Observable.add_observer`` - """ - pass class Pickleable(object): def _getstate(self): @@ -121,6 +84,50 @@ class Nameable(Parentable): if self.has_parent(): self._direct_parent_._name_changed(self, from_name) + +class Parameterizable(Parentable): + def __init__(self, *args, **kwargs): + super(Parameterizable, self).__init__(*args, **kwargs) + from GPy.core.parameterization.array_core import ParamList + _parameters_ = ParamList() + + def parameter_names(self, add_name=False): + if add_name: + return [adjust_name_for_printing(self.name) + "." + xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] + return [xi for x in self._parameters_ for xi in x.parameter_names(add_name=True)] + + def _collect_gradient(self, target): + import itertools + [p._collect_gradient(target[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], update=update) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + self.parameters_changed() + + def parameters_changed(self): + """ + This method gets called when parameters have changed. + Another way of listening to param changes is to + add self as a listener to the param, such that + updates get passed through. See :py:function:``GPy.core.param.Observable.add_observer`` + """ + pass + + def _notify_parameters_changed(self): + self.parameters_changed() + if self.has_parent(): + self._direct_parent_._notify_parameters_changed() + + class Gradcheckable(Parentable): #=========================================================================== # Gradchecking @@ -322,7 +329,7 @@ class Constrainable(Nameable, Indexable, Parameterizable): print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) which.add(transform, self._raveled_index()) if update: - self._highest_parent_.parameters_changed() + self._notify_parameters_changed() def _remove_from_index_operations(self, which, transforms): if len(transforms) == 0: From 345e5b3e7c324c59bf4148aeeb56fe746574270d Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 15:05:38 +0000 Subject: [PATCH 31/40] some updates for params changes and likelihood fixes --- GPy/core/parameterization/array_core.py | 22 ++++++++++++++++--- GPy/core/parameterization/param.py | 2 +- .../latent_function_inference/laplace.py | 2 ++ 3 files changed, 22 insertions(+), 4 deletions(-) diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index 78b9e8f2..9f2c7ae6 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -37,10 +37,26 @@ class ObservableArray(np.ndarray, Observable): if obj is None: return self._observers_ = getattr(obj, '_observers_', None) + def _s_not_empty(self, s): + # this checks whether there is something picked by this slice. + return True + # TODO: disarmed, for performance increase, + if not isinstance(s, (list,tuple,np.ndarray)): + return True + if isinstance(s, (list,tuple)): + return len(s)!=0 + if isinstance(s, np.ndarray): + if s.dtype is bool: + return np.all(s) + else: + return s.size != 0 + def __setitem__(self, s, val, update=True): - super(ObservableArray, self).__setitem__(s, val) - if update: - self._notify_observers() + if self._s_not_empty(s): + super(ObservableArray, self).__setitem__(s, val) + if update: + self._notify_observers() + def __getslice__(self, start, stop): return self.__getitem__(slice(start, stop)) def __setslice__(self, start, stop, val): diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 46fa3377..37e710a8 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -146,7 +146,7 @@ class Param(ObservableArray, Constrainable, Gradcheckable, Indexable, Parameteri def __setitem__(self, s, val, update=True): super(Param, self).__setitem__(s, val, update=update) #self._notify_tied_parameters() - if update: + if update and self._s_not_empty(s): self._notify_parameters_changed() #=========================================================================== diff --git a/GPy/inference/latent_function_inference/laplace.py b/GPy/inference/latent_function_inference/laplace.py index 1c524d71..906a7867 100644 --- a/GPy/inference/latent_function_inference/laplace.py +++ b/GPy/inference/latent_function_inference/laplace.py @@ -216,6 +216,8 @@ class Laplace(object): """ if not log_concave: #print "Under 1e-10: {}".format(np.sum(W < 1e-6)) + # W[W<1e-6] = 1e-6 + # NOTE: when setting a parameter inside parameters_changed it will allways come to closed update circles!!! W.__setitem__(W < 1e-6, 1e-6, update=False) # FIXME-HACK: This is a hack since GPy can't handle negative variances which can occur # If the likelihood is non-log-concave. We wan't to say that there is a negative variance # To cause the posterior to become less certain than the prior and likelihood, From a4212e19049f36cda4a2edf4327e2eb07d03af0f Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Fri, 14 Feb 2014 15:23:41 +0000 Subject: [PATCH 32/40] final prior computation issues killed --- GPy/core/parameterization/parameter_core.py | 3 +-- GPy/core/parameterization/priors.py | 2 ++ 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 9741b4ef..275198b2 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -228,10 +228,9 @@ class Constrainable(Nameable, Indexable, Parameterizable): def log_prior(self): """evaluate the prior""" - import numpy as np if self.priors.size > 0: x = self._get_params() - return np.sum([p.lnpdf(x[ind]) for p, ind in self.priors.iteritems()]) + return reduce(lambda a,b: a+b, [p.lnpdf(x[ind]).sum() for p, ind in self.priors.iteritems()], 0) return 0. def _log_prior_gradients(self): diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py index 3155be64..906fe003 100644 --- a/GPy/core/parameterization/priors.py +++ b/GPy/core/parameterization/priors.py @@ -21,6 +21,8 @@ class Prior(object): from ...plotting.matplot_dep import priors_plots priors_plots.univariate_plot(self) + def __repr__(self, *args, **kwargs): + return self.__str__() class Gaussian(Prior): """ From 6defb2dbde3630fc1015488b2e2825c20d35acc3 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 17 Feb 2014 08:45:23 +0000 Subject: [PATCH 33/40] added index testing --- GPy/testing/index_operations_tests.py | 60 +++++++++++++++++++++++++++ 1 file changed, 60 insertions(+) create mode 100644 GPy/testing/index_operations_tests.py diff --git a/GPy/testing/index_operations_tests.py b/GPy/testing/index_operations_tests.py new file mode 100644 index 00000000..d5ef7007 --- /dev/null +++ b/GPy/testing/index_operations_tests.py @@ -0,0 +1,60 @@ +''' +Created on 12 Feb 2014 + +@author: maxz +''' +import unittest +import numpy as np +from GPy.core.parameterization.index_operations import ParameterIndexOperations,\ + ParameterIndexOperationsView + +one, two, three = 'one', 'two', 'three' + +class Test(unittest.TestCase): + + def setUp(self): + self.param_index = ParameterIndexOperations() + self.param_index.add(one, [3]) + self.param_index.add(two, [0,5]) + self.param_index.add(three, [2,4,7]) + + def test_remove(self): + self.param_index.remove(three, np.r_[3:10]) + self.assertListEqual(self.param_index[three].tolist(), [2]) + self.param_index.remove(one, [1]) + self.assertListEqual(self.param_index[one].tolist(), [3]) + + def test_index_view(self): + #======================================================================= + # 0 1 2 3 4 5 6 7 8 9 + # one + # two two + # three three three + # view: [0 1 2 3 4 5 ] + #======================================================================= + view = ParameterIndexOperationsView(self.param_index, 2, 6) + self.assertSetEqual(set(view.properties()), set([one, two, three])) + for v,p in zip(view.properties_for(np.r_[:6]), self.param_index.properties_for(np.r_[2:2+6])): + self.assertEqual(v, p) + self.assertSetEqual(set(view[two]), set([3])) + self.assertSetEqual(set(self.param_index[two]), set([0, 5])) + view.add(two, np.array([0])) + self.assertSetEqual(set(view[two]), set([0,3])) + self.assertSetEqual(set(self.param_index[two]), set([0, 2, 5])) + view.clear() + for v,p in zip(view.properties_for(np.r_[:6]), self.param_index.properties_for(np.r_[2:2+6])): + self.assertEqual(v, p) + self.assertEqual(v, []) + param_index = ParameterIndexOperations() + param_index.add(one, [3]) + param_index.add(two, [0,5]) + param_index.add(three, [2,4,7]) + view2 = ParameterIndexOperationsView(param_index, 2, 6) + view.update(view2) + for [i,v],[i2,v2] in zip(sorted(param_index.items()), sorted(self.param_index.items())): + self.assertEqual(i, i2) + self.assertTrue(np.all(v == v2)) + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'Test.test_index_view'] + unittest.main() \ No newline at end of file From 20f749ff0df806419e30cc68e891e3f2d4d22898 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 17 Feb 2014 08:45:34 +0000 Subject: [PATCH 34/40] caching changes --- 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 374f9600..8806e6d3 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -29,7 +29,7 @@ class Cacher(object): return self.cached_outputs[-1] def on_cache_changed(self, X): - print id(X) + #print id(X) i = self.cached_inputs.index(X) self.inputs_changed[i] = True From 825d3c2154b3f90818c8374ffdf998dd428bff9c Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 17 Feb 2014 09:03:44 +0000 Subject: [PATCH 35/40] beginning of bgplvm with missing data --- GPy/core/sparse_gp.py | 14 ++++++++------ GPy/models/bayesian_gplvm.py | 37 ++++++++++++++++++++++++++++++++++-- 2 files changed, 43 insertions(+), 8 deletions(-) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 09e84c12..dd0f6c2c 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -54,19 +54,21 @@ class SparseGP(GP): self.add_parameter(self.Z, index=0) self.parameters_changed() - - def parameters_changed(self): - 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) - - #The derivative of the bound wrt the inducing inputs Z ( unless they're all fixed) + def _update_gradients_Z(self, add=False): + #The derivative of the bound wrt the inducing inputs Z ( unless they're all fixed) if not self.Z.is_fixed: - self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) + if add: self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) + else: self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) if self.X_variance is None: self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) else: self.Z.gradient += self.kern.dpsi1_dZ(self.grad_dict['dL_dpsi1'], self.Z, self.X, self.X_variance) self.Z.gradient += self.kern.dpsi2_dZ(self.grad_dict['dL_dpsi2'], self.Z, self.X, self.X_variance) + def parameters_changed(self): + 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._update_gradients_Z(add=False) + def _raw_predict(self, Xnew, X_variance_new=None, which_parts='all', full_cov=False): """ Make a prediction for the latent function values diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index acb04350..914ca4ae 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -72,9 +72,10 @@ class BayesianGPLVM(SparseGP, GPLVM): return 0.5 * (var_mean + var_S) - 0.5 * self.input_dim * self.num_data def parameters_changed(self): - super(BayesianGPLVM, self).parameters_changed() - self._log_marginal_likelihood -= self.KL_divergence() + 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._update_gradients_Z(add=False) + self._log_marginal_likelihood -= self.KL_divergence() dL_dmu, dL_dS = self.dL_dmuS() # dL: @@ -161,6 +162,38 @@ class BayesianGPLVM(SparseGP, GPLVM): 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.q.mean.gradient = dL_dmu + self.q.variance.gradient = dL_dS + + # dKL: + self.q.mean.gradient -= self.X + self.q.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): """ objective function for fitting the latent variables for test points From 619278a1c804fb59a139852ef37a7da02715fe2d Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 17 Feb 2014 11:06:08 +0000 Subject: [PATCH 36/40] missing sys --- GPy/inference/optimization/optimization.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/GPy/inference/optimization/optimization.py b/GPy/inference/optimization/optimization.py index d9be46ce..45586a1d 100644 --- a/GPy/inference/optimization/optimization.py +++ b/GPy/inference/optimization/optimization.py @@ -59,8 +59,9 @@ class Optimizer(): """ See GPy.plotting.matplot_dep.inference_plots """ + import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." - from ..plotting.matplot_dep import inference_plots + from ...plotting.matplot_dep import inference_plots inference_plots.plot_optimizer(self) From 07737b7ba732a566d0c3ed6671f47f58eb89b0e4 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 17 Feb 2014 12:04:40 +0000 Subject: [PATCH 37/40] rbf andl inear fixes --- GPy/kern/parts/bias.py | 3 ++- GPy/kern/parts/linear.py | 2 +- GPy/kern/parts/rbf.py | 4 ++-- GPy/util/__init__.py | 1 + 4 files changed, 6 insertions(+), 4 deletions(-) diff --git a/GPy/kern/parts/bias.py b/GPy/kern/parts/bias.py index 303074a8..d2301bcd 100644 --- a/GPy/kern/parts/bias.py +++ b/GPy/kern/parts/bias.py @@ -14,7 +14,8 @@ class Bias(Kernpart): :type variance: float """ super(Bias, self).__init__(input_dim, name) - self.variance = Param("variance", variance) + from ...core.parameterization.transformations import Logexp + self.variance = Param("variance", variance, Logexp()) self.add_parameter(self.variance) def K(self,X,X2,target): diff --git a/GPy/kern/parts/linear.py b/GPy/kern/parts/linear.py index 1b805219..db7c9834 100644 --- a/GPy/kern/parts/linear.py +++ b/GPy/kern/parts/linear.py @@ -61,7 +61,7 @@ class Linear(Kernpart): def update_gradients_full(self, dL_dK, X): #self.variances.gradient[:] = 0 - self._param_grad_helper(dL_dK, X, self.variances.gradient) + self._param_grad_helper(dL_dK, X, None, self.variances.gradient) def update_gradients_sparse(self, dL_dKmm, dL_dKnm, dL_dKdiag, X, Z): tmp = dL_dKdiag[:, None] * X ** 2 diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index 4a8ebd4d..027aa382 100644 --- a/GPy/kern/parts/rbf.py +++ b/GPy/kern/parts/rbf.py @@ -264,7 +264,7 @@ class RBF(Kernpart): } """ num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim - X, dvardLdK = param_to_array(X, dvardLdK) + X, dvardLdK, var_len3 = param_to_array(X, dvardLdK, var_len3) weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options) else: code = """ @@ -281,7 +281,7 @@ class RBF(Kernpart): } """ num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim - X, X2, dvardLdK = param_to_array(X, X2, dvardLdK) + X, X2, dvardLdK, var_len3 = param_to_array(X, X2, dvardLdK, var_len3) weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options) return target diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py index 5a335027..c10fea4c 100644 --- a/GPy/util/__init__.py +++ b/GPy/util/__init__.py @@ -10,6 +10,7 @@ import datasets import mocap import decorators import classification +import subarray_and_sorting import caching try: From 3e4013ecf3ce4538887ea7aa10d49f4ab841df02 Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 17 Feb 2014 12:05:20 +0000 Subject: [PATCH 38/40] fixed Observable-weave clash in rbf --- GPy/kern/parts/rbf.py | 1 + 1 file changed, 1 insertion(+) diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index 4a8ebd4d..986d48aa 100644 --- a/GPy/kern/parts/rbf.py +++ b/GPy/kern/parts/rbf.py @@ -247,6 +247,7 @@ class RBF(Kernpart): target = np.zeros(self.input_dim) dvardLdK = self._K_dvar * dL_dK var_len3 = self.variance / np.power(self.lengthscale, 3) + var_len3 = param_to_array(var_len3) if X2 is None: # save computation for the symmetrical case dvardLdK = dvardLdK + dvardLdK.T From afcffcf039484bd5d6ac0be496658d2393ed2a4f Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 17 Feb 2014 12:06:19 +0000 Subject: [PATCH 39/40] bad git merge --- GPy/kern/parts/rbf.py | 1 - 1 file changed, 1 deletion(-) diff --git a/GPy/kern/parts/rbf.py b/GPy/kern/parts/rbf.py index e549a14a..027aa382 100644 --- a/GPy/kern/parts/rbf.py +++ b/GPy/kern/parts/rbf.py @@ -247,7 +247,6 @@ class RBF(Kernpart): target = np.zeros(self.input_dim) dvardLdK = self._K_dvar * dL_dK var_len3 = self.variance / np.power(self.lengthscale, 3) - var_len3 = param_to_array(var_len3) if X2 is None: # save computation for the symmetrical case dvardLdK = dvardLdK + dvardLdK.T From e96d40b4cccd5c176a5cc54f4ae5ae49f40dc48e Mon Sep 17 00:00:00 2001 From: James Hensman Date: Mon, 17 Feb 2014 12:10:41 +0000 Subject: [PATCH 40/40] removed sampling keyword (sampling is a silly thing to have as an option --- GPy/plotting/matplot_dep/models_plots.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index 177b9a95..4e6589b4 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -134,7 +134,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', m, _ = model._raw_predict(Xgrid, which_parts=which_parts) Y = model.likelihood.Y else: - m, _, _, _ = model.predict(Xgrid, which_parts=which_parts,sampling=False) + m, _, _, _ = model.predict(Xgrid, which_parts=which_parts) Y = model.likelihood.data for d in which_data_ycols: m_d = m[:,d].reshape(resolution, resolution).T