merging

GPy/core/gp.py

@@ -33,13 +33,13 @@ class GP(Model):
 
         assert X.ndim == 2
         if isinstance(X, (ObsAr, VariationalPosterior)):
-            self.X = X
+            self.X = X.copy()
-        else: self.X = ObsAr(X)
+        else: self.X = ObsAr(X.copy())
 
         self.num_data, self.input_dim = self.X.shape
 
         assert Y.ndim == 2
-        self.Y = ObsAr(Y)
+        self.Y = ObsAr(Y.copy())
         assert Y.shape[0] == self.num_data
         _, self.output_dim = self.Y.shape
 
@@ -276,5 +276,9 @@ class GP(Model):
         TODO: valid args
         """
         self.inference_method.on_optimization_start()
-        super(GP, self).optimize(optimizer, start, **kwargs)
+        try:
-        self.inference_method.on_optimization_end()
+            super(GP, self).optimize(optimizer, start, **kwargs)
+        except KeyboardInterrupt:
+            print "KeyboardInterrupt caught, calling on_optimization_end() to round things up"
+            self.inference_method.on_optimization_end()
+            raise

GPy/core/model.py

@@ -20,7 +20,7 @@ class Model(Parameterized):
         super(Model, self).__init__(name)  # Parameterized.__init__(self)
         self.optimization_runs = []
         self.sampling_runs = []
-        self.preferred_optimizer = 'scg'
+        self.preferred_optimizer = 'bfgs'
 
     def log_likelihood(self):
         raise NotImplementedError, "this needs to be implemented to use the model class"
@@ -61,7 +61,7 @@ class Model(Parameterized):
         on the current machine.
 
         """
-        initial_parameters = self._get_params_transformed()
+        initial_parameters = self.optimizer_array
 
         if parallel:
             try:
@@ -124,13 +124,15 @@ class Model(Parameterized):
 
         For probabilistic models this is the negative log_likelihood
         (including the MAP prior), so we return it here. If your model is not 
-        probabilistic, just return your objective here!
+        probabilistic, just return your objective to minimize here!
         """
         return -float(self.log_likelihood()) - self.log_prior()
 
     def objective_function_gradients(self):
         """
         The gradients for the objective function for the given algorithm.
+        The gradients are w.r.t. the *negative* objective function, as 
+        this framework works with *negative* log-likelihoods as a default.
 
         You can find the gradient for the parameters in self.gradient at all times.
         This is the place, where gradients get stored for parameters.
@@ -141,7 +143,7 @@ class Model(Parameterized):
 
         For probabilistic models this is the gradient of the negative log_likelihood
         (including the MAP prior), so we return it here. If your model is not 
-        probabilistic, just return your gradient here!
+        probabilistic, just return your *negative* gradient here!
         """
         return -(self._log_likelihood_gradients() + self._log_prior_gradients())
 
@@ -157,7 +159,8 @@ class Model(Parameterized):
         :type x: np.array
         """
         try:
-            self._set_params_transformed(x)
+            # self._set_params_transformed(x)
+            self.optimizer_array = x
             obj_grads = self._transform_gradients(self.objective_function_gradients())
             self._fail_count = 0
         except (LinAlgError, ZeroDivisionError, ValueError):
@@ -180,7 +183,7 @@ class Model(Parameterized):
         :parameter type: np.array
         """
         try:
-            self._set_params_transformed(x)
+            self.optimizer_array = x
             obj = self.objective_function()
             self._fail_count = 0
         except (LinAlgError, ZeroDivisionError, ValueError):
@@ -192,7 +195,7 @@ class Model(Parameterized):
 
     def _objective_grads(self, x):
         try:
-            self._set_params_transformed(x)
+            self.optimizer_array = x
             obj_f, obj_grads = self.objective_function(), self._transform_gradients(self.objective_function_gradients())
             self._fail_count = 0
         except (LinAlgError, ZeroDivisionError, ValueError):
@@ -226,7 +229,7 @@ class Model(Parameterized):
             optimizer = self.preferred_optimizer
 
         if start == None:
-            start = self._get_params_transformed()
+            start = self.optimizer_array
 
         optimizer = optimization.get_optimizer(optimizer)
         opt = optimizer(start, model=self, **kwargs)
@@ -235,7 +238,7 @@ class Model(Parameterized):
 
         self.optimization_runs.append(opt)
 
-        self._set_params_transformed(opt.x_opt)
+        self.optimizer_array = opt.x_opt
 
     def optimize_SGD(self, momentum=0.1, learning_rate=0.01, iterations=20, **kwargs):
         # assert self.Y.shape[1] > 1, "SGD only works with D > 1"
@@ -260,7 +263,7 @@ class Model(Parameterized):
            The gradient is considered correct if the ratio of the analytical
            and numerical gradients is within <tolerance> of unity.
         """
-        x = self._get_params_transformed().copy()
+        x = self.optimizer_array.copy()
 
         if not verbose:
             # make sure only to test the selected parameters
@@ -270,8 +273,8 @@ class Model(Parameterized):
                 transformed_index = self._raveled_index_for(target_param)
                 if self._has_fixes():
                     indices = np.r_[:self.size]
-                    which = (transformed_index[:,None]==indices[self._fixes_][None,:]).nonzero()
+                    which = (transformed_index[:, None] == indices[self._fixes_][None, :]).nonzero()
-                    transformed_index = (indices-(~self._fixes_).cumsum())[transformed_index[which[0]]]
+                    transformed_index = (indices - (~self._fixes_).cumsum())[transformed_index[which[0]]]
 
                 if transformed_index.size == 0:
                     print "No free parameters to check"
@@ -290,7 +293,7 @@ class Model(Parameterized):
             gradient = gradient[transformed_index]
 
             denominator = (2 * np.dot(dx, gradient))
-            global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator)
+            global_ratio = (f1 - f2) / np.where(denominator == 0., 1e-32, denominator)
             global_diff = np.abs(f1 - f2) < tolerance and np.allclose(gradient, 0, atol=tolerance)
             if global_ratio is np.nan:
                 global_ratio = 0
@@ -319,10 +322,10 @@ class Model(Parameterized):
                 param_index = self._raveled_index_for(target_param)
                 if self._has_fixes():
                     indices = np.r_[:self.size]
-                    which = (param_index[:,None]==indices[self._fixes_][None,:]).nonzero()
+                    which = (param_index[:, None] == indices[self._fixes_][None, :]).nonzero()
                     param_index = param_index[which[0]]
-                    transformed_index = (indices-(~self._fixes_).cumsum())[param_index]
+                    transformed_index = (indices - (~self._fixes_).cumsum())[param_index]
-                    #print param_index, transformed_index
+                    # print param_index, transformed_index
                 else:
                     transformed_index = param_index
 
@@ -340,7 +343,7 @@ class Model(Parameterized):
                 xx[xind] -= 2.*step
                 f2 = self._objective(xx)
                 numerical_gradient = (f1 - f2) / (2 * step)
-                if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind]
+                if np.all(gradient[xind] == 0): ratio = (f1 - f2) == gradient[xind]
                 else: ratio = (f1 - f2) / (2 * step * gradient[xind])
                 difference = np.abs((f1 - f2) / 2 / step - gradient[xind])
 
@@ -358,7 +361,7 @@ class Model(Parameterized):
                 grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4])
                 print grad_string
 
-            self._set_params_transformed(x)
+            self.optimizer_array = x
             return ret
 
 

GPy/core/parameterization/lists_and_dicts.py

@@ -77,8 +77,18 @@ class ObserverList(object):
             self._poc.insert(ins, (priority, weakref.ref(observer), callble))
 
     def __str__(self):
+        from . import ObsAr, Param
+        from parameter_core import Parameterizable
         ret = []
         curr_p = None
+        
+        def frmt(o):
+            if isinstance(o, ObsAr):
+                return 'ObsArr <{}>'.format(hex(id(o)))
+            elif isinstance(o, (Param,Parameterizable)):
+                return '{}'.format(o.hierarchy_name())
+            else:
+                return repr(o)                
         for p, o, c in self:
             curr = ''
             if curr_p != p:
@@ -87,8 +97,9 @@ class ObserverList(object):
             else: curr_pre = " "*len(pre)
             curr_p = p
             curr += curr_pre
-            ret.append(curr + ", ".join(map(repr, [o,c])))
+            
-        return '\n'.join(ret)
+            ret.append(curr + ", ".join([frmt(o), str(c)]))
+            return '\n'.join(ret)
 
     def flush(self):
         """

GPy/core/parameterization/observable_array.py

@@ -30,16 +30,22 @@ class ObsAr(np.ndarray, Pickleable, Observable):
     def __array_wrap__(self, out_arr, context=None):
         return out_arr.view(np.ndarray)
 
+    def _setup_observers(self):
+        # do not setup anything, as observable arrays do not have default observers
+        pass
+    
     def copy(self):
+        from lists_and_dicts import ObserverList
         memo = {}
         memo[id(self)] = self
+        memo[id(self.observers)] = ObserverList()
         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))
+        Pickleable.__setstate__(s, copy.deepcopy(self.__getstate__(), memo))
         return s
 
     def __reduce__(self):

GPy/core/parameterization/param.py

@@ -4,7 +4,7 @@
 import itertools
 import numpy
 np = numpy
-from parameter_core import Parameterizable, adjust_name_for_printing
+from parameter_core import Parameterizable, adjust_name_for_printing, Pickleable
 from observable_array import ObsAr
 
 ###### printing
@@ -173,36 +173,6 @@ class Param(Parameterizable, ObsAr):
     def _ensure_fixes(self):
         if not self._has_fixes(): self._fixes_ = numpy.ones(self._realsize_, dtype=bool)
 
-    #===========================================================================
-    # parameterizable
-    #===========================================================================
-    def traverse(self, visit, *args, **kwargs):
-        """
-        Traverse the hierarchy performing visit(self, *args, **kwargs) at every node passed by.
-        See "visitor pattern" in literature. This is implemented in pre-order fashion.
-
-        This will function will just call visit on self, as Param are leaf nodes.
-        """
-        self.__visited = True
-        visit(self, *args, **kwargs)
-        self.__visited = False
-
-    def traverse_parents(self, visit, *args, **kwargs):
-        """
-        Traverse the hierarchy upwards, visiting all parents and their children, except self.
-        See "visitor pattern" in literature. This is implemented in pre-order fashion.
-
-        Example:
-
-        parents = []
-        self.traverse_parents(parents.append)
-        print parents
-        """
-        if self.has_parent():
-            self.__visited = True
-            self._parent_._traverse_parents(visit, *args, **kwargs)
-            self.__visited = False
-
     #===========================================================================
     # Convenience
     #===========================================================================
@@ -217,13 +187,23 @@ class Param(Parameterizable, ObsAr):
     #===========================================================================
     # Pickling and copying
     #===========================================================================
+    def copy(self):
+        return Parameterizable.copy(self, which=self)
+    
     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))
+        Pickleable.__setstate__(s, copy.deepcopy(self.__getstate__(), memo))
         return s
+    def _setup_observers(self):
+        """
+        Setup the default observers
         
+        1: pass through to parent, if present
+        """
+        if self.has_parent():
+            self.add_observer(self._parent_, self._parent_._pass_through_notify_observers, -np.inf)
     
     #===========================================================================
     # Printing -> done

GPy/core/parameterization/parameter_core.py

@@ -16,8 +16,9 @@ Observable Pattern for patameterization
 from transformations import Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED
 import numpy as np
 import re
+import logging
 
-__updated__ = '2014-05-20'
+__updated__ = '2014-05-21'
 
 class HierarchyError(Exception):
     """
@@ -49,7 +50,6 @@ class Observable(object):
     as an observer. Every time the observable changes, it sends a notification with
     self as only argument to all its observers.
     """
-    _updated = True
     _updates = True
     def __init__(self, *args, **kwargs):
         super(Observable, self).__init__()
@@ -58,13 +58,19 @@ class Observable(object):
 
     @property
     def updates(self):
-        self._updates = self._highest_parent_._updates
+        p = getattr(self, '_highest_parent_', None)
+        if p is not None:
+            self._updates = p._updates
         return self._updates
 
     @updates.setter
     def updates(self, ups):
         assert isinstance(ups, bool), "updates are either on (True) or off (False)"
-        self._highest_parent_._updates = ups
+        p = getattr(self, '_highest_parent_', None)
+        if p is not None:
+            p._updates = ups
+        else:
+            self._updates = ups
         if ups:
             self._trigger_params_changed()
 
@@ -86,7 +92,7 @@ class Observable(object):
         for poc in self.observers:
             _, obs, clble = poc
             if callble is not None:
-                if (obs == observer) and (callble == clble):
+                if (obs is observer) and (callble == clble):
                     to_remove.append(poc)
             else:
                 if obs is observer:
@@ -172,6 +178,7 @@ class Pickleable(object):
     """
     def __init__(self, *a, **kw):
         super(Pickleable, self).__init__()
+
     #===========================================================================
     # Pickling operations
     #===========================================================================
@@ -192,37 +199,46 @@ class Pickleable(object):
     #===========================================================================
     # copy and pickling
     #===========================================================================
-    def copy(self):
+    def copy(self, memo=None, which=None):
         """
         Returns a (deep) copy of the current parameter handle. 
 
         All connections to parents of the copy will be cut.
+        
+        :param dict memo: memo for deepcopy
+        :param Parameterized which: parameterized object which started the copy process [default: self] 
         """
         #raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy"
+        if memo is None:
+            memo = {}
         import copy
-        memo = {}
         # the next part makes sure that we do not include parents in any form:
         parents = []
-        self.traverse_parents(parents.append) # collect parents
+        if which is None:
+            which = self
+        which.traverse_parents(parents.append) # collect parents
         for p in parents:
-            memo[id(p)] = None # set all parents to be None, so they will not be copied
+            if not memo.has_key(id(p)):memo[id(p)] = None # set all parents to be None, so they will not be copied
-        memo[id(self.gradient)] = None # reset the gradient
+        if not memo.has_key(id(self.gradient)):memo[id(self.gradient)] = None # reset the gradient
-        memo[id(self.param_array)] = None # and param_array
+        if not memo.has_key(id(self._fixes_)):memo[id(self._fixes_)] = None # fixes have to be reset, as this is now highest parent
-        memo[id(self._fixes_)] = None # fixes have to be reset, as this is now highest parent
+        copy = copy.deepcopy(self, memo) # and start the copy
-        c = copy.deepcopy(self, memo) # and start the copy
+        copy._parent_index_ = None
-        c._parent_index_ = None
+        copy._trigger_params_changed()
-        return c
+        return copy
 
     def __deepcopy__(self, memo):
         s = self.__new__(self.__class__) # fresh instance
         memo[id(self)] = s # be sure to break all cycles --> self is already done
         import copy
-        s.__dict__.update(copy.deepcopy(self.__dict__, memo)) # standard copy
+        s.__setstate__(copy.deepcopy(self.__getstate__(), memo)) # standard copy
         return s
 
     def __getstate__(self):
         ignore_list = ['_param_array_', # parameters get set from bottom to top
                        '_gradient_array_', # as well as gradients
+                       '_optimizer_copy_',
+                       'logger',
+                       'observers',
                        '_fixes_', # and fixes
                        '_Cacher_wrap__cachers', # never pickle cachers
                        ]
@@ -234,7 +250,11 @@ class Pickleable(object):
  
     def __setstate__(self, state):
         self.__dict__.update(state)
-        return self
+        from lists_and_dicts import ObserverList
+        self.observers = ObserverList()
+        self._setup_observers()
+        self._optimizer_copy_transformed = False
+
 
 class Gradcheckable(Pickleable, Parentable):
     """
@@ -387,7 +407,6 @@ class Indexable(Nameable, Observable):
         if value is not None:
             self[:] = value
 
-        #index = self._raveled_index()
         index = self.unconstrain()
         index = self._add_to_index_operations(self.constraints, index, __fixed__, warning)
         self._highest_parent_._set_fixed(self, index)
@@ -612,40 +631,104 @@ class OptimizationHandlable(Indexable):
     """
     This enables optimization handles on an Object as done in GPy 0.4.
 
-    `..._transformed`: make sure the transformations and constraints etc are handled
+    `..._optimizer_copy_transformed`: make sure the transformations and constraints etc are handled
     """
     def __init__(self, name, default_constraint=None, *a, **kw):
         super(OptimizationHandlable, self).__init__(name, default_constraint=default_constraint, *a, **kw)
+        self._optimizer_copy_ = None
+        self._optimizer_copy_transformed = False
 
-    def _get_params_transformed(self):
+    #===========================================================================
-        # transformed parameters (apply un-transformation rules)
+    # Optimizer copy
-        p = self.param_array.copy()
+    #===========================================================================    
-        [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
+    @property
-        if self.has_parent() and self.constraints[__fixed__].size != 0:
+    def optimizer_array(self):
-            fixes = np.ones(self.size).astype(bool)
+        """
-            fixes[self.constraints[__fixed__]] = FIXED
+        Array for the optimizer to work on.
-            return p[fixes]
+        This array always lives in the space for the optimizer.
-        elif self._has_fixes():
+        Thus, it is untransformed, going from Transformations.
-            return p[self._fixes_]
-        return p
         
-    def _set_params_transformed(self, p):
+        Setting this array, will make sure the transformed parameters for this model
+        will be set accordingly. It has to be set with an array, retrieved from
+        this method, as e.g. fixing will resize the array.  
+        
+        The optimizer should only interfere with this array, such that transofrmations
+        are secured.
         """
-        Set parameters p, but make sure they get transformed before setting.
+        if self.__dict__.get('_optimizer_copy_', None) is None or self.size != self._optimizer_copy_.size:
-        This means, the optimizer sees p, whereas the model sees transformed(p), 
+            self._optimizer_copy_ = np.empty(self.size)
-        such that, the parameters the model sees are in the right domain.
+        
-        """
+        if not self._optimizer_copy_transformed:
-        if not(p is self.param_array):
+            self._optimizer_copy_.flat = self.param_array.flat
+            [np.put(self._optimizer_copy_, ind, c.finv(self.param_array[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)
                 fixes[self.constraints[__fixed__]] = FIXED
-                self.param_array.flat[fixes] = p
+                return self._optimizer_copy_[fixes]
-            elif self._has_fixes(): self.param_array.flat[self._fixes_] = p
+            elif self._has_fixes():
-            else: self.param_array.flat = p
+                return self._optimizer_copy_[self._fixes_]
-        [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) 
+            self._optimizer_copy_transformed = True
-         for c, ind in self.constraints.iteritems() if c != __fixed__]
+        
+        return self._optimizer_copy_
+    
+    @optimizer_array.setter
+    def optimizer_array(self, p):
+        """
+        Make sure the optimizer copy does not get touched, thus, we only want to 
+        set the values *inside* not the array itself.
+        
+        Also we want to update param_array in here.
+        """
+        f = None
+        if self.has_parent() and self.constraints[__fixed__].size != 0:
+            f = np.ones(self.size).astype(bool)
+            f[self.constraints[__fixed__]] = FIXED
+        elif self._has_fixes(): 
+            f = self._fixes_
+        if f is None:
+            self.param_array.flat = p
+            [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) 
+             for c, ind in self.constraints.iteritems() if c != __fixed__]
+        else:
+            self.param_array.flat[f] = p
+            [np.put(self.param_array, ind[f[ind]], c.f(self.param_array.flat[ind[f[ind]]])) 
+             for c, ind in self.constraints.iteritems() if c != __fixed__]
+    
+        self._optimizer_copy_transformed = False
         self._trigger_params_changed()
 
+    def _get_params_transformed(self):
+        raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!"
+#         # transformed parameters (apply un-transformation rules)
+#         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)
+#             fixes[self.constraints[__fixed__]] = FIXED
+#             return p[fixes]
+#         elif self._has_fixes():
+#             return p[self._fixes_]
+#         return p
+# 
+    def _set_params_transformed(self, p):
+        raise DeprecationWarning, "_get|set_params{_optimizer_copy_transformed} is deprecated, use self.optimizer array insetad!"
+
+#         """
+#         Set parameters p, but make sure they get transformed before setting.
+#         This means, the optimizer sees p, whereas the model sees transformed(p), 
+#         such that, the parameters the model sees are in the right domain.
+#         """
+#         if not(p is self.param_array):
+#             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
+#         [np.put(self.param_array, ind, c.f(self.param_array.flat[ind])) 
+#          for c, ind in self.constraints.iteritems() if c != __fixed__]
+#         self._trigger_params_changed()
+
     def _trigger_params_changed(self, trigger_parent=True):
         """
         First tell all children to update,
@@ -724,7 +807,7 @@ class OptimizationHandlable(Indexable):
         x = rand_gen(loc=loc, scale=scale, size=self._size_transformed(), *args, **kwargs)
         # now draw from prior where possible
         [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None]
-        self._set_params_transformed(x)  # makes sure all of the tied parameters get the same init (since there's only one prior object...)
+        self.optimizer_array = x  # makes sure all of the tied parameters get the same init (since there's only one prior object...)
 
     #===========================================================================
     # For shared memory arrays. This does nothing in Param, but sets the memory
@@ -783,6 +866,7 @@ class Parameterizable(OptimizationHandlable):
         self.parameters = ArrayList()
         self._param_array_ = None
         self._added_names_ = set()
+        self.logger = logging.getLogger(self.__class__.__name__)
         self.__visited = False # for traversing in reverse order we need to know if we were here already
 
     @property
@@ -893,6 +977,11 @@ class Parameterizable(OptimizationHandlable):
         self._remove_parameter_name(None, old_name)
         self._add_parameter_name(param)
 
+    def __setstate__(self, state):
+        super(Parameterizable, self).__setstate__(state)
+        self.logger = logging.getLogger(self.__class__.__name__)
+        return self
+
     #===========================================================================
     # notification system
     #===========================================================================
@@ -900,7 +989,16 @@ class Parameterizable(OptimizationHandlable):
         self.parameters_changed()
     def _pass_through_notify_observers(self, me, which=None):
         self.notify_observers(which=which)
+    def _setup_observers(self):
+        """
+        Setup the default observers
         
+        1: parameters_changed_notify
+        2: pass through to parent, if present
+        """
+        self.add_observer(self, self._parameters_changed_notification, -100)
+        if self.has_parent():
+            self.add_observer(self._parent_, self._parent_._pass_through_notify_observers, -np.inf)
     #===========================================================================
     # From being parentable, we have to define the parent_change notification
     #===========================================================================

GPy/core/parameterization/parameterized.py

@@ -292,12 +292,16 @@ class Parameterized(Parameterizable):
         except Exception as e:
             print "WARNING: caught exception {!s}, trying to continue".format(e)
 
-    def copy(self):
+    def copy(self, memo=None):
-        c = super(Parameterized, self).copy()
+        if memo is None:
-        c._connect_parameters()
+            memo = {}
-        c._connect_fixes()
+        memo[id(self.optimizer_array)] = None # and param_array
-        c._notify_parent_change()
+        memo[id(self.param_array)] = None # and param_array
-        return c
+        copy = super(Parameterized, self).copy(memo)
+        copy._connect_parameters()
+        copy._connect_fixes()
+        copy._notify_parent_change()
+        return copy
 
     #===========================================================================
     # Printing:

GPy/defaults.cfg

@@ -1,9 +1,14 @@
-# This is the configuration file for GPy
+# This is the default configuration file for GPy
 
+# Do note edit this file.
+
+# For machine specific changes (i.e. those specific to a given installation) edit GPy/installation.cfg
+
+# For user specific changes edit $HOME/.gpy_user.cfg
 [parallel]
 # Enable openmp support. This speeds up some computations, depending on the number
 # of cores available. Setting up a compiler with openmp support can be difficult on
-# some platforms, hence this option.
+# some platforms, hence by default it is off.
 openmp=False
 
 [datasets]

GPy/examples/dimensionality_reduction.py

@@ -99,7 +99,7 @@ def sparse_gplvm_oil(optimize=True, verbose=0, plot=True, N=100, Q=6, num_induci
         m.kern.plot_ARD()
     return m
 
-def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4, sigma=.2):
+def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=25, Q=4, sigma=.2):
     import GPy
     from GPy.util.datasets import swiss_roll_generated
     from GPy.models import BayesianGPLVM
@@ -144,16 +144,15 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4
     m = BayesianGPLVM(Y, Q, X=X, X_variance=S, num_inducing=num_inducing, Z=Z, kernel=kernel)
     m.data_colors = c
     m.data_t = t
-    m['noise_variance'] = Y.var() / 100.
     
     if optimize:
-        m.optimize('scg', messages=verbose, max_iters=2e3)
+        m.optimize('bfgs', messages=verbose, max_iters=2e3)
 
     if plot:
         fig = plt.figure('fitted')
         ax = fig.add_subplot(111)
         s = m.input_sensitivity().argsort()[::-1][:2]
-        ax.scatter(*m.X.T[s], c=c)
+        ax.scatter(*m.X.mean.T[s], c=c)
 
     return m
 
@@ -172,14 +171,14 @@ def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40,
     m.data_labels = data['Y'][:N].argmax(axis=1)
     
     if optimize:
-        m.optimize('scg', messages=verbose, max_iters=max_iters, gtol=.05)
+        m.optimize('bfgs', messages=verbose, max_iters=max_iters, gtol=.05)
 
     if plot:
         fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
         m.plot_latent(ax=latent_axes, labels=m.data_labels)
         data_show = GPy.plotting.matplot_dep.visualize.vector_show((m.Y[0,:]))
         lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(param_to_array(m.X.mean)[0:1,:], # @UnusedVariable
-            m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
+            m, data_show, latent_axes=latent_axes, sense_axes=sense_axes, labels=m.data_labels)
         raw_input('Press enter to finish')
         plt.close(fig)
     return m
@@ -386,18 +385,17 @@ def brendan_faces(optimize=True, verbose=True, plot=True):
     Yn = Y - Y.mean()
     Yn /= Yn.std()
 
-    m = GPy.models.GPLVM(Yn, Q)
+    m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20)
 
     # optimize
-    m.constrain('rbf|noise|white', GPy.transformations.LogexpClipped())
 
-    if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
+    if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
 
     if plot:
         ax = m.plot_latent(which_indices=(0, 1))
-        y = m.likelihood.Y[0, :]
+        y = m.Y[0, :]
         data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(20, 28), transpose=True, order='F', invert=False, scale=False)
-        GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
+        lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
         raw_input('Press enter to finish')
 
     return m
@@ -411,13 +409,14 @@ def olivetti_faces(optimize=True, verbose=True, plot=True):
     Yn = Y - Y.mean()
     Yn /= Yn.std()
 
-    m = GPy.models.GPLVM(Yn, Q)
+    m = GPy.models.BayesianGPLVM(Yn, Q, num_inducing=20)
-    if optimize: m.optimize('scg', messages=verbose, max_iters=1000)
+    
+    if optimize: m.optimize('bfgs', messages=verbose, max_iters=1000)
     if plot:
         ax = m.plot_latent(which_indices=(0, 1))
         y = m.likelihood.Y[0, :]
         data_show = GPy.plotting.matplot_dep.visualize.image_show(y[None, :], dimensions=(112, 92), transpose=False, invert=False, scale=False)
-        GPy.plotting.matplot_dep.visualize.lvm(m.X[0, :].copy(), m, data_show, ax)
+        lvm = GPy.plotting.matplot_dep.visualize.lvm(m.X.mean[0, :].copy(), m, data_show, ax)
         raw_input('Press enter to finish')
 
     return m
@@ -500,9 +499,8 @@ def robot_wireless(optimize=True, verbose=True, plot=True):
 
     data = GPy.util.datasets.robot_wireless()
     # optimize
-    m = GPy.models.GPLVM(data['Y'], 2)
+    m = GPy.models.BayesianGPLVM(data['Y'], 4, num_inducing=25)
     if optimize: m.optimize(messages=verbose, max_f_eval=10000)
-    m._set_params(m._get_params())
     if plot:
         m.plot_latent()
 
@@ -523,7 +521,11 @@ def stick_bgplvm(model=None, optimize=True, verbose=True, plot=True):
     m.likelihood.variance = 0.001
 
     # optimize
-    if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10)
+    try:
+        if optimize: m.optimize('bfgs', messages=verbose, max_iters=5e3, bfgs_factor=10)
+    except KeyboardInterrupt:
+        print "Keyboard interrupt, continuing to plot and return"
+
     if plot:
         fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
         plt.sca(latent_axes)

GPy/inference/latent_function_inference/expectation_propagation.py

@@ -32,7 +32,7 @@ class EP(LatentFunctionInference):
         pass
 
     def inference(self, kern, X, likelihood, Y, Y_metadata=None, Z=None):
-        num_data, output_dim = X.shape
+        num_data, output_dim = Y.shape
         assert output_dim ==1, "ep in 1D only (for now!)"
 
         K = kern.K(X)

GPy/installation.cfg

@@ -0,0 +1,2 @@
+# This is the local configuration file for GPy
+

GPy/kern/__init__.py

@@ -14,6 +14,7 @@ from _src.ODE_UYC import ODE_UYC
 from _src.ODE_st import ODE_st
 from _src.ODE_t import ODE_t
 from _src.poly import Poly
+from _src.splitKern import SplitKern
 
 # TODO: put this in an init file somewhere
 #I'm commenting this out because the files were not added. JH. Remember to add the files before commiting

GPy/kern/_src/independent_outputs.py

@@ -22,7 +22,7 @@ def index_to_slices(index):
     """
 
     #contruct the return structure
-    ind = np.asarray(index,dtype=np.int64)
+    ind = np.asarray(index,dtype=np.int)
     ret = [[] for i in range(ind.max()+1)]
 
     #find the switchpoints

GPy/mappings/additive.py

@@ -0,0 +1,61 @@
+# Copyright (c) 2013, GPy authors (see AUTHORS.txt).
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+import numpy as np
+from ..core.mapping import Mapping
+import GPy
+
+class Additive(Mapping):
+    """
+    Mapping based on adding two existing mappings together.
+
+    .. math::
+
+       f(\mathbf{x}*) = f_1(\mathbf{x}*) + f_2(\mathbf(x)*)
+
+    :param mapping1: first mapping to add together.
+    :type mapping1: GPy.mappings.Mapping
+    :param mapping2: second mapping to add together.
+    :type mapping2: GPy.mappings.Mapping
+    :param tensor: whether or not to use the tensor product of input spaces
+    :type tensor: bool
+
+    """
+
+    def __init__(self, mapping1, mapping2, tensor=False):
+        if tensor:
+            input_dim = mapping1.input_dim + mapping2.input_dim
+        else:
+            input_dim = mapping1.input_dim
+            assert(mapping1.input_dim==mapping2.input_dim)
+        assert(mapping1.output_dim==mapping2.output_dim)
+        output_dim = mapping1.output_dim
+        Mapping.__init__(self, input_dim=input_dim, output_dim=output_dim)
+        self.mapping1 = mapping1
+        self.mapping2 = mapping2
+        self.num_params = self.mapping1.num_params + self.mapping2.num_params
+        self.name = self.mapping1.name + '+' + self.mapping2.name
+    def _get_param_names(self):
+        return self.mapping1._get_param_names + self.mapping2._get_param_names
+
+    def _get_params(self):
+        return np.hstack((self.mapping1._get_params() self.mapping2._get_params()))
+
+    def _set_params(self, x):
+        self.mapping1._set_params(x[:self.mapping1.num_params])
+        self.mapping2._set_params(x[self.mapping1.num_params:])
+        
+    def randomize(self):
+        self.mapping1._randomize()
+        self.mapping2._randomize()
+
+    def f(self, X):
+        return self.mapping1.f(X) + self.mapping2.f(X)
+
+    def df_dtheta(self, dL_df, X):
+        self._df_dA = (dL_df[:, :, None]*self.kern.K(X, self.X)[:, None, :]).sum(0).T
+        self._df_dbias = (dL_df.sum(0))
+        return np.hstack((self._df_dA.flatten(), self._df_dbias))
+
+    def df_dX(self, dL_df, X):
+        return self.kern.dK_dX((dL_df[:, None, :]*self.A[None, :, :]).sum(2), X, self.X) 

GPy/models/bayesian_gplvm.py

@@ -10,6 +10,7 @@ from ..util import linalg
 from ..core.parameterization.variational import NormalPosterior, NormalPrior, VariationalPosterior
 from ..inference.latent_function_inference.var_dtc_parallel import update_gradients
 from ..inference.latent_function_inference.var_dtc_gpu import VarDTC_GPU
+import logging
 
 class BayesianGPLVM(SparseGP):
     """
@@ -25,8 +26,10 @@ class BayesianGPLVM(SparseGP):
     """
     def __init__(self, Y, input_dim, X=None, X_variance=None, init='PCA', num_inducing=10,
                  Z=None, kernel=None, inference_method=None, likelihood=None, name='bayesian gplvm', **kwargs):
+        self.logger = logging.getLogger(self.__class__.__name__)
         if X == None:
             from ..util.initialization import initialize_latent
+            self.logger.info("initializing latent space X with method {}".format(init))
             X, fracs = initialize_latent(init, input_dim, Y)
         else:
             fracs = np.ones(input_dim)
@@ -36,7 +39,6 @@ class BayesianGPLVM(SparseGP):
         if X_variance is None:
             X_variance = np.random.uniform(0,.1,X.shape)
 
-
         if Z is None:
             Z = np.random.permutation(X.copy())[:num_inducing]
         assert Z.shape[1] == X.shape[1]
@@ -52,11 +54,14 @@ class BayesianGPLVM(SparseGP):
         X = NormalPosterior(X, X_variance)
 
         if inference_method is None:
-            if np.any(np.isnan(Y)):
+            inan = np.isnan(Y)
+            if np.any(inan):
                 from ..inference.latent_function_inference.var_dtc import VarDTCMissingData
-                inference_method = VarDTCMissingData()
+                self.logger.debug("creating inference_method with var_dtc missing data")
+                inference_method = VarDTCMissingData(inan=inan)
             else:
                 from ..inference.latent_function_inference.var_dtc import VarDTC
+                self.logger.debug("creating inference_method var_dtc")
                 inference_method = VarDTC()
 
         SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method, name, **kwargs)
@@ -100,36 +105,41 @@ class BayesianGPLVM(SparseGP):
         Notes:
         This will only work with a univariate Gaussian likelihood (for now)
         """
-        assert not self.likelihood.is_heteroscedastic
         N_test = Y.shape[0]
         input_dim = self.Z.shape[1]
+        
         means = np.zeros((N_test, input_dim))
         covars = np.zeros((N_test, input_dim))
         
-        dpsi0 = -0.5 * self.input_dim * self.likelihood.precision
+        dpsi0 = -0.5 * self.input_dim / self.likelihood.variance
-        dpsi2 = self.dL_dpsi2[0][None, :, :] # TODO: this may change if we ignore het. likelihoods
+        dpsi2 = self.grad_dict['dL_dpsi2'][0][None, :, :] # TODO: this may change if we ignore het. likelihoods
-        V = self.likelihood.precision * Y
+        V = Y/self.likelihood.variance
 
         #compute CPsi1V
-        if self.Cpsi1V is None:
+        #if self.Cpsi1V is None:
-            psi1V = np.dot(self.psi1.T, self.likelihood.V)
+        #    psi1V = np.dot(self.psi1.T, self.likelihood.V)
-            tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
+        #    tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
-            tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1)
+        #    tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1)
-            self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1)
+        #    self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1)
 
-        dpsi1 = np.dot(self.Cpsi1V, V.T)
+        dpsi1 = np.dot(self.posterior.woodbury_vector, V.T)
 
-        start = np.zeros(self.input_dim * 2)
+        #start = np.zeros(self.input_dim * 2)
+        
+
+        from scipy.optimize import minimize
 
         for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]):
-            args = (self.kern, self.Z, dpsi0, dpsi1_n.T, dpsi2)
+            args = (input_dim, self.kern.copy(), self.Z, dpsi0, dpsi1_n.T, dpsi2)
-            xopt, fopt, neval, status = SCG(f=latent_cost, gradf=latent_grad, x=start, optargs=args, display=False)
+            res = minimize(latent_cost_and_grad, jac=True, x0=np.hstack((means[n], covars[n])), args=args, method='BFGS')
-
+            xopt = res.x
             mu, log_S = xopt.reshape(2, 1, -1)
             means[n] = mu[0].copy()
             covars[n] = np.exp(log_S[0]).copy()
 
-        return means, covars
+        X = NormalPosterior(means, covars)
+
+        return X 
 
     def dmu_dX(self, Xnew):
         """
@@ -161,57 +171,26 @@ class BayesianGPLVM(SparseGP):
         return dim_reduction_plots.plot_steepest_gradient_map(self,*args,**kwargs)
     
 
-def latent_cost_and_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
+def latent_cost_and_grad(mu_S, input_dim, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
     """
     objective function for fitting the latent variables for test points
     (negative log-likelihood: should be minimised!)
     """
-    mu, log_S = mu_S.reshape(2, 1, -1)
+    mu = mu_S[:input_dim][None]
+    log_S = mu_S[input_dim:][None]
     S = np.exp(log_S)
 
-    psi0 = kern.psi0(Z, mu, S)
+    X = NormalPosterior(mu, S)
-    psi1 = kern.psi1(Z, mu, S)
-    psi2 = kern.psi2(Z, mu, S)
 
-    lik = dL_dpsi0 * psi0 + np.dot(dL_dpsi1.flatten(), psi1.flatten()) + np.dot(dL_dpsi2.flatten(), psi2.flatten()) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
+    psi0 = kern.psi0(Z, X)
+    psi1 = kern.psi1(Z, X)
+    psi2 = kern.psi2(Z, X)
 
-    mu0, S0 = kern.dpsi0_dmuS(dL_dpsi0, Z, mu, S)
+    lik = dL_dpsi0 * psi0.sum() + np.einsum('ij,kj->...', dL_dpsi1, psi1) + np.einsum('ijk,lkj->...', dL_dpsi2, psi2) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
-    mu1, S1 = kern.dpsi1_dmuS(dL_dpsi1, Z, mu, S)
-    mu2, S2 = kern.dpsi2_dmuS(dL_dpsi2, Z, mu, S)
     
-    dmu = mu0 + mu1 + mu2 - mu
+    dLdmu, dLdS = kern.gradients_qX_expectations(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, X)    
+    dmu = dLdmu - mu
     # dS = S0 + S1 + S2 -0.5 + .5/S
-    dlnS = S * (S0 + S1 + S2 - 0.5) + .5
+    dlnS = S * (dLdS - 0.5) + .5
+    
     return -lik, -np.hstack((dmu.flatten(), dlnS.flatten()))
-
-def latent_cost(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
-    """
-    objective function for fitting the latent variables (negative log-likelihood: should be minimised!)
-    This is the same as latent_cost_and_grad but only for the objective
-    """
-    mu, log_S = mu_S.reshape(2, 1, -1)
-    S = np.exp(log_S)
-
-    psi0 = kern.psi0(Z, mu, S)
-    psi1 = kern.psi1(Z, mu, S)
-    psi2 = kern.psi2(Z, mu, S)
-
-    lik = dL_dpsi0 * psi0 + np.dot(dL_dpsi1.flatten(), psi1.flatten()) + np.dot(dL_dpsi2.flatten(), psi2.flatten()) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
-    return -float(lik)
-
-def latent_grad(mu_S, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
-    """
-    This is the same as latent_cost_and_grad but only for the grad
-    """
-    mu, log_S = mu_S.reshape(2, 1, -1)
-    S = np.exp(log_S)
-
-    mu0, S0 = kern.dpsi0_dmuS(dL_dpsi0, Z, mu, S)
-    mu1, S1 = kern.dpsi1_dmuS(dL_dpsi1, Z, mu, S)
-    mu2, S2 = kern.dpsi2_dmuS(dL_dpsi2, Z, mu, S)
-
-    dmu = mu0 + mu1 + mu2 - mu
-    # dS = S0 + S1 + S2 -0.5 + .5/S
-    dlnS = S * (S0 + S1 + S2 - 0.5) + .5
-
-    return -np.hstack((dmu.flatten(), dlnS.flatten()))

GPy/models/mrd.py

@@ -2,10 +2,8 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-import itertools
+import itertools, logging
-import pylab
 
-from ..core import Model
 from ..kern import Kern
 from ..core.parameterization.variational import NormalPosterior, NormalPrior
 from ..core.parameterization import Param, Parameterized
@@ -61,15 +59,18 @@ class MRD(SparseGP):
                  inference_method=None, likelihoods=None, name='mrd', Ynames=None):
         super(GP, self).__init__(name)
 
+        self.logger = logging.getLogger(self.__class__.__name__)
         self.input_dim = input_dim
         self.num_inducing = num_inducing
 
         if isinstance(Ylist, dict):
             Ynames, Ylist = zip(*Ylist.items())
 
+        self.logger.debug("creating observable arrays")
         self.Ylist = [ObsAr(Y) for Y in Ylist]
 
         if Ynames is None:
+            self.logger.debug("creating Ynames")
             Ynames = ['Y{}'.format(i) for i in range(len(Ylist))]
         self.names = Ynames
         assert len(self.names) == len(self.Ylist), "one name per dataset, or None if Ylist is a dict"
@@ -81,13 +82,15 @@ class MRD(SparseGP):
                 inan = np.isnan(y)
                 if np.any(inan):
                     if not warned:
-                        print "WARING: NaN values detected, make sure initx method can cope with NaN values or provide starting latent space X"
+                        self.logger.warn("WARNING: NaN values detected, make sure initx method can cope with NaN values or provide starting latent space X")
                         warned = True
                     self.inference_method.append(VarDTCMissingData(limit=1, inan=inan))
                 else:
                     self.inference_method.append(VarDTC(limit=1))
+                self.logger.debug("created inference method <{}>".format(hex(id(self.inference_method[-1]))))
         else:
             if not isinstance(inference_method, InferenceMethodList):
+                self.logger.debug("making inference_method an InferenceMethodList")
                 inference_method = InferenceMethodList(inference_method)
             self.inference_method = inference_method
 
@@ -101,18 +104,19 @@ class MRD(SparseGP):
         self.num_inducing = self.Z.shape[0] # ensure M==N if M>N
 
         # sort out the kernels
+        self.logger.info("building kernels")
         if kernel is None:
             from ..kern import RBF
-            self.kernels = [RBF(input_dim, ARD=1, lengthscale=fracs[i]) for i in range(len(Ylist))]
+            kernels = [RBF(input_dim, ARD=1, lengthscale=fracs[i]) for i in range(len(Ylist))]
         elif isinstance(kernel, Kern):
-            self.kernels = []
+            kernels = []
             for i in range(len(Ylist)):
                 k = kernel.copy()
-                self.kernels.append(k)
+                kernels.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.kernels = kernel
+            kernels = kernel
 
         if X_variance is None:
             X_variance = np.random.uniform(0.1, 0.2, X.shape)
@@ -121,17 +125,17 @@ class MRD(SparseGP):
         self.X = NormalPosterior(X, X_variance)
 
         if likelihoods is None:
-            self.likelihoods = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))]
+            likelihoods = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))]
-        else: self.likelihoods = likelihoods
+        else: likelihoods = likelihoods
 
+        self.logger.info("adding X and Z")
         self.add_parameters(self.X, self.Z)
 
         self.bgplvms = []
         self.num_data = Ylist[0].shape[0]
 
-        for i, n, k, l, Y in itertools.izip(itertools.count(), Ynames, self.kernels, self.likelihoods, self.Ylist):
+        for i, n, k, l, Y in itertools.izip(itertools.count(), Ynames, kernels, likelihoods, Ylist):
             assert Y.shape[0] == self.num_data, "All datasets need to share the number of datapoints, and those have to correspond to one another"
-
             p = Parameterized(name=n)
             p.add_parameter(k)
             p.kern = k
@@ -141,6 +145,7 @@ class MRD(SparseGP):
             self.bgplvms.append(p)
 
         self.posterior = None
+        self.logger.info("init done")
         self._in_init_ = False
 
     def parameters_changed(self):
@@ -148,8 +153,9 @@ class MRD(SparseGP):
         self.posteriors = []
         self.Z.gradient[:] = 0.
         self.X.gradient[:] = 0.
-
+        for y, b, i in itertools.izip(self.Ylist, self.bgplvms, self.inference_method):
-        for y, k, l, i in itertools.izip(self.Ylist, self.kernels, self.likelihoods, self.inference_method):
+            self.logger.info('working on im <{}>'.format(hex(id(i))))
+            k, l = b.kern, b.likelihood
             posterior, lml, grad_dict = i.inference(k, self.X, self.Z, l, y)
 
             self.posteriors.append(posterior)
@@ -177,11 +183,11 @@ class MRD(SparseGP):
             self.X.mean.gradient += dL_dmean
             self.X.variance.gradient += dL_dS
 
-        # update for the KL divergence
         self.posterior = self.posteriors[0]
-        self.kern = self.kernels[0]
+        self.kern = self.bgplvms[0].kern
-        self.likelihood = self.likelihoods[0]
+        self.likelihood = self.bgplvms[0].likelihood
 
+        # update for the KL divergence
         self.variational_prior.update_gradients_KL(self.X)
         self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X)
 
@@ -219,8 +225,9 @@ class MRD(SparseGP):
         return Z
 
     def _handle_plotting(self, fignum, axes, plotf, sharex=False, sharey=False):
+        import matplotlib.pyplot as plt
         if axes is None:
-            fig = pylab.figure(num=fignum)
+            fig = plt.figure(num=fignum)
         sharex_ax = None
         sharey_ax = None
         plots = []
@@ -242,8 +249,8 @@ class MRD(SparseGP):
                 raise ValueError("Need one axes per latent dimension input_dim")
             plots.append(plotf(i, g, ax))
             if sharey_ax is not None:
-                pylab.setp(ax.get_yticklabels(), visible=False)
+                plt.setp(ax.get_yticklabels(), visible=False)
-        pylab.draw()
+        plt.draw()
         if axes is None:
             try:
                 fig.tight_layout()
@@ -257,8 +264,8 @@ class MRD(SparseGP):
         This predicts the output mean and variance for the dataset given in Ylist[Yindex]
         """
         self.posterior = self.posteriors[Yindex]
-        self.kern = self.kernels[Yindex]
+        self.kern = self.bgplvms[0].kern
-        self.likelihood = self.likelihoods[Yindex]
+        self.likelihood = self.bgplvms[0].likelihood
         return super(MRD, self).predict(Xnew, full_cov, Y_metadata, kern)
 
     #===============================================================================
@@ -300,11 +307,12 @@ class MRD(SparseGP):
         """
         import sys
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
+        from matplotlib import pyplot as plt
         from ..plotting.matplot_dep import dim_reduction_plots
         if "Yindex" not in predict_kwargs:
             predict_kwargs['Yindex'] = 0
         if ax is None:
-            fig = pylab.figure(num=fignum)
+            fig = plt.figure(num=fignum)
             ax = fig.add_subplot(111)
         else:
             fig = ax.figure
@@ -321,10 +329,7 @@ class MRD(SparseGP):
         return plot
 
     def __getstate__(self):
-        # TODO:
+        state = super(MRD, self).__getstate__()
-        import copy
-        state = copy.copy(self.__dict__)
-        del state['kernels']
         del state['kern']
         del state['likelihood']
         return state
@@ -332,7 +337,6 @@ class MRD(SparseGP):
     def __setstate__(self, state):
         # TODO:
         super(MRD, self).__setstate__(state)
-        self.kernels = [p.kern for p in self.bgplvms]
+        self.kern = self.bgplvms[0].kern
-        self.kern = self.kernels[0]
+        self.likelihood = self.bgplvms[0].likelihood
-        self.likelihood = self.likelihoods[0]
         self.parameters_changed()

GPy/plotting/matplot_dep/visualize.py

@@ -98,9 +98,9 @@ class lvm(matplotlib_show):
         """
         if vals is None:
             if isinstance(model.X, VariationalPosterior):
-                vals = param_to_array(model.X.mean)
+                vals = model.X.mean.values
             else:
-                vals = param_to_array(model.X)
+                vals = model.X.values
         if len(vals.shape)==1:
             vals = vals[None,:]
         matplotlib_show.__init__(self, vals, axes=latent_axes)
@@ -133,7 +133,7 @@ class lvm(matplotlib_show):
 
     def modify(self, vals):
         """When latent values are modified update the latent representation and ulso update the output visualization."""
-        self.vals = vals.copy()
+        self.vals = vals.view(np.ndarray).copy()
         y = self.model.predict(self.vals)[0]
         self.data_visualize.modify(y)
         self.latent_handle.set_data(self.vals[0,self.latent_index[0]], self.vals[0,self.latent_index[1]])
@@ -218,6 +218,7 @@ class lvm_dimselect(lvm):
         self.labels = labels
         lvm.__init__(self,vals,model,data_visualize,latent_axes,sense_axes,latent_index)
         self.show_sensitivities()
+        print self.latent_values
         print "use left and right mouse buttons to select dimensions"
 
 
@@ -247,6 +248,7 @@ class lvm_dimselect(lvm):
 
 
     def on_leave(self,event):
+        print type(self.latent_values)
         latent_values = self.latent_values.copy()
         y = self.model.predict(latent_values[None,:])[0]
         self.data_visualize.modify(y)

GPy/testing/model_tests.py

@@ -94,22 +94,18 @@ class MiscTests(unittest.TestCase):
         np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
 
         m.kern.lengthscale.randomize()
-        m._trigger_params_changed()
         m2.kern.lengthscale = m.kern.lengthscale
         np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
 
         m.kern.lengthscale.randomize()
-        m._trigger_params_changed()
         m2['.*lengthscale'] = m.kern.lengthscale
         np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
 
         m.kern.lengthscale.randomize()
-        m._trigger_params_changed()
         m2['.*lengthscale'] = m.kern['.*lengthscale']
         np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
 
         m.kern.lengthscale.randomize()
-        m._trigger_params_changed()
         m2.kern.lengthscale = m.kern['.*lengthscale']
         np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
 
@@ -130,6 +126,23 @@ class MiscTests(unittest.TestCase):
         m2.kern[:] = m.kern[''].values()
         np.testing.assert_equal(m.log_likelihood(), m2.log_likelihood())
 
+    def test_big_model(self):
+        m = GPy.examples.dimensionality_reduction.mrd_simulation(optimize=0, plot=0, plot_sim=0)
+        m.X.fix()
+        print m
+        m.unfix()
+        m.checkgrad()
+        print m
+        m.fix()
+        print m
+        m.inducing_inputs.unfix()
+        print m
+        m.checkgrad()
+        m.unfix()
+        m.checkgrad()
+        m.checkgrad()
+        print m
+
     def test_model_set_params(self):
         m = GPy.models.GPRegression(self.X, self.Y)
         lengthscale = np.random.uniform()

GPy/testing/pickle_tests.py

@@ -16,8 +16,7 @@ from GPy.core.parameterization.priors import Gaussian
 from GPy.kern._src.rbf import RBF
 from GPy.kern._src.linear import Linear
 from GPy.kern._src.static import Bias, White
-from GPy.examples.dimensionality_reduction import mrd_simulation,\
+from GPy.examples.dimensionality_reduction import mrd_simulation
-    bgplvm_simulation
 from GPy.examples.regression import toy_rbf_1d_50
 from GPy.core.parameterization.variational import NormalPosterior
 from GPy.models.gp_regression import GPRegression
@@ -90,6 +89,7 @@ class Test(ListDictTestCase):
         self.assertIs(pcopy.constraints, pcopy.rbf.lengthscale.constraints._param_index_ops)
         self.assertIs(pcopy.constraints, pcopy.linear.constraints._param_index_ops)
         self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
+        pcopy.gradient = 10 # gradient does not get copied anymore
         self.assertListEqual(par.gradient_full.tolist(), pcopy.gradient_full.tolist())
         self.assertSequenceEqual(str(par), str(pcopy))
         self.assertIsNot(par.param_array, pcopy.param_array)
@@ -126,8 +126,8 @@ class Test(ListDictTestCase):
     def test_modelrecreation(self):
         par = toy_rbf_1d_50(optimize=0, plot=0)
         pcopy = GPRegression(par.X.copy(), par.Y.copy(), kernel=par.kern.copy())
-        self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
+        np.testing.assert_allclose(par.param_array, pcopy.param_array)
-        self.assertListEqual(par.gradient_full.tolist(), pcopy.gradient_full.tolist())
+        np.testing.assert_allclose(par.gradient_full, pcopy.gradient_full)
         self.assertSequenceEqual(str(par), str(pcopy))
         self.assertIsNot(par.param_array, pcopy.param_array)
         self.assertIsNot(par.gradient_full, pcopy.gradient_full)
@@ -140,7 +140,7 @@ class Test(ListDictTestCase):
             par.pickle(f)
             f.seek(0)
             pcopy = pickle.load(f)
-        self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
+        np.testing.assert_allclose(par.param_array, pcopy.param_array)
         np.testing.assert_allclose(par.gradient_full, pcopy.gradient_full)
         self.assertSequenceEqual(str(par), str(pcopy))
         self.assert_(pcopy.checkgrad())
@@ -151,6 +151,7 @@ class Test(ListDictTestCase):
         par = NormalPosterior(X,Xv)
         par.gradient = 10
         pcopy = par.copy()
+        pcopy.gradient = 10
         self.assertListEqual(par.param_array.tolist(), pcopy.param_array.tolist())
         self.assertListEqual(par.gradient_full.tolist(), pcopy.gradient_full.tolist())
         self.assertSequenceEqual(str(par), str(pcopy))
@@ -175,6 +176,7 @@ class Test(ListDictTestCase):
         self.assertSequenceEqual(str(par), str(pcopy))
         self.assertIsNot(par.param_array, pcopy.param_array)
         self.assertIsNot(par.gradient_full, pcopy.gradient_full)
+        self.assertTrue(par.checkgrad())
         self.assertTrue(pcopy.checkgrad())
         self.assert_(np.any(pcopy.gradient!=0.0))
         with tempfile.TemporaryFile('w+b') as f:

GPy/util/caching.py

@@ -1,9 +1,7 @@
 from ..core.parameterization.parameter_core import Observable
-import itertools, collections, weakref
+import collections, weakref, logging
 
 class Cacher(object):
-
-
     def __init__(self, operation, limit=5, ignore_args=(), force_kwargs=()):
         """
         Parameters:
@@ -12,13 +10,15 @@ class Cacher(object):
         :param int limit: depth of cacher
         :param [int] ignore_args: list of indices, pointing at arguments to ignore in *args of operation(*args). This includes self!
         :param [str] force_kwargs: list of kwarg names (strings). If a kwarg with that name is given, the cacher will force recompute and wont cache anything.
+        :param int verbose: verbosity level. 0: no print outs, 1: casual print outs, 2: debug level print outs
         """
         self.limit = int(limit)
         self.ignore_args = ignore_args
         self.force_kwargs = force_kwargs
-        self.operation=operation
+        self.operation = operation
         self.order = collections.deque()
-        self.cached_inputs = {} # point from cache_ids to a list of [ind_ids], which where used in cache cache_id
+        self.cached_inputs = {}  # point from cache_ids to a list of [ind_ids], which where used in cache cache_id
+        self.logger = logging.getLogger("cache")
 
         #=======================================================================
         # point from each ind_id to [ref(obj), cache_ids]
@@ -27,53 +27,75 @@ class Cacher(object):
         self.cached_input_ids = {} 
         #=======================================================================
 
-        self.cached_outputs = {} # point from cache_ids to outputs
+        self.cached_outputs = {}  # point from cache_ids to outputs
-        self.inputs_changed = {} # point from cache_ids to bools
+        self.inputs_changed = {}  # point from cache_ids to bools
+
+    def id(self, obj):
+        """returns the self.id of an object, to be used in caching individual self.ids"""
+        return hex(id(obj))
 
     def combine_inputs(self, args, kw):
         "Combines the args and kw in a unique way, such that ordering of kwargs does not lead to recompute"
+        self.logger.debug("combining args and kw")
         return args + tuple(c[1] for c in sorted(kw.items(), key=lambda x: x[0]))
 
     def prepare_cache_id(self, combined_args_kw, ignore_args):
-        "get the cacheid (conc. string of argument ids in order) ignoring ignore_args"
+        "get the cacheid (conc. string of argument self.ids in order) ignoring ignore_args"
-        return "".join(str(id(a)) for i,a in enumerate(combined_args_kw) if i not in ignore_args)
+        cache_id = "".join(self.id(a) for i, a in enumerate(combined_args_kw) if i not in ignore_args)
+        self.logger.debug("cache_id={} was created".format(cache_id))
+        return cache_id
 
     def ensure_cache_length(self, cache_id):
         "Ensures the cache is within its limits and has one place free"
+        self.logger.debug("cache length gets ensured")
         if len(self.order) == self.limit:
+            self.logger.debug("cache limit of l={} was reached".format(self.limit))
             # we have reached the limit, so lets release one element
             cache_id = self.order.popleft()
+            self.logger.debug("cach_id '{}' gets removed".format(cache_id))
             combined_args_kw = self.cached_inputs[cache_id]
             for ind in combined_args_kw:
-                ind_id = id(ind)
+                if ind is not None:
-                ref, cache_ids = self.cached_input_ids[ind_id]
+                    ind_id = self.id(ind)
-                if len(cache_ids) == 1 and ref() is not None:
+                    tmp = self.cached_input_ids.get(ind_id, None)
-                    ref().remove_observer(self, self.on_cache_changed)
+                    if tmp is not None:
-                    del self.cached_input_ids[ind_id]
+                        ref, cache_ids = tmp
-                else:
+                        if len(cache_ids) == 1 and ref() is not None:
-                    cache_ids.remove(cache_id)
+                            ref().remove_observer(self, self.on_cache_changed)
-                    self.cached_input_ids[ind_id] = [ref, cache_ids]
+                            del self.cached_input_ids[ind_id]
+                        else:
+                            cache_ids.remove(cache_id)
+                            self.cached_input_ids[ind_id] = [ref, cache_ids]
+            self.logger.debug("removing caches")
             del self.cached_outputs[cache_id]
             del self.inputs_changed[cache_id]
             del self.cached_inputs[cache_id]
 
-    def add_to_cache(self, cache_id, combined_args_kw, output):
+    def add_to_cache(self, cache_id, inputs, output):
+        """This adds cache_id to the cache, with inputs and output"""
         self.inputs_changed[cache_id] = False
         self.cached_outputs[cache_id] = output
         self.order.append(cache_id)
-        self.cached_inputs[cache_id] = combined_args_kw
+        self.cached_inputs[cache_id] = inputs
-        for a in combined_args_kw:
+        for a in inputs:
-            ind_id = id(a)
+            if a is not None:
-            v = self.cached_input_ids.get(ind_id, [weakref.ref(a), []])
+                ind_id = self.id(a)
-            v[1].append(cache_id)
+                v = self.cached_input_ids.get(ind_id, [weakref.ref(a), []])
-            if len(v[1]) == 1:
+                self.logger.debug("cache_id '{}' gets stored".format(cache_id))
-                a.add_observer(self, self.on_cache_changed)
+                v[1].append(cache_id)
-            self.cached_input_ids[ind_id] = v
+                if len(v[1]) == 1:
+                    self.logger.debug("adding observer to object {}".format(repr(a)))
+                    a.add_observer(self, self.on_cache_changed)
+                self.cached_input_ids[ind_id] = v
 
     def __call__(self, *args, **kw):
         """
         A wrapper function for self.operation,
         """
+        #=======================================================================
+        # !WARNING CACHE OFFSWITCH!
+        # return self.operation(*args, **kw)
+        #=======================================================================
 
         # 1: Check whether we have forced recompute arguments:
         if len(self.force_kwargs) != 0:
@@ -81,27 +103,33 @@ class Cacher(object):
                 if k in kw and kw[k] is not None:
                     return self.operation(*args, **kw)
 
-        # 2: prepare_cache_id and get the unique id string for this call
+        # 2: prepare_cache_id and get the unique self.id string for this call
         inputs = self.combine_inputs(args, kw)
         cache_id = self.prepare_cache_id(inputs, self.ignore_args)
-    
         # 2: if anything is not cachable, we will just return the operation, without caching
-        if reduce(lambda a,b: a or (not isinstance(b, Observable)), inputs, False):
+        if reduce(lambda a, b: a or (not (isinstance(b, Observable) or b is None)), inputs, False):
+            self.logger.info("some inputs are not observable: returning without caching")
+            self.logger.debug(str(map(lambda x: isinstance(x, Observable) or x is None, inputs)))
+            self.logger.debug(str(map(repr, inputs)))
             return self.operation(*args, **kw)
         # 3&4: check whether this cache_id has been cached, then has it changed?
         try:
             if(self.inputs_changed[cache_id]):
-                # 4: This happens, when elements have changed for this cache id
+                self.logger.debug("{} already seen, but inputs changed. refreshing cacher".format(cache_id))
+                # 4: This happens, when elements have changed for this cache self.id
                 self.inputs_changed[cache_id] = False
                 self.cached_outputs[cache_id] = self.operation(*args, **kw)
         except KeyError:
+            self.logger.info("{} never seen, creating cache entry".format(cache_id))
             # 3: This is when we never saw this chache_id:
             self.ensure_cache_length(cache_id)
             self.add_to_cache(cache_id, inputs, self.operation(*args, **kw))
         except:
+            self.logger.error("an error occurred while trying to run caching for {}, resetting".format(cache_id))
             self.reset()
             raise
         # 5: We have seen this cache_id and it is cached:
+        self.logger.info("returning cache {}".format(cache_id))
         return self.cached_outputs[cache_id]
 
     def on_cache_changed(self, direct, which=None):
@@ -110,10 +138,13 @@ class Cacher(object):
 
         this function gets 'hooked up' to the inputs when we cache them, and upon their elements being changed we update here.
         """
-        for ind_id in [id(direct), id(which)]:
+        for what in [direct, which]:
-            _, cache_ids = self.cached_input_ids.get(ind_id, [None, []])
+            if what is not None:
-            for cache_id in cache_ids:
+                ind_id = self.id(what)
-                self.inputs_changed[cache_id] = True
+                _, cache_ids = self.cached_input_ids.get(ind_id, [None, []])
+                for cache_id in cache_ids:
+                    self.logger.info("callback from {} changed inputs from {}".format(ind_id, self.inputs_changed[cache_id]))
+                    self.inputs_changed[cache_id] = True
 
     def reset(self):
         """
@@ -150,7 +181,7 @@ class Cacher_wrap(object):
         return partial(self, obj)
     def __call__(self, *args, **kwargs):
         obj = args[0]
-        #import ipdb;ipdb.set_trace()
+        # import ipdb;ipdb.set_trace()
         try:
             caches = obj.__cachers
         except AttributeError:

GPy/util/config.py

@@ -5,18 +5,19 @@ import ConfigParser
 import os
 config = ConfigParser.ConfigParser()
 
-home = os.getenv('HOME') or os.getenv('USERPROFILE')
+# This is the default configuration file that always needs to be present.
-user_file = os.path.join(home,'.gpy_config.cfg')
+default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'defaults.cfg'))
-default_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'gpy_config.cfg'))
-# print user_file, os.path.isfile(user_file)
-# print default_file, os.path.isfile(default_file)
 
-# 1. check if the user has a ~/.gpy_config.cfg
+# These files are optional 
-if os.path.isfile(user_file):
+# This specifies configurations that are typically specific to the machine (it is found alongside the GPy installation).
-    config.read(user_file)
+local_file = os.path.abspath(os.path.join(os.path.dirname( __file__ ), '..', 'machine.cfg'))
-elif os.path.isfile(default_file):
+
-    # 2. if not, use the default one
+# This specifies configurations specific to the user (it is found in the user home directory)
-    config.read(default_file)
+home = os.getenv('HOME') or os.getenv('USERPROFILE')
-else:
+user_file = os.path.join(home,'.gpy_user.cfg')
-    #3. panic
+
-    raise ValueError, "no configuration file found"
+# Read in the given files.
+config.readfp(open(default_file))
+config.read([local_file, user_file])
+if not config:
+    raise ValueError, "No configuration file found at either " + user_file + " or " + local_file + " or " + default_file + "."

GPy/util/data_resources.json

@@ -57,6 +57,20 @@
             "http://www.cs.nyu.edu/~roweis/data/"
         ]
     },
+    "cifar-10": {
+        "citation": "Learning Multiple Layers of Features from Tiny Images, Alex Krizhevsky, 2009, Tech report available here: http://www.cs.toronto.edu/~kriz/learning-features-2009-TR.pdf",
+        "details": "The CIFAR-10 and CIFAR-100 are labeled subsets of the 80 million tiny images dataset. They were collected by Alex Krizhevsky, Vinod Nair, and Geoffrey Hinton. Details are available on this webpage: http://www.cs.toronto.edu/~kriz/cifar.html. The CIFAR-10 dataset consists of 60000 32x32 colour images in 10 classes, with 6000 images per class. There are 50000 training images and 10000 test images.",
+        "files": [
+            [
+                "cifar-10-python.tar.gz"
+            ]
+        ],
+        "license": null,
+        "size": 0,
+        "urls": [
+            "http://www.cs.toronto.edu/~kriz/"
+        ]
+    },
     "cmu_mocap_full": {
         "citation": "Please include this in your acknowledgements: The data used in this project was obtained from mocap.cs.cmu.edu.\\nThe database was created with funding from NSF EIA-0196217.",
         "details": "CMU Motion Capture data base. Captured by a Vicon motion capture system consisting of 12 infrared MX-40 cameras, each of which is capable of recording at 120 Hz with images of 4 megapixel resolution. Motions are captured in a working volume of approximately 3m x 8m. The capture subject wears 41 markers and a stylish black garment.",
@@ -123,11 +137,39 @@
             ]
         ],
         "license": null,
-        "size": 0,
+        "size": 20258,
         "urls": [
             "http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/drosophila_protein/"
         ]
     },
+    "spellman_yeast": {
+        "citation": "Paul T. Spellman, Gavin Sherlock, Michael Q. Zhang, Vishwanath R. Iyer, Kirk Anders, Michael B. Eisen, Patrick O. Brown, David Botstein, and Bruce Futcher 'Comprehensive Identification of Cell Cycle-regulated Genes of the Yeast Saccharomyces cerevisiae by Microarray Hybridization.'  Molecular Biology of the Cell 9, 3273-3297",
+        "details": "Two colour spotted cDNA array data set of a series of experiments to identify which genes in Yeast are cell cycle regulated.",
+        "files": [
+            [
+                "combined.txt"
+            ]
+        ],
+        "license": null,
+        "size": 2510955,
+        "urls": [
+            "http://genome-www.stanford.edu/cellcycle/data/rawdata/"
+        ]
+    },
+    "lee_yeast_ChIP": {
+        "citation": "Tong Ihn Lee, Nicola J. Rinaldi, Francois Robert, Duncan T. Odom, Ziv Bar-Joseph, Georg K. Gerber, Nancy M. Hannett, Christopher T. Harbison, Craig M. Thompson, Itamar Simon, Julia Zeitlinger, Ezra G. Jennings, Heather L. Murray, D. Benjamin Gordon, Bing Ren, John J. Wyrick, Jean-Bosco Tagne, Thomas L. Volkert, Ernest Fraenkel, David K. Gifford, Richard A. Young 'Transcriptional Regulatory Networks in Saccharomyces cerevisiae' Science 298 (5594) pg 799--804. DOI: 10.1126/science.1075090",
+        "details": "Binding location analysis for 106 regulators in yeast. The data consists of p-values for binding of regulators to genes derived from ChIP-chip experiments.",
+        "files": [
+            [
+                "binding_by_gene.tsv"
+            ]
+        ],
+        "license": null,
+        "size": 1674161,
+        "urls": [
+            "http://jura.wi.mit.edu/young_public/regulatory_network/"
+        ]
+    },
     "epomeo_gpx": {
         "citation": "",
         "details": "Five different GPS traces of the same run up Mount Epomeo in Ischia. The traces are from different sources. endomondo_1 and endomondo_2 are traces from the mobile phone app Endomondo, with a split in the middle. garmin_watch_via_endomondo is the trace from a Garmin watch, with a segment missing about 4 kilometers in. viewranger_phone and viewranger_tablet are traces from a phone and a tablet through the viewranger app. The viewranger_phone data comes from the same mobile phone as the Endomondo data (i.e. there are 3 GPS devices, but one device recorded two traces).",

GPy/util/datasets.py

@@ -240,7 +240,7 @@ def cmu_urls_files(subj_motions, messages = True):
         if not os.path.exists(cur_skel_file):
             # Current skel file doesn't exist.
             if not os.path.isdir(skel_dir):
-                os.mkdir(skel_dir)
+                os.makedirs(skel_dir)
             # Add skel file to list.
             url_required = True
             file_download.append(subjects[i] + '.asf')
@@ -367,8 +367,8 @@ def sod1_mouse(data_set='sod1_mouse'):
     if not data_available(data_set):
         download_data(data_set)
     from pandas import read_csv
-    dirpath = os.path.join(data_path, data_set)
+    dir_path = os.path.join(data_path, data_set)
-    filename = os.path.join(dirpath, 'sod1_C57_129_exprs.csv')
+    filename = os.path.join(dir_path, 'sod1_C57_129_exprs.csv')
     Y = read_csv(filename, header=0, index_col=0)
     num_repeats=4
     num_time=4
@@ -376,12 +376,48 @@ def sod1_mouse(data_set='sod1_mouse'):
     X = 1
     return data_details_return({'X': X, 'Y': Y}, data_set)
 
+def spellman_yeast(data_set='spellman_yeast'):
+    if not data_available(data_set):
+        download_data(data_set)
+    from pandas import read_csv
+    dir_path = os.path.join(data_path, data_set)
+    filename = os.path.join(dir_path, 'combined.txt')
+    Y = read_csv(filename, header=0, index_col=0, sep='\t')
+    return data_details_return({'Y': Y}, data_set)
+
+def spellman_yeast_cdc(data_set='spellman_yeast'):
+    if not data_available(data_set):
+        download_data(data_set)
+    from pandas import read_csv
+    dir_path = os.path.join(data_path, data_set)
+    filename = os.path.join(dir_path, 'combined.txt')
+    Y = read_csv(filename, header=0, index_col=0, sep='\t')
+    t = np.asarray([10, 30, 50, 70, 80, 90, 100, 110, 120, 130, 140, 150, 170, 180, 190, 200, 210, 220, 230, 240, 250, 270, 290])
+    times = ['cdc15_'+str(time) for time in t]
+    Y = Y[times].T
+    t = t[:, None]
+    return data_details_return({'Y' : Y, 't': t, 'info': 'Time series of synchronized yeast cells from the CDC-15 experiment of Spellman et al (1998).'}, data_set)
+
+def lee_yeast_ChIP(data_set='lee_yeast_ChIP'):
+    if not data_available(data_set):
+        download_data(data_set)
+    from pandas import read_csv
+    import zipfile
+    dir_path = os.path.join(data_path, data_set)
+    filename = os.path.join(dir_path, 'binding_by_gene.tsv')
+    X = read_csv(filename, header=1, index_col=0, sep='\t')
+    transcription_factors = [col for col in X.columns if col[:7] != 'Unnamed']    
+    annotations = X[['Unnamed: 1', 'Unnamed: 2', 'Unnamed: 3']]
+    X = X[transcription_factors]
+    return data_details_return({'annotations' : annotations, 'X' : X, 'transcription_factors': transcription_factors}, data_set)
+    
+
 def fruitfly_tomancak(data_set='fruitfly_tomancak', gene_number=None):
     if not data_available(data_set):
         download_data(data_set)
     from pandas import read_csv
-    dirpath = os.path.join(data_path, data_set)
+    dir_path = os.path.join(data_path, data_set)
-    filename = os.path.join(dirpath, 'tomancak_exprs.csv')
+    filename = os.path.join(dir_path, 'tomancak_exprs.csv')
     Y = read_csv(filename, header=0, index_col=0).T
     num_repeats = 3
     num_time = 12
@@ -395,8 +431,8 @@ def drosophila_protein(data_set='drosophila_protein'):
     if not data_available(data_set):
         download_data(data_set)
     from pandas import read_csv
-    dirpath = os.path.join(data_path, data_set)
+    dir_path = os.path.join(data_path, data_set)
-    filename = os.path.join(dirpath, 'becker_et_al.csv')
+    filename = os.path.join(dir_path, 'becker_et_al.csv')
     Y = read_csv(filename, header=0)
     return data_details_return({'Y': Y}, data_set)
 
@@ -404,8 +440,8 @@ def drosophila_knirps(data_set='drosophila_protein'):
     if not data_available(data_set):
         download_data(data_set)
     from pandas import read_csv
-    dirpath = os.path.join(data_path, data_set)
+    dir_path = os.path.join(data_path, data_set)
-    filename = os.path.join(dirpath, 'becker_et_al.csv')
+    filename = os.path.join(dir_path, 'becker_et_al.csv')
     # in the csv file we have facts_kni and ext_kni. We treat facts_kni as protein and ext_kni as mRNA
     df = read_csv(filename, header=0)
     t = df['t'][:,None]
@@ -426,31 +462,59 @@ def drosophila_knirps(data_set='drosophila_protein'):
     return data_details_return({'Y': Y, 'X': X}, data_set)
 
 # This will be for downloading google trends data.
-def google_trends(query_terms=['big data', 'machine learning', 'data science'], data_set='google_trends'):
+def google_trends(query_terms=['big data', 'machine learning', 'data science'], data_set='google_trends', refresh_data=False):
-    """Data downloaded from Google trends for given query terms. Warning, if you use this function multiple times in a row you get blocked due to terms of service violations."""
+    """Data downloaded from Google trends for given query terms. Warning, if you use this function multiple times in a row you get blocked due to terms of service violations. The function will cache the result of your query, if you wish to refresh an old query set refresh_data to True. The function is inspired by this notebook: http://nbviewer.ipython.org/github/sahuguet/notebooks/blob/master/GoogleTrends%20meet%20Notebook.ipynb"""
-    # Inspired by this notebook:
+    query_terms.sort()
-    # http://nbviewer.ipython.org/github/sahuguet/notebooks/blob/master/GoogleTrends%20meet%20Notebook.ipynb
+    import pandas
         
-    # quote the query terms.
+    # Create directory name for data
-    for i, element in enumerate(query_terms):
+    dir_path = os.path.join(data_path,'google_trends')
-        query_terms[i] = urllib2.quote(element)
+    if not os.path.isdir(dir_path):
-    query = 'http://www.google.com/trends/fetchComponent?q=%s&cid=TIMESERIES_GRAPH_0&export=3' % ",".join(query_terms)
+        os.makedirs(dir_path)
+    dir_name = '-'.join(query_terms)
+    dir_name = dir_name.replace(' ', '_')
+    dir_path = os.path.join(dir_path,dir_name)
+    file = 'data.csv'
+    file_name = os.path.join(dir_path,file)
+    if not os.path.exists(file_name) or refresh_data:
+        print "Accessing Google trends to acquire the data. Note that repeated accesses will result in a block due to a google terms of service violation. Failure at this point may be due to such blocks."
+        # quote the query terms.
+        quoted_terms = []
+        for term in query_terms:
+            quoted_terms.append(urllib2.quote(term))
+        print "Query terms: ", ', '.join(query_terms)
 
-    data = urllib2.urlopen(query).read()
+        print "Fetching query:"
+        query = 'http://www.google.com/trends/fetchComponent?q=%s&cid=TIMESERIES_GRAPH_0&export=3' % ",".join(quoted_terms)
 
-    # In the notebook they did some data cleaning: remove Javascript header+footer, and translate new Date(....,..,..) into YYYY-MM-DD.
+        data = urllib2.urlopen(query).read()
-    header = """// Data table response\ngoogle.visualization.Query.setResponse("""
+        print "Done."
-    data = data[len(header):-2]
+        # In the notebook they did some data cleaning: remove Javascript header+footer, and translate new Date(....,..,..) into YYYY-MM-DD.
-    data = re.sub('new Date\((\d+),(\d+),(\d+)\)', (lambda m: '"%s-%02d-%02d"' % (m.group(1).strip(), 1+int(m.group(2)), int(m.group(3)))), data)
+        header = """// Data table response\ngoogle.visualization.Query.setResponse("""
-    timeseries = json.loads(data)
+        data = data[len(header):-2]
-    #import pandas as pd
+        data = re.sub('new Date\((\d+),(\d+),(\d+)\)', (lambda m: '"%s-%02d-%02d"' % (m.group(1).strip(), 1+int(m.group(2)), int(m.group(3)))), data)
-    columns = [k['label'] for k in timeseries['table']['cols']]
+        timeseries = json.loads(data)
-    rows = map(lambda x: [k['v'] for k in x['c']], timeseries['table']['rows'])
+        columns = [k['label'] for k in timeseries['table']['cols']]
-    terms = len(columns)-1
+        rows = map(lambda x: [k['v'] for k in x['c']], timeseries['table']['rows'])
-    X = np.asarray([(pb.datestr2num(row[0]), i) for i in range(terms) for row in rows ])
+        df = pandas.DataFrame(rows, columns=columns)
-    Y = np.asarray([[row[i+1]] for i in range(terms) for row in rows ])
+        if not os.path.isdir(dir_path):
+            os.makedirs(dir_path)
+
+        df.to_csv(file_name)
+    else:
+        print "Reading cached data for google trends. To refresh the cache set 'refresh_data=True' when calling this function."
+        print "Query terms: ", ', '.join(query_terms)
+
+        df = pandas.read_csv(file_name, parse_dates=[0])
+
+    columns = df.columns
+    terms = len(query_terms)
+    import datetime
+    X = np.asarray([(row, i) for i in range(terms) for row in df.index])
+    Y = np.asarray([[df.ix[row][query_terms[i]]] for i in range(terms) for row in df.index ])
     output_info = columns[1:]
-    return data_details_return({'X': X, 'Y': Y, 'query_terms': output_info, 'info': "Data downloaded from google trends with query terms: " + ', '.join(output_info) + '.'}, data_set)
+        
+    return data_details_return({'data frame' : df, 'X': X, 'Y': Y, 'query_terms': output_info, 'info': "Data downloaded from google trends with query terms: " + ', '.join(output_info) + '.'}, data_set)
     
 # The data sets
 def oil(data_set='three_phase_oil_flow'):
@@ -594,7 +658,7 @@ def ripley_synth(data_set='ripley_prnn_data'):
     ytest = test[:, 2:3]
     return data_details_return({'X': X, 'Y': y, 'Xtest': Xtest, 'Ytest': ytest, 'info': 'Synthetic data generated by Ripley for a two class classification problem.'}, data_set)
 
-def mauna_loa(data_set='mauna_loa', num_train=543, refresh_data=False):
+def mauna_loa(data_set='mauna_loa', num_train=545, refresh_data=False):
     path = os.path.join(data_path, data_set)
     if data_available(data_set) and not refresh_data:
         print 'Using cached version of the data set, to use latest version set refresh_data to True'
@@ -635,7 +699,7 @@ def osu_run1(data_set='osu_run1', sample_every=4):
     return data_details_return({'Y': Y, 'connect' : connect}, data_set)
 
 def swiss_roll_generated(num_samples=1000, sigma=0.0):
-    with open(os.path.join(data_path, 'swiss_roll.pickle')) as f:
+    with open(os.path.join(os.path.dirname(__file__), 'datasets', 'swiss_roll.pickle')) as f:
         data = pickle.load(f)
     Na = data['Y'].shape[0]
     perm = np.random.permutation(np.r_[:Na])[:num_samples]
@@ -688,15 +752,15 @@ def hapmap3(data_set='hapmap3'):
     except ImportError as i:
         raise i, "Need pandas for hapmap dataset, make sure to install pandas (http://pandas.pydata.org/) before loading the hapmap dataset"
 
-    dirpath = os.path.join(data_path,'hapmap3')
+    dir_path = os.path.join(data_path,'hapmap3')
     hapmap_file_name = 'hapmap3_r2_b36_fwd.consensus.qc.poly'
-    unpacked_files = [os.path.join(dirpath, hapmap_file_name+ending) for ending in ['.ped', '.map']]
+    unpacked_files = [os.path.join(dir_path, hapmap_file_name+ending) for ending in ['.ped', '.map']]
     unpacked_files_exist = reduce(lambda a, b:a and b, map(os.path.exists, unpacked_files))
 
     if not unpacked_files_exist and not data_available(data_set):
         download_data(data_set)
 
-    preprocessed_data_paths = [os.path.join(dirpath,hapmap_file_name + file_name) for file_name in \
+    preprocessed_data_paths = [os.path.join(dir_path,hapmap_file_name + file_name) for file_name in \
                                ['.snps.pickle',
                                 '.info.pickle',
                                 '.nan.pickle']]
@@ -739,7 +803,7 @@ def hapmap3(data_set='hapmap3'):
         mapnp = np.loadtxt(unpacked_files[1], dtype=str)
         status=write_status('reading relationships.txt...', 42, status)
         # and metainfo:
-        infodf = DataFrame.from_csv(os.path.join(dirpath,'./relationships_w_pops_121708.txt'), header=0, sep='\t')
+        infodf = DataFrame.from_csv(os.path.join(dir_path,'./relationships_w_pops_121708.txt'), header=0, sep='\t')
         infodf.set_index('IID', inplace=1)
         status=write_status('filtering nan...', 45, status)
         snpstr = snpstrnp[:,6:].astype('S1').reshape(snpstrnp.shape[0], -1, 2)
@@ -804,12 +868,12 @@ def singlecell(data_set='singlecell'):
         download_data(data_set)
 
     from pandas import read_csv
-    dirpath = os.path.join(data_path, data_set)
+    dir_path = os.path.join(data_path, data_set)
-    filename = os.path.join(dirpath, 'singlecell.csv')
+    filename = os.path.join(dir_path, 'singlecell.csv')
     Y = read_csv(filename, header=0, index_col=0)
     genes = Y.columns
     labels = Y.index
-    # data = np.loadtxt(os.path.join(dirpath, 'singlecell.csv'), delimiter=",", dtype=str)
+    # data = np.loadtxt(os.path.join(dir_path, 'singlecell.csv'), delimiter=",", dtype=str)
     return data_details_return({'Y': Y, 'info' : "qPCR singlecell experiment in Mouse, measuring 48 gene expressions in 1-64 cell states. The labels have been created as in Guo et al. [2010]",
                                 'genes': genes, 'labels':labels,
                                 }, data_set)
@@ -1109,6 +1173,30 @@ def creep_data(data_set='creep_rupture'):
     X = all_data[:, features].copy()
     return data_details_return({'X': X, 'y': y}, data_set)
 
+def cifar10_patches(data_set='cifar-10'):
+    """The Candian Institute for Advanced Research 10 image data set. Code for loading in this data is taken from this Boris Babenko's blog post, original code available here: http://bbabenko.tumblr.com/post/86756017649/learning-low-level-vision-feautres-in-10-lines-of-code"""
+    dir_path = os.path.join(data_path, data_set)
+    filename = os.path.join(dir_path, 'cifar-10-python.tar.gz')
+    if not data_available(data_set):
+        download_data(data_set)
+        import tarfile
+        # This code is from Boris Babenko's blog post.
+        # http://bbabenko.tumblr.com/post/86756017649/learning-low-level-vision-feautres-in-10-lines-of-code
+        tfile = tarfile.open(filename, 'r:gz')
+        tfile.extractall(dir_path)
+
+    with open(os.path.join(dir_path, 'cifar-10-batches-py','data_batch_1'),'rb') as f:
+        data = pickle.load(f)
+
+    images = data['data'].reshape((-1,3,32,32)).astype('float32')/255
+    images = np.rollaxis(images, 1, 4)
+    patches = np.zeros((0,5,5,3))
+    for x in range(0,32-5,5):
+        for y in range(0,32-5,5):
+            patches = np.concatenate((patches, images[:,x:x+5,y:y+5,:]), axis=0)
+    patches = patches.reshape((patches.shape[0],-1))    
+    return data_details_return({'Y': patches, "info" : "32x32 pixel patches extracted from the CIFAR-10 data by Boris Babenko to demonstrate k-means features."}, data_set)
+
 def cmu_mocap_49_balance(data_set='cmu_mocap'):
     """Load CMU subject 49's one legged balancing motion that was used by Alvarez, Luengo and Lawrence at AISTATS 2009."""
     train_motions = ['18', '19']

GPy/util/subarray_and_sorting.py

@@ -4,9 +4,9 @@
 .. moduleauthor:: Max Zwiessele <ibinbei@gmail.com>
 
 '''
-__updated__ = '2014-05-20'
+__updated__ = '2014-05-21'
 
-import numpy as np
+import numpy as np, logging
 
 def common_subarrays(X, axis=0):
     """
@@ -48,7 +48,17 @@ def common_subarrays(X, axis=0):
     assert X.ndim == 2 and axis in (0,1), "Only implemented for 2D arrays"
     subarrays = defaultdict(list)
     cnt = count()
-    np.apply_along_axis(lambda x: iadd(subarrays[tuple(x)], [cnt.next()]), 1-axis, X)
+    logger = logging.getLogger("common_subarrays")
+    def accumulate(x, s, c):
+        logger.debug("creating tuple")
+        t = tuple(x)
+        logger.debug("tuple done")
+        col = c.next()
+        iadd(s[t], [col])
+        logger.debug("added col {}".format(col))
+        return None
+    if axis == 0: [accumulate(x, subarrays, cnt) for x in X]
+    else: [accumulate(x, subarrays, cnt) for x in X.T]
     return subarrays
 
 if __name__ == '__main__':