From 56d749ded8434d6e09d91cfba20f2291c320f2b2 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 4 Mar 2014 17:31:10 +0000 Subject: [PATCH 01/11] indentation... --- GPy/util/linalg.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/GPy/util/linalg.py b/GPy/util/linalg.py index f2b372be..4745c4aa 100644 --- a/GPy/util/linalg.py +++ b/GPy/util/linalg.py @@ -78,24 +78,24 @@ def force_F_ordered(A): # 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." + 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." From 0df263956fc0e30890e9a03465c93acecbc28881 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Tue, 4 Mar 2014 17:32:46 +0000 Subject: [PATCH 02/11] Dont call parameters_changed ever yourself anymore and parameters are now inplace once in memory --- GPy/core/gp.py | 2 - GPy/core/model.py | 19 +- GPy/core/parameterization/array_core.py | 52 +-- GPy/core/parameterization/index_operations.py | 3 +- GPy/core/parameterization/param.py | 54 +-- GPy/core/parameterization/parameter_core.py | 429 ++++++++++++++---- GPy/core/parameterization/parameterized.py | 72 ++- GPy/core/parameterization/transformations.py | 4 +- GPy/core/sparse_gp.py | 15 +- GPy/examples/__init__.py | 1 + .../latent_function_inference/var_dtc.py | 17 +- GPy/kern/_src/coregionalize.py | 2 - GPy/kern/_src/periodic.py | 1 - GPy/kern/_src/sympykern.py | 3 - GPy/models/bayesian_gplvm.py | 1 - GPy/models/gp_regression.py | 1 - GPy/models/mrd.py | 159 +++++-- GPy/plotting/matplot_dep/kernel_plots.py | 2 +- GPy/testing/likelihood_tests.py | 10 +- GPy/testing/observable_tests.py | 21 +- GPy/testing/parameterized_tests.py | 17 +- 21 files changed, 601 insertions(+), 284 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 7cc39dde..1add8268 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -65,8 +65,6 @@ class GP(Model): self.add_parameter(self.kern) self.add_parameter(self.likelihood) - 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, Y_metadata=self.Y_metadata) diff --git a/GPy/core/model.py b/GPy/core/model.py index bf8915c6..38b5eb29 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -15,6 +15,7 @@ import itertools class Model(Parameterized): _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) self.optimization_runs = [] @@ -25,14 +26,8 @@ class Model(Parameterized): raise NotImplementedError, "this needs to be implemented to use the model class" def _log_likelihood_gradients(self): - g = np.zeros(self.size) - try: - [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) - return g - raise NotImplementedError, "this needs to be implemented to use the model class" - + return self.gradient + def _getstate(self): """ Get the current state of the class. @@ -208,8 +203,8 @@ class Model(Parameterized): try: self._set_params_transformed(x) obj_f = -float(self.log_likelihood()) - self.log_prior() - self._fail_count = 0 obj_grads = -self._transform_gradients(self._log_likelihood_gradients() + self._log_prior_gradients()) + self._fail_count = 0 except (LinAlgError, ZeroDivisionError, ValueError) as e: if self._fail_count >= self._allowed_failures: raise e @@ -275,9 +270,8 @@ class Model(Parameterized): The gradient is considered correct if the ratio of the analytical and numerical gradients is within of unity. """ - - x = self._get_params_transformed() - + x = self._get_params_transformed().copy() + if not verbose: # make sure only to test the selected parameters if target_param is None: @@ -368,6 +362,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 diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index cf353ead..27801e23 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -49,7 +49,7 @@ class ObservableArray(np.ndarray, Observable): def __setitem__(self, s, val): if self._s_not_empty(s): super(ObservableArray, self).__setitem__(s, val) - self._notify_observers(self[s]) + self.notify_observers(self[s]) def __getslice__(self, start, stop): return self.__getitem__(slice(start, stop)) @@ -65,149 +65,149 @@ class ObservableArray(np.ndarray, Observable): def __ilshift__(self, *args, **kwargs): r = np.ndarray.__ilshift__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __irshift__(self, *args, **kwargs): r = np.ndarray.__irshift__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __ixor__(self, *args, **kwargs): r = np.ndarray.__ixor__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __ipow__(self, *args, **kwargs): r = np.ndarray.__ipow__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __ifloordiv__(self, *args, **kwargs): r = np.ndarray.__ifloordiv__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __isub__(self, *args, **kwargs): r = np.ndarray.__isub__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __ior__(self, *args, **kwargs): r = np.ndarray.__ior__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __itruediv__(self, *args, **kwargs): r = np.ndarray.__itruediv__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __idiv__(self, *args, **kwargs): r = np.ndarray.__idiv__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __iand__(self, *args, **kwargs): r = np.ndarray.__iand__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __imod__(self, *args, **kwargs): r = np.ndarray.__imod__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __iadd__(self, *args, **kwargs): r = np.ndarray.__iadd__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r def __imul__(self, *args, **kwargs): r = np.ndarray.__imul__(self, *args, **kwargs) - self._notify_observers() + self.notify_observers() return r # def __rrshift__(self, *args, **kwargs): # r = np.ndarray.__rrshift__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __ror__(self, *args, **kwargs): # r = np.ndarray.__ror__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rxor__(self, *args, **kwargs): # r = np.ndarray.__rxor__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rdivmod__(self, *args, **kwargs): # r = np.ndarray.__rdivmod__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __radd__(self, *args, **kwargs): # r = np.ndarray.__radd__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rdiv__(self, *args, **kwargs): # r = np.ndarray.__rdiv__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rtruediv__(self, *args, **kwargs): # r = np.ndarray.__rtruediv__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rshift__(self, *args, **kwargs): # r = np.ndarray.__rshift__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rmul__(self, *args, **kwargs): # r = np.ndarray.__rmul__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rpow__(self, *args, **kwargs): # r = np.ndarray.__rpow__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rsub__(self, *args, **kwargs): # r = np.ndarray.__rsub__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r # def __rfloordiv__(self, *args, **kwargs): # r = np.ndarray.__rfloordiv__(self, *args, **kwargs) -# self._notify_observers() +# self.notify_observers() # return r diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index a9f3768e..f8f6ab5b 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -62,6 +62,7 @@ class ParameterIndexOperations(object): def clear(self): self._properties.clear() + @property def size(self): return reduce(lambda a,b: a+b.size, self.iterindices(), 0) @@ -165,7 +166,7 @@ class ParameterIndexOperationsView(object): for i, ind in self.items(): self._param_index_ops.remove(i, ind+self._offset) - + @property def size(self): return reduce(lambda a,b: a+b.size, self.iterindices(), 0) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 8f4dbaf3..2917cad7 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -54,7 +54,7 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): obj._tied_to_me_ = SetDict() obj._tied_to_ = [] obj._original_ = True - obj._gradient_ = None + obj._gradient_array_ = numpy.zeros(obj.shape, dtype=numpy.float64) return obj def __init__(self, name, input_array, default_constraint=None, *a, **kw): @@ -77,7 +77,7 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): # see InfoArray.__array_finalize__ for comments if obj is None: return super(Param, self).__array_finalize__(obj) - self._direct_parent_ = getattr(obj, '_direct_parent_', None) + self._parent_ = getattr(obj, '_parent_', None) self._parent_index_ = getattr(obj, '_parent_index_', None) self._default_constraint_ = getattr(obj, '_default_constraint_', None) self._current_slice_ = getattr(obj, '_current_slice_', None) @@ -89,16 +89,18 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): self._updated_ = getattr(obj, '_updated_', None) self._original_ = getattr(obj, '_original_', None) self._name = getattr(obj, 'name', None) - self._gradient_ = getattr(obj, '_gradient_', None) + self._gradient_array_ = getattr(obj, '_gradient_array_', None) self.constraints = getattr(obj, 'constraints', None) self.priors = getattr(obj, 'priors', None) - + @property + def _param_array_(self): + return self + @property def gradient(self): - if self._gradient_ is None: - self._gradient_ = numpy.zeros(self._realshape_) - return self._gradient_[self._current_slice_] + return self._gradient_array_[self._current_slice_] + @gradient.setter def gradient(self, val): self.gradient[:] = val @@ -110,7 +112,7 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): func, args, state = super(Param, self).__reduce__() return func, args, (state, (self.name, - self._direct_parent_, + self._parent_, self._parent_index_, self._default_constraint_, self._current_slice_, @@ -135,7 +137,7 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): self._current_slice_ = state.pop() self._default_constraint_ = state.pop() self._parent_index_ = state.pop() - self._direct_parent_ = state.pop() + self._parent_ = state.pop() self.name = state.pop() def copy(self, *args): @@ -148,20 +150,20 @@ class Param(OptimizationHandlable, ObservableArray, Gradcheckable): #=========================================================================== # get/set parameters #=========================================================================== - def _set_params(self, param, trigger_parent=True): - self.flat = param - if trigger_parent: min_priority = None - else: min_priority = -numpy.inf - self._notify_observers(None, min_priority) - - def _get_params(self): - return self.flat - - def _collect_gradient(self, target): - target += self.gradient.flat - - def _set_gradient(self, g): - self.gradient = g.reshape(self._realshape_) +# def _set_params(self, param, trigger_parent=True): +# self.flat = param +# if trigger_parent: min_priority = None +# else: min_priority = -numpy.inf +# self.notify_observers(None, min_priority) +# +# def _get_params(self): +# return self.flat +# +# def _collect_gradient(self, target): +# target += self.gradient.flat +# +# def _set_gradient(self, g): +# self.gradient = g.reshape(self._realshape_) #=========================================================================== # Array operations -> done @@ -362,7 +364,7 @@ class ParamConcatenation(object): parents = dict() for p in self.params: if p.has_parent(): - parent = p._direct_parent_ + parent = p._parent_ level = 0 while parent is not None: if parent in parents: @@ -370,7 +372,7 @@ class ParamConcatenation(object): else: parents[parent] = level level += 1 - parent = parent._direct_parent_ + parent = parent._parent_ import operator self.parents = map(lambda x: x[0], sorted(parents.iteritems(), key=operator.itemgetter(1))) #=========================================================================== @@ -397,7 +399,7 @@ class ParamConcatenation(object): #=========================================================================== def update_all_params(self): for par in self.parents: - par._notify_observers(-numpy.inf) + par.notify_observers(-numpy.inf) def constrain(self, constraint, warning=True): [param.constrain(constraint, trigger_parent=False) for param in self.params] diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 109df6e9..c9372b58 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -1,26 +1,51 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) +""" +Core module for parameterization. +This module implements all parameterization techniques, split up in modular bits. + +HierarchyError: +raised when an error with the hierarchy occurs (circles etc.) + +Observable: +Observable Pattern for patameterization + + +""" from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED import numpy as np +import itertools __updated__ = '2013-12-16' class HierarchyError(Exception): """ - Gets thrown when something is wrong with the parameter hierarchy + Gets thrown when something is wrong with the parameter hierarchy. """ def adjust_name_for_printing(name): + """ + Make sure a name can be printed, alongside used as a variable name. + """ if name is not None: - return name.replace(" ", "_").replace(".", "_").replace("-", "").replace("+", "").replace("!", "").replace("*", "").replace("/", "") + return name.replace(" ", "_").replace(".", "_").replace("-", "_m_").replace("+", "_p_").replace("!", "_I_").replace("**", "_xx_").replace("*", "_x_").replace("/", "_l_").replace("@",'_at_') return '' class Observable(object): + """ + Observable pattern for parameterization. + + This Object allows for observers to register with self and a (bound!) function + as an observer. Every time the observable changes, it sends a notification with + self as only argument to all its observers. + """ _updated = True def __init__(self, *args, **kwargs): self._observer_callables_ = [] - + def __del__(self, *args, **kwargs): + del self._observer_callables_ + def add_observer(self, observer, callble, priority=0): self._insert_sorted(priority, observer, callble) @@ -35,8 +60,8 @@ class Observable(object): to_remove.append((p, obs, clble)) for r in to_remove: self._observer_callables_.remove(r) - - def _notify_observers(self, which=None, min_priority=None): + + def notify_observers(self, which=None, min_priority=None): """ Notifies all observers. Which is the element, which kicked off this notification loop. @@ -67,6 +92,41 @@ class Observable(object): self._observer_callables_.insert(ins, (p, o, c)) class Pickleable(object): + """ + Make an object pickleable (See python doc 'pickling'). + + This class allows for pickling support by Memento pattern. + _getstate returns a memento of the class, which gets pickled. + _setstate() (re-)sets the state of the class to the memento + """ + #=========================================================================== + # Pickling operations + #=========================================================================== + def pickle(self, f, protocol=-1): + """ + :param f: either filename or open file object to write to. + if it is an open buffer, you have to make sure to close + it properly. + :param protocol: pickling protocol to use, python-pickle for details. + """ + import cPickle + if isinstance(f, str): + with open(f, 'w') as f: + cPickle.dump(self, f, protocol) + else: + cPickle.dump(self, f, protocol) + def __getstate__(self): + if self._has_get_set_state(): + return self._getstate() + return self.__dict__ + def __setstate__(self, state): + if self._has_get_set_state(): + self._setstate(state) + # TODO: maybe parameters_changed() here? + return + self.__dict__ = state + def _has_get_set_state(self): + return '_getstate' in vars(self.__class__) and '_setstate' in vars(self.__class__) def _getstate(self): """ Returns the state of this class in a memento pattern. @@ -93,70 +153,145 @@ class Pickleable(object): #=============================================================================== class Parentable(object): - _direct_parent_ = None + """ + Enable an Object to have a parent. + + Additionally this adds the parent_index, which is the index for the parent + to look for in its parameter list. + """ + _parent_ = None _parent_index_ = None 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) + """ + Return whether this parentable object currently has a parent. + """ + return self._parent_ is not None def _parent_changed(self): + """ + Gets called, when the parent changed, so we can adjust our + inner attributes according to the new parent. + """ raise NotImplementedError, "shouldnt happen, Parentable objects need to be able to change their parent" + def _disconnect_parent(self, *args, **kw): + """ + Disconnect this object from its parent + """ + raise NotImplementedError, "Abstaract superclass" + @property def _highest_parent_(self): - if self._direct_parent_ is None: + """ + Gets the highest parent by traversing up to the root node of the hierarchy. + """ + if self._parent_ is None: return self - return self._direct_parent_._highest_parent_ + return self._parent_._highest_parent_ - def _notify_parameters_changed(self): - raise NotImplementedError, "shouldnt happen, abstract superclass" - + def _notify_parent_change(self): + """ + Dont do anything if in leaf node + """ + pass class Nameable(Parentable): + """ + Make an object nameable inside the hierarchy. + """ def __init__(self, name, *a, **kw): super(Nameable, self).__init__(*a, **kw) self._name = name or self.__class__.__name__ @property def name(self): + """ + The name of this object + """ return self._name @name.setter def name(self, name): + """ + Set the name of this object. + Tell the parent if the name has changed. + """ from_name = self.name assert isinstance(name, str) self._name = name if self.has_parent(): - self._direct_parent_._name_changed(self, from_name) + self._parent_._name_changed(self, from_name) def hierarchy_name(self, adjust_for_printing=True): + """ + return the name for this object with the parents names attached by dots. + + :param bool adjust_for_printing: whether to call :func:`~adjust_for_printing()` + on the names, recursively + """ if adjust_for_printing: adjust = lambda x: adjust_name_for_printing(x) else: adjust = lambda x: x if self.has_parent(): - return self._direct_parent_.hierarchy_name() + "." + adjust(self.name) + return self._parent_.hierarchy_name() + "." + adjust(self.name) return adjust(self.name) class Gradcheckable(Parentable): + """ + Adds the functionality for an object to be gradcheckable. + It is just a thin wrapper of a call to the highest parent for now. + TODO: Can be done better, by only changing parameters of the current parameter handle, + such that object hierarchy only has to change for those. + """ def __init__(self, *a, **kw): super(Gradcheckable, self).__init__(*a, **kw) def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): + """ + Check the gradient of this parameter with respect to the highest parent's + objective function. + This is a three point estimate of the gradient, wiggling at the parameters + with a stepsize step. + The check passes if either the ratio or the difference between numerical and + analytical gradient is smaller then tolerance. + + :param bool verbose: whether each parameter shall be checked individually. + :param float step: the stepsize for the numerical three point gradient estimate. + :param flaot tolerance: the tolerance for the gradient ratio or difference. + """ if self.has_parent(): return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance) def _checkgrad(self, param): + """ + Perform the checkgrad on the model. + TODO: this can be done more efficiently, when doing it inside here + """ raise NotImplementedError, "Need log likelihood to check gradient against" class Indexable(object): + """ + Enable enraveled indexes and offsets for this object. + The raveled index of an object is the index for its parameters in a flattened int array. + """ def _raveled_index(self): + """ + Flattened array of ints, specifying the index of this object. + This has to account for shaped parameters! + """ raise NotImplementedError, "Need to be able to get the raveled Index" def _internal_offset(self): + """ + The offset for this parameter inside its parent. + This has to account for shaped parameters! + """ return 0 def _offset_for(self, param): + """ + Return the offset of the param inside this parameterized object. + This does not need to account for shaped parameters, as it + basically just sums up the parameter sizes which come before param. + """ raise NotImplementedError, "shouldnt happen, offset required from non parameterization object?" def _raveled_index_for(self, param): @@ -169,6 +304,15 @@ class Indexable(object): class Constrainable(Nameable, Indexable): + """ + Make an object constrainable with Priors and Transformations. + TODO: Mappings!! + Adding a constraint to a Parameter means to tell the highest parent that + the constraint was added and making sure that all parameters covered + by this object are indeed conforming to the constraint. + + :func:`constrain()` and :func:`unconstrain()` are main methods here + """ def __init__(self, name, default_constraint=None, *a, **kw): super(Constrainable, self).__init__(name=name, *a, **kw) self._default_constraint_ = default_constraint @@ -178,12 +322,16 @@ class Constrainable(Nameable, Indexable): if self._default_constraint_ is not None: self.constrain(self._default_constraint_) - def _disconnect_parent(self, constr=None): + def _disconnect_parent(self, constr=None, *args, **kw): + """ + From Parentable: + disconnect the parent and set the new constraints to constr + """ if constr is None: constr = self.constraints.copy() self.constraints.clear() self.constraints = constr - self._direct_parent_ = None + self._parent_ = None self._parent_index_ = None self._connect_fixes() self._notify_parent_change() @@ -193,7 +341,7 @@ class Constrainable(Nameable, Indexable): #=========================================================================== def constrain_fixed(self, value=None, warning=True, trigger_parent=True): """ - Constrain this paramter to be fixed to the current value it carries. + Constrain this parameter to be fixed to the current value it carries. :param warning: print a warning for overwriting constraints. """ @@ -237,11 +385,20 @@ class Constrainable(Nameable, Indexable): #=========================================================================== # Prior Operations #=========================================================================== - def set_prior(self, prior, warning=True, trigger_parent=True): + def set_prior(self, prior, warning=True): + """ + Set the prior for this object to prior. + :param :class:`~GPy.priors.Prior` prior: a prior to set for this parameter + :param bool warning: whether to warn if another prior was set for this parameter + """ repriorized = self.unset_priors() self._add_to_index_operations(self.priors, repriorized, prior, warning) def unset_priors(self, *priors): + """ + Un-set all priors given from this parameter handle. + + """ return self._remove_from_index_operations(self.priors, priors) def log_prior(self): @@ -274,7 +431,7 @@ class Constrainable(Nameable, Indexable): :py:class:`GPy.core.transformations.Transformation`. """ if isinstance(transform, Transformation): - self._set_params(transform.initialize(self._get_params()), trigger_parent=trigger_parent) + self._param_array_[:] = transform.initialize(self._param_array_) reconstrained = self.unconstrain() self._add_to_index_operations(self.constraints, reconstrained, transform, warning) @@ -333,6 +490,10 @@ class Constrainable(Nameable, Indexable): self.unconstrain(Logistic(lower, upper)) def _parent_changed(self, parent): + """ + From Parentable: + Called when the parent changed + """ 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) @@ -340,14 +501,25 @@ class Constrainable(Nameable, Indexable): for p in self._parameters_: p._parent_changed(parent) - def _add_to_index_operations(self, which, reconstrained, transform, warning): + def _add_to_index_operations(self, which, reconstrained, what, warning): + """ + Helper preventing copy code. + This addes the given what (transformation, prior etc) to parameter index operations which. + revonstrained are reconstrained indices. + warn when reconstraining parameters if warning is True. + TODO: find out which parameters have changed specifically + """ if warning and reconstrained.size > 0: # TODO: figure out which parameters have changed and only print those print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) - which.add(transform, self._raveled_index()) + which.add(what, self._raveled_index()) - def _remove_from_index_operations(self, which, transforms): - if len(transforms) == 0: + def _remove_from_index_operations(self, which, what): + """ + Helper preventing copy code. + Remove given what (transform prior etc) from which param index ops. + """ + if len(what) == 0: transforms = which.properties() removed = np.empty((0,), dtype=int) for t in transforms: @@ -359,36 +531,65 @@ class Constrainable(Nameable, Indexable): return removed class OptimizationHandlable(Constrainable, Observable): + """ + This enables optimization handles on an Object as done in GPy 0.4. + + transformed: make sure the transformations and constraints etc are handled + """ + def transform(self): + [np.put(self._param_array_, ind, c.finv(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + + def untransform(self): + [np.put(self._param_array_, ind, c.f(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + def _get_params_transformed(self): # transformed parameters (apply transformation rules) - p = self._get_params() + p = self._param_array_.copy() [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] if self._has_fixes(): return p[self._fixes_] return p - + def _set_params_transformed(self, p): - # inverse apply transformations for parameters and set the resulting parameters - self._set_params(self._untransform_params(p)) + if self._has_fixes(): self._param_array_[self._fixes_] = p.copy() + else: self._param_array_[:] = p.copy() + self.untransform() + self._trigger_params_changed() + + def _trigger_params_changed(self, trigger_parent=True): + [p._trigger_params_changed(trigger_parent=False) for p in self._parameters_] + if trigger_parent: min_priority = None + else: min_priority = -np.inf + self.notify_observers(None, min_priority) def _size_transformed(self): return self.size - self.constraints[__fixed__].size - - def _untransform_params(self, p): - p = p.copy() - if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp - [np.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - return p - - def _get_params(self): +# +# def _untransform_params(self, p): +# # inverse apply transformations for parameters +# #p = p.copy() +# if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp +# [np.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] +# return p +# +# def _get_params(self): +# """ +# get all parameters +# """ +# return self._param_array_ +# p = np.empty(self.size, dtype=np.float64) +# if self.size == 0: +# return p +# [np.put(p, ind, par._get_params()) for ind, par in itertools.izip(self._param)] +# return p + +# def _set_params(self, params, trigger_parent=True): +# self._param_array_.flat = params +# if trigger_parent: min_priority = None +# else: min_priority = -np.inf +# self.notify_observers(None, min_priority) # 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, trigger_parent=True): - # don't overwrite this anymore! - raise NotImplementedError, "This needs to be implemented in Param and Parametrizable" + #raise NotImplementedError, "Abstract superclass: This needs to be implemented in Param and Parameterizable" #=========================================================================== # Optimization handles: @@ -396,6 +597,7 @@ class OptimizationHandlable(Constrainable, Observable): def _get_param_names(self): n = np.array([p.hierarchy_name() + '[' + str(i) + ']' for p in self.flattened_parameters for i in p._indices()]) return n + def _get_param_names_transformed(self): n = self._get_param_names() if self._has_fixes(): @@ -405,19 +607,16 @@ class OptimizationHandlable(Constrainable, Observable): #=========================================================================== # Randomizeable #=========================================================================== - def randomize(self): + def randomize(self, rand_gen=np.random.normal, loc=0, scale=1, *args, **kwargs): """ Randomize the model. - Make this draw from the prior if one exists, else draw from N(0,1) + Make this draw from the prior if one exists, else draw from given random generator """ # first take care of all parameters (from N(0,1)) - # x = self._get_params_transformed() - x = np.random.randn(self._size_transformed()) - x = self._untransform_params(x) + x = rand_gen(loc=loc, scale=scale, size=self._size_transformed(), *args, **kwargs) # now draw from prior where possible [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None] - self._set_params(x) - # self._set_params_transformed(self._get_params_transformed()) # makes sure all of the tied parameters get the same init (since there's only one prior object...) + self._set_params_transformed(x) # makes sure all of the tied parameters get the same init (since there's only one prior object...) class Parameterizable(OptimizationHandlable): def __init__(self, *args, **kwargs): @@ -427,6 +626,13 @@ class Parameterizable(OptimizationHandlable): self._added_names_ = set() def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): + """ + Get the names of all parameters of this model. + + :param bool add_self: whether to add the own name in front of names + :param bool adjust_for_printing: whether to call `adjust_name_for_printing` on names + :param bool recursive: whether to traverse through hierarchy and append leaf node names + """ if adjust_for_printing: adjust = lambda x: adjust_name_for_printing(x) else: adjust = lambda x: x if recursive: names = [xi for x in self._parameters_ for xi in x.parameter_names(add_self=True, adjust_for_printing=adjust_for_printing)] @@ -438,8 +644,11 @@ class Parameterizable(OptimizationHandlable): def num_params(self): return len(self._parameters_) - def _add_parameter_name(self, param): + def _add_parameter_name(self, param, ignore_added_names=False): pname = adjust_name_for_printing(param.name) + if ignore_added_names: + self.__dict__[pname] = param + return # and makes sure to not delete programmatically added parameters if pname in self.__dict__: if not (param is self.__dict__[pname]): @@ -461,28 +670,42 @@ class Parameterizable(OptimizationHandlable): def _name_changed(self, param, old_name): self._remove_parameter_name(None, old_name) self._add_parameter_name(param) - - def _collect_gradient(self, target): - import itertools - [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + + #========================================================================= + # Gradient handling + #========================================================================= + @property + def gradient(self): + return self._gradient_array_ + + @gradient.setter + def gradient(self, val): + self._gradient_array_[:] = val + #=========================================================================== + # def _collect_gradient(self, target): + # [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + #=========================================================================== - def _set_params(self, params, trigger_parent=True): - import itertools - [p._set_params(params[s], trigger_parent=False) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - if trigger_parent: min_priority = None - else: min_priority = -np.inf - self._notify_observers(None, min_priority) + #=========================================================================== + # def _set_params(self, params, trigger_parent=True): + # [p._set_params(params[s], trigger_parent=False) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + # if trigger_parent: min_priority = None + # else: min_priority = -np.inf + # self.notify_observers(None, min_priority) + #=========================================================================== - def _set_gradient(self, g): - import itertools - [p._set_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + #=========================================================================== + # def _set_gradient(self, g): + # [p._set_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + #=========================================================================== - def add_parameter(self, param, index=None): + def add_parameter(self, param, index=None, _ignore_added_names=False): """ :param parameters: the parameters to add :type parameters: list of or one :py:class:`GPy.core.param.Param` :param [index]: index of where to put parameters - + + :param bool _ignore_added_names: whether the name of the parameter overrides a possibly existing field Add all parameters to this param class, you can insert parameters at any given index using the :func:`list.insert` syntax @@ -494,12 +717,12 @@ class Parameterizable(OptimizationHandlable): self.add_parameter(param, index) elif param not in self._parameters_: if param.has_parent(): - parent = param._direct_parent_ + parent = param._parent_ while parent is not None: if parent is self: - raise HierarchyError, "You cannot add a parameter twice into the hirarchy" - parent = parent._direct_parent_ - param._direct_parent_.remove_parameter(param) + raise HierarchyError, "You cannot add a parameter twice into the hierarchy" + parent = parent._parent_ + param._parent_.remove_parameter(param) # make sure the size is set if index is None: self.constraints.update(param.constraints, self.size) @@ -517,7 +740,7 @@ class Parameterizable(OptimizationHandlable): self.size += param.size - self._connect_parameters() + self._connect_parameters(ignore_added_names=_ignore_added_names) self._notify_parent_change() self._connect_fixes() else: @@ -551,14 +774,14 @@ class Parameterizable(OptimizationHandlable): self._connect_parameters() self._notify_parent_change() - parent = self._direct_parent_ + parent = self._parent_ while parent is not None: parent._connect_fixes() parent._connect_parameters() parent._notify_parent_change() - parent = parent._direct_parent_ + parent = parent._parent_ - def _connect_parameters(self): + def _connect_parameters(self, ignore_added_names=False): # connect parameterlist to this parameterized object # This just sets up the right connection for the params objects # to be used as parameters @@ -567,14 +790,35 @@ class Parameterizable(OptimizationHandlable): if not hasattr(self, "_parameters_") or len(self._parameters_) < 1: # no parameters for this class return - sizes = [0] + old_size = 0 + self._param_array_ = np.empty(self.size, dtype=np.float64) + self._gradient_array_ = np.empty(self.size, dtype=np.float64) + self._param_slices_ = [] + for i, p in enumerate(self._parameters_): - p._direct_parent_ = self + p._parent_ = self p._parent_index_ = i - sizes.append(p.size + sizes[-1]) - self._param_slices_.append(slice(sizes[-2], sizes[-1])) - self._add_parameter_name(p) + + pslice = slice(old_size, old_size+p.size) + pi_old_size = old_size + for pi in p.flattened_parameters: + pislice = slice(pi_old_size, pi_old_size+pi.size) + + self._param_array_[pislice] = pi._param_array_.flat + self._gradient_array_[pislice] = pi._gradient_array_.flat + + pi._param_array_.data = self._param_array_[pislice].data + pi._gradient_array_.data = self._gradient_array_[pislice].data + + pi_old_size += pi.size + + p._param_array_.data = self._param_array_[pslice].data + p._gradient_array_.data = self._gradient_array_[pslice].data + + self._param_slices_.append(pslice) + self._add_parameter_name(p, ignore_added_names=ignore_added_names) + old_size += p.size #=========================================================================== # notification system @@ -582,7 +826,7 @@ class Parameterizable(OptimizationHandlable): def _parameters_changed_notification(self, which): self.parameters_changed() def _pass_through_notify_observers(self, which): - self._notify_observers(which) + self.notify_observers(which) #=========================================================================== # TODO: not working yet @@ -595,7 +839,7 @@ class Parameterizable(OptimizationHandlable): dc = dict() for k, v in self.__dict__.iteritems(): - if k not in ['_direct_parent_', '_parameters_', '_parent_index_', '_observer_callables_'] + self.parameter_names(): + if k not in ['_parent_', '_parameters_', '_parent_index_', '_observer_callables_'] + self.parameter_names(recursive=False): if isinstance(v, (Constrainable, ParameterIndexOperations, ParameterIndexOperationsView)): dc[k] = v.copy() else: @@ -603,7 +847,7 @@ class Parameterizable(OptimizationHandlable): if k == '_parameters_': params = [p.copy() for p in v] - dc['_direct_parent_'] = None + dc['_parent_'] = None dc['_parent_index_'] = None dc['_observer_callables_'] = [] dc['_parameters_'] = ArrayList() @@ -615,11 +859,20 @@ class Parameterizable(OptimizationHandlable): s.__dict__ = dc for p in params: - import ipdb;ipdb.set_trace() - s.add_parameter(p) + s.add_parameter(p, _ignore_added_names=True) return s - + + #=========================================================================== + # From being parentable, we have to define the parent_change notification + #=========================================================================== + def _notify_parent_change(self): + """ + Notify all parameters that the parent has changed + """ + for p in self._parameters_: + p._parent_changed(self) + def parameters_changed(self): """ This method gets called when parameters have changed. diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 3ef99a35..50dd7b3b 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -11,6 +11,12 @@ from parameter_core import Pickleable, Parameterizable, adjust_name_for_printing from transformations import __fixed__ from lists_and_dicts import ArrayList +class ParametersChangedMeta(type): + def __call__(self, *args, **kw): + instance = super(ParametersChangedMeta, self).__call__(*args, **kw) + instance.parameters_changed() + return instance + class Parameterized(Parameterizable, Pickleable, Gradcheckable): """ Parameterized class @@ -53,6 +59,12 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): If you want to operate on all parameters use m[''] to wildcard select all paramters and concatenate them. Printing m[''] will result in printing of all parameters in detail. """ + #=========================================================================== + # Metaclass for parameters changed after init. + # This makes sure, that parameters changed will always be called after __init__ + # **Never** call parameters_changed() yourself + __metaclass__ = ParametersChangedMeta + #=========================================================================== def __init__(self, name=None, *a, **kw): super(Parameterized, self).__init__(name=name, parent=None, parent_index=None, *a, **kw) self._in_init_ = True @@ -88,39 +100,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): return G return node - #=========================================================================== - # Pickling operations - #=========================================================================== - def pickle(self, f, protocol=-1): - """ - :param f: either filename or open file object to write to. - if it is an open buffer, you have to make sure to close - it properly. - :param protocol: pickling protocol to use, python-pickle for details. - """ - if isinstance(f, str): - with open(f, 'w') as f: - cPickle.dump(self, f, protocol) - else: - cPickle.dump(self, f, protocol) - - def copy(self): - c = super(Parameterized, self).copy() - c.add_observer(c, c._parameters_changed_notification, -100) - return c - def __getstate__(self): - if self._has_get_set_state(): - return self._getstate() - 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 - return - self.__dict__ = state - def _has_get_set_state(self): - return '_getstate' in vars(self.__class__) and '_setstate' in vars(self.__class__) def _getstate(self): """ Get the current state of the class, @@ -149,25 +129,33 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self._connect_parameters() self.parameters_changed() #=========================================================================== + # Override copy to handle programmatically added observers + #=========================================================================== + def copy(self): + c = super(Pickleable, self).copy() + c.add_observer(c, c._parameters_changed_notification, -100) + return c + + #=========================================================================== # Gradient control #=========================================================================== def _transform_gradients(self, g): 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__] - for p in self.flattened_parameters: - for t, i in p._tied_to_me_.iteritems(): - g[self._offset_for(p) + numpy.array(list(i))] += g[self._raveled_index_for(t)] + [numpy.put(g, i, g[i] * c.gradfactor(self._param_array_[i])) for c, i in self.constraints.iteritems() if c != __fixed__] if self._has_fixes(): return g[self._fixes_] return g + + #=========================================================================== + # Indexable + #=========================================================================== def _offset_for(self, param): # get the offset in the parameterized index array for param if param.has_parent(): - if param._direct_parent_._get_original(param) in self._parameters_: - 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) + if param._parent_._get_original(param) in self._parameters_: + return self._param_slices_[param._parent_._get_original(param)._parent_index_].start + return self._offset_for(param._parent_) + param._parent_._offset_for(param) return 0 def _raveled_index_for(self, param): @@ -229,8 +217,8 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): return ParamConcatenation(paramlist) def __setitem__(self, name, value, paramlist=None): - if isinstance(name, slice): - self[''][name] = value + if isinstance(name, (slice, tuple, np.ndarray)): + self._param_array_[name] = value else: try: param = self.__getitem__(name, paramlist) except AttributeError as a: raise a diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index 60fcc469..5cda8d46 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -10,7 +10,7 @@ import sys #_lim_val = -np.log(sys.float_info.epsilon) _exp_lim_val = np.finfo(np.float64).max -_lim_val = np.log(_exp_lim_val)# +_lim_val = np.log(_exp_lim_val) #=============================================================================== # Fixing constants @@ -57,7 +57,7 @@ class Logexp(Transformation): return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -_lim_val, _lim_val)))) #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x))) def finv(self, f): - return np.where(f>_lim_val, f, np.log(np.exp(f) - 1.)) + return np.where(f>_lim_val, f, np.log(np.exp(f+1e-20) - 1.)) def gradfactor(self, f): return np.where(f>_lim_val, 1., 1 - np.exp(-f)) def initialize(self, f): diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 4230d8ba..f4f34a5e 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -48,7 +48,6 @@ class SparseGP(GP): GP.__init__(self, X, Y, kernel, likelihood, inference_method=inference_method, name=name) self.add_parameter(self.Z, index=0) - self.parameters_changed() def has_uncertain_inputs(self): return isinstance(self.X, VariationalPosterior) @@ -60,11 +59,9 @@ class SparseGP(GP): #gradients wrt kernel dL_dKmm = self.grad_dict.pop('dL_dKmm') self.kern.update_gradients_full(dL_dKmm, self.Z, None) - target = np.zeros(self.kern.size) - self.kern._collect_gradient(target) + target = self.kern.gradient.copy() self.kern.update_gradients_expectations(variational_posterior=self.X, Z=self.Z, **self.grad_dict) - self.kern._collect_gradient(target) - self.kern._set_gradient(target) + self.kern.gradient += target #gradients wrt Z self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) @@ -72,14 +69,12 @@ class SparseGP(GP): self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) else: #gradients wrt kernel - target = np.zeros(self.kern.size) self.kern.update_gradients_diag(self.grad_dict['dL_dKdiag'], self.X) - self.kern._collect_gradient(target) + target = self.kern.gradient.copy() self.kern.update_gradients_full(self.grad_dict['dL_dKnm'], self.X, self.Z) - self.kern._collect_gradient(target) + target += self.kern.gradient self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None) - self.kern._collect_gradient(target) - self.kern._set_gradient(target) + self.kern.gradient += target #gradients wrt Z self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) diff --git a/GPy/examples/__init__.py b/GPy/examples/__init__.py index 2f74858a..87548553 100644 --- a/GPy/examples/__init__.py +++ b/GPy/examples/__init__.py @@ -6,3 +6,4 @@ import regression import dimensionality_reduction import tutorials import stochastic +import non_gaussian \ No newline at end of file diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 64707298..1d998fcb 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -19,12 +19,16 @@ class VarDTC(object): """ const_jitter = 1e-6 - def __init__(self): + def __init__(self, limit=1): #self._YYTfactor_cache = caching.cache() from ...util.caching import Cacher - self.get_trYYT = Cacher(self._get_trYYT, 1) - self.get_YYTfactor = Cacher(self._get_YYTfactor, 1) + self.get_trYYT = Cacher(self._get_trYYT, limit) + self.get_YYTfactor = Cacher(self._get_YYTfactor, limit) + def set_limit(self, limit): + self.get_trYYT.limit = limit + self.get_YYTfactor.limit = limit + def _get_trYYT(self, Y): return param_to_array(np.sum(np.square(Y))) @@ -175,11 +179,14 @@ class VarDTC(object): return post, log_marginal, grad_dict class VarDTCMissingData(object): - def __init__(self): + def __init__(self, limit=1): from ...util.caching import Cacher - self._Y = Cacher(self._subarray_computations, 1) + self._Y = Cacher(self._subarray_computations, limit) pass + def set_limit(self, limit): + self._Y.limit = limit + def _subarray_computations(self, Y): inan = np.isnan(Y) has_none = inan.any() diff --git a/GPy/kern/_src/coregionalize.py b/GPy/kern/_src/coregionalize.py index 6679eba4..1381b611 100644 --- a/GPy/kern/_src/coregionalize.py +++ b/GPy/kern/_src/coregionalize.py @@ -51,8 +51,6 @@ class Coregionalize(Kern): assert kappa.shape==(self.output_dim, ) self.kappa = Param('kappa', kappa, Logexp()) self.add_parameters(self.W, self.kappa) - self.parameters_changed() - def parameters_changed(self): self.B = np.dot(self.W, self.W.T) + np.diag(self.kappa) diff --git a/GPy/kern/_src/periodic.py b/GPy/kern/_src/periodic.py index e4e659a2..36ff3527 100644 --- a/GPy/kern/_src/periodic.py +++ b/GPy/kern/_src/periodic.py @@ -34,7 +34,6 @@ class Periodic(Kern): self.lengthscale = Param('lengthscale', np.float64(lengthscale), Logexp()) self.period = Param('period', np.float64(period), Logexp()) self.add_parameters(self.variance, self.lengthscale, self.period) - self.parameters_changed() def _cos(self, alpha, omega, phase): def f(x): diff --git a/GPy/kern/_src/sympykern.py b/GPy/kern/_src/sympykern.py index 0688682a..920f47f3 100644 --- a/GPy/kern/_src/sympykern.py +++ b/GPy/kern/_src/sympykern.py @@ -124,9 +124,6 @@ class Sympykern(Kern): # generate the code for the covariance functions self._gen_code() - self.parameters_changed() # initializes caches - - def __add__(self,other): return spkern(self._sp_k+other._sp_k) diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index 8763426a..3617e260 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -49,7 +49,6 @@ class BayesianGPLVM(SparseGP): SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, name, **kwargs) self.add_parameter(self.X, index=0) - self.parameters_changed() def _getstate(self): """ diff --git a/GPy/models/gp_regression.py b/GPy/models/gp_regression.py index f8957906..5e83db09 100644 --- a/GPy/models/gp_regression.py +++ b/GPy/models/gp_regression.py @@ -28,7 +28,6 @@ class GPRegression(GP): likelihood = likelihoods.Gaussian() super(GPRegression, self).__init__(X, Y, kernel, likelihood, name='GP regression') - self.parameters_changed() def _getstate(self): return GP._getstate(self) diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 0423aecd..dd1c44ba 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -1,14 +1,17 @@ # ## Copyright (c) 2013, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -from GPy.core import Model -from GPy.core import SparseGP -from GPy.util.linalg import PCA -import numpy +import numpy as np import itertools import pylab -from GPy.kern import Kern -from GPy.models.bayesian_gplvm import BayesianGPLVM + +from ..core import Model, SparseGP +from ..util.linalg import PCA +from ..kern import Kern +from bayesian_gplvm import BayesianGPLVM +from ..core.parameterization.variational import NormalPosterior, NormalPrior +from ..inference.latent_function_inference.var_dtc import VarDTCMissingData +from ..likelihoods.gaussian import Gaussian class MRD2(Model): """ @@ -20,11 +23,101 @@ class MRD2(Model): to match up, whereas the dimensionality p_d can differ. :param [array-like] Ylist: List of datasets to apply MRD on - :param array-like q_mean: mean of starting latent space q in [n x q] - :param array-like q_variance: variance of starting latent space q in [n x q] - :param :class:`~GPy.inference.latent_function_inference + :param input_dim: latent dimensionality + :type input_dim: int + :param array-like X: mean of starting latent space q in [n x q] + :param array-like X_variance: variance of starting latent space q in [n x q] + :param initx: initialisation method for the latent space : + + * 'concat' - PCA on concatenation of all datasets + * 'single' - Concatenation of PCA on datasets, respectively + * 'random' - Random draw from a Normal(0,1) + + :type initx: ['concat'|'single'|'random'] + :param initz: initialisation method for inducing inputs + :type initz: 'permute'|'random' + :param num_inducing: number of inducing inputs to use + :param Z: initial inducing inputs + :param kernel: list of kernels or kernel to copy for each output + :type kernel: [GPy.kern.kern] | GPy.kern.kern | None (default) + :param :class:`~GPy.inference.latent_function_inference inference_method: the inference method to use + :param :class:`~GPy.likelihoods.likelihood.Likelihood` likelihood: the likelihood to use + :param str name: the name of this model """ + def __init__(self, Ylist, input_dim, X=None, X_variance=None, + initx = 'PCA', initz = 'permute', + num_inducing=10, Z=None, kernel=None, + inference_method=None, likelihood=None, name='mrd'): + super(MRD2, self).__init__(name) + + # sort out the kernels + if kernel is None: + from ..kern import RBF + self.kern = [RBF(input_dim, ARD=1, name='Y_{}'.format(i)) for i in range(len(Ylist))] + elif isinstance(kernel, Kern): + self.kern = [kernel.copy(name='Y_{}'.format(i)) for i in range(len(Ylist))] + else: + assert len(kernel) == len(Ylist), "need one kernel per output" + assert all([isinstance(k, Kern) for k in kernel]), "invalid kernel object detected!" + + self.input_dim = input_dim + self.num_inducing = num_inducing + + self._in_init_ = True + X = self._init_X(initx, Ylist) + self.Z = self._init_Z(initz, X) + self.num_inducing = self.Z.shape[0] # ensure M==N if M>N + + if X_variance is None: + X_variance = np.random.uniform(0,.2,X.shape) + + self.variational_prior = NormalPrior() + self.X = NormalPosterior(X, X_variance) + + if likelihood is None: + likelihood = Gaussian() + + if inference_method is None: + if any(np.any(np.isnan(y)) for y in Ylist): + self.inference_method = VarDTCMissingData(limit=len(Ylist)) + + self.Ylist = Ylist + + def parameters_changed(self): + for y in self.Ylist: + pass + + def _init_X(self, init='PCA', likelihood_list=None): + if likelihood_list is None: + likelihood_list = self.likelihood_list + Ylist = [] + for likelihood_or_Y in likelihood_list: + if type(likelihood_or_Y) is np.ndarray: + Ylist.append(likelihood_or_Y) + else: + Ylist.append(likelihood_or_Y.Y) + del likelihood_list + if init in "PCA_concat": + X = PCA(np.hstack(Ylist), self.input_dim)[0] + elif init in "PCA_single": + X = np.zeros((Ylist[0].shape[0], self.input_dim)) + for qs, Y in itertools.izip(np.array_split(np.arange(self.input_dim), len(Ylist)), Ylist): + X[:, qs] = PCA(Y, len(qs))[0] + else: # init == 'random': + X = np.random.randn(Ylist[0].shape[0], self.input_dim) + self.X = X + return X + + def _init_Z(self, init="permute", X=None): + if X is None: + X = self.X + if init in "permute": + Z = np.random.permutation(X.copy())[:self.num_inducing] + elif init in "random": + Z = np.random.randn(self.num_inducing, self.input_dim) * X.var() + self.Z = Z + return Z class MRD(Model): """ @@ -84,7 +177,7 @@ class MRD(Model): del self._in_init_ self.gref = self.bgplvms[0] - nparams = numpy.array([0] + [SparseGP._get_params(g).size - g.Z.size for g in self.bgplvms]) + nparams = np.array([0] + [SparseGP._get_params(g).size - g.Z.size for g in self.bgplvms]) self.nparams = nparams.cumsum() self.num_data = self.gref.num_data @@ -216,7 +309,7 @@ class MRD(Model): X_var = self.gref.X_variance.ravel() Z = self.gref.Z.ravel() thetas = [SparseGP._get_params(g)[g.Z.size:] for g in self.bgplvms] - params = numpy.hstack([X, X_var, Z, numpy.hstack(thetas)]) + params = np.hstack([X, X_var, Z, np.hstack(thetas)]) return params # def _set_var_params(self, g, X, X_var, Z): @@ -239,13 +332,13 @@ class MRD(Model): # set params for all: for g, s, e in itertools.izip(self.bgplvms, self.nparams, self.nparams[1:]): - g._set_params(numpy.hstack([X, X_var, Z, thetas[s:e]])) + g._set_params(np.hstack([X, X_var, Z, thetas[s:e]])) # self._set_var_params(g, X, X_var, Z) # self._set_kern_params(g, thetas[s:e].copy()) # g._compute_kernel_matrices() # if self.auto_scale_factor: -# g.scale_factor = numpy.sqrt(g.psi2.sum(0).mean() * g.likelihood.precision) -# # self.scale_factor = numpy.sqrt(self.psi2.sum(0).mean() * self.likelihood.precision) +# g.scale_factor = np.sqrt(g.psi2.sum(0).mean() * g.likelihood.precision) +# # self.scale_factor = np.sqrt(self.psi2.sum(0).mean() * self.likelihood.precision) # g._computations() @@ -264,48 +357,18 @@ class MRD(Model): dKLmu, dKLdS = self.gref.dKL_dmuS() dLdmu -= dKLmu dLdS -= dKLdS - dLdmuS = numpy.hstack((dLdmu.flatten(), dLdS.flatten())).flatten() + dLdmuS = np.hstack((dLdmu.flatten(), dLdS.flatten())).flatten() dldzt1 = reduce(lambda a, b: a + b, (SparseGP._log_likelihood_gradients(g)[:self.MQ] for g in self.bgplvms)) - return numpy.hstack((dLdmuS, + return np.hstack((dLdmuS, dldzt1, - numpy.hstack([numpy.hstack([g.dL_dtheta(), + np.hstack([np.hstack([g.dL_dtheta(), g.likelihood._gradients(\ partial=g.partial_for_likelihood)]) \ for g in self.bgplvms]))) - def _init_X(self, init='PCA', likelihood_list=None): - if likelihood_list is None: - likelihood_list = self.likelihood_list - Ylist = [] - for likelihood_or_Y in likelihood_list: - if type(likelihood_or_Y) is numpy.ndarray: - Ylist.append(likelihood_or_Y) - else: - Ylist.append(likelihood_or_Y.Y) - del likelihood_list - if init in "PCA_concat": - X = PCA(numpy.hstack(Ylist), self.input_dim)[0] - elif init in "PCA_single": - X = numpy.zeros((Ylist[0].shape[0], self.input_dim)) - for qs, Y in itertools.izip(numpy.array_split(numpy.arange(self.input_dim), len(Ylist)), Ylist): - X[:, qs] = PCA(Y, len(qs))[0] - else: # init == 'random': - X = numpy.random.randn(Ylist[0].shape[0], self.input_dim) - self.X = X - return X - def _init_Z(self, init="permute", X=None): - if X is None: - X = self.X - if init in "permute": - Z = numpy.random.permutation(X.copy())[:self.num_inducing] - elif init in "random": - Z = numpy.random.randn(self.num_inducing, self.input_dim) * X.var() - self.Z = Z - return Z - def _handle_plotting(self, fignum, axes, plotf, sharex=False, sharey=False): if axes is None: fig = pylab.figure(num=fignum) @@ -358,7 +421,7 @@ class MRD(Model): """ if titles is None: titles = [r'${}$'.format(name) for name in self.names] - ymax = reduce(max, [numpy.ceil(max(g.input_sensitivity())) for g in self.bgplvms]) + ymax = reduce(max, [np.ceil(max(g.input_sensitivity())) for g in self.bgplvms]) def plotf(i, g, ax): ax.set_ylim([0,ymax]) g.kern.plot_ARD(ax=ax, title=titles[i], *args, **kwargs) diff --git a/GPy/plotting/matplot_dep/kernel_plots.py b/GPy/plotting/matplot_dep/kernel_plots.py index b3227e43..2b990611 100644 --- a/GPy/plotting/matplot_dep/kernel_plots.py +++ b/GPy/plotting/matplot_dep/kernel_plots.py @@ -23,7 +23,7 @@ def add_bar_labels(fig, ax, bars, bottom=0): xi = patch.get_x() + patch.get_width() / 2. va = 'top' c = 'w' - t = TextPath((0, 0), "${xi}$".format(xi=xi), rotation=0, usetex=True, ha='center') + t = TextPath((0, 0), "${xi}$".format(xi=xi), rotation=0, ha='center') transform = transOffset if patch.get_extents().height <= t.get_extents().height + 3: va = 'bottom' diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index d4105e3c..d2a236dd 100644 --- a/GPy/testing/likelihood_tests.py +++ b/GPy/testing/likelihood_tests.py @@ -10,7 +10,7 @@ from functools import partial #np.random.seed(300) #np.random.seed(7) -np.seterr(divide='raise') +#np.seterr(divide='raise') def dparam_partial(inst_func, *args): """ If we have a instance method that needs to be called but that doesn't @@ -350,7 +350,7 @@ class TestNoiseModels(object): def t_logpdf(self, model, Y, f): print "\n{}".format(inspect.stack()[0][3]) print model - print model._get_params() + #print model._get_params() np.testing.assert_almost_equal( model.pdf(f.copy(), Y.copy()), np.exp(model.logpdf(f.copy(), Y.copy())) @@ -664,7 +664,8 @@ class LaplaceTests(unittest.TestCase): print m1 print m2 - m2._set_params(m1._get_params()) + m2.parameters_changed() + #m2._set_params(m1._get_params()) #Predict for training points to get posterior mean and variance post_mean, post_var, _, _ = m1.predict(X) @@ -700,7 +701,8 @@ class LaplaceTests(unittest.TestCase): np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2) #Check marginals are the same with random m1.randomize() - m2._set_params(m1._get_params()) + #m2._set_params(m1._get_params()) + m2.parameters_changed() np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2) #Check they are checkgradding diff --git a/GPy/testing/observable_tests.py b/GPy/testing/observable_tests.py index 4123f58a..ebda1630 100644 --- a/GPy/testing/observable_tests.py +++ b/GPy/testing/observable_tests.py @@ -8,14 +8,17 @@ from GPy.core.parameterization.parameterized import Parameterized from GPy.core.parameterization.param import Param import numpy +# One trigger in init +_trigger_start = -1 class ParamTestParent(Parameterized): - parent_changed_count = 0 + parent_changed_count = _trigger_start def parameters_changed(self): self.parent_changed_count += 1 class ParameterizedTest(Parameterized): - params_changed_count = 0 + # One trigger after initialization + params_changed_count = _trigger_start def parameters_changed(self): self.params_changed_count += 1 def _set_params(self, params, trigger_parent=True): @@ -92,29 +95,31 @@ class Test(unittest.TestCase): def test_set_params(self): self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') - self.par._set_params(numpy.ones(self.par.size)) + self.par._param_array_[:] = 1 + self.par._trigger_params_changed() self.assertEqual(self.par.params_changed_count, 1, 'now params changed') self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) - self.parent._set_params(numpy.ones(self.parent.size) * 2) + self.par._param_array_[:] = 2 + self.par._trigger_params_changed() self.assertEqual(self.par.params_changed_count, 2, 'now params changed') self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) def test_priority_notify(self): self.assertEqual(self.par.params_changed_count, 0) - self.par._notify_observers(0, None) + self.par.notify_observers(0, None) self.assertEqual(self.par.params_changed_count, 1) self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) - self.par._notify_observers(0, -numpy.inf) + self.par.notify_observers(0, -numpy.inf) self.assertEqual(self.par.params_changed_count, 2) self.assertEqual(self.parent.parent_changed_count, 1) def test_priority(self): self.par.add_observer(self, self._trigger, -1) self.par.add_observer(self, self._trigger_priority, 0) - self.par._notify_observers(0) + self.par.notify_observers(0) self.assertEqual(self._first, self._trigger_priority, 'priority should be first') self.assertEqual(self._second, self._trigger, 'priority should be first') @@ -123,7 +128,7 @@ class Test(unittest.TestCase): self.par.add_observer(self, self._trigger, 1) self.par.add_observer(self, self._trigger_priority, 0) - self.par._notify_observers(0) + self.par.notify_observers(0) self.assertEqual(self._first, self._trigger, 'priority should be second') self.assertEqual(self._second, self._trigger_priority, 'priority should be second') diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index 0da3f3ae..b2f57144 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -22,6 +22,10 @@ class Test(unittest.TestCase): self.test1.add_parameter(self.rbf, 0) self.test1.add_parameter(self.param) + x = np.linspace(-2,6,4)[:,None] + y = np.sin(x) + self.testmodel = GPy.models.GPRegression(x,y) + def test_add_parameter(self): self.assertEquals(self.rbf._parent_index_, 0) self.assertEquals(self.white._parent_index_, 1) @@ -38,7 +42,6 @@ class Test(unittest.TestCase): 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__, Logexp self.white.fix() @@ -89,6 +92,18 @@ class Test(unittest.TestCase): self.assertEqual(self.rbf.constraints._offset, 0) self.assertEqual(self.param.constraints._offset, 3) + def test_fixing_randomize(self): + self.white.fix(warning=False) + val = float(self.test1.white.variance) + self.test1.randomize() + self.assertEqual(val, self.white.variance) + + def test_fixing_optimize(self): + self.testmodel.kern.lengthscale.fix() + val = float(self.testmodel.kern.lengthscale) + self.testmodel.randomize() + self.assertEqual(val, self.testmodel.kern.lengthscale) + if __name__ == "__main__": #import sys;sys.argv = ['', 'Test.test_add_parameter'] unittest.main() \ No newline at end of file From 0a82427e295fd8c1e73918753fc3e9bf150957dd Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 5 Mar 2014 10:45:35 +0000 Subject: [PATCH 03/11] parameters once in memory --- GPy/core/parameterization/param.py | 2 +- GPy/core/parameterization/parameter_core.py | 84 ++++++++++++--------- GPy/core/parameterization/parameterized.py | 4 +- GPy/core/parameterization/variational.py | 2 +- 4 files changed, 51 insertions(+), 41 deletions(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 2917cad7..4d867487 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -15,7 +15,7 @@ __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision __print_threshold__ = 5 ###### -class Param(OptimizationHandlable, ObservableArray, Gradcheckable): +class Param(OptimizationHandlable, ObservableArray): """ Parameter object for GPy models. diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index c9372b58..a78cf02d 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -195,7 +195,41 @@ class Parentable(object): Dont do anything if in leaf node """ pass -class Nameable(Parentable): + +class Gradcheckable(Parentable): + """ + Adds the functionality for an object to be gradcheckable. + It is just a thin wrapper of a call to the highest parent for now. + TODO: Can be done better, by only changing parameters of the current parameter handle, + such that object hierarchy only has to change for those. + """ + def __init__(self, *a, **kw): + super(Gradcheckable, self).__init__(*a, **kw) + def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): + """ + Check the gradient of this parameter with respect to the highest parent's + objective function. + This is a three point estimate of the gradient, wiggling at the parameters + with a stepsize step. + The check passes if either the ratio or the difference between numerical and + analytical gradient is smaller then tolerance. + + :param bool verbose: whether each parameter shall be checked individually. + :param float step: the stepsize for the numerical three point gradient estimate. + :param flaot tolerance: the tolerance for the gradient ratio or difference. + """ + if self.has_parent(): + return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) + return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance) + def _checkgrad(self, param): + """ + Perform the checkgrad on the model. + TODO: this can be done more efficiently, when doing it inside here + """ + raise NotImplementedError, "Need log likelihood to check gradient against" + + +class Nameable(Gradcheckable): """ Make an object nameable inside the hierarchy. """ @@ -233,40 +267,6 @@ class Nameable(Parentable): return self._parent_.hierarchy_name() + "." + adjust(self.name) return adjust(self.name) - -class Gradcheckable(Parentable): - """ - Adds the functionality for an object to be gradcheckable. - It is just a thin wrapper of a call to the highest parent for now. - TODO: Can be done better, by only changing parameters of the current parameter handle, - such that object hierarchy only has to change for those. - """ - def __init__(self, *a, **kw): - super(Gradcheckable, self).__init__(*a, **kw) - def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): - """ - Check the gradient of this parameter with respect to the highest parent's - objective function. - This is a three point estimate of the gradient, wiggling at the parameters - with a stepsize step. - The check passes if either the ratio or the difference between numerical and - analytical gradient is smaller then tolerance. - - :param bool verbose: whether each parameter shall be checked individually. - :param float step: the stepsize for the numerical three point gradient estimate. - :param flaot tolerance: the tolerance for the gradient ratio or difference. - """ - if self.has_parent(): - return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) - return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance) - def _checkgrad(self, param): - """ - Perform the checkgrad on the model. - TODO: this can be done more efficiently, when doing it inside here - """ - raise NotImplementedError, "Need log likelihood to check gradient against" - - class Indexable(object): """ Enable enraveled indexes and offsets for this object. @@ -551,8 +551,10 @@ class OptimizationHandlable(Constrainable, Observable): return p def _set_params_transformed(self, p): - if self._has_fixes(): self._param_array_[self._fixes_] = p.copy() - else: self._param_array_[:] = p.copy() + if p is self._param_array_: + p = p.copy() + if self._has_fixes(): self._param_array_[self._fixes_] = p + else: self._param_array_[:] = p self.untransform() self._trigger_params_changed() @@ -611,6 +613,11 @@ class OptimizationHandlable(Constrainable, Observable): """ Randomize the model. Make this draw from the prior if one exists, else draw from given random generator + + :param rand_gen: numpy random number generator which takes args and kwargs + :param flaot loc: loc parameter for random number generator + :param float scale: scale parameter for random number generator + :param args, kwargs: will be passed through to random number generator """ # first take care of all parameters (from N(0,1)) x = rand_gen(loc=loc, scale=scale, size=self._size_transformed(), *args, **kwargs) @@ -623,6 +630,9 @@ class Parameterizable(OptimizationHandlable): super(Parameterizable, self).__init__(*args, **kwargs) from GPy.core.parameterization.lists_and_dicts import ArrayList _parameters_ = ArrayList() + self.size = 0 + self._param_array_ = np.empty(self.size, dtype=np.float64) + self._gradient_array_ = np.empty(self.size, dtype=np.float64) self._added_names_ = set() def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 50dd7b3b..6d06018a 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -7,7 +7,7 @@ import cPickle import itertools from re import compile, _pattern_type from param import ParamConcatenation -from parameter_core import Pickleable, Parameterizable, adjust_name_for_printing, Gradcheckable +from parameter_core import Pickleable, Parameterizable, adjust_name_for_printing from transformations import __fixed__ from lists_and_dicts import ArrayList @@ -17,7 +17,7 @@ class ParametersChangedMeta(type): instance.parameters_changed() return instance -class Parameterized(Parameterizable, Pickleable, Gradcheckable): +class Parameterized(Parameterizable, Pickleable): """ Parameterized class diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py index a209cb39..71921ab1 100644 --- a/GPy/core/parameterization/variational.py +++ b/GPy/core/parameterization/variational.py @@ -10,7 +10,7 @@ from param import Param from transformations import Logexp, Logistic class VariationalPrior(Parameterized): - def __init__(self, name=None, **kw): + def __init__(self, name='latent space', **kw): super(VariationalPrior, self).__init__(name=name, **kw) def KL_divergence(self, variational_posterior): From cde8722e1bfa79a7769c4b34494f579a1f8a2132 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 5 Mar 2014 11:19:46 +0000 Subject: [PATCH 04/11] printing for older numpy versions --- GPy/core/parameterization/param.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 4d867487..f8beb014 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -250,7 +250,7 @@ class Param(OptimizationHandlable, ObservableArray): @property def _description_str(self): if self.size <= 1: - return [str(numpy.take(self, 0))] + return [str(self.view(numpy.ndarray)[0])] else: return [str(self.shape)] def parameter_names(self, add_self=False, adjust_for_printing=False): if adjust_for_printing: From 8e22373a00cb1daf5782fde8148c425e87388db3 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 5 Mar 2014 11:32:54 +0000 Subject: [PATCH 05/11] some missing .Ks --- GPy/kern/_src/prod.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index bb809356..b324eaa7 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -39,17 +39,17 @@ class Prod(Kern): return self.k1.Kdiag(X[:,self.slice1]) * self.k2.Kdiag(X[:,self.slice2]) def update_gradients_full(self, dL_dK, X): - self.k1.update_gradients_full(dL_dK*self.k2(X[:,self.slice2]), X[:,self.slice1]) - self.k2.update_gradients_full(dL_dK*self.k1(X[:,self.slice1]), X[:,self.slice2]) + self.k1.update_gradients_full(dL_dK*self.k2.K(X[:,self.slice2]), X[:,self.slice1]) + self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2]) def gradients_X(self, dL_dK, X, X2=None): target = np.zeros(X.shape) if X2 is None: - target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2(X[:,self.slice2]), X[:,self.slice1], None) - target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1(X[:,self.slice1]), X[:,self.slice2], None) + target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2.K(X[:,self.slice2]), X[:,self.slice1], None) + target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2], None) else: - target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2(X[:,self.slice2], X2[:,self.slice2]), X[:,self.slice1], X2[:,self.slice1]) - target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1(X[:,self.slice1], X2[:,self.slice1]), X[:,self.slice2], X2[:,self.slice2]) + target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2.K(X[:,self.slice2], X2[:,self.slice2]), X[:,self.slice1], X2[:,self.slice1]) + target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1.K(X[:,self.slice1], X2[:,self.slice1]), X[:,self.slice2], X2[:,self.slice2]) return target def gradients_X_diag(self, dL_dKdiag, X): From f5e3e97794ea3049ee87f25799b616d64aa9d888 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 5 Mar 2014 12:44:53 +0000 Subject: [PATCH 06/11] gradcheck global diff --- GPy/core/model.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 38b5eb29..343d5f08 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -298,9 +298,11 @@ class Model(Parameterized): dx = dx[transformed_index] gradient = gradient[transformed_index] - + global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) - return (np.abs(1. - global_ratio) < tolerance) + num_grad =(np.abs((f1-f2)/-(2*dx)*np.where(gradient == 0, 1e-32, gradient))).mean() + + return (np.abs(1. - global_ratio) < tolerance) or (num_grad < tolerance) else: # check the gradient of each parameter individually, and do some pretty printing try: From b1ebeea9121b32e6648bf6ce20c5159c681a056c Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 5 Mar 2014 12:46:13 +0000 Subject: [PATCH 07/11] param concat fix --- GPy/core/parameterization/param.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index f8beb014..3ebeb566 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -393,7 +393,7 @@ class ParamConcatenation(object): if update: self.update_all_params() def _vals(self): - return numpy.hstack([p._get_params() for p in self.params]) + return numpy.hstack([p._param_array_ for p in self.params]) #=========================================================================== # parameter operations: #=========================================================================== From 988bad88a31434f5821a2b40e3f03e49425ecde2 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 5 Mar 2014 12:50:46 +0000 Subject: [PATCH 08/11] numerical global diff in gradcheck --- GPy/core/model.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index 343d5f08..e7924b2b 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -300,9 +300,9 @@ class Model(Parameterized): gradient = gradient[transformed_index] global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) - num_grad =(np.abs((f1-f2)/-(2*dx)*np.where(gradient == 0, 1e-32, gradient))).mean() + gloabl_diff = (f1 - f2) - (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) - return (np.abs(1. - global_ratio) < tolerance) or (num_grad < tolerance) + return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gloabl_diff) < tolerance) else: # check the gradient of each parameter individually, and do some pretty printing try: From 19c87c9f77ea43e02ed34ddd7b7306d0fb94af61 Mon Sep 17 00:00:00 2001 From: Ricardo Date: Wed, 5 Mar 2014 12:52:56 +0000 Subject: [PATCH 09/11] name added as a parameter of Prod --- GPy/kern/_src/kern.py | 4 ++-- GPy/kern/_src/prod.py | 8 +++++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index f632783b..47166156 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -129,7 +129,7 @@ class Kern(Parameterized): """ return self.prod(other, tensor=True) - def prod(self, other, tensor=False): + def prod(self, other, tensor=False, name=None): """ Multiply two kernels (either on the same space, or on the tensor product of the input space). @@ -142,4 +142,4 @@ class Kern(Parameterized): """ assert isinstance(other, Kern), "only kernels can be added to kernels..." from prod import Prod - return Prod(self, other, tensor) + return Prod(self, other, tensor, name) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index b324eaa7..51490687 100644 --- a/GPy/kern/_src/prod.py +++ b/GPy/kern/_src/prod.py @@ -15,14 +15,16 @@ class Prod(Kern): :rtype: kernel object """ - def __init__(self, k1, k2, tensor=False): + def __init__(self, k1, k2, tensor=False,name=None): if tensor: - super(Prod, self).__init__(k1.input_dim + k2.input_dim, k1.name + '_xx_' + k2.name) + name = k1.name + '_xx_' + k2.name if name is None else name + super(Prod, self).__init__(k1.input_dim + k2.input_dim, name) self.slice1 = slice(0,k1.input_dim) self.slice2 = slice(k1.input_dim,k1.input_dim+k2.input_dim) else: assert k1.input_dim == k2.input_dim, "Error: The input spaces of the kernels to multiply don't have the same dimension." - super(Prod, self).__init__(k1.input_dim, k1.name + '_x_' + k2.name) + name = k1.name + '_x_' + k2.name if name is None else name + super(Prod, self).__init__(k1.input_dim, name) self.slice1 = slice(0, self.input_dim) self.slice2 = slice(0, self.input_dim) self.k1 = k1 From cf5c6bf227ed39e90b9be013a6b184575f4523bb Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Wed, 5 Mar 2014 13:01:41 +0000 Subject: [PATCH 10/11] checkgrad divide by zero catches --- GPy/core/model.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/GPy/core/model.py b/GPy/core/model.py index e7924b2b..a858a62d 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -299,8 +299,9 @@ class Model(Parameterized): dx = dx[transformed_index] gradient = gradient[transformed_index] - global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) - gloabl_diff = (f1 - f2) - (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) + denominator = (2 * np.dot(dx, gradient)) + global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator) + gloabl_diff = (f1 - f2) - denominator return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gloabl_diff) < tolerance) else: @@ -348,7 +349,8 @@ class Model(Parameterized): xx[xind] -= 2.*step f2 = self.objective_function(xx) numerical_gradient = (f1 - f2) / (2 * step) - ratio = (f1 - f2) / (2 * step * gradient[xind]) + if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind] + else: 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: From 2f5d5dd3bfd5338f1707f921882682bcafedec74 Mon Sep 17 00:00:00 2001 From: Alan Saul Date: Wed, 5 Mar 2014 14:16:53 +0000 Subject: [PATCH 11/11] Made sampling default for non-gaussian likelihoods as a quick fix to allow plotting again for likelihoods without predictive values --- GPy/examples/regression.py | 2 +- GPy/kern/_src/sympykern.py | 29 ++++++++++++++--------------- GPy/likelihoods/likelihood.py | 2 +- 3 files changed, 16 insertions(+), 17 deletions(-) diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index aa6bbbf9..cc23410a 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -284,7 +284,7 @@ def toy_poisson_rbf_1d_laplace(optimize=True, plot=True): kern = GPy.kern.RBF(1) poisson_lik = GPy.likelihoods.Poisson() - laplace_inf = GPy.inference.latent_function_inference.LaplaceInference() + laplace_inf = GPy.inference.latent_function_inference.Laplace() # create simple GP Model m = GPy.core.GP(X, Y, kernel=kern, likelihood=poisson_lik, inference_method=laplace_inf) diff --git a/GPy/kern/_src/sympykern.py b/GPy/kern/_src/sympykern.py index 920f47f3..9878ec68 100644 --- a/GPy/kern/_src/sympykern.py +++ b/GPy/kern/_src/sympykern.py @@ -1,11 +1,10 @@ # Check Matthew Rocklin's blog post. -try: +try: import sympy as sp sympy_available=True from sympy.utilities.lambdify import lambdify except ImportError: sympy_available=False - exit() import numpy as np from kern import Kern @@ -36,7 +35,7 @@ class Sympykern(Kern): super(Sympykern, self).__init__(input_dim, name) self._sp_k = k - + # pull the variable names out of the symbolic covariance function. sp_vars = [e for e in k.atoms() if e.is_Symbol] self._sp_x= sorted([e for e in sp_vars if e.name[0:2]=='x_'],key=lambda x:int(x.name[2:])) @@ -51,7 +50,7 @@ class Sympykern(Kern): self._sp_kdiag = k for x, z in zip(self._sp_x, self._sp_z): self._sp_kdiag = self._sp_kdiag.subs(z, x) - + # If it is a multi-output covariance, add an input for indexing the outputs. self._real_input_dim = x_dim # Check input dim is number of xs + 1 if output_dim is >1 @@ -73,7 +72,7 @@ class Sympykern(Kern): # Extract names of shared parameters (those without a subscript) self._sp_theta = [theta for theta in thetas if theta not in self._sp_theta_i and theta not in self._sp_theta_j] - + self.num_split_params = len(self._sp_theta_i) self._split_theta_names = ["%s"%theta.name[:-2] for theta in self._sp_theta_i] # Add split parameters to the model. @@ -82,11 +81,11 @@ class Sympykern(Kern): setattr(self, theta, Param(theta, np.ones(self.output_dim), None)) self.add_parameter(getattr(self, theta)) - + self.num_shared_params = len(self._sp_theta) for theta_i, theta_j in zip(self._sp_theta_i, self._sp_theta_j): self._sp_kdiag = self._sp_kdiag.subs(theta_j, theta_i) - + else: self.num_split_params = 0 self._split_theta_names = [] @@ -107,10 +106,10 @@ class Sympykern(Kern): derivative_arguments = self._sp_x + self._sp_theta if self.output_dim > 1: derivative_arguments += self._sp_theta_i - + self.derivatives = {theta.name : sp.diff(self._sp_k,theta).simplify() for theta in derivative_arguments} self.diag_derivatives = {theta.name : sp.diff(self._sp_kdiag,theta).simplify() for theta in derivative_arguments} - + # This gives the parameters for the arg list. self.arg_list = self._sp_x + self._sp_z + self._sp_theta self.diag_arg_list = self._sp_x + self._sp_theta @@ -137,7 +136,7 @@ class Sympykern(Kern): for key in self.derivatives.keys(): setattr(self, '_Kdiag_diff_' + key, lambdify(self.diag_arg_list, self.diag_derivatives[key], 'numpy')) - def K(self,X,X2=None): + def K(self,X,X2=None): self._K_computations(X, X2) return self._K_function(**self._arguments) @@ -145,11 +144,11 @@ class Sympykern(Kern): def Kdiag(self,X): self._K_computations(X) return self._Kdiag_function(**self._diag_arguments) - + def _param_grad_helper(self,partial,X,Z,target): pass - + def gradients_X(self, dL_dK, X, X2=None): #if self._X is None or X.base is not self._X.base or X2 is not None: self._K_computations(X, X2) @@ -168,7 +167,7 @@ class Sympykern(Kern): gf = getattr(self, '_Kdiag_diff_' + x.name) dX[:, i] = gf(**self._diag_arguments)*dL_dK return dX - + def update_gradients_full(self, dL_dK, X, X2=None): # Need to extract parameters to local variables first self._K_computations(X, X2) @@ -193,7 +192,7 @@ class Sympykern(Kern): gradient += np.asarray([A[np.where(self._output_ind2==i)].T.sum() for i in np.arange(self.output_dim)]) setattr(parameter, 'gradient', gradient) - + def update_gradients_diag(self, dL_dKdiag, X): self._K_computations(X) @@ -209,7 +208,7 @@ class Sympykern(Kern): setattr(parameter, 'gradient', np.asarray([a[np.where(self._output_ind==i)].sum() for i in np.arange(self.output_dim)])) - + def _K_computations(self, X, X2=None): """Set up argument lists for the derivatives.""" # Could check if this needs doing or not, there could diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py index 701a5a2f..aff55533 100644 --- a/GPy/likelihoods/likelihood.py +++ b/GPy/likelihoods/likelihood.py @@ -358,7 +358,7 @@ class Likelihood(Parameterized): return dlogpdf_dtheta, dlogpdf_df_dtheta, d2logpdf_df2_dtheta - def predictive_values(self, mu, var, full_cov=False, sampling=False, num_samples=10000): + def predictive_values(self, mu, var, full_cov=False, sampling=True, num_samples=10000): """ Compute mean, variance and conficence interval (percentiles 5 and 95) of the prediction.