From f3b74fa85ff628de673e3a652a04d30880646487 Mon Sep 17 00:00:00 2001 From: Max Zwiessele Date: Mon, 31 Mar 2014 12:45:09 +0100 Subject: [PATCH] pickling and caching --- GPy/core/gp.py | 25 -- GPy/core/model.py | 32 +- GPy/core/parameterization/array_core.py | 105 +---- GPy/core/parameterization/index_operations.py | 13 +- GPy/core/parameterization/lists_and_dicts.py | 73 ++++ GPy/core/parameterization/param.py | 125 ++---- GPy/core/parameterization/parameter_core.py | 389 ++++++++---------- GPy/core/parameterization/parameterized.py | 31 +- GPy/core/parameterization/variational.py | 3 +- GPy/core/sparse_gp.py | 12 - GPy/core/svigp.py | 102 ++--- GPy/examples/dimensionality_reduction.py | 9 +- .../latent_function_inference/var_dtc.py | 12 +- GPy/kern/_src/kern.py | 19 +- GPy/kern/_src/kernel_slice_operations.py | 2 +- GPy/models/bayesian_gplvm.py | 13 - GPy/models/gp_regression.py | 5 - GPy/models/gplvm.py | 6 - GPy/models/mrd.py | 13 +- GPy/models/sparse_gp_classification.py | 8 - GPy/models/sparse_gp_regression.py | 8 - GPy/models/sparse_gplvm.py | 8 - GPy/models/svigp_regression.py | 7 - GPy/models/warped_gp.py | 8 - GPy/testing/kernel_tests.py | 122 +++--- GPy/testing/observable_tests.py | 4 +- GPy/testing/parameterized_tests.py | 4 +- GPy/util/caching.py | 9 + 28 files changed, 481 insertions(+), 686 deletions(-) diff --git a/GPy/core/gp.py b/GPy/core/gp.py index 6fc127ea..490bcc72 100644 --- a/GPy/core/gp.py +++ b/GPy/core/gp.py @@ -214,28 +214,3 @@ class GP(Model): """ return self.kern.input_sensitivity() - def _getstate(self): - """ - - Get the current state of the class, here we return everything that is - needed to recompute the model. - - """ - - return []#Model._getstate(self) + [self.X, -# self.num_data, -# self.input_dim, -# self.kern, -# self.likelihood, -# self.output_dim, -# ] - - def _setstate(self, state): - return - self.output_dim = state.pop() - self.likelihood = state.pop() - self.kern = state.pop() - self.input_dim = state.pop() - self.num_data = state.pop() - self.X = state.pop() - Model._setstate(self, state) diff --git a/GPy/core/model.py b/GPy/core/model.py index e04993cb..a39eceda 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -27,33 +27,6 @@ class Model(Parameterized): def _log_likelihood_gradients(self): return self.gradient - def _getstate(self): - """ - Get the current state of the class. - Inherited from Parameterized, so add those parameters to the state - - :return: list of states from the model. - - """ - return Parameterized._getstate(self) + \ - [self.priors, self.optimization_runs, - self.sampling_runs, self.preferred_optimizer] - - def _setstate(self, state): - """ - set state from previous call to _getstate - call Parameterized with the rest of the state - - :param state: the state of the model. - :type state: list as returned from _getstate. - - """ - self.preferred_optimizer = state.pop() - self.sampling_runs = state.pop() - self.optimization_runs = state.pop() - self.priors = state.pop() - Parameterized._setstate(self, state) - def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs): """ Perform random restarts of the model, and set the model to the best @@ -318,7 +291,10 @@ class Model(Parameterized): denominator = (2 * np.dot(dx, gradient)) global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator) - return np.abs(1. - global_ratio) < tolerance or np.abs(f1-f2).sum() + np.abs((2 * np.dot(dx, gradient))).sum() < tolerance + global_diff = np.abs(f1 - f2) < tolerance and np.allclose(gradient, 0, atol=tolerance) + if global_ratio is np.nan: + global_ratio = 0 + return np.abs(1. - global_ratio) < tolerance or np.abs(f1-f2).sum() + np.abs((2 * np.dot(dx, gradient))).sum() < tolerance or global_diff else: # check the gradient of each parameter individually, and do some pretty printing try: diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index ab8214f2..fc9d6cf2 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -1,12 +1,12 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) -__updated__ = '2014-03-24' +__updated__ = '2014-03-31' import numpy as np -from parameter_core import Observable +from parameter_core import Observable, Pickleable -class ObsAr(np.ndarray, Observable): +class ObsAr(np.ndarray, Pickleable, Observable): """ An ndarray which reports changes to its observers. The observers can add themselves with a callable, which @@ -30,13 +30,25 @@ class ObsAr(np.ndarray, Observable): def __array_wrap__(self, out_arr, context=None): return out_arr.view(np.ndarray) + def copy(self): + memo = {} + memo[id(self)] = self + return self.__deepcopy__(memo) + + def __deepcopy__(self, memo): + s = self.__new__(self.__class__, input_array=self.view(np.ndarray).copy()) + memo[id(self)] = s + import copy + s.__dict__.update(copy.deepcopy(self.__dict__, memo)) + return s + def __reduce__(self): - func, args, state = np.ndarray.__reduce__(self) - return func, args, (state, Observable._getstate(self)) + func, args, state = super(ObsAr, self).__reduce__() + return func, args, (state, Pickleable.__getstate__(self)) def __setstate__(self, state): np.ndarray.__setstate__(self, state[0]) - Observable._setstate(self, state[1]) + Pickleable.__setstate__(self, state[1]) def __setitem__(self, s, val): super(ObsAr, self).__setitem__(s, val) @@ -48,12 +60,6 @@ class ObsAr(np.ndarray, Observable): def __setslice__(self, start, stop, val): return self.__setitem__(slice(start, stop), val) - def __copy__(self, *args): - return ObsAr(self.view(np.ndarray).copy()) - - def copy(self, *args): - return self.__copy__(*args) - def __ilshift__(self, *args, **kwargs): r = np.ndarray.__ilshift__(self, *args, **kwargs) self.notify_observers() @@ -128,77 +134,4 @@ class ObsAr(np.ndarray, Observable): def __imul__(self, *args, **kwargs): r = np.ndarray.__imul__(self, *args, **kwargs) self.notify_observers() - return r - - -# def __rrshift__(self, *args, **kwargs): -# r = np.ndarray.__rrshift__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __ror__(self, *args, **kwargs): -# r = np.ndarray.__ror__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rxor__(self, *args, **kwargs): -# r = np.ndarray.__rxor__(self, *args, **kwargs) -# self.notify_observers() -# return r - - - -# def __rdivmod__(self, *args, **kwargs): -# r = np.ndarray.__rdivmod__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __radd__(self, *args, **kwargs): -# r = np.ndarray.__radd__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rdiv__(self, *args, **kwargs): -# r = np.ndarray.__rdiv__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rtruediv__(self, *args, **kwargs): -# r = np.ndarray.__rtruediv__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rshift__(self, *args, **kwargs): -# r = np.ndarray.__rshift__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rmul__(self, *args, **kwargs): -# r = np.ndarray.__rmul__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rpow__(self, *args, **kwargs): -# r = np.ndarray.__rpow__(self, *args, **kwargs) -# self.notify_observers() -# return r - - -# def __rsub__(self, *args, **kwargs): -# r = np.ndarray.__rsub__(self, *args, **kwargs) -# self.notify_observers() -# return r - -# def __rfloordiv__(self, *args, **kwargs): -# r = np.ndarray.__rfloordiv__(self, *args, **kwargs) -# self.notify_observers() -# return r - + return r \ No newline at end of file diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index e2a041f7..ebfe2904 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -24,12 +24,6 @@ class ParameterIndexOperations(object): for t, i in constraints.iteritems(): self.add(t, i) - def __getstate__(self): - return self._properties - - def __setstate__(self, state): - self._properties = state - def iteritems(self): return self._properties.iteritems() @@ -92,8 +86,10 @@ class ParameterIndexOperations(object): for i, v in parameter_index_view.iteritems(): self.add(i, v+offset) - def copy(self): + return self.__deepcopy__(None) + + def __deepcopy__(self, memo): return ParameterIndexOperations(dict(self.iteritems())) def __getitem__(self, prop): @@ -203,6 +199,9 @@ class ParameterIndexOperationsView(object): def copy(self): + return self.__deepcopy__(None) + + def __deepcopy__(self, memo): return ParameterIndexOperations(dict(self.iteritems())) pass diff --git a/GPy/core/parameterization/lists_and_dicts.py b/GPy/core/parameterization/lists_and_dicts.py index 31235952..6902c249 100644 --- a/GPy/core/parameterization/lists_and_dicts.py +++ b/GPy/core/parameterization/lists_and_dicts.py @@ -36,3 +36,76 @@ class ArrayList(list): index += 1 raise ValueError, "{} is not in list".format(item) pass + +class ObservablesList(object): + def __init__(self): + self._poc = [] + + def remove(self, value): + return self._poc.remove(value) + + + def __delitem__(self, ind): + return self._poc.__delitem__(ind) + + + def __setitem__(self, ind, item): + return self._poc.__setitem__(ind, item) + + + def __getitem__(self, ind): + return self._poc.__getitem__(ind) + + + def __repr__(self): + return self._poc.__repr__() + + + def append(self, obj): + return self._poc.append(obj) + + + def index(self, value): + return self._poc.index(value) + + + def extend(self, iterable): + return self._poc.extend(iterable) + + + def __str__(self): + return self._poc.__str__() + + + def __iter__(self): + return self._poc.__iter__() + + + def insert(self, index, obj): + return self._poc.insert(index, obj) + + + def __len__(self): + return self._poc.__len__() + + def __deepcopy__(self, memo): + s = ObservablesList() + import copy + s._poc = copy.deepcopy(self._poc, memo) + return s + + def __getstate__(self): + from ...util.caching import Cacher + obs = [] + for p, o, c in self: + if (getattr(o, c.__name__, None) is not None + and not isinstance(o, Cacher)): + obs.append((p,o,c.__name__)) + return obs + + def __setstate__(self, state): + self._poc = [] + for p, o, c in state: + self._poc.append((p,o,getattr(o, c))) + + pass diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index de16a1a0..f89b09df 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -43,14 +43,13 @@ class Param(OptimizationHandlable, ObsAr): _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, name=name, default_constraint=default_constraint)) + 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 obj._realndim_ = obj.ndim obj._original_ = True - obj._gradient_array_ = numpy.zeros(obj.shape, dtype=numpy.float64) return obj def __init__(self, name, input_array, default_constraint=None, *a, **kw): @@ -87,74 +86,30 @@ class Param(OptimizationHandlable, ObsAr): self.priors = getattr(obj, 'priors', None) @property - def _param_array_(self): + def param_array(self): return self + @property + def current_slice(self): + if self._current_slice_ is None: + return slice(0, self.shape[0], 1) + return self._current_slice_ + @property def gradient(self): + """ + Return a view on the gradient, which is in the same shape as this parameter is. + Note: this is not the real gradient array, it is just a view on it. + + To work on the real gradient array use: self.full_gradient + """ + if getattr(self, '_gradient_array_', None) is None: + self._gradient_array_ = numpy.empty(self._realshape_, dtype=numpy.float64) return self._gradient_array_[self._current_slice_] @gradient.setter def gradient(self, val): - self.gradient[:] = val - - #=========================================================================== - # Pickling operations - #=========================================================================== - def __reduce__(self): - func, args, state = super(Param, self).__reduce__() - return func, args, (state, - (self._name, - self._parent_, - self._parent_index_, - self._default_constraint_, - self._current_slice_, - self._realshape_, - self._realsize_, - self._realndim_, - self.constraints, - self.priors - ) - ) - - def __setstate__(self, state): - super(Param, self).__setstate__(state[0]) - state = list(state[1]) - self.priors = state.pop() - self.constraints = state.pop() - self._realndim_ = state.pop() - self._realsize_ = state.pop() - self._realshape_ = state.pop() - self._current_slice_ = state.pop() - self._default_constraint_ = state.pop() - self._parent_index_ = state.pop() - self._parent_ = state.pop() - self._name = state.pop() - - def copy(self, *args): - constr = self.constraints.copy() - priors = self.priors.copy() - p = Param(self.name, self.view(numpy.ndarray).copy(), self._default_constraint_) - p.constraints = constr - p.priors = priors - return p - #=========================================================================== - # 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_) + self._gradient_array_[self._current_slice_] = val #=========================================================================== # Array operations -> done @@ -172,24 +127,6 @@ class Param(OptimizationHandlable, ObsAr): def __setitem__(self, s, val): super(Param, self).__setitem__(s, val) - #=========================================================================== - # Index Operations: - #=========================================================================== - #def _internal_offset(self): - # internal_offset = 0 - # extended_realshape = numpy.cumprod((1,) + self._realshape_[:0:-1])[::-1] - # for i, si in enumerate(self._current_slice_[:self._realndim_]): - # if numpy.all(si == Ellipsis): - # continue - # if isinstance(si, slice): - # a = si.indices(self._realshape_[i])[0] - # elif isinstance(si, (list,numpy.ndarray,tuple)): - # a = si[0] - # else: a = si - # if a < 0: - # a = self._realshape_[i] + a - # internal_offset += a * extended_realshape[i] - # return internal_offset def _raveled_index(self, slice_index=None): # return an index array on the raveled array, which is formed by the current_slice @@ -235,13 +172,21 @@ class Param(OptimizationHandlable, ObsAr): def is_fixed(self): from transformations import __fixed__ return self.constraints[__fixed__].size == self.size - #def round(self, decimals=0, out=None): - # view = super(Param, self).round(decimals, out).view(Param) - # view.__array_finalize__(self) - # return view - #round.__doc__ = numpy.round.__doc__ + def _get_original(self, param): return self + + #=========================================================================== + # Pickling and copying + #=========================================================================== + def __deepcopy__(self, memo): + s = self.__new__(self.__class__, name=self.name, input_array=self.view(numpy.ndarray).copy()) + memo[id(self)] = s + import copy + s.__dict__.update(copy.deepcopy(self.__dict__, memo)) + return s + + #=========================================================================== # Printing -> done #=========================================================================== @@ -250,7 +195,8 @@ class Param(OptimizationHandlable, ObsAr): if self.size <= 1: return [str(self.view(numpy.ndarray)[0])] else: return [str(self.shape)] - def parameter_names(self, add_self=False, adjust_for_printing=False): + def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): + # this is just overwrighting the parameterized calls to parameter names, in order to maintain OOP if adjust_for_printing: return [adjust_name_for_printing(self.name)] return [self.name] @@ -261,6 +207,9 @@ class Param(OptimizationHandlable, ObsAr): def parameter_shapes(self): return [self.shape] @property + def num_params(self): + return 0 + @property 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 @@ -368,7 +317,7 @@ class ParamConcatenation(object): #=========================================================================== def __getitem__(self, s): ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; - params = [p._param_array_[ind[ps]] for p,ps in zip(self.params, self._param_slices_) if numpy.any(p._param_array_[ind[ps]])] + params = [p.param_array[ind[ps]] for p,ps in zip(self.params, self._param_slices_) if numpy.any(p.param_array[ind[ps]])] if len(params)==1: return params[0] return ParamConcatenation(params) def __setitem__(self, s, val, update=True): @@ -381,7 +330,7 @@ class ParamConcatenation(object): if update: self.update_all_params() def values(self): - return numpy.hstack([p._param_array_ for p in self.params]) + return numpy.hstack([p.param_array for p in self.params]) #=========================================================================== # parameter operations: #=========================================================================== diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index b804a61a..a60b8b38 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -13,10 +13,10 @@ Observable Pattern for patameterization """ -from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED +from transformations import Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED import numpy as np -__updated__ = '2014-03-24' +__updated__ = '2014-03-31' class HierarchyError(Exception): """ @@ -31,71 +31,8 @@ def adjust_name_for_printing(name): return name.replace(" ", "_").replace(".", "_").replace("-", "_m_").replace("+", "_p_").replace("!", "_I_").replace("**", "_xx_").replace("*", "_x_").replace("/", "_l_").replace("@", '_at_') return '' -class InterfacePickleFunctions(object): - def __init__(self, *a, **kw): - super(InterfacePickleFunctions, self).__init__() - def _getstate(self): - """ - Returns the state of this class in a memento pattern. - The state must be a list-like structure of all the fields - this class needs to run. - - See python doc "pickling" (`__getstate__` and `__setstate__`) for details. - """ - raise NotImplementedError, "To be able to use pickling you need to implement this method" - def _setstate(self, state): - """ - Set the state (memento pattern) of this class to the given state. - Usually this is just the counterpart to _getstate, such that - an object is a copy of another when calling - - copy = .__new__(*args,**kw)._setstate(._getstate()) - - See python doc "pickling" (`__getstate__` and `__setstate__`) for details. - """ - raise NotImplementedError, "To be able to use pickling you need to implement this method" - -class Pickleable(InterfacePickleFunctions): - """ - 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 - """ - def __init__(self, *a, **kw): - super(Pickleable, self).__init__() - #=========================================================================== - # 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__) - -class Observable(Pickleable): +class Observable(object): """ Observable pattern for parameterization. @@ -105,8 +42,9 @@ class Observable(Pickleable): """ _updated = True def __init__(self, *args, **kwargs): - super(Observable, self).__init__(*args, **kwargs) - self._observer_callables_ = [] + super(Observable, self).__init__() + from lists_and_dicts import ObservablesList + self._observer_callables_ = ObservablesList() def add_observer(self, observer, callble, priority=0): self._insert_sorted(priority, observer, callble) @@ -151,17 +89,11 @@ class Observable(Pickleable): ins += 1 self._observer_callables_.insert(ins, (p, o, c)) - def _getstate(self): - return [self._observer_callables_] - - def _setstate(self, state): - self._observer_callables_ = state.pop() - #=============================================================================== # Foundation framework for parameterized and param objects: #=============================================================================== -class Parentable(Observable): +class Parentable(object): """ Enable an Object to have a parent. @@ -171,7 +103,7 @@ class Parentable(Observable): _parent_ = None _parent_index_ = None def __init__(self, *args, **kwargs): - super(Parentable, self).__init__(*args, **kwargs) + super(Parentable, self).__init__() def has_parent(self): """ @@ -207,7 +139,84 @@ class Parentable(Observable): """ pass -class Gradcheckable(Parentable): +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 + """ + def __init__(self, *a, **kw): + super(Pickleable, self).__init__() + #=========================================================================== + # 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 as pickle + import pickle #TODO: cPickle + if isinstance(f, str): + with open(f, 'w') as f: + pickle.dump(self, f, protocol) + else: + pickle.dump(self, f, protocol) + + #=========================================================================== + # copy and pickling + #=========================================================================== + def copy(self): + """Returns a (deep) copy of the current model""" + #raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy" + import copy + memo = {} + memo[id(self._parent_)] = None + memo[id(self._parent_index_)] = None + memo[id(self.gradient)] = None + memo[id(self.param_array)] = None + memo[id(self._fixes_)] = None + c = copy.deepcopy(self, memo) + return c + + def __deepcopy__(self, memo): + s = self.__new__(self.__class__) + memo[id(self)] = s + import copy + s.__dict__.update(copy.deepcopy(self.__dict__, memo)) + return s + + def __getstate__(self): + ignore_list = ([#'_parent_', '_parent_index_', + #'_observer_callables_', + '_param_array_', '_gradient_array_', '_fixes_', + '_Cacher_wrap__cachers'] + #+ self.parameter_names(recursive=False) + ) + dc = dict() + for k,v in self.__dict__.iteritems(): + if k not in ignore_list: + #if hasattr(v, "__getstate__"): + #dc[k] = v.__getstate__() + #else: + dc[k] = v + return dc + + def __setstate__(self, state): + self.__dict__.update(state) + return self + + #def __getstate__(self, memo): + # raise NotImplementedError, "get state must be implemented to be able to pickle objects" + + #def __setstate__(self, memo): + # raise NotImplementedError, "set state must be implemented to be able to pickle objects" + +class Gradcheckable(Pickleable, 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. @@ -312,7 +321,7 @@ class Indexable(object): raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization object?" -class Constrainable(Nameable, Indexable): +class Constrainable(Nameable, Indexable, Observable): """ Make an object constrainable with Priors and Transformations. TODO: Mappings!! @@ -429,14 +438,14 @@ class Constrainable(Nameable, Indexable): def log_prior(self): """evaluate the prior""" if self.priors.size > 0: - x = self._param_array_ + x = self.param_array 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): """evaluate the gradients of the priors""" if self.priors.size > 0: - x = self._param_array_ + x = self.param_array ret = np.zeros(x.size) [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.iteritems()] return ret @@ -455,7 +464,7 @@ class Constrainable(Nameable, Indexable): Constrain the parameter to the given :py:class:`GPy.core.transformations.Transformation`. """ - self._param_array_[:] = transform.initialize(self._param_array_) + self.param_array[:] = transform.initialize(self.param_array) reconstrained = self.unconstrain() self._add_to_index_operations(self.constraints, reconstrained, transform, warning) self.notify_observers(self, None if trigger_parent else -np.inf) @@ -565,14 +574,14 @@ class OptimizationHandlable(Constrainable): super(OptimizationHandlable, self).__init__(name, default_constraint=default_constraint, *a, **kw) def transform(self): - [np.put(self._param_array_, ind, c.finv(self._param_array_.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + [np.put(self.param_array, ind, c.finv(self.param_array.flat[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_.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] def _get_params_transformed(self): # transformed parameters (apply transformation rules) - p = self._param_array_.copy() + 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_parent() and self.constraints[__fixed__].size != 0: fixes = np.ones(self.size).astype(bool) @@ -583,14 +592,14 @@ class OptimizationHandlable(Constrainable): return p def _set_params_transformed(self, p): - if p is self._param_array_: + if p is self.param_array: p = p.copy() if self.has_parent() and self.constraints[__fixed__].size != 0: fixes = np.ones(self.size).astype(bool) fixes[self.constraints[__fixed__]] = FIXED - self._param_array_.flat[fixes] = p - elif self._has_fixes(): self._param_array_.flat[self._fixes_] = p - else: self._param_array_.flat = p + self.param_array.flat[fixes] = p + elif self._has_fixes(): self.param_array.flat[self._fixes_] = p + else: self.param_array.flat = p self.untransform() self._trigger_params_changed() @@ -600,36 +609,29 @@ class OptimizationHandlable(Constrainable): def _size_transformed(self): return self.size - self.constraints[__fixed__].size -# -# 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! - # raise NotImplementedError, "Abstract superclass: This needs to be implemented in Param and Parameterizable" + @property + def num_params(self): + """ + Return the number of parameters of this parameter_handle. + Param objects will allways return 0. + """ + raise NotImplemented, "Abstract, please implement in respective classes" - #=========================================================================== - # Optimization handles: - #=========================================================================== + 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)] + else: names = [adjust(x.name) for x in self._parameters_] + if add_self: names = map(lambda x: adjust(self.name) + "." + x, names) + return names 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 @@ -663,16 +665,30 @@ class OptimizationHandlable(Constrainable): # For shared memory arrays. This does nothing in Param, but sets the memory # for all parameterized objects #=========================================================================== + @property + def full_gradient(self): + """ + Note to users: + This does not return the gradient in the right shape! Use self.gradient + for the right gradient array. + + To work on the gradient array, use this as the gradient handle. + This method exists for in memory use of parameters. + When trying to access the true gradient array, use this. + """ + self.gradient # <<< ensure _gradient_array_ + return self._gradient_array_ + def _propagate_param_grad(self, parray, garray): pi_old_size = 0 for pi in self._parameters_: pislice = slice(pi_old_size, pi_old_size + pi.size) - self._param_array_[pislice] = pi._param_array_.flat # , requirements=['C', 'W']).flat - self._gradient_array_[pislice] = pi._gradient_array_.flat # , requirements=['C', 'W']).flat + self.param_array[pislice] = pi.param_array.flat # , requirements=['C', 'W']).flat + self.full_gradient[pislice] = pi.full_gradient.flat # , requirements=['C', 'W']).flat - pi._param_array_.data = parray[pislice].data - pi._gradient_array_.data = garray[pislice].data + pi.param_array.data = parray[pislice].data + pi.full_gradient.data = garray[pislice].data pi._propagate_param_grad(parray[pislice], garray[pislice]) pi_old_size += pi.size @@ -681,26 +697,32 @@ class Parameterizable(OptimizationHandlable): def __init__(self, *args, **kwargs): super(Parameterizable, self).__init__(*args, **kwargs) from GPy.core.parameterization.lists_and_dicts import ArrayList - _parameters_ = ArrayList() + self._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): - """ - Get the names of all parameters of this model. + @property + def param_array(self): + if not hasattr(self, '_param_array_'): + self._param_array_ = np.empty(self.size, dtype=np.float64) + return self._param_array_ - :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)] - else: names = [adjust(x.name) for x in self._parameters_] - if add_self: names = map(lambda x: adjust(self.name) + "." + x, names) - return names + @param_array.setter + def param_array(self, arr): + self._param_array_ = arr + + #========================================================================= + # Gradient handling + #========================================================================= + @property + def gradient(self): + if not hasattr(self, '_gradient_array_'): + self._gradient_array_ = np.empty(self.size, dtype=np.float64) + return self._gradient_array_ + + @gradient.setter + def gradient(self, val): + self._gradient_array_[:] = val @property def num_params(self): @@ -737,34 +759,6 @@ class Parameterizable(OptimizationHandlable): self._remove_parameter_name(None, old_name) self._add_parameter_name(param) - #========================================================================= - # 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): - # [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): - # [p._set_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - #=========================================================================== - def add_parameter(self, param, index=None, _ignore_added_names=False): """ :param parameters: the parameters to add @@ -864,7 +858,7 @@ class Parameterizable(OptimizationHandlable): # no parameters for this class return old_size = 0 - self._param_array_ = np.empty(self.size, dtype=np.float64) + self.param_array = np.empty(self.size, dtype=np.float64) self._gradient_array_ = np.empty(self.size, dtype=np.float64) self._param_slices_ = [] @@ -874,15 +868,16 @@ class Parameterizable(OptimizationHandlable): pslice = slice(old_size, old_size + p.size) # first connect all children - p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice]) + p._propagate_param_grad(self.param_array[pslice], self.full_gradient[pslice]) # then connect children to self - self._param_array_[pslice] = p._param_array_.flat # , requirements=['C', 'W']).ravel(order='C') - self._gradient_array_[pslice] = p._gradient_array_.flat # , requirements=['C', 'W']).ravel(order='C') + self.param_array[pslice] = p.param_array.flat # , requirements=['C', 'W']).ravel(order='C') + self.full_gradient[pslice] = p.full_gradient.flat # , requirements=['C', 'W']).ravel(order='C') - if not p._param_array_.flags['C_CONTIGUOUS']: + if not p.param_array.flags['C_CONTIGUOUS']: + raise ValueError, "This should not happen! Please write an email to the developers with the code, which reproduces this error. All parameter arrays must be C_CONTIGUOUS" import ipdb;ipdb.set_trace() - p._param_array_.data = self._param_array_[pslice].data - p._gradient_array_.data = self._gradient_array_[pslice].data + p.param_array.data = self.param_array[pslice].data + p.full_gradient.data = self.full_gradient[pslice].data self._param_slices_.append(pslice) @@ -898,46 +893,22 @@ class Parameterizable(OptimizationHandlable): self.notify_observers(which=which) #=========================================================================== - # TODO: not working yet + # Pickling #=========================================================================== + def __setstate__(self, state): + super(Parameterizable, self).__setstate__(state) + self._connect_parameters() + self._connect_fixes() + self._notify_parent_change() + + self.parameters_changed() + def copy(self): - """Returns a (deep) copy of the current model""" - #raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy" - import copy - from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView - from .lists_and_dicts import ArrayList - - param_mapping = [[] for _ in range(self.num_params)] - - dc = dict() - for k, v in self.__dict__.iteritems(): - if k not in ['_parent_', '_parameters_', '_parent_index_', '_observer_callables_'] + self.parameter_names(recursive=False): - if v in self._parameters_: - param_mapping[self._parameters_.index(v)] += [k] - elif isinstance(v, (Constrainable, ParameterIndexOperations, ParameterIndexOperationsView)): - dc[k] = v.copy() - else: - dc[k] = copy.deepcopy(v) - if k == '_parameters_': - params = [p.copy() for p in v] - - dc['_parent_'] = None - dc['_parent_index_'] = None - dc['_observer_callables_'] = [] - dc['_parameters_'] = ArrayList() - dc['constraints'].clear() - dc['priors'].clear() - dc['size'] = 0 - - s = self.__new__(self.__class__) - s.__dict__ = dc - - for p, mlist in zip(params, param_mapping): - s.add_parameter(p, _ignore_added_names=True) - for m in mlist: - setattr(s, m, p) - return s - + c = super(Parameterizable, self).copy() + c._connect_parameters() + c._connect_fixes() + c._notify_parent_change() + return c #=========================================================================== # From being parentable, we have to define the parent_change notification #=========================================================================== diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index 529d3733..0760f8c6 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -101,34 +101,13 @@ class Parameterized(Parameterizable, Pickleable): return G return node - def _getstate(self): - """ - Get the current state of the class, - here just all the indices, rest can get recomputed - For inheriting from Parameterized: - - Allways append the state of the inherited object - and call down to the inherited object in _setstate!! - """ - return [] - - def _setstate(self, state): - self.parameters_changed() - #=========================================================================== - # Override copy to handle programmatically added observers - #=========================================================================== - def copy(self): - c = super(Parameterized, 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 - [numpy.put(g, i, g[i] * c.gradfactor(self._param_array_[i])) for c, i in self.constraints.iteritems() if c != __fixed__] + [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 @@ -160,7 +139,7 @@ class Parameterized(Parameterizable, Pickleable): this is not in the global view of things! """ return numpy.r_[:self.size] - + #=========================================================================== # Convenience for fixed, tied checking of param: #=========================================================================== @@ -175,7 +154,7 @@ class Parameterized(Parameterizable, Pickleable): # you can retrieve the original param through this method, by passing # the copy here return self._parameters_[param._parent_index_] - + #=========================================================================== # Get/set parameters: #=========================================================================== @@ -192,7 +171,7 @@ class Parameterized(Parameterizable, Pickleable): def __getitem__(self, name, paramlist=None): if isinstance(name, (int, slice, tuple, np.ndarray)): - return self._param_array_[name] + return self.param_array[name] else: if paramlist is None: paramlist = self.grep_param_names(name) @@ -208,7 +187,7 @@ class Parameterized(Parameterizable, Pickleable): def __setitem__(self, name, value, paramlist=None): if isinstance(name, (slice, tuple, np.ndarray)): try: - self._param_array_[name] = value + self.param_array[name] = value except: raise ValueError, "Setting by slice or index only allowed with array-like" self._trigger_params_changed() diff --git a/GPy/core/parameterization/variational.py b/GPy/core/parameterization/variational.py index ce39e2c9..f8fd165f 100644 --- a/GPy/core/parameterization/variational.py +++ b/GPy/core/parameterization/variational.py @@ -61,7 +61,7 @@ class SpikeAndSlabPrior(VariationalPrior): self.pi.gradient = (gamma/self.pi - (1.-gamma)/(1.-self.pi)).sum(axis=0) class VariationalPosterior(Parameterized): - def __init__(self, means=None, variances=None, name=None, *a, **kw): + def __init__(self, means=None, variances=None, name='latent space', *a, **kw): super(VariationalPosterior, self).__init__(name=name, *a, **kw) self.mean = Param("mean", means) self.variance = Param("variance", variances, Logexp()) @@ -119,6 +119,7 @@ class NormalPosterior(VariationalPosterior): import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." from ...plotting.matplot_dep import variational_plots + import matplotlib return variational_plots.plot(self,*args) class SpikeAndSlabPosterior(VariationalPosterior): diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 7bf0ca2a..7552b8ac 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -106,15 +106,3 @@ class SparseGP(GP): return mu, var - def _getstate(self): - """ - Get the current state of the class, - """ - return GP._getstate(self) + [ - self.Z, - self.num_inducing] - - def _setstate(self, state): - self.num_inducing = state.pop() - self.Z = state.pop() - GP._setstate(self, state) diff --git a/GPy/core/svigp.py b/GPy/core/svigp.py index a2c7acee..60e8371c 100644 --- a/GPy/core/svigp.py +++ b/GPy/core/svigp.py @@ -89,57 +89,57 @@ class SVIGP(GP): self._param_steplength_trace = [] self._vb_steplength_trace = [] - def _getstate(self): - steplength_params = [self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength] - return GP._getstate(self) + \ - [self.get_vb_param(), - self.Z, - self.num_inducing, - self.has_uncertain_inputs, - self.X_variance, - self.X_batch, - self.X_variance_batch, - steplength_params, - self.batchcounter, - self.batchsize, - self.epochs, - self.momentum, - self.data_prop, - self._param_trace, - self._param_steplength_trace, - self._vb_steplength_trace, - self._ll_trace, - self._grad_trace, - self.Y, - self._permutation, - self.iterations - ] - - def _setstate(self, state): - self.iterations = state.pop() - self._permutation = state.pop() - self.Y = state.pop() - self._grad_trace = state.pop() - self._ll_trace = state.pop() - self._vb_steplength_trace = state.pop() - self._param_steplength_trace = state.pop() - self._param_trace = state.pop() - self.data_prop = state.pop() - self.momentum = state.pop() - self.epochs = state.pop() - self.batchsize = state.pop() - self.batchcounter = state.pop() - steplength_params = state.pop() - (self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength) = steplength_params - self.X_variance_batch = state.pop() - self.X_batch = state.pop() - self.X_variance = state.pop() - self.has_uncertain_inputs = state.pop() - self.num_inducing = state.pop() - self.Z = state.pop() - vb_param = state.pop() - GP._setstate(self, state) - self.set_vb_param(vb_param) +# def _getstate(self): +# steplength_params = [self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength] +# return GP._getstate(self) + \ +# [self.get_vb_param(), +# self.Z, +# self.num_inducing, +# self.has_uncertain_inputs, +# self.X_variance, +# self.X_batch, +# self.X_variance_batch, +# steplength_params, +# self.batchcounter, +# self.batchsize, +# self.epochs, +# self.momentum, +# self.data_prop, +# self._param_trace, +# self._param_steplength_trace, +# self._vb_steplength_trace, +# self._ll_trace, +# self._grad_trace, +# self.Y, +# self._permutation, +# self.iterations +# ] +# +# def _setstate(self, state): +# self.iterations = state.pop() +# self._permutation = state.pop() +# self.Y = state.pop() +# self._grad_trace = state.pop() +# self._ll_trace = state.pop() +# self._vb_steplength_trace = state.pop() +# self._param_steplength_trace = state.pop() +# self._param_trace = state.pop() +# self.data_prop = state.pop() +# self.momentum = state.pop() +# self.epochs = state.pop() +# self.batchsize = state.pop() +# self.batchcounter = state.pop() +# steplength_params = state.pop() +# (self.hbar_t, self.tau_t, self.gbar_t, self.gbar_t1, self.gbar_t2, self.hbar_tp, self.tau_tp, self.gbar_tp, self.adapt_param_steplength, self.adapt_vb_steplength, self.vb_steplength, self.param_steplength) = steplength_params +# self.X_variance_batch = state.pop() +# self.X_batch = state.pop() +# self.X_variance = state.pop() +# self.has_uncertain_inputs = state.pop() +# self.num_inducing = state.pop() +# self.Z = state.pop() +# vb_param = state.pop() +# GP._setstate(self, state) +# self.set_vb_param(vb_param) def _compute_kernel_matrices(self): # kernel computations, using BGPLVM notation diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 8171a032..07623d6b 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -324,18 +324,15 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, def mrd_simulation(optimize=True, verbose=True, plot=True, plot_sim=True, **kw): from GPy import kern from GPy.models import MRD - from GPy.likelihoods import Gaussian D1, D2, D3, N, num_inducing, Q = 60, 20, 36, 60, 6, 5 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) #Ylist = [Ylist[0]] - k = [kern.Linear(Q, ARD=True) for _ in range(len(Ylist))] - m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="", initz='permute', **kw) + k = kern.Linear(Q, ARD=True) + m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="PCA_concat", initz='permute', **kw) - m['.*noise'] = [Y.var()/500. for Y in Ylist] - #for i, Y in enumerate(Ylist): - # m['.*Y_{}.*Gaussian.*noise'.format(i)] = Y.var(1) / 500. + m['.*noise'] = [Y.var()/40. for Y in Ylist] if optimize: print "Optimizing Model:" diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index 53f12722..0e10a175 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -23,6 +23,7 @@ class VarDTC(object): def __init__(self, limit=1): #self._YYTfactor_cache = caching.cache() from ...util.caching import Cacher + self.limit = limit self.get_trYYT = Cacher(self._get_trYYT, limit) self.get_YYTfactor = Cacher(self._get_YYTfactor, limit) @@ -33,6 +34,15 @@ class VarDTC(object): def _get_trYYT(self, Y): return param_to_array(np.sum(np.square(Y))) + def __getstate__(self): + return self.limit + + def __setstate__(self, state): + self.limit = state + from ...util.caching import Cacher + self.get_trYYT = Cacher(self._get_trYYT, self.limit) + self.get_YYTfactor = Cacher(self._get_YYTfactor, self.limit) + def _get_YYTfactor(self, Y): """ find a matrix L which satisfies LLT = YYT. @@ -126,7 +136,7 @@ class VarDTC(object): delit += output_dim * np.eye(num_inducing) # Compute dL_dKmm dL_dKmm = backsub_both_sides(Lm, delit) - + # derivatives of L w.r.t. psi dL_dpsi0, dL_dpsi1, dL_dpsi2 = _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 488745c5..de99bddb 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -15,7 +15,6 @@ class Kern(Parameterized): # found in kernel_slice_operations __metaclass__ = KernCallsViaSlicerMeta #=========================================================================== - _debug=False def __init__(self, input_dim, active_dims, name, *a, **kw): """ The base class for a kernel: a positive definite function @@ -175,22 +174,6 @@ class Kern(Parameterized): #else: kernels.append(other) return Prod([self, other], name) - def _getstate(self): - """ - Get the current state of the class, - here just all the indices, rest can get recomputed - """ - return super(Kern, self)._getstate() + [ - self.active_dims, - self.input_dim, - self._sliced_X] - - def _setstate(self, state): - self._sliced_X = state.pop() - self.input_dim = state.pop() - self.active_dims = state.pop() - super(Kern, self)._setstate(state) - class CombinationKernel(Kern): """ Abstract super class for combination kernels. @@ -220,7 +203,7 @@ class CombinationKernel(Kern): def get_input_dim_active_dims(self, kernels, extra_dims = None): active_dims = reduce(np.union1d, (np.r_[x.active_dims] for x in kernels), np.array([], dtype=int)) - input_dim = active_dims.max()+1 + (len(extra_dims) if extra_dims is not None else 0) + input_dim = active_dims.max()+1 + (len(np.r_[extra_dims]) if extra_dims is not None else 0) active_dims = slice(0, input_dim, 1) return input_dim, active_dims diff --git a/GPy/kern/_src/kernel_slice_operations.py b/GPy/kern/_src/kernel_slice_operations.py index ea5d2b0a..a4bb8f62 100644 --- a/GPy/kern/_src/kernel_slice_operations.py +++ b/GPy/kern/_src/kernel_slice_operations.py @@ -36,7 +36,7 @@ class _Slice_wrap(object): if self.k._sliced_X == 0: assert X.shape[1] > max(np.r_[self.k.active_dims]), "At least {} dimensional X needed".format(max(np.r_[self.k.active_dims])) self.X = self.k._slice_X(X) - self.X2 = self.k._slice_X(X2) if X2 is not None else None + self.X2 = self.k._slice_X(X2) if X2 is not None else X2 self.ret = True else: self.X = X diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py index ef3462f6..d623c8f1 100644 --- a/GPy/models/bayesian_gplvm.py +++ b/GPy/models/bayesian_gplvm.py @@ -2,7 +2,6 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) import numpy as np -from gplvm import GPLVM from .. import kern from ..core import SparseGP from ..likelihoods import Gaussian @@ -61,18 +60,6 @@ class BayesianGPLVM(SparseGP): SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, name, **kwargs) self.add_parameter(self.X, index=0) - def _getstate(self): - """ - Get the current state of the class, - here just all the indices, rest can get recomputed - """ - return SparseGP._getstate(self) + [self.init] - - def _setstate(self, state): - self._const_jitter = None - self.init = state.pop() - SparseGP._setstate(self, state) - def set_X_gradients(self, X, X_grad): """Set the gradients of the posterior distribution of X in its specific form.""" X.mean.gradient, X.variance.gradient = X_grad diff --git a/GPy/models/gp_regression.py b/GPy/models/gp_regression.py index 86e64a54..d56e72b9 100644 --- a/GPy/models/gp_regression.py +++ b/GPy/models/gp_regression.py @@ -29,8 +29,3 @@ class GPRegression(GP): super(GPRegression, self).__init__(X, Y, kernel, likelihood, name='GP regression', Y_metadata=Y_metadata) - def _getstate(self): - return GP._getstate(self) - - def _setstate(self, state): - return GP._setstate(self, state) diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py index fb7d93e7..542dcd31 100644 --- a/GPy/models/gplvm.py +++ b/GPy/models/gplvm.py @@ -44,12 +44,6 @@ class GPLVM(GP): super(GPLVM, self).parameters_changed() self.X.gradient = self.kern.gradients_X(self.grad_dict['dL_dK'], self.X, None) - def _getstate(self): - return GP._getstate(self) - - def _setstate(self, state): - GP._setstate(self, state) - def jacobian(self,X): target = np.zeros((X.shape[0],X.shape[1],self.output_dim)) for i in range(self.output_dim): diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index 36088e35..458a70a1 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -65,14 +65,17 @@ class MRD(Model): from ..kern import RBF self.kern = [RBF(input_dim, ARD=1, lengthscale=fracs[i], name='rbf'.format(i)) for i in range(len(Ylist))] elif isinstance(kernel, Kern): - self.kern = [kernel.copy(name='{}'.format(kernel.name, i)) for i in range(len(Ylist))] + self.kern = [] + for i in range(len(Ylist)): + k = kernel.copy() + self.kern.append(k) 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.kern = kernel if X_variance is None: - X_variance = np.random.uniform(0, .1, X.shape) + X_variance = np.random.uniform(0.1, 0.2, X.shape) self.variational_prior = NormalPrior() self.X = NormalPosterior(X, X_variance) @@ -108,8 +111,8 @@ class MRD(Model): def parameters_changed(self): self._log_marginal_likelihood = 0 self.posteriors = [] - self.Z.gradient = 0. - self.X.gradient = 0. + self.Z.gradient[:] = 0. + self.X.gradient[:] = 0. for y, k, l, i in itertools.izip(self.Ylist, self.kern, self.likelihood, self.inference_method): posterior, lml, grad_dict = i.inference(k, self.X, self.Z, l, y) @@ -160,6 +163,8 @@ class MRD(Model): X = np.random.randn(Ylist[0].shape[0], self.input_dim) fracs = X.var(0) fracs = [fracs]*self.input_dim + X -= X.mean() + X /= X.std() return X, fracs def _init_Z(self, init="permute", X=None): diff --git a/GPy/models/sparse_gp_classification.py b/GPy/models/sparse_gp_classification.py index 96f7ac5a..e2c77d95 100644 --- a/GPy/models/sparse_gp_classification.py +++ b/GPy/models/sparse_gp_classification.py @@ -46,11 +46,3 @@ class SparseGPClassification(SparseGP): SparseGP.__init__(self, X, likelihood, kernel, Z=Z, normalize_X=normalize_X) self.ensure_default_constraints() - def _getstate(self): - return SparseGP._getstate(self) - - - def _setstate(self, state): - return SparseGP._setstate(self, state) - - pass diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py index 7edb93e4..f4d5513e 100644 --- a/GPy/models/sparse_gp_regression.py +++ b/GPy/models/sparse_gp_regression.py @@ -51,14 +51,6 @@ class SparseGPRegression(SparseGP): SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC()) - def _getstate(self): - return SparseGP._getstate(self) - - def _setstate(self, state): - return SparseGP._setstate(self, state) - - - class SparseGPRegressionUncertainInput(SparseGP): """ Gaussian Process model for regression with Gaussian variance on the inputs (X_variance) diff --git a/GPy/models/sparse_gplvm.py b/GPy/models/sparse_gplvm.py index 5c10d0b8..638da63e 100644 --- a/GPy/models/sparse_gplvm.py +++ b/GPy/models/sparse_gplvm.py @@ -28,14 +28,6 @@ class SparseGPLVM(SparseGPRegression, GPLVM): SparseGPRegression.__init__(self, X, Y, kernel=kernel, num_inducing=num_inducing) self.ensure_default_constraints() - def _getstate(self): - return SparseGPRegression._getstate(self) - - - def _setstate(self, state): - return SparseGPRegression._setstate(self, state) - - def _get_param_names(self): return (sum([['X_%i_%i' % (n, q) for q in range(self.input_dim)] for n in range(self.num_data)], []) + SparseGPRegression._get_param_names(self)) diff --git a/GPy/models/svigp_regression.py b/GPy/models/svigp_regression.py index 3faa1cab..3397e31e 100644 --- a/GPy/models/svigp_regression.py +++ b/GPy/models/svigp_regression.py @@ -43,10 +43,3 @@ class SVIGPRegression(SVIGP): SVIGP.__init__(self, X, likelihood, kernel, Z, q_u=q_u, batchsize=batchsize) self.load_batch() - def _getstate(self): - return GPBase._getstate(self) - - - def _setstate(self, state): - return GPBase._setstate(self, state) - diff --git a/GPy/models/warped_gp.py b/GPy/models/warped_gp.py index d78f31df..4b982ed2 100644 --- a/GPy/models/warped_gp.py +++ b/GPy/models/warped_gp.py @@ -30,14 +30,6 @@ class WarpedGP(GP): GP.__init__(self, X, likelihood, kernel, normalize_X=normalize_X) self._set_params(self._get_params()) - def _getstate(self): - return GP._getstate(self) - - - def _setstate(self, state): - return GP._setstate(self, state) - - def _scale_data(self, Y): self._Ymax = Y.max() self._Ymin = Y.min() diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 944a054f..bda64f8a 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -5,6 +5,7 @@ import unittest import numpy as np import GPy import sys +from GPy.core.parameterization.param import Param verbose = 0 @@ -30,7 +31,7 @@ class Kern_check_model(GPy.core.Model): dL_dK = np.ones((X.shape[0], X2.shape[0])) self.kernel = kernel - self.X = GPy.core.parameterization.Param('X',X) + self.X = X self.X2 = X2 self.dL_dK = dL_dK @@ -77,10 +78,11 @@ class Kern_check_dK_dX(Kern_check_model): """This class allows gradient checks for the gradient of a kernel with respect to X. """ def __init__(self, kernel=None, dL_dK=None, X=None, X2=None): Kern_check_model.__init__(self,kernel=kernel,dL_dK=dL_dK, X=X, X2=X2) + self.X = Param('X',X) self.add_parameter(self.X) def parameters_changed(self): - self.X.gradient = self.kernel.gradients_X(self.dL_dK, self.X, self.X2) + self.X.gradient[:] = self.kernel.gradients_X(self.dL_dK, self.X, self.X2) class Kern_check_dKdiag_dX(Kern_check_dK_dX): """This class allows gradient checks for the gradient of a kernel diagonal with respect to X. """ @@ -91,7 +93,7 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX): return (np.diag(self.dL_dK)*self.kernel.Kdiag(self.X)).sum() def parameters_changed(self): - self.X.gradient = self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X) + self.X.gradient[:] = self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X) @@ -127,6 +129,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb if not result: print("Positive definite check failed for " + kern.name + " covariance function.") pass_checks = False + assert(result) return False if verbose: @@ -138,6 +141,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb print("Gradient of K(X, X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:") Kern_check_dK_dtheta(kern, X=X, X2=None).checkgrad(verbose=True) pass_checks = False + assert(result) return False if verbose: @@ -149,6 +153,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb print("Gradient of K(X, X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:") Kern_check_dK_dtheta(kern, X=X, X2=X2).checkgrad(verbose=True) pass_checks = False + assert(result) return False if verbose: @@ -165,6 +170,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb print("Gradient of Kdiag(X) wrt theta failed for " + kern.name + " covariance function. Gradient values as follows:") Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=True) pass_checks = False + assert(result) return False if verbose: @@ -183,6 +189,8 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb if not result: print("Gradient of K(X, X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:") testmodel.checkgrad(verbose=True) + import ipdb;ipdb.set_trace() + assert(result) pass_checks = False return False @@ -202,6 +210,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb if not result: print("Gradient of K(X, X2) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:") testmodel.checkgrad(verbose=True) + assert(result) pass_checks = False return False @@ -219,6 +228,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb print("Gradient of Kdiag(X) wrt X failed for " + kern.name + " covariance function. Gradient values as follows:") Kern_check_dKdiag_dX(kern, X=X).checkgrad(verbose=True) pass_checks = False + assert(result) return False return pass_checks @@ -227,7 +237,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb class KernelGradientTestsContinuous(unittest.TestCase): def setUp(self): - self.N, self.D = 100, 5 + self.N, self.D = 10, 5 self.X = np.random.randn(self.N,self.D) self.X2 = np.random.randn(self.N+10,self.D) @@ -339,59 +349,59 @@ class KernelTestsMiscellaneous(unittest.TestCase): self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=[self.linear, self.rbf]), self.linear.K(self.X)+self.rbf.K(self.X))) self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=self.sumkern.parts[0]), self.rbf.K(self.X))) -class KernelTestsNonContinuous(unittest.TestCase): - def setUp(self): - N0 = 3 - N1 = 9 - N2 = 4 - N = N0+N1+N2 - self.D = 3 - self.X = np.random.randn(N, self.D+1) - indices = np.random.random_integers(0, 2, size=N) - self.X[indices==0, -1] = 0 - self.X[indices==1, -1] = 1 - self.X[indices==2, -1] = 2 - #self.X = self.X[self.X[:, -1].argsort(), :] - self.X2 = np.random.randn((N0+N1)*2, self.D+1) - self.X2[:(N0*2), -1] = 0 - self.X2[(N0*2):, -1] = 1 - - def test_IndependentOutputs(self): - k = GPy.kern.RBF(self.D) - kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single') - self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1)) - k = [GPy.kern.RBF(1, active_dims=[1], name='rbf1'), GPy.kern.RBF(self.D, name='rbf012'), GPy.kern.RBF(2, active_dims=[0,2], name='rbf02')] - kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split') - self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1)) - - def test_ODE_UY(self): - kern = GPy.kern.ODE_UY(2, active_dims=[0, self.D]) - X = self.X[self.X[:,-1]!=2] - X2 = self.X2[self.X2[:,-1]!=2] - self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) +# class KernelTestsNonContinuous(unittest.TestCase): +# def setUp(self): +# N0 = 3 +# N1 = 9 +# N2 = 4 +# N = N0+N1+N2 +# self.D = 3 +# self.X = np.random.randn(N, self.D+1) +# indices = np.random.random_integers(0, 2, size=N) +# self.X[indices==0, -1] = 0 +# self.X[indices==1, -1] = 1 +# self.X[indices==2, -1] = 2 +# #self.X = self.X[self.X[:, -1].argsort(), :] +# self.X2 = np.random.randn((N0+N1)*2, self.D+1) +# self.X2[:(N0*2), -1] = 0 +# self.X2[(N0*2):, -1] = 1 +# +# def test_IndependentOutputs(self): +# k = GPy.kern.RBF(self.D) +# kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single') +# self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1)) +# k = [GPy.kern.RBF(1, active_dims=[1], name='rbf1'), GPy.kern.RBF(self.D, name='rbf012'), GPy.kern.RBF(2, active_dims=[0,2], name='rbf02')] +# kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split') +# self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1)) +# +# def test_ODE_UY(self): +# kern = GPy.kern.ODE_UY(2, active_dims=[0, self.D]) +# X = self.X[self.X[:,-1]!=2] +# X2 = self.X2[self.X2[:,-1]!=2] +# self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) if __name__ == "__main__": print "Running unit tests, please be (very) patient..." - #unittest.main() - np.random.seed(0) - N0 = 3 - N1 = 9 - N2 = 4 - N = N0+N1+N2 - D = 3 - X = np.random.randn(N, D+1) - indices = np.random.random_integers(0, 2, size=N) - X[indices==0, -1] = 0 - X[indices==1, -1] = 1 - X[indices==2, -1] = 2 - #X = X[X[:, -1].argsort(), :] - X2 = np.random.randn((N0+N1)*2, D+1) - X2[:(N0*2), -1] = 0 - X2[(N0*2):, -1] = 1 - k = [GPy.kern.RBF(1, active_dims=[1], name='rbf1'), GPy.kern.RBF(D, name='rbf012'), GPy.kern.RBF(2, active_dims=[0,2], name='rbf02')] - kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split') - assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) - k = GPy.kern.RBF(D) - kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single') - assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) + unittest.main() +# np.random.seed(0) +# N0 = 3 +# N1 = 9 +# N2 = 4 +# N = N0+N1+N2 +# D = 3 +# X = np.random.randn(N, D+1) +# indices = np.random.random_integers(0, 2, size=N) +# X[indices==0, -1] = 0 +# X[indices==1, -1] = 1 +# X[indices==2, -1] = 2 +# #X = X[X[:, -1].argsort(), :] +# X2 = np.random.randn((N0+N1)*2, D+1) +# X2[:(N0*2), -1] = 0 +# X2[(N0*2):, -1] = 1 +# k = [GPy.kern.RBF(1, active_dims=[1], name='rbf1'), GPy.kern.RBF(D, name='rbf012'), GPy.kern.RBF(2, active_dims=[0,2], name='rbf02')] +# kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split') +# assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) +# k = GPy.kern.RBF(D) +# kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single') +# assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1)) diff --git a/GPy/testing/observable_tests.py b/GPy/testing/observable_tests.py index 90623703..05794dc3 100644 --- a/GPy/testing/observable_tests.py +++ b/GPy/testing/observable_tests.py @@ -93,12 +93,12 @@ class Test(unittest.TestCase): def test_set_params(self): self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') - self.par._param_array_[:] = 1 + 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.par._param_array_[:] = 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) diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index dc59449f..911cde0b 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -86,9 +86,9 @@ class ParameterizedTest(unittest.TestCase): self.assertListEqual(self.test1.constraints[Logexp()].tolist(), range(self.param.size, self.param.size+self.rbf.size)) def test_remove_parameter_param_array_grad_array(self): - val = self.test1.kern._param_array_.copy() + val = self.test1.kern.param_array.copy() self.test1.kern.remove_parameter(self.white) - self.assertListEqual(self.test1.kern._param_array_.tolist(), val[:2].tolist()) + self.assertListEqual(self.test1.kern.param_array.tolist(), val[:2].tolist()) def test_add_parameter_already_in_hirarchy(self): self.assertRaises(HierarchyError, self.test1.add_parameter, self.white._parameters_[0]) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 0886d0c6..ced56727 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -97,6 +97,15 @@ class Cacher(object): self.cached_outputs = [] self.inputs_changed = [] + def __deepcopy__(self, memo=None): + return Cacher(self.operation, self.limit, self.ignore_args, self.force_kwargs) + + def __getstate__(self, memo=None): + raise NotImplementedError, "Trying to pickle Cacher object with function {}, pickling functions not possible.".format(str(self.operation)) + + def __setstate__(self, memo=None): + raise NotImplementedError, "Trying to pickle Cacher object with function {}, pickling functions not possible.".format(str(self.operation)) + @property def __name__(self): return self.operation.__name__