Merge branch 'params' of github.com:SheffieldML/GPy into params

GPy/core/model.py

@@ -253,7 +253,7 @@ class Model(Parameterized):
         sgd.run()
         self.optimization_runs.append(sgd)
 
-    def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3, _debug=False):
+    def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3):
         """
         Check the gradient of the ,odel by comparing to a numerical
         estimate.  If the verbose flag is passed, invividual
@@ -349,13 +349,6 @@ class Model(Parameterized):
                 xx[xind] -= 2.*step
                 f2 = self.objective_function(xx)
                 numerical_gradient = (f1 - f2) / (2 * step)
-                if _debug:
-                    for p in self.kern.flattened_parameters:
-                        p._parent_._debug=True
-                    self.gradient[xind] = numerical_gradient
-                    self._set_params_transformed(x)
-                    for p in self.kern.flattened_parameters:
-                        p._parent_._debug=False
                 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])

GPy/core/parameterization/param.py

@@ -446,8 +446,8 @@ class ParamConcatenation(object):
     def untie(self, *ties):
         [param.untie(*ties) for param in self.params]
 
-    def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3, _debug=False):
+    def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3):
-        return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance, _debug=_debug)
+        return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance)
     #checkgrad.__doc__ = Gradcheckable.checkgrad.__doc__
 
     __lt__ = lambda self, val: self._vals() < val

GPy/core/parameterization/parameter_core.py

@@ -41,7 +41,7 @@ class Observable(object):
     """
     _updated = True
     def __init__(self, *args, **kwargs):
-        super(Observable, self).__init__(*args, **kwargs)
+        super(Observable, self).__init__()
         self._observer_callables_ = []
 
     def add_observer(self, observer, callble, priority=0):
@@ -160,7 +160,7 @@ class Parentable(object):
     _parent_ = None
     _parent_index_ = None
     def __init__(self, *args, **kwargs):
-        super(Parentable, self).__init__(*args, **kwargs)
+        super(Parentable, self).__init__()
 
     def has_parent(self):
         """
@@ -206,7 +206,7 @@ class Gradcheckable(Parentable):
     def __init__(self, *a, **kw):
         super(Gradcheckable, self).__init__(*a, **kw)
 
-    def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3, _debug=False):
+    def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3):
         """
         Check the gradient of this parameter with respect to the highest parent's
         objective function.
@@ -220,10 +220,10 @@ class Gradcheckable(Parentable):
         :param flaot tolerance: the tolerance for the gradient ratio or difference.
         """
         if self.has_parent():
-            return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance, _debug=_debug)
+            return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance)
-        return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance, _debug=_debug)
+        return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance)
 
-    def _checkgrad(self, param, verbose=0, step=1e-6, tolerance=1e-3, _debug=False):
+    def _checkgrad(self, param, verbose=0, step=1e-6, tolerance=1e-3):
         """
         Perform the checkgrad on the model.
         TODO: this can be done more efficiently, when doing it inside here
@@ -275,7 +275,7 @@ class Indexable(object):
     The raveled index of an object is the index for its parameters in a flattened int array.
     """
     def __init__(self, *a, **kw):
-        super(Indexable, self).__init__(*a, **kw)
+        super(Indexable, self).__init__()
 
     def _raveled_index(self):
         """
@@ -308,7 +308,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!!
@@ -352,9 +352,11 @@ class Constrainable(Nameable, Indexable):
         """
         if value is not None:
             self[:] = value
-        self.constrain(__fixed__, warning=warning, trigger_parent=trigger_parent)
+        reconstrained = self.unconstrain()
+        self._add_to_index_operations(self.constraints, reconstrained, __fixed__, warning)
         rav_i = self._highest_parent_._raveled_index_for(self)
         self._highest_parent_._set_fixed(rav_i)
+        self.notify_observers(self, None if trigger_parent else -np.inf)
     fix = constrain_fixed
 
     def unconstrain_fixed(self):
@@ -435,10 +437,10 @@ class Constrainable(Nameable, Indexable):
         Constrain the parameter to the given
         :py:class:`GPy.core.transformations.Transformation`.
         """
-        if isinstance(transform, 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)
 
     def unconstrain(self, *transforms):
         """
@@ -535,7 +537,7 @@ class Constrainable(Nameable, Indexable):
 
         return removed
 
-class OptimizationHandlable(Constrainable, Observable):
+class OptimizationHandlable(Constrainable):
     """
     This enables optimization handles on an Object as done in GPy 0.4.
 
@@ -568,9 +570,7 @@ class OptimizationHandlable(Constrainable, Observable):
 
     def _trigger_params_changed(self, trigger_parent=True):
         [p._trigger_params_changed(trigger_parent=False) for p in self._parameters_]
-        if trigger_parent: min_priority = None
+        self.notify_observers(None, None if trigger_parent else -np.inf)
-        else: min_priority = -np.inf
-        self.notify_observers(None, min_priority)
 
     def _size_transformed(self):
         return self.size - self.constraints[__fixed__].size
@@ -694,6 +694,10 @@ class Parameterizable(OptimizationHandlable):
         elif pname not in dir(self):
             self.__dict__[pname] = param
             self._added_names_.add(pname)
+        else:
+            print "WARNING: added a parameter with formatted name {}, which is already a member of {} object. Trying to change the parameter name to\n   {}".format(pname, self.__class__, param.name+"_")
+            param.name += "_"
+            self._add_parameter_name(param, ignore_added_names)
 
     def _remove_parameter_name(self, param=None, pname=None):
         assert param is None or pname is None, "can only delete either param by name, or the name of a param"

GPy/examples/coreg_example.py

@@ -0,0 +1,30 @@
+import numpy as np
+import pylab as pb
+import GPy
+pb.ion()
+
+X1 = 100 * np.random.rand(100)[:,None]
+X2 = 100 * np.random.rand(100)[:,None]
+#X1.sort()
+#X2.sort()
+
+Y1 = np.sin(X1/10.) + np.random.rand(100)[:,None]
+Y2 = np.cos(X2/10.) + np.random.rand(100)[:,None]
+
+
+
+
+Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")]
+kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=12,kernels_list=Mlist,name='H')
+
+
+m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern)
+m.optimize()
+
+fig = pb.figure()
+ax0 = fig.add_subplot(211)
+ax1 = fig.add_subplot(212)
+slices = GPy.util.multioutput.get_slices([Y1,Y2])
+m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
+m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
+

GPy/kern/_src/add.py

@@ -3,7 +3,6 @@
 
 import numpy as np
 import itertools
-from ...core.parameterization import Parameterized
 from ...util.caching import Cache_this
 from kern import CombinationKernel
 

GPy/kern/_src/kern.py

@@ -156,7 +156,7 @@ class Kern(Parameterized):
         other.active_dims += self.input_dim
         return self.prod(other)
 
-    def prod(self, other, name=None):
+    def prod(self, other, name='mul'):
         """
         Multiply two kernels (either on the same space, or on the tensor
         product of the input space).
@@ -169,12 +169,12 @@ class Kern(Parameterized):
         """
         assert isinstance(other, Kern), "only kernels can be added to kernels..."
         from prod import Prod
-        kernels = []
+        #kernels = []
-        if isinstance(self, Prod): kernels.extend(self._parameters_)
+        #if isinstance(self, Prod): kernels.extend(self._parameters_)
-        else: kernels.append(self)
+        #else: kernels.append(self)
-        if isinstance(other, Prod): kernels.extend(other._parameters_)
+        #if isinstance(other, Prod): kernels.extend(other._parameters_)
-        else: kernels.append(other)
+        #else: kernels.append(other)
-        return Prod(self, other, name)
+        return Prod([self, other], name)
 
     def _getstate(self):
         """
@@ -195,8 +195,10 @@ class Kern(Parameterized):
 class CombinationKernel(Kern):
     def __init__(self, kernels, name):
         assert all([isinstance(k, Kern) for k in kernels])
+        # make sure the active dimensions of all underlying kernels are covered:
         ma = reduce(lambda a,b: max(a, max(b)), (x.active_dims for x in kernels), 0)
         input_dim = np.r_[0:ma+1]
+        # initialize the kernel with the full input_dim
         super(CombinationKernel, self).__init__(input_dim, name)
         self.add_parameters(*kernels)
 

GPy/kern/_src/kernel_slice_operations.py

@@ -9,17 +9,17 @@ class KernCallsViaSlicerMeta(ParametersChangedMeta):
     def __call__(self, *args, **kw):
         instance = super(ParametersChangedMeta, self).__call__(*args, **kw)
         instance.K = _slice_wrapper(instance, instance.K)
-        instance.Kdiag = _slice_wrapper(instance, instance.Kdiag, True)
+        instance.Kdiag = _slice_wrapper(instance, instance.Kdiag, diag=True)
-        instance.update_gradients_full = _slice_wrapper(instance, instance.update_gradients_full, False, True)
+        instance.update_gradients_full = _slice_wrapper(instance, instance.update_gradients_full, diag=False, derivative=True)
-        instance.update_gradients_diag = _slice_wrapper(instance, instance.update_gradients_diag, True, True)
+        instance.update_gradients_diag = _slice_wrapper(instance, instance.update_gradients_diag, diag=True, derivative=True)
-        instance.gradients_X = _slice_wrapper(instance, instance.gradients_X, False, True)
+        instance.gradients_X = _slice_wrapper(instance, instance.gradients_X, diag=False, derivative=True)
-        instance.gradients_X_diag = _slice_wrapper(instance, instance.gradients_X_diag, True, True)
+        instance.gradients_X_diag = _slice_wrapper(instance, instance.gradients_X_diag, diag=True, derivative=True)
-        instance.psi0 = _slice_wrapper(instance, instance.psi0, False, False)
+        instance.psi0 = _slice_wrapper(instance, instance.psi0, diag=False, derivative=False)
-        instance.psi1 = _slice_wrapper(instance, instance.psi1, False, False)
+        instance.psi1 = _slice_wrapper(instance, instance.psi1, diag=False, derivative=False)
-        instance.psi2 = _slice_wrapper(instance, instance.psi2, False, False)
+        instance.psi2 = _slice_wrapper(instance, instance.psi2, diag=False, derivative=False)
-        instance.update_gradients_expectations = _slice_wrapper(instance, instance.update_gradients_expectations, psi_stat=True)
+        instance.update_gradients_expectations = _slice_wrapper(instance, instance.update_gradients_expectations, derivative=True, psi_stat=True)
-        instance.gradients_Z_expectations = _slice_wrapper(instance, instance.gradients_Z_expectations, psi_stat_Z=True)
+        instance.gradients_Z_expectations = _slice_wrapper(instance, instance.gradients_Z_expectations, derivative=True, psi_stat_Z=True)
-        instance.gradients_qX_expectations = _slice_wrapper(instance, instance.gradients_qX_expectations, psi_stat=True)
+        instance.gradients_qX_expectations = _slice_wrapper(instance, instance.gradients_qX_expectations, derivative=True, psi_stat=True)
         instance.parameters_changed()
         return instance
 
@@ -44,6 +44,28 @@ def _slice_wrapper(kern, operation, diag=False, derivative=False, psi_stat=False
                 finally:
                     kern._sliced_X -= 1
                 return ret
+        elif psi_stat:
+            def x_slice_wrapper(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
+                Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior
+                kern._sliced_X += 1
+                try:
+                    ret = operation(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior)
+                except:
+                    raise
+                finally:
+                    kern._sliced_X -= 1
+                return ret
+        elif psi_stat_Z:
+            def x_slice_wrapper(dL_dpsi1, dL_dpsi2, Z, variational_posterior):
+                Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior
+                kern._sliced_X += 1
+                try:
+                    ret = operation(dL_dpsi1, dL_dpsi2, Z, variational_posterior)
+                except:
+                    raise
+                finally:
+                    kern._sliced_X -= 1
+                return ret
         else:
             def x_slice_wrapper(dL_dK, X, X2=None):
                 X, X2 = kern._slice_X(X) if not kern._sliced_X else X, kern._slice_X(X2) if X2 is not None and not kern._sliced_X else X2
@@ -55,28 +77,6 @@ def _slice_wrapper(kern, operation, diag=False, derivative=False, psi_stat=False
                 finally:
                     kern._sliced_X -= 1
                 return ret
-    elif psi_stat:
-        def x_slice_wrapper(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
-            Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior
-            kern._sliced_X += 1
-            try:
-                ret = operation(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior)
-            except:
-                raise
-            finally:
-                kern._sliced_X -= 1
-            return ret
-    elif psi_stat_Z:
-        def x_slice_wrapper(dL_dpsi1, dL_dpsi2, Z, variational_posterior):
-            Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior
-            kern._sliced_X += 1
-            try:
-                ret = operation(dL_dpsi1, dL_dpsi2, Z, variational_posterior)
-            except:
-                raise
-            finally:
-                kern._sliced_X -= 1
-            return ret
     else:
         if diag:
             def x_slice_wrapper(X, *args, **kw):

GPy/kern/_src/mlp.py

@@ -96,12 +96,12 @@ class MLP(Kern):
             vec = (X*X).sum(1)*self.weight_variance+self.bias_variance + 1.
             return 2*four_over_tau*self.weight_variance*self.variance*((X[None, :, :]/denom[:, :, None] - vec[None, :, None]*X[:, None, :]*(numer/denom3)[:, :, None])*(dL_dK/np.sqrt(1-arg*arg))[:, :, None]).sum(1)
 
-    def dKdiag_dX(self, dL_dKdiag, X, target):
+    def gradients_X_diag(self, dL_dKdiag, X):
         """Gradient of diagonal of covariance with respect to X"""
         self._K_diag_computations(X)
         arg = self._K_diag_asin_arg
         denom = self._K_diag_denom
-        numer = self._K_diag_numer
+        #numer = self._K_diag_numer
         return four_over_tau*2.*self.weight_variance*self.variance*X*(1./denom*(1. - arg)*dL_dKdiag/(np.sqrt(1-arg*arg)))[:, None]
 
 

GPy/kern/_src/periodic.py

@@ -85,8 +85,9 @@ class PeriodicExponential(Periodic):
         self.b = [1]
 
         self.basis_alpha = np.ones((self.n_basis,))
-        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))[:,0]
+        self.basis_omega = (2*np.pi*np.arange(1,self.n_freq+1)/self.period).repeat(2)
-        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
+        self.basis_phi =   np.zeros(self.n_freq * 2)
+        self.basis_phi[::2] = -np.pi/2
 
         self.G = self.Gram_matrix()
         self.Gi = np.linalg.inv(self.G)
@@ -100,7 +101,6 @@ class PeriodicExponential(Periodic):
         Flower = np.array(self._cos(self.basis_alpha,self.basis_omega,self.basis_phi)(self.lower))[:,None]
         return(self.lengthscale/(2*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T))
 
-    #@silence_errors
     def update_gradients_full(self, dL_dK, X, X2=None):
         """derivative of the covariance matrix with respect to the parameters (shape is N x num_inducing x num_params)"""
         if X2 is None: X2 = X
@@ -194,8 +194,9 @@ class PeriodicMatern32(Periodic):
         self.b = [1,self.lengthscale**2/3]
 
         self.basis_alpha = np.ones((self.n_basis,))
-        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))
+        self.basis_omega = (2*np.pi*np.arange(1,self.n_freq+1)/self.period).repeat(2)
-        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
+        self.basis_phi =   np.zeros(self.n_freq * 2)
+        self.basis_phi[::2] = -np.pi/2
 
         self.G = self.Gram_matrix()
         self.Gi = np.linalg.inv(self.G)
@@ -212,8 +213,8 @@ class PeriodicMatern32(Periodic):
         return(self.lengthscale**3/(12*np.sqrt(3)*self.variance) * Gint + 1./self.variance*np.dot(Flower,Flower.T) + self.lengthscale**2/(3.*self.variance)*np.dot(F1lower,F1lower.T))
 
 
-    @silence_errors
+    #@silence_errors
-    def update_gradients_full(self,dL_dK,X,X2,target):
+    def update_gradients_full(self,dL_dK,X,X2):
         """derivative of the covariance matrix with respect to the parameters (shape is num_data x num_inducing x num_params)"""
         if X2 is None: X2 = X
         FX  = self._cos(self.basis_alpha[None,:],self.basis_omega[None,:],self.basis_phi[None,:])(X)
@@ -307,8 +308,9 @@ class PeriodicMatern52(Periodic):
         self.b  = [9./8, 9*self.lengthscale**4/200., 3*self.lengthscale**2/5., 3*self.lengthscale**2/(5*8.), 3*self.lengthscale**2/(5*8.)]
 
         self.basis_alpha = np.ones((2*self.n_freq,))
-        self.basis_omega = np.array(sum([[i*2*np.pi/self.period]*2 for i in  range(1,self.n_freq+1)],[]))
+        self.basis_omega = (2*np.pi*np.arange(1,self.n_freq+1)/self.period).repeat(2)
-        self.basis_phi =   np.array(sum([[-np.pi/2, 0.]  for i in range(1,self.n_freq+1)],[]))
+        self.basis_phi =   np.zeros(self.n_freq * 2)
+        self.basis_phi[::2] = -np.pi/2
 
         self.G = self.Gram_matrix()
         self.Gi = np.linalg.inv(self.G)

GPy/kern/_src/prod.py

@@ -17,7 +17,8 @@ class Prod(CombinationKernel):
     :rtype: kernel object
 
     """
-    def __init__(self, kernels, name='prod'):
+    def __init__(self, kernels, name='mul'):
+        assert len(kernels) == 2, 'only implemented for two kernels as of yet'
         super(Prod, self).__init__(kernels, name)
 
     @Cache_this(limit=2, force_kwargs=['which_parts'])
@@ -37,26 +38,28 @@ class Prod(CombinationKernel):
             which_parts = self.parts
         return reduce(np.multiply, (p.Kdiag(X) for p in which_parts))
 
-    def update_gradients_full(self, dL_dK, X):
+    def update_gradients_full(self, dL_dK, X, X2=None):
         for k1,k2 in itertools.combinations(self.parts, 2):
-            k1._sliced_X = k1._sliced_X2 = k2._sliced_X = k2._sliced_X2 = True
+            k1.update_gradients_full(dL_dK*k2.K(X, X2), X, X2)
-            k1.update_gradients_full(dL_dK*k2.K(X, X))
+            k2.update_gradients_full(dL_dK*k1.K(X, X2), X, X2)
-            self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2])
+
+    def update_gradients_diag(self, dL_dKdiag, X):
+        for k1,k2 in itertools.combinations(self.parts, 2):
+            k1.update_gradients_diag(dL_dKdiag*k2.Kdiag(X), X)
+            k2.update_gradients_diag(dL_dKdiag*k1.Kdiag(X), X)
 
     def gradients_X(self, dL_dK, X, X2=None):
         target = np.zeros(X.shape)
-        if X2 is None:
+        for k1,k2 in itertools.combinations(self.parts, 2):
-            target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2.K(X[:,self.slice2]), X[:,self.slice1], None)
+            target[:,k1.active_dims] += k1.gradients_X(dL_dK*k2.K(X, X2), X, X2)
-            target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2], None)
+            target[:,k2.active_dims] += k2.gradients_X(dL_dK*k1.K(X, X2), X, X2)
-        else:
-            target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2.K(X[:,self.slice2], X2[:,self.slice2]), X[:,self.slice1], X2[:,self.slice1])
-            target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1.K(X[:,self.slice1], X2[:,self.slice1]), X[:,self.slice2], X2[:,self.slice2])
         return target
 
     def gradients_X_diag(self, dL_dKdiag, X):
         target = np.zeros(X.shape)
-        target[:,self.slice1] = self.k1.gradients_X(dL_dKdiag*self.k2.Kdiag(X[:,self.slice2]), X[:,self.slice1])
+        for k1,k2 in itertools.combinations(self.parts, 2):
-        target[:,self.slice2] += self.k2.gradients_X(dL_dKdiag*self.k1.Kdiag(X[:,self.slice1]), X[:,self.slice2])
+            target[:,k1.active_dims] += k1.gradients_X(dL_dKdiag*k2.Kdiag(X), X)
+            target[:,k2.active_dims] += k2.gradients_X(dL_dKdiag*k1.Kdiag(X), X)
         return target
 
 

GPy/kern/_src/sympykern.py

@@ -116,6 +116,7 @@ class Sympykern(Kern):
         if self.output_dim > 1:
             self.arg_list += self._sp_theta_i + self._sp_theta_j
             self.diag_arg_list += self._sp_theta_i
+
         # psi_stats aren't yet implemented.
         if False:
             self.compute_psi_stats()

GPy/models/mrd.py

@@ -45,7 +45,6 @@ class MRD(Model):
     :param str name: the name of this model
     :param [str] Ynames: the names for the datasets given, must be of equal length as Ylist or None
     """
-    
     def __init__(self, Ylist, input_dim, X=None, X_variance=None,
                  initx = 'PCA', initz = 'permute',
                  num_inducing=10, Z=None, kernel=None,

GPy/testing/bgplvm_tests.py

@@ -1,85 +0,0 @@
-# Copyright (c) 2012, Nicolo Fusi
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-
-import unittest
-import numpy as np
-import GPy
-from ..models import BayesianGPLVM
-
-class BGPLVMTests(unittest.TestCase):
-    def test_bias_kern(self):
-        N, num_inducing, input_dim, D = 10, 3, 2, 4
-        X = np.random.rand(N, input_dim)
-        k = GPy.kern.RBF(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        K = k.K(X)
-        Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
-        Y -= Y.mean(axis=0)
-        k = GPy.kern.bias(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
-        m.randomize()
-        self.assertTrue(m.checkgrad())
-
-    def test_linear_kern(self):
-        N, num_inducing, input_dim, D = 10, 3, 2, 4
-        X = np.random.rand(N, input_dim)
-        k = GPy.kern.RBF(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        K = k.K(X)
-        Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
-        Y -= Y.mean(axis=0)
-        k = GPy.kern.Linear(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
-        m.randomize()
-        self.assertTrue(m.checkgrad())
-
-    def test_rbf_kern(self):
-        N, num_inducing, input_dim, D = 10, 3, 2, 4
-        X = np.random.rand(N, input_dim)
-        k = GPy.kern.RBF(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        K = k.K(X)
-        Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
-        Y -= Y.mean(axis=0)
-        k = GPy.kern.RBF(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
-        m.randomize()
-        self.assertTrue(m.checkgrad())
-
-    def test_rbf_bias_kern(self):
-        N, num_inducing, input_dim, D = 10, 3, 2, 4
-        X = np.random.rand(N, input_dim)
-        k = GPy.kern.RBF(input_dim) +  GPy.kern.Bias(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        K = k.K(X)
-        Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
-        Y -= Y.mean(axis=0)
-        k = GPy.kern.RBF(input_dim) + GPy.kern.Bias(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
-        m.randomize()
-        self.assertTrue(m.checkgrad())
-
-    def test_rbf_line_kern(self):
-        N, num_inducing, input_dim, D = 10, 3, 2, 4
-        X = np.random.rand(N, input_dim)
-        k = GPy.kern.RBF(input_dim) +  GPy.kern.Linear(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        K = k.K(X)
-        Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
-        Y -= Y.mean(axis=0)
-        k = GPy.kern.RBF(input_dim) + GPy.kern.Bias(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
-        m.randomize()
-        self.assertTrue(m.checkgrad())
-
-    def test_linear_bias_kern(self):
-        N, num_inducing, input_dim, D = 30, 5, 4, 30
-        X = np.random.rand(N, input_dim)
-        k = GPy.kern.Linear(input_dim) +  GPy.kern.Bias(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        K = k.K(X)
-        Y = np.random.multivariate_normal(np.zeros(N),K,input_dim).T
-        Y -= Y.mean(axis=0)
-        k = GPy.kern.Linear(input_dim) + GPy.kern.Bias(input_dim) + GPy.kern.White(input_dim, 0.00001)
-        m = BayesianGPLVM(Y, input_dim, kernel=k, num_inducing=num_inducing)
-        m.randomize()
-        self.assertTrue(m.checkgrad())
-
-
-if __name__ == "__main__":
-    print "Running unit tests, please be (very) patient..."
-    unittest.main()

GPy/testing/kernel_tests.py

@@ -33,9 +33,10 @@ class Kern_check_model(GPy.core.Model):
         self.X2 = X2
         self.dL_dK = dL_dK
 
-    def is_positive_definite(self):
+    def is_positive_semi_definite(self):
         v = np.linalg.eig(self.kernel.K(self.X))[0]
-        if any(v<-10*sys.float_info.epsilon):
+        if any(v.real<=-1e-10):
+            print v.real.min()
             return False
         else:
             return True
@@ -89,7 +90,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, self.X)
+        self.X.gradient =  self.kernel.gradients_X_diag(self.dL_dK.diagonal(), self.X)
 
 
 
@@ -119,7 +120,7 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb
 
     if verbose:
         print("Checking covariance function is positive definite.")
-    result = Kern_check_model(kern, X=X).is_positive_definite()
+    result = Kern_check_model(kern, X=X).is_positive_semi_definite()
     if result and verbose:
         print("Check passed.")
     if not result:
@@ -214,21 +215,67 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb
 
 class KernelGradientTestsContinuous(unittest.TestCase):
     def setUp(self):
-        self.X = np.random.randn(100,2)
+        self.N, self.D = 100, 5
-        self.X2 = np.random.randn(110,2)
+        self.X = np.random.randn(self.N,self.D)
+        self.X2 = np.random.randn(self.N+10,self.D)
 
         continuous_kerns = ['RBF', 'Linear']
         self.kernclasses = [getattr(GPy.kern, s) for s in continuous_kerns]
 
     def test_Matern32(self):
-        k = GPy.kern.Matern32(2)
+        k = GPy.kern.Matern32(self.D)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+    def test_Prod(self):
+        k = GPy.kern.Matern32([2,3]) * GPy.kern.RBF([0,4]) + GPy.kern.Linear(self.D)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+    def test_Add(self):
+        k = GPy.kern.Matern32([2,3]) + GPy.kern.RBF([0,4]) + GPy.kern.Linear(self.D)
+        k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
     def test_Matern52(self):
-        k = GPy.kern.Matern52(2)
+        k = GPy.kern.Matern52(self.D)
+        k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
-    #TODO: turn off grad checkingwrt X for indexed kernels liek coregionalize
+    def test_RBF(self):
+        k = GPy.kern.RBF(self.D)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+    def test_Linear(self):
+        k = GPy.kern.Linear(self.D)
+        k.randomize()
+        self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+
+#TODO: turn off grad checkingwrt X for indexed kernels liek coregionalize
+# class KernelGradientTestsContinuous1D(unittest.TestCase):
+#     def setUp(self):
+#         self.N, self.D = 100, 1
+#         self.X = np.random.randn(self.N,self.D)
+#         self.X2 = np.random.randn(self.N+10,self.D)
+# 
+#         continuous_kerns = ['RBF', 'Linear']
+#         self.kernclasses = [getattr(GPy.kern, s) for s in continuous_kerns]
+# 
+#     def test_PeriodicExponential(self):
+#         k = GPy.kern.PeriodicExponential(self.D)
+#         k.randomize()
+#         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+# 
+#     def test_PeriodicMatern32(self):
+#         k = GPy.kern.PeriodicMatern32(self.D)
+#         k.randomize()
+#         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
+# 
+#     def test_PeriodicMatern52(self):
+#         k = GPy.kern.PeriodicMatern52(self.D)
+#         k.randomize()
+#         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
 
 class KernelTestsMiscellaneous(unittest.TestCase):
@@ -237,7 +284,7 @@ class KernelTestsMiscellaneous(unittest.TestCase):
         N, D = 100, 10
         self.X = np.linspace(-np.pi, +np.pi, N)[:,None] * np.ones(D)
         self.rbf = GPy.kern.RBF(range(2))
-        self.linear = GPy.kern.Linear((3,5,6))
+        self.linear = GPy.kern.Linear((3,6))
         self.matern = GPy.kern.Matern32(np.array([2,4,7]))
         self.sumkern = self.rbf + self.linear
         self.sumkern += self.matern
@@ -251,6 +298,7 @@ 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)))
 
+
 if __name__ == "__main__":
     print "Running unit tests, please be (very) patient..."
     unittest.main()

GPy/testing/likelihood_tests.py

@@ -541,7 +541,8 @@ class TestNoiseModels(object):
             #import ipdb; ipdb.set_trace()
             #NOTE this test appears to be stochastic for some likelihoods (student t?)
             # appears to all be working in test mode right now...
-
+        #if isinstance(model, GPy.likelihoods.StudentT):
+        #    import ipdb;ipdb.set_trace()
         assert m.checkgrad(step=step)
 
     ###########
@@ -700,7 +701,6 @@ class LaplaceTests(unittest.TestCase):
         np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2)
         #Check marginals are the same with random
         m1.randomize()
-        import ipdb;ipdb.set_trace()
         m2[:] = m1[:]
 
         np.testing.assert_almost_equal(m1.log_likelihood(), m2.log_likelihood(), decimal=2)

GPy/testing/mrd_tests.py

@@ -1,32 +0,0 @@
-# Copyright (c) 2013, Max Zwiessele
-# Licensed under the BSD 3-clause license (see LICENSE.txt)
-'''
-Created on 10 Apr 2013
-
-@author: maxz
-'''
-
-import unittest
-import numpy as np
-import GPy
-
-class MRDTests(unittest.TestCase):
-
-    def test_gradients(self):
-        num_m = 3
-        N, num_inducing, input_dim, D = 20, 8, 6, 20
-        X = np.random.rand(N, input_dim)
-
-        k = GPy.kern.linear(input_dim) + GPy.kern.bias(input_dim) + GPy.kern.white(input_dim)
-        K = k.K(X)
-
-        Ylist = [np.random.multivariate_normal(np.zeros(N), K, input_dim).T for _ in range(num_m)]
-        likelihood_list = [GPy.likelihoods.Gaussian(Y) for Y in Ylist]
-
-        m = GPy.models.MRD(likelihood_list, input_dim=input_dim, kernels=k, num_inducing=num_inducing)
-
-        self.assertTrue(m.checkgrad())
-
-if __name__ == "__main__":
-    print "Running unit tests, please be (very) patient..."
-    unittest.main()

GPy/testing/unit_tests.py

@@ -34,7 +34,7 @@ class GradientTests(unittest.TestCase):
         model_fit = getattr(GPy.models, model_type)
 
         # noise = GPy.kern.White(dimension)
-        kern = kern #  + noise
+        kern = kern  #  + noise
         if uncertain_inputs:
             m = model_fit(X, Y, kernel=kern, X_variance=np.random.rand(X.shape[0], X.shape[1]))
         else:
@@ -60,13 +60,14 @@ class GradientTests(unittest.TestCase):
 
     def test_GPRegression_mlp_1d(self):
         ''' Testing the GP regression with mlp kernel with white kernel on 1d data '''
-        mlp = GPy.kern.mlp(1)
+        mlp = GPy.kern.MLP(1)
         self.check_model(mlp, model_type='GPRegression', dimension=1)
 
-    def test_GPRegression_poly_1d(self):
+    #TODO:
-        ''' Testing the GP regression with polynomial kernel with white kernel on 1d data '''
+    #def test_GPRegression_poly_1d(self):
-        mlp = GPy.kern.Poly(1, degree=5)
+    #    ''' Testing the GP regression with polynomial kernel with white kernel on 1d data '''
-        self.check_model(mlp, model_type='GPRegression', dimension=1)
+    #    mlp = GPy.kern.Poly(1, degree=5)
+    #    self.check_model(mlp, model_type='GPRegression', dimension=1)
 
     def test_GPRegression_matern52_1D(self):
         ''' Testing the GP regression with matern52 kernel on 1d data '''
@@ -163,14 +164,14 @@ class GradientTests(unittest.TestCase):
         rbflin = GPy.kern.RBF(2) + GPy.kern.Linear(2)
         self.check_model(rbflin, model_type='SparseGPRegression', dimension=2)
 
-    #@unittest.expectedFailure
+    # @unittest.expectedFailure
     def test_SparseGPRegression_rbf_linear_white_kern_2D_uncertain_inputs(self):
         ''' Testing the sparse GP regression with rbf, linear kernel on 2d data with uncertain inputs'''
         rbflin = GPy.kern.RBF(2) + GPy.kern.Linear(2)
         raise unittest.SkipTest("This is not implemented yet!")
         self.check_model(rbflin, model_type='SparseGPRegression', dimension=2, uncertain_inputs=1)
 
-    #@unittest.expectedFailure
+    # @unittest.expectedFailure
     def test_SparseGPRegression_rbf_linear_white_kern_1D_uncertain_inputs(self):
         ''' Testing the sparse GP regression with rbf, linear kernel on 1d data with uncertain inputs'''
         rbflin = GPy.kern.RBF(1) + GPy.kern.Linear(1)
@@ -202,7 +203,7 @@ class GradientTests(unittest.TestCase):
         X = np.hstack([np.random.normal(5, 2, N / 2), np.random.normal(10, 2, N / 2)])[:, None]
         Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None]
         kernel = GPy.kern.RBF(1)
-        m = GPy.models.GPClassification(X,Y,kernel=kernel)
+        m = GPy.models.GPClassification(X, Y, kernel=kernel)
         m.update_likelihood_approximation()
         self.assertTrue(m.checkgrad())
 
@@ -212,11 +213,11 @@ class GradientTests(unittest.TestCase):
         Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None]
         Z = np.linspace(0, 15, 4)[:, None]
         kernel = GPy.kern.RBF(1)
-        m = GPy.models.SparseGPClassification(X,Y,kernel=kernel,Z=Z)
+        m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, Z=Z)
-        #distribution = GPy.likelihoods.likelihood_functions.Bernoulli()
+        # distribution = GPy.likelihoods.likelihood_functions.Bernoulli()
-        #likelihood = GPy.likelihoods.EP(Y, distribution)
+        # likelihood = GPy.likelihoods.EP(Y, distribution)
-        #m = GPy.core.SparseGP(X, likelihood, kernel, Z)
+        # m = GPy.core.SparseGP(X, likelihood, kernel, Z)
-        #m.ensure_default_constraints()
+        # m.ensure_default_constraints()
         m.update_likelihood_approximation()
         self.assertTrue(m.checkgrad())
 
@@ -224,8 +225,8 @@ class GradientTests(unittest.TestCase):
         N = 20
         X = np.hstack([np.random.rand(N / 2) + 1, np.random.rand(N / 2) - 1])[:, None]
         k = GPy.kern.RBF(1) + GPy.kern.White(1)
-        Y = np.hstack([np.ones(N/2),np.zeros(N/2)])[:,None]
+        Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None]
-        m = GPy.models.FITCClassification(X, Y, kernel = k)
+        m = GPy.models.FITCClassification(X, Y, kernel=k)
         m.update_likelihood_approximation()
         self.assertTrue(m.checkgrad())
 
@@ -238,7 +239,7 @@ class GradientTests(unittest.TestCase):
         Y = np.vstack((Y1, Y2))
 
         k1 = GPy.kern.RBF(1)
-        m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
+        m = GPy.models.GPMultioutputRegression(X_list=[X1, X2], Y_list=[Y1, Y2], kernel_list=[k1])
         m.constrain_fixed('.*rbf_var', 1.)
         self.assertTrue(m.checkgrad())
 
@@ -251,7 +252,7 @@ class GradientTests(unittest.TestCase):
         Y = np.vstack((Y1, Y2))
 
         k1 = GPy.kern.RBF(1)
-        m = GPy.models.SparseGPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
+        m = GPy.models.SparseGPMultioutputRegression(X_list=[X1, X2], Y_list=[Y1, Y2], kernel_list=[k1])
         m.constrain_fixed('.*rbf_var', 1.)
         self.assertTrue(m.checkgrad())
 

GPy/util/caching.py

@@ -52,29 +52,33 @@ class Cacher(object):
 
         #if the result is cached, return the cached computation
         state = [all(a is b for a, b in itertools.izip_longest(args, cached_i)) for cached_i in self.cached_inputs]
-        if any(state):
+        try:
-            i = state.index(True)
+            if any(state):
-            if self.inputs_changed[i]:
+                i = state.index(True)
-                #(elements of) the args have changed since we last computed: update
+                if self.inputs_changed[i]:
-                self.cached_outputs[i] = self.operation(*args, **kw)
+                    #(elements of) the args have changed since we last computed: update
-                self.inputs_changed[i] = False
+                    self.cached_outputs[i] = self.operation(*args, **kw)
-            return self.cached_outputs[i]
+                    self.inputs_changed[i] = False
-        else:
+                return self.cached_outputs[i]
-            #first time we've seen these arguments: compute
+            else:
+                #first time we've seen these arguments: compute
 
-            #first make sure the depth limit isn't exceeded
+                #first make sure the depth limit isn't exceeded
-            if len(self.cached_inputs) == self.limit:
+                if len(self.cached_inputs) == self.limit:
-                args_ = self.cached_inputs.pop(0)
+                    args_ = self.cached_inputs.pop(0)
-                [a.remove_observer(self, self.on_cache_changed) for a in args_ if a is not None]
+                    [a.remove_observer(self, self.on_cache_changed) for a in args_ if a is not None]
-                self.inputs_changed.pop(0)
+                    self.inputs_changed.pop(0)
-                self.cached_outputs.pop(0)
+                    self.cached_outputs.pop(0)
-
+                #compute
-            #compute
+                self.cached_inputs.append(oa_all)
-            self.cached_inputs.append(oa_all)
+                self.cached_outputs.append(self.operation(*args, **kw))
-            self.cached_outputs.append(self.operation(*args, **kw))
+                self.inputs_changed.append(False)
-            self.inputs_changed.append(False)
+                [a.add_observer(self, self.on_cache_changed) for a in observable_args]
-            [a.add_observer(self, self.on_cache_changed) for a in observable_args]
+                return self.cached_outputs[-1]#return
-            return self.cached_outputs[-1]#return
+        except:
+            raise
+        finally:
+            self.reset()
 
     def on_cache_changed(self, arg):
         """
@@ -84,7 +88,7 @@ class Cacher(object):
         """
         self.inputs_changed = [any([a is arg for a in args]) or old_ic for args, old_ic in zip(self.cached_inputs, self.inputs_changed)]
 
-    def reset(self, obj):
+    def reset(self):
         """
         Totally reset the cache
         """

GPy/util/data_resources.json

@@ -32,6 +32,33 @@
       "details":"Artificially generated data of silhouettes given poses. Note that the data does not display a left/right ambiguity because across the entire data set one of the arms sticks out more the the other, disambiguating the pose as to which way the individual is facing.",
       "size":1
    },
+   "football_data":{
+      "files":[
+         [
+	     "E0.csv", "E1.csv", "E2.csv", "E3.csv"
+         ]
+      ],
+      "citation":"",
+      "license":null,
+      "urls":[
+          "http://www.football-data.co.uk/mmz4281/"
+      ],
+      "details":"Results of English football matches since 1993/94 season.",
+      "size":1
+   },
+   "google_trends":{
+      "files":[
+         [
+         ]
+      ],
+      "citation":"",
+      "license":null,
+      "urls":[
+          "http://www.google.com/trends/"
+      ],
+      "details":"Google trends results.",
+      "size":0
+   },
    "osu_accad":{
       "files":[
          [

GPy/util/datasets.py

@@ -1,5 +1,8 @@
+import csv
 import os
+import copy
 import numpy as np
+import pylab as pb
 import GPy
 import scipy.io
 import cPickle as pickle
@@ -7,6 +10,8 @@ import zipfile
 import tarfile
 import datetime
 import json
+import re
+
 ipython_available=True
 try:
     import IPython
@@ -32,11 +37,18 @@ neil_url = 'http://staffwww.dcs.shef.ac.uk/people/N.Lawrence/dataset_mirror/'
 # Read data resources from json file.
 # Don't do this when ReadTheDocs is scanning as it breaks things
 on_rtd = os.environ.get('READTHEDOCS', None) == 'True' #Checks if RTD is scanning
+
 if not (on_rtd):
     path = os.path.join(os.path.dirname(__file__), 'data_resources.json')
     json_data=open(path).read()
     data_resources = json.loads(json_data)
 
+if not (on_rtd):
+    path = os.path.join(os.path.dirname(__file__), 'football_teams.json')
+    json_data=open(path).read()
+    football_dict = json.loads(json_data)
+
+    
 
 def prompt_user(prompt):
     """Ask user for agreeing to data set licenses."""
@@ -276,6 +288,74 @@ def della_gatta_TRP63_gene_expression(data_set='della_gatta', gene_number=None):
 
     
 
+def football_data(season='1314', data_set='football_data'):
+    """Football data from English games since 1993. This downloads data from football-data.co.uk for the given season. """
+    def league2num(string):
+        league_dict = {'E0':0, 'E1':1, 'E2': 2, 'E3': 3, 'EC':4}
+        return league_dict[string]
+
+    def football2num(string):
+        if football_dict.has_key(string):
+            return football_dict[string]
+        else:
+            football_dict[string] = len(football_dict)+1
+            return len(football_dict)+1
+
+    data_set_season = data_set + '_' + season
+    data_resources[data_set_season] = copy.deepcopy(data_resources[data_set])
+    data_resources[data_set_season]['urls'][0]+=season + '/'
+    start_year = int(year[0:2])
+    end_year = int(year[2:4])
+    files = ['E0.csv', 'E1.csv', 'E2.csv', 'E3.csv']
+    if start_year>4 and start_year < 93:
+        files += ['EC.csv']
+    data_resources[data_set_season]['files'] = [files]
+    if not data_available(data_set_season):
+        download_data(data_set_season)
+    for file in reversed(files):
+        filename = os.path.join(data_path, data_set_season, file)
+        # rewrite files removing blank rows.
+        writename = os.path.join(data_path, data_set_season, 'temp.csv')
+        input = open(filename, 'rb')
+        output = open(writename, 'wb')
+        writer = csv.writer(output)
+        for row in csv.reader(input):
+            if any(field.strip() for field in row):
+                writer.writerow(row)
+        input.close()
+        output.close()
+        table = np.loadtxt(writename,skiprows=1, usecols=(0, 1, 2, 3, 4, 5), converters = {0: league2num, 1: pb.datestr2num, 2:football2num, 3:football2num}, delimiter=',')
+        X = table[:, :4]
+        Y = table[:, 4:]
+    return data_details_return({'X': X, 'Y': Y}, 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'):
+    """Data downloaded from Google trends for given query terms."""
+    # Inspired by this notebook:
+    # http://nbviewer.ipython.org/github/sahuguet/notebooks/blob/master/GoogleTrends%20meet%20Notebook.ipynb
+
+    # quote the query terms.
+    for i, element in enumerate(query_terms):
+        query_terms[i] = urllib2.quote(element)
+    query = 'http://www.google.com/trends/fetchComponent?q=%s&cid=TIMESERIES_GRAPH_0&export=3' % ",".join(query_terms)
+
+    data = urllib2.urlopen(query).read()
+
+    # In the notebook they did some data cleaning: remove Javascript header+footer, and translate new Date(....,..,..) into YYYY-MM-DD.
+    header = """// Data table response\ngoogle.visualization.Query.setResponse("""
+    data = data[len(header):-2]
+    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)
+    timeseries = json.loads(data)
+    #import pandas as pd
+    columns = [k['label'] for k in timeseries['table']['cols']]
+    rows = map(lambda x: [k['v'] for k in x['c']], timeseries['table']['rows'])
+    terms = len(columns)-1
+    X = np.asarray([(pb.datestr2num(row[0]), i) for i in range(terms) for row in rows ])
+    Y = np.asarray([[row[i+1]] for i in range(terms) for row in rows ])
+    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)
+    
 # The data sets
 def oil(data_set='three_phase_oil_flow'):
     """The three phase oil data from Bishop and James (1993)."""

GPy/util/football_teams.json

@@ -0,0 +1 @@
+{"Canvey Island": 94, "Crewe": 21, "Fleetwood Town": 134, "Wrexham": 89, "Barnet": 69, "Ipswich": 29, "Rochdale": 84, "Bristol Rvs": 70, "Liverpool": 10, "Chelsea": 20, "York": 113, "Newcastle": 18, "QPR": 28, "Middlesboro": 116, "Tranmere": 68, "Bury": 72, "Luton": 24, "AFC Wimbledon": 126, "West Ham": 15, "Braintree Town": 135, "Bournemouth": 58, "Hayes & Yeading": 130, "Rushden & D": 81, "Weymouth": 120, "Chesterfield": 48, "Exeter": 104, "Barnsley": 45, "Aldershot": 95, "Gateshead": 129, "Hartlepool": 55, "Newport County": 132, "Crystal Palace": 23, "Ebbsfleet": 123, "Wigan": 19, "Shrewsbury": 83, "Hereford": 105, "Stevenage": 111, "Grimsby": 73, "Crawley Town": 114, "Morecambe": 109, "Oldham": 61, "Aston Villa": 1, "Bristol City": 51, "Gravesend": 103, "Huddersfield": 60, "Reading": 33, "Nuneaton Town": 140, "AFC Telford United": 137, "Wycombe": 91, "Leeds": 43, "Colchester": 54, "Rotherham": 63, "Southport": 100, "Southampton": 37, "Darlington": 82, "Blackburn": 16, "Bath City": 133, "Yeovil": 62, "Leyton Orient": 75, "Forest Green": 101, "Chester": 80, "Halifax": 110, "Portsmouth": 11, "Woking": 108, "Histon": 125, "Man City": 7, "Northampton": 78, "Arsenal": 17, "Charlton": 14, "Middlesbrough": 9, "Watford": 41, "Nott'm Forest": 59, "Eastbourne Borough": 131, "Hull": 27, "Barrow": 127, "Doncaster": 52, "Carlisle": 92, "Gillingham": 53, "Accrington": 93, "Dartford": 139, "Altrincham": 112, "Scarborough": 106, "Northwich": 117, "Farsley": 124, "Tamworth": 96, "St. Albans": 119, "Alfreton Town": 136, "Mansfield": 86, "Macclesfield": 76, "Torquay": 87, "Brighton": 26, "Bradford": 56, "Lincoln": 77, "Brentford": 49, "Everton": 3, "Cambridge": 102, "Sheffield United": 35, "Stockport": 85, "Bolton": 2, "Southend": 65, "Cheltenham": 71, "Walsall": 64, "Preston": 42, "Peterboro": 79, "Birmingham": 6, "Boston": 90, "Burton": 97, "West Brom": 8, "Man United": 4, "Stafford Rangers": 118, "Wimbledon": 115, "Scunthorpe": 50, "Kidderminster": 107, "Millwall": 44, "Swansea": 67, "Norwich": 31, "Burnley": 22, "Sunderland": 13, "Sheffield Weds": 40, "Fulham": 5, "Dag and Red": 99, "Oxford": 74, "Stoke": 39, "Tottenham": 12, "Kettering Town": 128, "Coventry": 32, "Wolves": 38, "Port Vale": 66, "Milton Keynes Dons": 57, "Plymouth": 34, "Derby": 25, "Notts County": 88, "Leicester": 36, "Droylsden": 121, "Blackpool": 47, "Salisbury": 122, "Cardiff": 30, "Grays": 98, "Swindon": 46, "Hyde United": 138}

GPy/util/multioutput.py

@@ -54,8 +54,8 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'):
         kernel.input_dim = input_dim
         warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
 
-    #K = kernel.prod(GPy.kern.Coregionalize(input_dim, num_outputs,W_rank,W,kappa,name='B'),tensor=True,name=name)
+    K = kernel.prod(GPy.kern.Coregionalize([input_dim], num_outputs,W_rank,W,kappa,name='B'),name=name)
-    K = kernel.prod(GPy.kern.Coregionalize(input_dim, num_outputs,W_rank,W,kappa,name='B'),name=name)
+    #K = kernel ** GPy.kern.Coregionalize(input_dim, num_outputs,W_rank,W,kappa, name= 'B')
     K['.*variance'] = 1.
     K['.*variance'].fix()
     return K