pickling and caching

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)

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:

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)
+        func, args, state = super(ObsAr, self).__reduce__()
-        return func, args, (state, Observable._getstate(self))
+        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
+        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
-

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
 

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

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
+        self._gradient_array_[self._current_slice_] = 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_)
 
     #===========================================================================
     # 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:
     #===========================================================================

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):
+class Observable(object):
-        """
-        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 = <classname>.__new__(*args,**kw)._setstate(<to_be_copied>._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(<memento>) (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):
     """
     Observable pattern for parameterization.
 
@@ -105,8 +42,9 @@ class Observable(Pickleable):
     """
     _updated = True
     def __init__(self, *args, **kwargs):
-        super(Observable, self).__init__(*args, **kwargs)
+        super(Observable, self).__init__()
-        self._observer_callables_ = []
+        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(<memento>) (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
+            self.param_array.flat[fixes] = p
-        elif self._has_fixes(): self._param_array_.flat[self._fixes_] = p
+        elif self._has_fixes(): self.param_array.flat[self._fixes_] = p
-        else: self._param_array_.flat = 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):
+    @property
-#         self._param_array_.flat = params
+    def num_params(self):
-#         if trigger_parent: min_priority = None
+        """
-#         else: min_priority = -np.inf
+        Return the number of parameters of this parameter_handle.
-#         self.notify_observers(None, min_priority)
+        Param objects will allways return 0.
-        # don't overwrite this anymore!
+        """
-        # raise NotImplementedError, "Abstract superclass: This needs to be implemented in Param and Parameterizable"
+        raise NotImplemented, "Abstract, please implement in respective classes"
 
-    #===========================================================================
+    def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True):
-    # Optimization handles:
+        """
-    #===========================================================================
+        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.param_array[pislice] = pi.param_array.flat  # , requirements=['C', 'W']).flat
-            self._gradient_array_[pislice] = pi._gradient_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.param_array.data = parray[pislice].data
-            pi._gradient_array_.data = garray[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):
+    @property
-        """
+    def param_array(self):
-        Get the names of all parameters of this model.
+        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_array.setter
-        :param bool adjust_for_printing: whether to call `adjust_name_for_printing` on names
+    def param_array(self, arr):
-        :param bool recursive: whether to traverse through hierarchy and append leaf node names
+        self._param_array_ = arr
-        """
+
-        if adjust_for_printing: adjust = lambda x: adjust_name_for_printing(x)
+    #=========================================================================
-        else: adjust = lambda x: x
+    # Gradient handling
-        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_]
+    @property
-        if add_self: names = map(lambda x: adjust(self.name) + "." + x, names)
+    def gradient(self):
-        return names
+        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.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.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.param_array.data = self.param_array[pslice].data
-            p._gradient_array_.data = self._gradient_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"""
+        c = super(Parameterizable, self).copy()
-        #raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy"
+        c._connect_parameters()
-        import copy
+        c._connect_fixes()
-        from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView
+        c._notify_parent_change()
-        from .lists_and_dicts import ArrayList
+        return c
-        
-        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
-
     #===========================================================================
     # From being parentable, we have to define the parent_change notification
     #===========================================================================

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()

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):

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)

GPy/core/svigp.py

@@ -89,57 +89,57 @@ class SVIGP(GP):
         self._param_steplength_trace = []
         self._vb_steplength_trace = []
 
-    def _getstate(self):
+#     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]
+#         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) + \
+#         return GP._getstate(self) + \
-            [self.get_vb_param(),
+#             [self.get_vb_param(),
-             self.Z,
+#              self.Z,
-             self.num_inducing,
+#              self.num_inducing,
-             self.has_uncertain_inputs,
+#              self.has_uncertain_inputs,
-             self.X_variance,
+#              self.X_variance,
-             self.X_batch,
+#              self.X_batch,
-             self.X_variance_batch,
+#              self.X_variance_batch,
-             steplength_params,
+#              steplength_params,
-             self.batchcounter,
+#              self.batchcounter,
-             self.batchsize,
+#              self.batchsize,
-             self.epochs,
+#              self.epochs,
-             self.momentum,
+#              self.momentum,
-             self.data_prop,
+#              self.data_prop,
-             self._param_trace,
+#              self._param_trace,
-             self._param_steplength_trace,
+#              self._param_steplength_trace,
-             self._vb_steplength_trace,
+#              self._vb_steplength_trace,
-             self._ll_trace,
+#              self._ll_trace,
-             self._grad_trace,
+#              self._grad_trace,
-             self.Y,
+#              self.Y,
-             self._permutation,
+#              self._permutation,
-             self.iterations
+#              self.iterations
-            ]
+#             ]
-
+# 
-    def _setstate(self, state):
+#     def _setstate(self, state):
-        self.iterations = state.pop()
+#         self.iterations = state.pop()
-        self._permutation = state.pop()
+#         self._permutation = state.pop()
-        self.Y = state.pop()
+#         self.Y = state.pop()
-        self._grad_trace = state.pop()
+#         self._grad_trace = state.pop()
-        self._ll_trace = state.pop()
+#         self._ll_trace = state.pop()
-        self._vb_steplength_trace = state.pop()
+#         self._vb_steplength_trace = state.pop()
-        self._param_steplength_trace = state.pop()
+#         self._param_steplength_trace = state.pop()
-        self._param_trace = state.pop()
+#         self._param_trace = state.pop()
-        self.data_prop = state.pop()
+#         self.data_prop = state.pop()
-        self.momentum = state.pop()
+#         self.momentum = state.pop()
-        self.epochs = state.pop()
+#         self.epochs = state.pop()
-        self.batchsize = state.pop()
+#         self.batchsize = state.pop()
-        self.batchcounter = state.pop()
+#         self.batchcounter = state.pop()
-        steplength_params = 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.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_variance_batch = state.pop()
-        self.X_batch = state.pop()
+#         self.X_batch = state.pop()
-        self.X_variance = state.pop()
+#         self.X_variance = state.pop()
-        self.has_uncertain_inputs = state.pop()
+#         self.has_uncertain_inputs = state.pop()
-        self.num_inducing = state.pop()
+#         self.num_inducing = state.pop()
-        self.Z = state.pop()
+#         self.Z = state.pop()
-        vb_param = state.pop()
+#         vb_param = state.pop()
-        GP._setstate(self, state)
+#         GP._setstate(self, state)
-        self.set_vb_param(vb_param)
+#         self.set_vb_param(vb_param)
 
     def _compute_kernel_matrices(self):
         # kernel computations, using BGPLVM notation

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))]
+    k = kern.Linear(Q, ARD=True)
-    m = MRD(Ylist, input_dim=Q, num_inducing=num_inducing, kernel=k, initx="", initz='permute', **kw)
+    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]
+    m['.*noise'] = [Y.var()/40. for Y in Ylist]
-    #for i, Y in enumerate(Ylist):
-    #    m['.*Y_{}.*Gaussian.*noise'.format(i)] = Y.var(1) / 500.
 
     if optimize:
         print "Optimizing Model:"

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,

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
 

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

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

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)

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):

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.Z.gradient[:] = 0.
-        self.X.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):

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

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)

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))

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)
-

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()

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):
+# class KernelTestsNonContinuous(unittest.TestCase):
-    def setUp(self):
+#     def setUp(self):
-        N0 = 3
+#         N0 = 3
-        N1 = 9
+#         N1 = 9
-        N2 = 4
+#         N2 = 4
-        N = N0+N1+N2
+#         N = N0+N1+N2
-        self.D = 3
+#         self.D = 3
-        self.X = np.random.randn(N, self.D+1)
+#         self.X = np.random.randn(N, self.D+1)
-        indices = np.random.random_integers(0, 2, size=N)
+#         indices = np.random.random_integers(0, 2, size=N)
-        self.X[indices==0, -1] = 0
+#         self.X[indices==0, -1] = 0
-        self.X[indices==1, -1] = 1
+#         self.X[indices==1, -1] = 1
-        self.X[indices==2, -1] = 2
+#         self.X[indices==2, -1] = 2
-        #self.X = self.X[self.X[:, -1].argsort(), :]
+#         #self.X = self.X[self.X[:, -1].argsort(), :]
-        self.X2 = np.random.randn((N0+N1)*2, self.D+1)
+#         self.X2 = np.random.randn((N0+N1)*2, self.D+1)
-        self.X2[:(N0*2), -1] = 0
+#         self.X2[:(N0*2), -1] = 0
-        self.X2[(N0*2):, -1] = 1
+#         self.X2[(N0*2):, -1] = 1
-
+#  
-    def test_IndependentOutputs(self):
+#     def test_IndependentOutputs(self):
-        k = GPy.kern.RBF(self.D)
+#         k = GPy.kern.RBF(self.D)
-        kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single')
+#         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))
+#         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')]
+#         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')
+#         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))
+#         self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1))
-
+#  
-    def test_ODE_UY(self):
+#     def test_ODE_UY(self):
-        kern = GPy.kern.ODE_UY(2, active_dims=[0, self.D])
+#         kern = GPy.kern.ODE_UY(2, active_dims=[0, self.D])
-        X = self.X[self.X[:,-1]!=2]
+#         X = self.X[self.X[:,-1]!=2]
-        X2 = self.X2[self.X2[:,-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))
+#         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()
+    unittest.main()
-    np.random.seed(0)
+#     np.random.seed(0)
-    N0 = 3
+#     N0 = 3
-    N1 = 9
+#     N1 = 9
-    N2 = 4
+#     N2 = 4
-    N = N0+N1+N2
+#     N = N0+N1+N2
-    D = 3
+#     D = 3
-    X = np.random.randn(N, D+1)
+#     X = np.random.randn(N, D+1)
-    indices = np.random.random_integers(0, 2, size=N)
+#     indices = np.random.random_integers(0, 2, size=N)
-    X[indices==0, -1] = 0
+#     X[indices==0, -1] = 0
-    X[indices==1, -1] = 1
+#     X[indices==1, -1] = 1
-    X[indices==2, -1] = 2
+#     X[indices==2, -1] = 2
-    #X = X[X[:, -1].argsort(), :]
+#     #X = X[X[:, -1].argsort(), :]
-    X2 = np.random.randn((N0+N1)*2, D+1)
+#     X2 = np.random.randn((N0+N1)*2, D+1)
-    X2[:(N0*2), -1] = 0
+#     X2[:(N0*2), -1] = 0
-    X2[(N0*2):, -1] = 1
+#     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')]
+#     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')
+#     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))
+#     assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
-    k = GPy.kern.RBF(D)
+#     k = GPy.kern.RBF(D)
-    kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single')
+#     kern = GPy.kern.IndependentOutputs(k, -1, 'ind_single')
-    assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
+#     assert(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))

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)

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])

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__