diff --git a/GPy/core/model.py b/GPy/core/model.py index 710c1b22..c2a9ed23 100644 --- a/GPy/core/model.py +++ b/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]) diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 8eb10608..cad20a8a 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -94,15 +94,15 @@ class Param(OptimizationHandlable, ObservableArray): @property def _param_array_(self): return self - + @property def gradient(self): return self._gradient_array_[self._current_slice_] - + @gradient.setter def gradient(self, val): self.gradient[:] = val - + #=========================================================================== # Pickling operations #=========================================================================== @@ -135,7 +135,7 @@ class Param(OptimizationHandlable, ObservableArray): self._parent_index_ = state.pop() self._parent_ = state.pop() self.name = state.pop() - + def copy(self, *args): constr = self.constraints.copy() priors = self.priors.copy() @@ -151,13 +151,13 @@ class Param(OptimizationHandlable, ObservableArray): # if trigger_parent: min_priority = None # else: min_priority = -numpy.inf # self.notify_observers(None, min_priority) -# +# # def _get_params(self): # return self.flat -# +# # def _collect_gradient(self, target): # target += self.gradient.flat -# +# # def _set_gradient(self, g): # self.gradient = g.reshape(self._realshape_) @@ -173,10 +173,10 @@ class Param(OptimizationHandlable, ObservableArray): try: new_arr._current_slice_ = s; new_arr._original_ = self.base is new_arr.base except AttributeError: pass # returning 0d array or float, double etc return new_arr - + def __setitem__(self, s, val): super(Param, self).__setitem__(s, val) - + #=========================================================================== # Index Operations: #=========================================================================== @@ -195,7 +195,7 @@ class Param(OptimizationHandlable, ObservableArray): a = self._realshape_[i] + a internal_offset += a * extended_realshape[i] return internal_offset - + def _raveled_index(self, slice_index=None): # return an index array on the raveled array, which is formed by the current_slice # of this object @@ -203,7 +203,7 @@ class Param(OptimizationHandlable, ObservableArray): ind = self._indices(slice_index) if ind.ndim < 2: ind = ind[:, None] return numpy.asarray(numpy.apply_along_axis(lambda x: numpy.sum(extended_realshape * x), 1, ind), dtype=int) - + def _expand_index(self, slice_index=None): # this calculates the full indexing arrays from the slicing objects given by get_item for _real..._ attributes # it basically translates slices to their respective index arrays and turns negative indices around @@ -245,7 +245,7 @@ class Param(OptimizationHandlable, ObservableArray): #=========================================================================== @property def _description_str(self): - if self.size <= 1: + if self.size <= 1: return [str(self.view(numpy.ndarray)[0])] else: return [str(self.shape)] def parameter_names(self, add_self=False, adjust_for_printing=False): @@ -356,7 +356,7 @@ class ParamConcatenation(object): self._param_sizes = [p.size for p in self.params] startstops = numpy.cumsum([0] + self._param_sizes) self._param_slices_ = [slice(start, stop) for start,stop in zip(startstops, startstops[1:])] - + parents = dict() for p in self.params: if p.has_parent(): @@ -396,7 +396,7 @@ class ParamConcatenation(object): def update_all_params(self): for par in self.parents: par.notify_observers(-numpy.inf) - + def constrain(self, constraint, warning=True): [param.constrain(constraint, trigger_parent=False) for param in self.params] self.update_all_params() @@ -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): - return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance, _debug=_debug) + def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): + return self.params[0]._highest_parent_._checkgrad(self, verbose, step, tolerance) #checkgrad.__doc__ = Gradcheckable.checkgrad.__doc__ __lt__ = lambda self, val: self._vals() < val diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 2a61c970..51b6cddf 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -1,7 +1,7 @@ # Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Licensed under the BSD 3-clause license (see LICENSE.txt) """ -Core module for parameterization. +Core module for parameterization. This module implements all parameterization techniques, split up in modular bits. HierarchyError: @@ -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): @@ -61,7 +61,7 @@ class Observable(object): def notify_observers(self, which=None, min_priority=None): """ - Notifies all observers. Which is the element, which kicked off this + Notifies all observers. Which is the element, which kicked off this notification loop. NOTE: notifies only observers with priority p > min_priority! @@ -91,11 +91,11 @@ class Observable(object): class Pickleable(object): """ - Make an object pickleable (See python doc 'pickling'). + Make an object pickleable (See python doc 'pickling'). This class allows for pickling support by Memento pattern. _getstate returns a memento of the class, which gets pickled. - _setstate() (re-)sets the state of the class to the memento + _setstate() (re-)sets the state of the class to the memento """ #=========================================================================== # Pickling operations @@ -112,14 +112,14 @@ class Pickleable(object): with open(f, 'w') as f: cPickle.dump(self, f, protocol) else: - cPickle.dump(self, f, protocol) + cPickle.dump(self, f, protocol) def __getstate__(self): if self._has_get_set_state(): return self._getstate() return self.__dict__ def __setstate__(self, state): if self._has_get_set_state(): - self._setstate(state) + self._setstate(state) # TODO: maybe parameters_changed() here? return self.__dict__ = state @@ -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): """ @@ -201,18 +201,18 @@ class Gradcheckable(Parentable): Adds the functionality for an object to be gradcheckable. It is just a thin wrapper of a call to the highest parent for now. TODO: Can be done better, by only changing parameters of the current parameter handle, - such that object hierarchy only has to change for those. + such that object hierarchy only has to change for those. """ def __init__(self, *a, **kw): super(Gradcheckable, self).__init__(*a, **kw) - def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3, _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 + Check the gradient of this parameter with respect to the highest parent's objective function. This is a three point estimate of the gradient, wiggling at the parameters with a stepsize step. - The check passes if either the ratio or the difference between numerical and + The check passes if either the ratio or the difference between numerical and analytical gradient is smaller then tolerance. :param bool verbose: whether each parameter shall be checked individually. @@ -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._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) - 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,22 +275,22 @@ 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): """ Flattened array of ints, specifying the index of this object. This has to account for shaped parameters! """ raise NotImplementedError, "Need to be able to get the raveled Index" - + def _internal_offset(self): """ - The offset for this parameter inside its parent. + The offset for this parameter inside its parent. This has to account for shaped parameters! """ return 0 - + def _offset_for(self, param): """ Return the offset of the param inside this parameterized object. @@ -298,24 +298,24 @@ class Indexable(object): basically just sums up the parameter sizes which come before param. """ raise NotImplementedError, "shouldnt happen, offset required from non parameterization object?" - + def _raveled_index_for(self, param): """ get the raveled index for a param that is an int array, containing the indexes for the flattened param inside this parameterized logic. """ - raise NotImplementedError, "shouldnt happen, raveld index transformation required from non parameterization 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!! Adding a constraint to a Parameter means to tell the highest parent that the constraint was added and making sure that all parameters covered by this object are indeed conforming to the constraint. - + :func:`constrain()` and :func:`unconstrain()` are main methods here """ def __init__(self, name, default_constraint=None, *a, **kw): @@ -326,7 +326,7 @@ class Constrainable(Nameable, Indexable): self.priors = ParameterIndexOperations() if self._default_constraint_ is not None: self.constrain(self._default_constraint_) - + def _disconnect_parent(self, constr=None, *args, **kw): """ From Parentable: @@ -340,7 +340,7 @@ class Constrainable(Nameable, Indexable): self._parent_index_ = None self._connect_fixes() self._notify_parent_change() - + #=========================================================================== # Fixing Parameters: #=========================================================================== @@ -352,24 +352,26 @@ 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): """ This parameter will no longer be fixed. """ unconstrained = self.unconstrain(__fixed__) - self._highest_parent_._set_unfixed(unconstrained) + self._highest_parent_._set_unfixed(unconstrained) unfix = unconstrain_fixed - + def _set_fixed(self, index): if not self._has_fixes(): self._fixes_ = np.ones(self.size, dtype=bool) self._fixes_[index] = FIXED if np.all(self._fixes_): self._fixes_ = None # ==UNFIXED - + def _set_unfixed(self, index): if not self._has_fixes(): self._fixes_ = np.ones(self.size, dtype=bool) # rav_i = self._raveled_index_for(param)[index] @@ -383,7 +385,7 @@ class Constrainable(Nameable, Indexable): self._fixes_[fixed_indices] = FIXED else: self._fixes_ = None - + def _has_fixes(self): return hasattr(self, "_fixes_") and self._fixes_ is not None @@ -398,21 +400,21 @@ class Constrainable(Nameable, Indexable): """ repriorized = self.unset_priors() self._add_to_index_operations(self.priors, repriorized, prior, warning) - + def unset_priors(self, *priors): """ Un-set all priors given from this parameter handle. - + """ return self._remove_from_index_operations(self.priors, priors) - + def log_prior(self): """evaluate the prior""" if self.priors.size > 0: x = self._get_params() return reduce(lambda a, b: a + b, [p.lnpdf(x[ind]).sum() for p, ind in self.priors.iteritems()], 0) return 0. - + def _log_prior_gradients(self): """evaluate the gradients of the priors""" if self.priors.size > 0: @@ -421,7 +423,7 @@ class Constrainable(Nameable, Indexable): [np.put(ret, ind, p.lnpdf_grad(x[ind])) for p, ind in self.priors.iteritems()] return ret return 0. - + #=========================================================================== # Constrain operations -> done #=========================================================================== @@ -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): """ @@ -448,7 +450,7 @@ class Constrainable(Nameable, Indexable): transformats of this parameter object. """ return self._remove_from_index_operations(self.constraints, transforms) - + def constrain_positive(self, warning=True, trigger_parent=True): """ :param warning: print a warning if re-constraining parameters. @@ -493,7 +495,7 @@ class Constrainable(Nameable, Indexable): Remove (lower, upper) bounded constrain from this parameter/ """ self.unconstrain(Logistic(lower, upper)) - + def _parent_changed(self, parent): """ From Parentable: @@ -522,7 +524,7 @@ class Constrainable(Nameable, Indexable): def _remove_from_index_operations(self, which, what): """ Helper preventing copy code. - Remove given what (transform prior etc) from which param index ops. + Remove given what (transform prior etc) from which param index ops. """ if len(what) == 0: transforms = which.properties() @@ -532,10 +534,10 @@ class Constrainable(Nameable, Indexable): removed = np.union1d(removed, unconstrained) if t is __fixed__: self._highest_parent_._set_unfixed(unconstrained) - + return removed -class OptimizationHandlable(Constrainable, Observable): +class OptimizationHandlable(Constrainable): """ This enables optimization handles on an Object as done in GPy 0.4. @@ -543,13 +545,13 @@ class OptimizationHandlable(Constrainable, Observable): """ def __init__(self, name, default_constraint=None, *a, **kw): super(OptimizationHandlable, self).__init__(name, default_constraint=default_constraint, *a, **kw) - + def transform(self): [np.put(self._param_array_, ind, c.finv(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - + def untransform(self): [np.put(self._param_array_, ind, c.f(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - + def _get_params_transformed(self): # transformed parameters (apply transformation rules) p = self._param_array_.copy() @@ -565,23 +567,21 @@ class OptimizationHandlable(Constrainable, Observable): else: self._param_array_[:] = p self.untransform() self._trigger_params_changed() - + def _trigger_params_changed(self, trigger_parent=True): [p._trigger_params_changed(trigger_parent=False) for p in self._parameters_] - if trigger_parent: min_priority = None - else: min_priority = -np.inf - self.notify_observers(None, min_priority) - + self.notify_observers(None, None if trigger_parent else -np.inf) + def _size_transformed(self): return self.size - self.constraints[__fixed__].size -# +# # def _untransform_params(self, p): # # inverse apply transformations for parameters # #p = p.copy() # if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp # [np.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] # return p -# +# # def _get_params(self): # """ # get all parameters @@ -592,7 +592,7 @@ class OptimizationHandlable(Constrainable, Observable): # return p # [np.put(p, ind, par._get_params()) for ind, par in itertools.izip(self._param)] # return p - + # def _set_params(self, params, trigger_parent=True): # self._param_array_.flat = params # if trigger_parent: min_priority = None @@ -600,14 +600,14 @@ class OptimizationHandlable(Constrainable, Observable): # self.notify_observers(None, min_priority) # don't overwrite this anymore! #raise NotImplementedError, "Abstract superclass: This needs to be implemented in Param and Parameterizable" - + #=========================================================================== # Optimization handles: #=========================================================================== def _get_param_names(self): n = np.array([p.hierarchy_name() + '[' + str(i) + ']' for p in self.flattened_parameters for i in p._indices()]) return n - + def _get_param_names_transformed(self): n = self._get_param_names() if self._has_fixes(): @@ -621,7 +621,7 @@ class OptimizationHandlable(Constrainable, Observable): """ Randomize the model. Make this draw from the prior if one exists, else draw from given random generator - + :param rand_gen: numpy random number generator which takes args and kwargs :param flaot loc: loc parameter for random number generator :param float scale: scale parameter for random number generator @@ -663,7 +663,7 @@ class Parameterizable(OptimizationHandlable): def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): """ - Get the names of all parameters of this model. + Get the names of all parameters of this model. :param bool add_self: whether to add the own name in front of names :param bool adjust_for_printing: whether to call `adjust_name_for_printing` on names @@ -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" @@ -712,7 +716,7 @@ class Parameterizable(OptimizationHandlable): #========================================================================= @property def gradient(self): - return self._gradient_array_ + return self._gradient_array_ @gradient.setter def gradient(self, val): @@ -821,8 +825,8 @@ class Parameterizable(OptimizationHandlable): # connect parameterlist to this parameterized object # This just sets up the right connection for the params objects # to be used as parameters - # it also sets the constraints for each parameter to the constraints - # of their respective parents + # it also sets the constraints for each parameter to the constraints + # of their respective parents if not hasattr(self, "_parameters_") or len(self._parameters_) < 1: # no parameters for this class return @@ -837,7 +841,7 @@ class Parameterizable(OptimizationHandlable): pslice = slice(old_size, old_size+p.size) # first connect all children - p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice]) + p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice]) # then connect children to self self._param_array_[pslice] = p._param_array_.ravel()#, requirements=['C', 'W']).ravel(order='C') self._gradient_array_[pslice] = p._gradient_array_.ravel()#, requirements=['C', 'W']).ravel(order='C') @@ -879,7 +883,7 @@ class Parameterizable(OptimizationHandlable): dc[k] = copy.deepcopy(v) if k == '_parameters_': params = [p.copy() for p in v] - + dc['_parent_'] = None dc['_parent_index_'] = None dc['_observer_callables_'] = [] @@ -890,12 +894,12 @@ class Parameterizable(OptimizationHandlable): s = self.__new__(self.__class__) s.__dict__ = dc - + for p in params: s.add_parameter(p, _ignore_added_names=True) - + return s - + #=========================================================================== # From being parentable, we have to define the parent_change notification #=========================================================================== diff --git a/GPy/examples/coreg_example.py b/GPy/examples/coreg_example.py new file mode 100644 index 00000000..967758c6 --- /dev/null +++ b/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) + diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py index fdebdfac..1e386c01 100644 --- a/GPy/kern/_src/add.py +++ b/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 diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index 0aa414ca..dc6eceb4 100644 --- a/GPy/kern/_src/kern.py +++ b/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 = [] - if isinstance(self, Prod): kernels.extend(self._parameters_) - else: kernels.append(self) - if isinstance(other, Prod): kernels.extend(other._parameters_) - else: kernels.append(other) - return Prod(self, other, name) + #kernels = [] + #if isinstance(self, Prod): kernels.extend(self._parameters_) + #else: kernels.append(self) + #if isinstance(other, Prod): kernels.extend(other._parameters_) + #else: kernels.append(other) + 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) diff --git a/GPy/kern/_src/kernel_slice_operations.py b/GPy/kern/_src/kernel_slice_operations.py index c1774a35..ff33cc24 100644 --- a/GPy/kern/_src/kernel_slice_operations.py +++ b/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.update_gradients_full = _slice_wrapper(instance, instance.update_gradients_full, False, True) - instance.update_gradients_diag = _slice_wrapper(instance, instance.update_gradients_diag, True, True) - instance.gradients_X = _slice_wrapper(instance, instance.gradients_X, False, True) - instance.gradients_X_diag = _slice_wrapper(instance, instance.gradients_X_diag, True, True) - instance.psi0 = _slice_wrapper(instance, instance.psi0, False, False) - instance.psi1 = _slice_wrapper(instance, instance.psi1, False, False) - instance.psi2 = _slice_wrapper(instance, instance.psi2, False, False) - instance.update_gradients_expectations = _slice_wrapper(instance, instance.update_gradients_expectations, psi_stat=True) - instance.gradients_Z_expectations = _slice_wrapper(instance, instance.gradients_Z_expectations, psi_stat_Z=True) - instance.gradients_qX_expectations = _slice_wrapper(instance, instance.gradients_qX_expectations, psi_stat=True) + instance.Kdiag = _slice_wrapper(instance, instance.Kdiag, diag=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, diag=True, derivative=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, diag=True, derivative=True) + instance.psi0 = _slice_wrapper(instance, instance.psi0, diag=False, derivative=False) + instance.psi1 = _slice_wrapper(instance, instance.psi1, diag=False, derivative=False) + instance.psi2 = _slice_wrapper(instance, instance.psi2, diag=False, derivative=False) + 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, derivative=True, psi_stat_Z=True) + instance.gradients_qX_expectations = _slice_wrapper(instance, instance.gradients_qX_expectations, derivative=True, psi_stat=True) instance.parameters_changed() return instance @@ -44,7 +44,29 @@ def _slice_wrapper(kern, operation, diag=False, derivative=False, psi_stat=False finally: kern._sliced_X -= 1 return ret - else: + 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 kern._sliced_X += 1 @@ -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): diff --git a/GPy/kern/_src/mlp.py b/GPy/kern/_src/mlp.py index 85792acd..ee15d967 100644 --- a/GPy/kern/_src/mlp.py +++ b/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] diff --git a/GPy/kern/_src/periodic.py b/GPy/kern/_src/periodic.py index 36ff3527..6b423a57 100644 --- a/GPy/kern/_src/periodic.py +++ b/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_phi = np.array(sum([[-np.pi/2, 0.] 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.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_phi = np.array(sum([[-np.pi/2, 0.] 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.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 - def update_gradients_full(self,dL_dK,X,X2,target): + #@silence_errors + 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_phi = np.array(sum([[-np.pi/2, 0.] 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.zeros(self.n_freq * 2) + self.basis_phi[::2] = -np.pi/2 self.G = self.Gram_matrix() self.Gi = np.linalg.inv(self.G) diff --git a/GPy/kern/_src/prod.py b/GPy/kern/_src/prod.py index bd664681..f3b2b50f 100644 --- a/GPy/kern/_src/prod.py +++ b/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, X)) - self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2]) + k1.update_gradients_full(dL_dK*k2.K(X, X2), X, X2) + k2.update_gradients_full(dL_dK*k1.K(X, X2), X, X2) + + 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: - target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2.K(X[:,self.slice2]), X[:,self.slice1], None) - target[:,self.slice2] += self.k2.gradients_X(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2], None) - else: - target[:,self.slice1] += self.k1.gradients_X(dL_dK*self.k2.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]) + for k1,k2 in itertools.combinations(self.parts, 2): + target[:,k1.active_dims] += k1.gradients_X(dL_dK*k2.K(X, X2), X, X2) + target[:,k2.active_dims] += k2.gradients_X(dL_dK*k1.K(X, X2), X, X2) 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]) - target[:,self.slice2] += self.k2.gradients_X(dL_dKdiag*self.k1.Kdiag(X[:,self.slice1]), X[:,self.slice2]) + for k1,k2 in itertools.combinations(self.parts, 2): + 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 diff --git a/GPy/kern/_src/sympykern.py b/GPy/kern/_src/sympykern.py index 9878ec68..3f6b5445 100644 --- a/GPy/kern/_src/sympykern.py +++ b/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() diff --git a/GPy/models/mrd.py b/GPy/models/mrd.py index b547f2d1..17949012 100644 --- a/GPy/models/mrd.py +++ b/GPy/models/mrd.py @@ -15,13 +15,13 @@ from ..likelihoods import Gaussian class MRD(Model): """ - Apply MRD to all given datasets Y in Ylist. - + Apply MRD to all given datasets Y in Ylist. + Y_i in [n x p_i] - - The samples n in the datasets need + + The samples n in the datasets need to match up, whereas the dimensionality p_d can differ. - + :param [array-like] Ylist: List of datasets to apply MRD on :param input_dim: latent dimensionality :type input_dim: int @@ -45,13 +45,12 @@ 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, + def __init__(self, Ylist, input_dim, X=None, X_variance=None, initx = 'PCA', initz = 'permute', - num_inducing=10, Z=None, kernel=None, + num_inducing=10, Z=None, kernel=None, inference_method=None, likelihood=None, name='mrd', Ynames=None): super(MRD, self).__init__(name) - + # sort out the kernels if kernel is None: from ..kern import RBF @@ -64,23 +63,23 @@ class MRD(Model): self.kern = kernel self.input_dim = input_dim self.num_inducing = num_inducing - + self.Ylist = Ylist self._in_init_ = True X = self._init_X(initx, Ylist) self.Z = Param('inducing inputs', self._init_Z(initz, X)) self.num_inducing = self.Z.shape[0] # ensure M==N if M>N - + if X_variance is None: X_variance = np.random.uniform(0, .2, X.shape) - + self.variational_prior = NormalPrior() self.X = NormalPosterior(X, X_variance) - + if likelihood is None: self.likelihood = [Gaussian(name='Gaussian_noise'.format(i)) for i in range(len(Ylist))] else: self.likelihood = likelihood - + if inference_method is None: self.inference_method= [] for y in Ylist: @@ -91,12 +90,12 @@ class MRD(Model): else: self.inference_method = inference_method self.inference_method.set_limit(len(Ylist)) - + self.add_parameters(self.X, self.Z) - + if Ynames is None: Ynames = ['Y{}'.format(i) for i in range(len(Ylist))] - + for i, n, k, l in itertools.izip(itertools.count(), Ynames, self.kern, self.likelihood): p = Parameterized(name=n) p.add_parameter(k) @@ -104,23 +103,23 @@ class MRD(Model): setattr(self, 'Y{}'.format(i), p) self.add_parameter(p) self._in_init_ = False - + def parameters_changed(self): self._log_marginal_likelihood = 0 self.posteriors = [] self.Z.gradient = 0. self.X.mean.gradient = 0. self.X.variance.gradient = 0. - + for y, k, l, i in itertools.izip(self.Ylist, self.kern, self.likelihood, self.inference_method): posterior, lml, grad_dict = i.inference(k, self.X, self.Z, l, y) - + self.posteriors.append(posterior) self._log_marginal_likelihood += lml - + # likelihood gradients l.update_gradients(grad_dict.pop('partial_for_likelihood')) - + #gradients wrt kernel dL_dKmm = grad_dict.pop('dL_dKmm') k.update_gradients_full(dL_dKmm, self.Z, None) @@ -132,7 +131,7 @@ class MRD(Model): self.Z.gradient += k.gradients_X(dL_dKmm, self.Z) self.Z.gradient += k.gradients_Z_expectations( grad_dict['dL_dpsi1'], grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) - + dL_dmean, dL_dS = k.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, **grad_dict) self.X.mean.gradient += dL_dmean self.X.variance.gradient += dL_dS diff --git a/GPy/testing/bgplvm_tests.py b/GPy/testing/bgplvm_tests.py deleted file mode 100644 index fd55d314..00000000 --- a/GPy/testing/bgplvm_tests.py +++ /dev/null @@ -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() diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 2789d1de..d54b3871 100644 --- a/GPy/testing/kernel_tests.py +++ b/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.X2 = np.random.randn(110,2) + self.N, self.D = 100, 5 + 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() diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py index b1db94a7..461a00c4 100644 --- a/GPy/testing/likelihood_tests.py +++ b/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) diff --git a/GPy/testing/mrd_tests.py b/GPy/testing/mrd_tests.py deleted file mode 100644 index 40fcb86a..00000000 --- a/GPy/testing/mrd_tests.py +++ /dev/null @@ -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() diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index b2f57144..6555b8f4 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -16,21 +16,21 @@ class Test(unittest.TestCase): from GPy.core.parameterization import Param from GPy.core.parameterization.transformations import Logistic self.param = Param('param', np.random.rand(25,2), Logistic(0, 1)) - + self.test1 = GPy.core.Parameterized("test model") self.test1.add_parameter(self.white) self.test1.add_parameter(self.rbf, 0) self.test1.add_parameter(self.param) - + x = np.linspace(-2,6,4)[:,None] y = np.sin(x) self.testmodel = GPy.models.GPRegression(x,y) - + def test_add_parameter(self): self.assertEquals(self.rbf._parent_index_, 0) self.assertEquals(self.white._parent_index_, 1) pass - + def test_fixes(self): self.white.fix(warning=False) self.test1.remove_parameter(self.test1.param) @@ -41,18 +41,18 @@ class Test(unittest.TestCase): self.test1.add_parameter(self.white, 0) self.assertListEqual(self.test1._fixes_.tolist(),[FIXED,UNFIXED,UNFIXED]) - + def test_remove_parameter(self): from GPy.core.parameterization.transformations import FIXED, UNFIXED, __fixed__, Logexp self.white.fix() self.test1.remove_parameter(self.white) self.assertIs(self.test1._fixes_,None) - + self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) self.assertEquals(self.white.constraints._offset, 0) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) - + self.test1.add_parameter(self.white, 0) self.assertIs(self.test1.constraints, self.white.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) @@ -60,17 +60,17 @@ class Test(unittest.TestCase): self.assertListEqual(self.test1.constraints[__fixed__].tolist(), [0]) self.assertIs(self.white._fixes_,None) self.assertListEqual(self.test1._fixes_.tolist(),[FIXED] + [UNFIXED] * 52) - + self.test1.remove_parameter(self.white) self.assertIs(self.test1._fixes_,None) self.assertListEqual(self.white._fixes_.tolist(), [FIXED]) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.param.constraints._param_index_ops) self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1]) - + def test_add_parameter_already_in_hirarchy(self): self.assertRaises(HierarchyError, self.test1.add_parameter, self.white._parameters_[0]) - + def test_default_constraints(self): self.assertIs(self.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops) self.assertIs(self.test1.constraints, self.rbf.constraints._param_index_ops) @@ -83,7 +83,7 @@ class Test(unittest.TestCase): self.rbf.constrain(GPy.transformations.Square(), False) self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), range(2)) self.assertListEqual(self.test1.constraints[GPy.transformations.Logexp()].tolist(), [2]) - + self.test1.remove_parameter(self.rbf) self.assertListEqual(self.test1.constraints[GPy.transformations.Square()].tolist(), []) diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py index 0cb4cd66..a7ebe6fe 100644 --- a/GPy/testing/unit_tests.py +++ b/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): - ''' Testing the GP regression with polynomial kernel with white kernel on 1d data ''' - mlp = GPy.kern.Poly(1, degree=5) - self.check_model(mlp, model_type='GPRegression', dimension=1) + #TODO: + #def test_GPRegression_poly_1d(self): + # ''' Testing the GP regression with polynomial kernel with white kernel on 1d data ''' + # 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) - #distribution = GPy.likelihoods.likelihood_functions.Bernoulli() - #likelihood = GPy.likelihoods.EP(Y, distribution) - #m = GPy.core.SparseGP(X, likelihood, kernel, Z) - #m.ensure_default_constraints() + m = GPy.models.SparseGPClassification(X, Y, kernel=kernel, Z=Z) + # distribution = GPy.likelihoods.likelihood_functions.Bernoulli() + # likelihood = GPy.likelihoods.EP(Y, distribution) + # m = GPy.core.SparseGP(X, likelihood, kernel, Z) + # 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] - m = GPy.models.FITCClassification(X, Y, kernel = k) + Y = np.hstack([np.ones(N / 2), np.zeros(N / 2)])[:, None] + 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()) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 96b4cee9..792d82e2 100644 --- a/GPy/util/caching.py +++ b/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): - i = state.index(True) - if self.inputs_changed[i]: - #(elements of) the args have changed since we last computed: update - self.cached_outputs[i] = self.operation(*args, **kw) - self.inputs_changed[i] = False - return self.cached_outputs[i] - else: - #first time we've seen these arguments: compute + try: + if any(state): + i = state.index(True) + if self.inputs_changed[i]: + #(elements of) the args have changed since we last computed: update + self.cached_outputs[i] = self.operation(*args, **kw) + self.inputs_changed[i] = False + return self.cached_outputs[i] + else: + #first time we've seen these arguments: compute - #first make sure the depth limit isn't exceeded - if len(self.cached_inputs) == self.limit: - args_ = self.cached_inputs.pop(0) - [a.remove_observer(self, self.on_cache_changed) for a in args_ if a is not None] - self.inputs_changed.pop(0) - self.cached_outputs.pop(0) - - #compute - self.cached_inputs.append(oa_all) - self.cached_outputs.append(self.operation(*args, **kw)) - self.inputs_changed.append(False) - [a.add_observer(self, self.on_cache_changed) for a in observable_args] - return self.cached_outputs[-1]#return + #first make sure the depth limit isn't exceeded + if len(self.cached_inputs) == self.limit: + args_ = self.cached_inputs.pop(0) + [a.remove_observer(self, self.on_cache_changed) for a in args_ if a is not None] + self.inputs_changed.pop(0) + self.cached_outputs.pop(0) + #compute + self.cached_inputs.append(oa_all) + self.cached_outputs.append(self.operation(*args, **kw)) + self.inputs_changed.append(False) + [a.add_observer(self, self.on_cache_changed) for a in observable_args] + 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 """ diff --git a/GPy/util/data_resources.json b/GPy/util/data_resources.json index ca15bf2d..57b79f10 100644 --- a/GPy/util/data_resources.json +++ b/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":[ [ diff --git a/GPy/util/datasets.py b/GPy/util/datasets.py index 23f5d0c8..54e42733 100644 --- a/GPy/util/datasets.py +++ b/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.""" @@ -274,8 +286,76 @@ def della_gatta_TRP63_gene_expression(data_set='della_gatta', gene_number=None): Y = Y[:, None] return data_details_return({'X': X, 'Y': Y, 'gene_number' : gene_number}, data_set) + +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).""" diff --git a/GPy/util/football_teams.json b/GPy/util/football_teams.json new file mode 100644 index 00000000..a4eb9c38 --- /dev/null +++ b/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} \ No newline at end of file diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py index bec0e490..79022a5f 100644 --- a/GPy/util/multioutput.py +++ b/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