diff --git a/GPy/core/model.py b/GPy/core/model.py index 6514d73a..d27cbc69 100644 --- a/GPy/core/model.py +++ b/GPy/core/model.py @@ -60,20 +60,6 @@ class Model(Parameterized): self.priors = state.pop() Parameterized._setstate(self, state) - def randomize(self): - """ - Randomize the model. - Make this draw from the prior if one exists, else draw from N(0,1) - """ - # first take care of all parameters (from N(0,1)) - # x = self._get_params_transformed() - x = np.random.randn(self.size_transformed) - x = self._untransform_params(x) - # now draw from prior where possible - [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None] - self._set_params(x) - # self._set_params_transformed(self._get_params_transformed()) # makes sure all of the tied parameters get the same init (since there's only one prior object...) - def optimize_restarts(self, num_restarts=10, robust=False, verbose=True, parallel=False, num_processes=None, **kwargs): """ Perform random restarts of the model, and set the model to the best @@ -240,6 +226,11 @@ class Model(Parameterized): TODO: valid args """ + if self.is_fixed: + raise RuntimeError, "Cannot optimize, when everything is fixed" + if self.size == 0: + raise RuntimeError, "Model without parameters cannot be minimized" + if optimizer is None: optimizer = self.preferred_optimizer @@ -279,7 +270,7 @@ class Model(Parameterized): and numerical gradients is within of unity. """ - x = self._get_params_transformed().copy() + x = self._get_params_transformed() if not verbose: # make sure only to test the selected parameters @@ -297,7 +288,7 @@ class Model(Parameterized): return # just check the global ratio - dx = np.zeros_like(x) + dx = np.zeros(x.shape) dx[transformed_index] = step * np.sign(np.random.uniform(-1, 1, transformed_index.size)) # evaulate around the point x @@ -308,9 +299,8 @@ class Model(Parameterized): dx = dx[transformed_index] gradient = gradient[transformed_index] - numerical_gradient = (f1 - f2) / (2 * dx) global_ratio = (f1 - f2) / (2 * np.dot(dx, np.where(gradient == 0, 1e-32, gradient))) - return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gradient - numerical_gradient).mean() < tolerance) + return (np.abs(1. - global_ratio) < tolerance) else: # check the gradient of each parameter individually, and do some pretty printing try: diff --git a/GPy/core/parameterization/array_core.py b/GPy/core/parameterization/array_core.py index a338ceed..cf353ead 100644 --- a/GPy/core/parameterization/array_core.py +++ b/GPy/core/parameterization/array_core.py @@ -6,19 +6,6 @@ __updated__ = '2013-12-16' import numpy as np from parameter_core import Observable -class ParamList(list): - """ - List to store ndarray-likes in. - It will look for 'is' instead of calling __eq__ on each element. - """ - def __contains__(self, other): - for el in self: - if el is other: - return True - return False - - pass - class ObservableArray(np.ndarray, Observable): """ An ndarray which reports changes to its observers. @@ -62,10 +49,11 @@ class ObservableArray(np.ndarray, Observable): def __setitem__(self, s, val): if self._s_not_empty(s): super(ObservableArray, self).__setitem__(s, val) - self._notify_observers() + self._notify_observers(self[s]) def __getslice__(self, start, stop): return self.__getitem__(slice(start, stop)) + def __setslice__(self, start, stop, val): return self.__setitem__(slice(start, stop), val) diff --git a/GPy/core/parameterization/index_operations.py b/GPy/core/parameterization/index_operations.py index b5399741..a9f3768e 100644 --- a/GPy/core/parameterization/index_operations.py +++ b/GPy/core/parameterization/index_operations.py @@ -5,47 +5,7 @@ Created on Oct 2, 2013 ''' import numpy from numpy.lib.function_base import vectorize -from param import Param -from collections import defaultdict - -class ParamDict(defaultdict): - def __init__(self): - """ - Default will be self._default, if not set otherwise - """ - defaultdict.__init__(self, self.default_factory) - - def __getitem__(self, key): - try: - return defaultdict.__getitem__(self, key) - except KeyError: - for a in self.iterkeys(): - if numpy.all(a==key) and a._parent_index_==key._parent_index_: - return defaultdict.__getitem__(self, a) - raise - - def __contains__(self, key): - if defaultdict.__contains__(self, key): - return True - for a in self.iterkeys(): - if numpy.all(a==key) and a._parent_index_==key._parent_index_: - return True - return False - - def __setitem__(self, key, value): - if isinstance(key, Param): - for a in self.iterkeys(): - if numpy.all(a==key) and a._parent_index_==key._parent_index_: - return super(ParamDict, self).__setitem__(a, value) - defaultdict.__setitem__(self, key, value) - -class SetDict(ParamDict): - def default_factory(self): - return set() - -class IntArrayDict(ParamDict): - def default_factory(self): - return numpy.int_([]) +from lists_and_dicts import IntArrayDict class ParameterIndexOperations(object): ''' @@ -194,9 +154,13 @@ class ParameterIndexOperationsView(object): def shift_right(self, start, size): - raise NotImplementedError, 'Shifting only supported in original ParamIndexOperations' - + self._param_index_ops.shift_right(start+self._offset, size) + def shift_left(self, start, size): + self._param_index_ops.shift_left(start+self._offset, size) + self._offset -= size + self._size -= size + def clear(self): for i, ind in self.items(): self._param_index_ops.remove(i, ind+self._offset) @@ -232,9 +196,7 @@ class ParameterIndexOperationsView(object): def __getitem__(self, prop): ind = self._filter_index(self._param_index_ops[prop]) - if ind.size > 0: - return ind - raise KeyError, prop + return ind def __str__(self, *args, **kwargs): import pprint diff --git a/GPy/core/parameterization/lists_and_dicts.py b/GPy/core/parameterization/lists_and_dicts.py new file mode 100644 index 00000000..5b13b3b5 --- /dev/null +++ b/GPy/core/parameterization/lists_and_dicts.py @@ -0,0 +1,35 @@ +''' +Created on 27 Feb 2014 + +@author: maxz +''' + +from collections import defaultdict +class DefaultArrayDict(defaultdict): + def __init__(self): + """ + Default will be self._default, if not set otherwise + """ + defaultdict.__init__(self, self.default_factory) + +class SetDict(DefaultArrayDict): + def default_factory(self): + return set() + +class IntArrayDict(DefaultArrayDict): + def default_factory(self): + import numpy as np + return np.int_([]) + +class ArrayList(list): + """ + List to store ndarray-likes in. + It will look for 'is' instead of calling __eq__ on each element. + """ + def __contains__(self, other): + for el in self: + if el is other: + return True + return False + + pass diff --git a/GPy/core/parameterization/param.py b/GPy/core/parameterization/param.py index 89d3a4e4..22610a70 100644 --- a/GPy/core/parameterization/param.py +++ b/GPy/core/parameterization/param.py @@ -3,8 +3,8 @@ import itertools import numpy -from parameter_core import Constrainable, Gradcheckable, Indexable, Parentable, adjust_name_for_printing -from array_core import ObservableArray, ParamList +from parameter_core import OptimizationHandlable, Gradcheckable, adjust_name_for_printing +from array_core import ObservableArray ###### printing __constraints_name__ = "Constraint" @@ -15,7 +15,7 @@ __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision __print_threshold__ = 5 ###### -class Param(Constrainable, ObservableArray, Gradcheckable): +class Param(OptimizationHandlable, ObservableArray, Gradcheckable): """ Parameter object for GPy models. @@ -50,7 +50,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable): obj._realsize_ = obj.size obj._realndim_ = obj.ndim obj._updated_ = False - from index_operations import SetDict + from lists_and_dicts import SetDict obj._tied_to_me_ = SetDict() obj._tied_to_ = [] obj._original_ = True @@ -148,8 +148,11 @@ class Param(Constrainable, ObservableArray, Gradcheckable): #=========================================================================== # get/set parameters #=========================================================================== - def _set_params(self, param, update=True): + def _set_params(self, param, trigger_parent=True): self.flat = param + if trigger_parent: min_priority = None + else: min_priority = -numpy.inf + self._notify_observers(None, min_priority) def _get_params(self): return self.flat @@ -172,11 +175,9 @@ class Param(Constrainable, ObservableArray, Gradcheckable): 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) - if self.has_parent(): - self._direct_parent_._notify_parameters_changed() - #self._notify_observers() #=========================================================================== # Index Operations: @@ -204,6 +205,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable): 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 @@ -230,7 +232,8 @@ class Param(Constrainable, ObservableArray, Gradcheckable): #=========================================================================== @property def is_fixed(self): - return self._highest_parent_._is_fixed(self) + from transformations import __fixed__ + return self.constraints[__fixed__].size == self.size #def round(self, decimals=0, out=None): # view = super(Param, self).round(decimals, out).view(Param) # view.__array_finalize__(self) @@ -267,7 +270,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable): return [t._short() for t in self._tied_to_] or [''] def __repr__(self, *args, **kwargs): name = "\033[1m{x:s}\033[0;0m:\n".format( - x=self.hirarchy_name()) + x=self.hierarchy_name()) return name + super(Param, self).__repr__(*args, **kwargs) def _ties_for(self, rav_index): # size = sum(p.size for p in self._tied_to_) @@ -301,12 +304,12 @@ class Param(Constrainable, ObservableArray, Gradcheckable): gen = map(lambda x: " ".join(map(str, x)), gen) return reduce(lambda a, b:max(a, len(b)), gen, len(header)) def _max_len_values(self): - return reduce(lambda a, b:max(a, len("{x:=.{0}g}".format(__precision__, x=b))), self.flat, len(self.hirarchy_name())) + return reduce(lambda a, b:max(a, len("{x:=.{0}g}".format(__precision__, x=b))), self.flat, len(self.hierarchy_name())) def _max_len_index(self, ind): return reduce(lambda a, b:max(a, len(str(b))), ind, len(__index_name__)) def _short(self): # short string to print - name = self.hirarchy_name() + name = self.hierarchy_name() if self._realsize_ < 2: return name ind = self._indices() @@ -329,8 +332,8 @@ class Param(Constrainable, ObservableArray, Gradcheckable): if lp is None: lp = self._max_len_names(prirs, __tie_name__) sep = '-' header_format = " {i:{5}^{2}s} | \033[1m{x:{5}^{1}s}\033[0;0m | {c:{5}^{0}s} | {p:{5}^{4}s} | {t:{5}^{3}s}" - if only_name: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hirarchy_name(), c=sep*lc, i=sep*li, t=sep*lt, p=sep*lp) # nice header for printing - else: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hirarchy_name(), c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # nice header for printing + if only_name: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hierarchy_name(), c=sep*lc, i=sep*li, t=sep*lt, p=sep*lp) # nice header for printing + else: header = header_format.format(lc, lx, li, lt, lp, ' ', x=self.hierarchy_name(), c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # nice header for printing if not ties: ties = itertools.cycle(['']) return "\n".join([header] + [" {i!s:^{3}s} | {x: >{1}.{2}g} | {c:^{0}s} | {p:^{5}s} | {t:^{4}s} ".format(lc, lx, __precision__, li, lt, lp, x=x, c=" ".join(map(str, c)), p=" ".join(map(str, p)), t=(t or ''), i=i) for i, x, c, t, p in itertools.izip(indices, vals, constr_matrix, ties, prirs)]) # return all the constraints with right indices # except: return super(Param, self).__str__() @@ -345,7 +348,8 @@ class ParamConcatenation(object): See :py:class:`GPy.core.parameter.Param` for more details on constraining. """ # self.params = params - self.params = ParamList([]) + from lists_and_dicts import ArrayList + self.params = ArrayList([]) for p in params: for p in p.flattened_parameters: if p not in self.params: @@ -353,6 +357,21 @@ 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(): + parent = p._direct_parent_ + level = 0 + while parent is not None: + if parent in parents: + parents[parent] = max(level, parents[parent]) + else: + parents[parent] = level + level += 1 + parent = parent._direct_parent_ + import operator + self.parents = map(lambda x: x[0], sorted(parents.iteritems(), key=operator.itemgetter(1))) #=========================================================================== # Get/set items, enable broadcasting #=========================================================================== @@ -366,24 +385,26 @@ class ParamConcatenation(object): val = val._vals() ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; vals = self._vals(); vals[s] = val; del val - [numpy.place(p, ind[ps], vals[ps]) and update and p._notify_observers() + [numpy.place(p, ind[ps], vals[ps]) for p, ps in zip(self.params, self._param_slices_)] + if update: + self.update_all_params() def _vals(self): return numpy.hstack([p._get_params() for p in self.params]) #=========================================================================== # parameter operations: #=========================================================================== def update_all_params(self): - for p in self.params: - p._notify_observers() - + for par in self.parents: + par._notify_observers(-numpy.inf) + def constrain(self, constraint, warning=True): - [param.constrain(constraint, update=False) for param in self.params] + [param.constrain(constraint, trigger_parent=False) for param in self.params] self.update_all_params() constrain.__doc__ = Param.constrain.__doc__ def constrain_positive(self, warning=True): - [param.constrain_positive(warning, update=False) for param in self.params] + [param.constrain_positive(warning, trigger_parent=False) for param in self.params] self.update_all_params() constrain_positive.__doc__ = Param.constrain_positive.__doc__ @@ -393,12 +414,12 @@ class ParamConcatenation(object): fix = constrain_fixed def constrain_negative(self, warning=True): - [param.constrain_negative(warning, update=False) for param in self.params] + [param.constrain_negative(warning, trigger_parent=False) for param in self.params] self.update_all_params() constrain_negative.__doc__ = Param.constrain_negative.__doc__ def constrain_bounded(self, lower, upper, warning=True): - [param.constrain_bounded(lower, upper, warning, update=False) for param in self.params] + [param.constrain_bounded(lower, upper, warning, trigger_parent=False) for param in self.params] self.update_all_params() constrain_bounded.__doc__ = Param.constrain_bounded.__doc__ diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py index 6afa94cb..e7344fa5 100644 --- a/GPy/core/parameterization/parameter_core.py +++ b/GPy/core/parameterization/parameter_core.py @@ -2,34 +2,58 @@ # Licensed under the BSD 3-clause license (see LICENSE.txt) from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED +import heapq __updated__ = '2013-12-16' +class HierarchyError(Exception): + """ + Gets thrown when something is wrong with the parameter hierarchy + """ + def adjust_name_for_printing(name): if name is not None: return name.replace(" ", "_").replace(".", "_").replace("-", "").replace("+", "").replace("!", "").replace("*", "").replace("/", "") return '' class Observable(object): + _updated = True def __init__(self, *args, **kwargs): - from collections import defaultdict - self._observer_callables_ = defaultdict(list) - - def add_observer(self, observer, callble): - self._observer_callables_[observer].append(callble) + self._observer_callables_ = [] + def add_observer(self, observer, callble, priority=0): + heapq.heappush(self._observer_callables_, (priority, observer, callble)) + def remove_observer(self, observer, callble=None): - if observer in self._observer_callables_: - if callble is None: - del self._observer_callables_[observer] - elif callble in self._observer_callables_[observer]: - self._observer_callables_[observer].remove(callble) - if len(self._observer_callables_[observer]) == 0: - self.remove_observer(observer) - - def _notify_observers(self): - [[callble(self) for callble in callables] - for callables in self._observer_callables_.itervalues()] + to_remove = [] + for p, obs, clble in self._observer_callables_: + if callble is not None: + if (obs == observer) and (callble == clble): + to_remove.append((p, obs, clble)) + else: + if obs is observer: + to_remove.append((p, obs, clble)) + for r in to_remove: + self._observer_callables_.remove(r) + + def _notify_observers(self, which=None, min_priority=None): + """ + Notifies all observers. Which is the element, which kicked off this + notification loop. + + NOTE: notifies only observers with priority p > min_priority! + ^^^^^^^^^^^^^^^^ + + :param which: object, which started this notification loop + :param min_priority: only notify observers with priority > min_priority + if min_priority is None, notify all observers in order + """ + if which is None: + which = self + if min_priority is None: + [callble(which) for _, _, callble in heapq.nlargest(len(self._observer_callables_), self._observer_callables_)] + else: + [callble(which) for p, _, callble in heapq.nlargest(len(self._observer_callables_), self._observer_callables_) if p > min_priority] class Pickleable(object): def _getstate(self): @@ -95,11 +119,11 @@ class Nameable(Parentable): self._name = name if self.has_parent(): self._direct_parent_._name_changed(self, from_name) - def hirarchy_name(self, adjust_for_printing=True): + def hierarchy_name(self, adjust_for_printing=True): if adjust_for_printing: adjust = lambda x: adjust_name_for_printing(x) else: adjust = lambda x: x if self.has_parent(): - return self._direct_parent_.hirarchy_name() + "." + adjust(self.name) + return self._direct_parent_.hierarchy_name() + "." + adjust(self.name) return adjust(self.name) @@ -156,7 +180,7 @@ class Constrainable(Nameable, Indexable): #=========================================================================== # Fixing Parameters: #=========================================================================== - def constrain_fixed(self, value=None, warning=True): + def constrain_fixed(self, value=None, warning=True, trigger_parent=True): """ Constrain this paramter to be fixed to the current value it carries. @@ -164,7 +188,7 @@ class Constrainable(Nameable, Indexable): """ if value is not None: self[:] = value - self.constrain(__fixed__, warning=warning) + self.constrain(__fixed__, warning=warning, trigger_parent=trigger_parent) rav_i = self._highest_parent_._raveled_index_for(self) self._highest_parent_._set_fixed(rav_i) fix = constrain_fixed @@ -205,9 +229,9 @@ class Constrainable(Nameable, Indexable): #=========================================================================== # Prior Operations #=========================================================================== - def set_prior(self, prior, warning=True, update=True): + def set_prior(self, prior, warning=True, trigger_parent=True): repriorized = self.unset_priors() - self._add_to_index_operations(self.priors, repriorized, prior, warning, update) + self._add_to_index_operations(self.priors, repriorized, prior, warning) def unset_priors(self, *priors): return self._remove_from_index_operations(self.priors, priors) @@ -233,7 +257,7 @@ class Constrainable(Nameable, Indexable): # Constrain operations -> done #=========================================================================== - def constrain(self, transform, warning=True, update=True): + def constrain(self, transform, warning=True, trigger_parent=True): """ :param transform: the :py:class:`GPy.core.transformations.Transformation` to constrain the this parameter to. @@ -243,9 +267,9 @@ class Constrainable(Nameable, Indexable): :py:class:`GPy.core.transformations.Transformation`. """ if isinstance(transform, Transformation): - self._set_params(transform.initialize(self._get_params()), update=False) + self._set_params(transform.initialize(self._get_params()), trigger_parent=trigger_parent) reconstrained = self.unconstrain() - self._add_to_index_operations(self.constraints, reconstrained, transform, warning, update) + self._add_to_index_operations(self.constraints, reconstrained, transform, warning) def unconstrain(self, *transforms): """ @@ -256,30 +280,30 @@ class Constrainable(Nameable, Indexable): """ return self._remove_from_index_operations(self.constraints, transforms) - def constrain_positive(self, warning=True, update=True): + def constrain_positive(self, warning=True, trigger_parent=True): """ :param warning: print a warning if re-constraining parameters. Constrain this parameter to the default positive constraint. """ - self.constrain(Logexp(), warning=warning, update=update) + self.constrain(Logexp(), warning=warning, trigger_parent=trigger_parent) - def constrain_negative(self, warning=True, update=True): + def constrain_negative(self, warning=True, trigger_parent=True): """ :param warning: print a warning if re-constraining parameters. Constrain this parameter to the default negative constraint. """ - self.constrain(NegativeLogexp(), warning=warning, update=update) + self.constrain(NegativeLogexp(), warning=warning, trigger_parent=trigger_parent) - def constrain_bounded(self, lower, upper, warning=True, update=True): + def constrain_bounded(self, lower, upper, warning=True, trigger_parent=True): """ :param lower, upper: the limits to bound this parameter to :param warning: print a warning if re-constraining parameters. Constrain this parameter to lie within the given range. """ - self.constrain(Logistic(lower, upper), warning=warning, update=update) + self.constrain(Logistic(lower, upper), warning=warning, trigger_parent=trigger_parent) def unconstrain_positive(self): """ @@ -309,12 +333,11 @@ class Constrainable(Nameable, Indexable): for p in self._parameters_: p._parent_changed(parent) - def _add_to_index_operations(self, which, reconstrained, transform, warning, update): + def _add_to_index_operations(self, which, reconstrained, transform, warning): if warning and reconstrained.size > 0: + # TODO: figure out which parameters have changed and only print those print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name) which.add(transform, self._raveled_index()) - if update: - self._notify_observers() def _remove_from_index_operations(self, which, transforms): if len(transforms) == 0: @@ -329,12 +352,76 @@ class Constrainable(Nameable, Indexable): return removed +class OptimizationHandlable(Constrainable, Observable): + def _get_params_transformed(self): + # transformed parameters (apply transformation rules) + p = self._get_params() + [np.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] + if self._has_fixes(): + return p[self._fixes_] + return p + + def _set_params_transformed(self, p): + # inverse apply transformations for parameters and set the resulting parameters + self._set_params(self._untransform_params(p)) + + def _size_transformed(self): + return self.size - self.constraints[__fixed__].size + + def _untransform_params(self, p): + 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): + # don't overwrite this anymore! + if not self.size: + return np.empty(shape=(0,), dtype=np.float64) + return np.hstack([x._get_params() for x in self._parameters_ if x.size > 0]) -class Parameterizable(Constrainable, Observable): + def _set_params(self, params, trigger_parent=True): + # don't overwrite this anymore! + raise NotImplementedError, "This needs to be implemented in Param and Parametrizable" + + #=========================================================================== + # 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(): + return n[self._fixes_] + return n + + #=========================================================================== + # Randomizeable + #=========================================================================== + def randomize(self): + """ + Randomize the model. + Make this draw from the prior if one exists, else draw from N(0,1) + """ + import numpy as np + # first take care of all parameters (from N(0,1)) + # x = self._get_params_transformed() + x = np.random.randn(self._size_transformed()) + x = self._untransform_params(x) + # now draw from prior where possible + [np.put(x, ind, p.rvs(ind.size)) for p, ind in self.priors.iteritems() if not p is None] + self._set_params(x) + # self._set_params_transformed(self._get_params_transformed()) # makes sure all of the tied parameters get the same init (since there's only one prior object...) + + +import numpy as np + +class Parameterizable(OptimizationHandlable): def __init__(self, *args, **kwargs): super(Parameterizable, self).__init__(*args, **kwargs) - from GPy.core.parameterization.array_core import ParamList - _parameters_ = ParamList() + from GPy.core.parameterization.lists_and_dicts import ArrayList + _parameters_ = ArrayList() self._added_names_ = set() def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): @@ -357,7 +444,7 @@ class Parameterizable(Constrainable, Observable): if pname in self._added_names_: del self.__dict__[pname] self._add_parameter_name(param) - else: + elif pname not in dir(self): self.__dict__[pname] = param self._added_names_.add(pname) @@ -377,28 +464,26 @@ class Parameterizable(Constrainable, Observable): import itertools [p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + def _set_params(self, params, trigger_parent=True): + import itertools + [p._set_params(params[s], trigger_parent=False) for p, s in itertools.izip(self._parameters_, self._param_slices_)] + if trigger_parent: min_priority = None + else: min_priority = -np.inf + self._notify_observers(None, min_priority) + def _set_gradient(self, g): import itertools [p._set_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - def _get_params(self): - import numpy as np - # don't overwrite this anymore! - if not self.size: - return np.empty(shape=(0,), dtype=np.float64) - return np.hstack([x._get_params() for x in self._parameters_ if x.size > 0]) - - def _set_params(self, params, update=True): - # don't overwrite this anymore! - import itertools - [p._set_params(params[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] - self._notify_parameters_changed() - + + #=========================================================================== + # TODO: not working yet + #=========================================================================== def copy(self): """Returns a (deep) copy of the current model""" import copy from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView - from .array_core import ParamList + from .lists_and_dicts import ArrayList dc = dict() for k, v in self.__dict__.iteritems(): @@ -412,7 +497,7 @@ class Parameterizable(Constrainable, Observable): dc['_direct_parent_'] = None dc['_parent_index_'] = None - dc['_parameters_'] = ParamList() + dc['_parameters_'] = ArrayList() dc['constraints'].clear() dc['priors'].clear() dc['size'] = 0 @@ -424,12 +509,6 @@ class Parameterizable(Constrainable, Observable): s.add_parameter(p) return s - - def _notify_parameters_changed(self): - self.parameters_changed() - self._notify_observers() - if self.has_parent(): - self._direct_parent_._notify_parameters_changed() def parameters_changed(self): """ diff --git a/GPy/core/parameterization/parameterized.py b/GPy/core/parameterization/parameterized.py index f5fcc6ad..6fd60442 100644 --- a/GPy/core/parameterization/parameterized.py +++ b/GPy/core/parameterization/parameterized.py @@ -7,9 +7,9 @@ import cPickle import itertools from re import compile, _pattern_type from param import ParamConcatenation -from parameter_core import Constrainable, Pickleable, Parentable, Observable, Parameterizable, adjust_name_for_printing, Gradcheckable +from parameter_core import Pickleable, Parameterizable, adjust_name_for_printing, Gradcheckable from transformations import __fixed__ -from array_core import ParamList +from lists_and_dicts import ArrayList class Parameterized(Parameterizable, Pickleable, Gradcheckable): """ @@ -56,8 +56,9 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): def __init__(self, name=None, *a, **kw): super(Parameterized, self).__init__(name=name, parent=None, parent_index=None, *a, **kw) self._in_init_ = True - self._parameters_ = ParamList() + self._parameters_ = ArrayList() self.size = sum(p.size for p in self._parameters_) + self.add_observer(self, self._parameters_changed_notification, -100) if not self._has_fixes(): self._fixes_ = None self._param_slices_ = [] @@ -65,7 +66,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): del self._in_init_ def build_pydot(self, G=None): - import pydot + import pydot # @UnresolvedImport iamroot = False if G is None: G = pydot.Dot(graph_type='digraph') @@ -104,6 +105,14 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.remove_parameter(param) self.add_parameter(param, index) elif param not in self._parameters_: + if param.has_parent(): + parent = param._direct_parent_ + while parent is not None: + if parent is self: + from GPy.core.parameterization.parameter_core import HierarchyError + raise HierarchyError, "You cannot add a parameter twice into the hirarchy" + parent = parent._direct_parent_ + param._direct_parent_.remove_parameter(param) # make sure the size is set if index is None: self.constraints.update(param.constraints, self.size) @@ -116,12 +125,16 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self.constraints.update(param.constraints, start) self.priors.update(param.priors, start) self._parameters_.insert(index, param) + + param.add_observer(self, self._pass_through_notify_observers, -np.inf) + self.size += param.size + + self._connect_parameters() + self._notify_parent_change() + self._connect_fixes() else: raise RuntimeError, """Parameter exists already added and no copy made""" - self._connect_parameters() - self._notify_parent_change() - self._connect_fixes() def add_parameters(self, *parameters): @@ -144,12 +157,19 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): del self._parameters_[param._parent_index_] param._disconnect_parent() - param.remove_observer(self, self._notify_parameters_changed) + param.remove_observer(self, self._pass_through_notify_observers) self.constraints.shift_left(start, param.size) + self._connect_fixes() self._connect_parameters() self._notify_parent_change() + parent = self._direct_parent_ + while parent is not None: + parent._connect_fixes() + parent._connect_parameters() + parent._notify_parent_change() + parent = parent._direct_parent_ def _connect_parameters(self): # connect parameterlist to this parameterized object @@ -170,6 +190,13 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): self._add_parameter_name(p) #=========================================================================== + # notification system + #=========================================================================== + def _parameters_changed_notification(self, which): + self.parameters_changed() + def _pass_through_notify_observers(self, which): + self._notify_observers(which) + #=========================================================================== # Pickling operations #=========================================================================== def pickle(self, f, protocol=-1): @@ -237,42 +264,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): g[self._offset_for(p) + numpy.array(list(i))] += g[self._raveled_index_for(t)] if self._has_fixes(): return g[self._fixes_] return g - #=========================================================================== - # Optimization handles: - #=========================================================================== - def _get_param_names(self): - n = numpy.array([p.hirarchy_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(): - return n[self._fixes_] - return n - def _get_params_transformed(self): - # transformed parameters (apply transformation rules) - p = self._get_params() - [numpy.put(p, ind, c.finv(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - if self._has_fixes(): - return p[self._fixes_] - return p - def _set_params_transformed(self, p): - # inverse apply transformations for parameters and set the resulting parameters - self._set_params(self._untransform_params(p)) - def _untransform_params(self, p): - p = p.copy() - if self._has_fixes(): tmp = self._get_params(); tmp[self._fixes_] = p; p = tmp; del tmp - [numpy.put(p, ind, c.f(p[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__] - return p - #=========================================================================== - # Indexable Handling - #=========================================================================== - def _backtranslate_index(self, param, ind): - # translate an index in parameterized indexing into the index of param - ind = ind - self._offset_for(param) - ind = ind[ind >= 0] - internal_offset = param._internal_offset() - ind = ind[ind < param.size + internal_offset] - return ind + def _offset_for(self, param): # get the offset in the parameterized index array for param if param.has_parent(): @@ -297,34 +289,22 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): this is not in the global view of things! """ return numpy.r_[:self.size] - #=========================================================================== - # Fixing parameters: - #=========================================================================== - def _fixes_for(self, param): - if self._has_fixes(): - return self._fixes_[self._raveled_index_for(param)] - return numpy.ones(self.size, dtype=bool)[self._raveled_index_for(param)] + #=========================================================================== # Convenience for fixed, tied checking of param: #=========================================================================== - def fixed_indices(self): - return np.array([x.is_fixed for x in self._parameters_]) - def _is_fixed(self, param): - # returns if the whole param is fixed - if not self._has_fixes(): - return False - return not self._fixes_[self._raveled_index_for(param)].any() - # return not self._fixes_[self._offset_for(param): self._offset_for(param)+param._realsize_].any() @property def is_fixed(self): for p in self._parameters_: if not p.is_fixed: return False return True + def _get_original(self, param): # if advanced indexing is activated it happens that the array is a copy # you can retrieve the original param through this method, by passing # the copy here return self._parameters_[param._parent_index_] + #=========================================================================== # Get/set parameters: #=========================================================================== @@ -365,7 +345,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): # Printing: #=========================================================================== def _short(self): - return self.hirarchy_name() + return self.hierarchy_name() @property def flattened_parameters(self): return [xi for x in self._parameters_ for xi in x.flattened_parameters] @@ -373,11 +353,6 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable): def _parameter_sizes_(self): return [x.size for x in self._parameters_] @property - def size_transformed(self): - if self._has_fixes(): - return sum(self._fixes_) - return self.size - @property def parameter_shapes(self): return [xi for x in self._parameters_ for xi in x.parameter_shapes] @property diff --git a/GPy/core/parameterization/priors.py b/GPy/core/parameterization/priors.py index 906fe003..29adc923 100644 --- a/GPy/core/parameterization/priors.py +++ b/GPy/core/parameterization/priors.py @@ -64,6 +64,36 @@ class Gaussian(Prior): return np.random.randn(n) * self.sigma + self.mu +class Uniform(Prior): + domain = _REAL + _instances = [] + def __new__(cls, lower, upper): # Singleton: + if cls._instances: + cls._instances[:] = [instance for instance in cls._instances if instance()] + for instance in cls._instances: + if instance().lower == lower and instance().upper == upper: + return instance() + o = super(Prior, cls).__new__(cls, lower, upper) + cls._instances.append(weakref.ref(o)) + return cls._instances[-1]() + + def __init__(self, lower, upper): + self.lower = float(lower) + self.upper = float(upper) + + def __str__(self): + return "[" + str(np.round(self.lower)) + ', ' + str(np.round(self.upper)) + ']' + + def lnpdf(self, x): + region = (x>=self.lower) * (x<=self.upper) + return region + + def lnpdf_grad(self, x): + return np.zeros(x.shape) + + def rvs(self, n): + return np.random.uniform(self.lower, self.upper, size=n) + class LogGaussian(Prior): """ Implementation of the univariate *log*-Gaussian probability function, coupled with random variables. diff --git a/GPy/core/parameterization/transformations.py b/GPy/core/parameterization/transformations.py index 36291ca3..60fcc469 100644 --- a/GPy/core/parameterization/transformations.py +++ b/GPy/core/parameterization/transformations.py @@ -6,8 +6,11 @@ import numpy as np from domains import _POSITIVE,_NEGATIVE, _BOUNDED import weakref +import sys +#_lim_val = -np.log(sys.float_info.epsilon) + _exp_lim_val = np.finfo(np.float64).max -_lim_val = np.log(_exp_lim_val)#-np.log(sys.float_info.epsilon) +_lim_val = np.log(_exp_lim_val)# #=============================================================================== # Fixing constants @@ -35,7 +38,6 @@ class Transformation(object): """ produce a sensible initial value for f(x)""" raise NotImplementedError def plot(self, xlabel=r'transformed $\theta$', ylabel=r'$\theta$', axes=None, *args,**kw): - import sys assert "matplotlib" in sys.modules, "matplotlib package has not been imported." import matplotlib.pyplot as plt from ...plotting.matplot_dep import base_plots @@ -52,7 +54,7 @@ class Transformation(object): class Logexp(Transformation): domain = _POSITIVE def f(self, x): - return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -np.inf, _lim_val)))) + return np.where(x>_lim_val, x, np.log(1. + np.exp(np.clip(x, -_lim_val, _lim_val)))) #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x))) def finv(self, f): return np.where(f>_lim_val, f, np.log(np.exp(f) - 1.)) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 751a295d..b06ffbc7 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -85,11 +85,11 @@ class SparseGP(GP): self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) - def _raw_predict(self, Xnew, X_variance_new=None, full_cov=False): + def _raw_predict(self, Xnew, full_cov=False): """ Make a prediction for the latent function values """ - if X_variance_new is None: + if not isinstance(Xnew, VariationalPosterior): Kx = self.kern.K(self.Z, Xnew) mu = np.dot(Kx.T, self.posterior.woodbury_vector) if full_cov: @@ -100,13 +100,13 @@ class SparseGP(GP): Kxx = self.kern.Kdiag(Xnew) var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T else: - Kx = self.kern.psi1(self.Z, Xnew, X_variance_new) - mu = np.dot(Kx, self.Cpsi1V) + Kx = self.kern.psi1(self.Z, Xnew) + mu = np.dot(Kx, self.posterior.woodbury_vector) if full_cov: raise NotImplementedError, "TODO" else: - Kxx = self.kern.psi0(self.Z, Xnew, X_variance_new) - psi2 = self.kern.psi2(self.Z, Xnew, X_variance_new) + Kxx = self.kern.psi0(self.Z, Xnew) + psi2 = self.kern.psi2(self.Z, Xnew) var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1) return mu, var diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py index 2044f08d..9ebb54a2 100644 --- a/GPy/examples/dimensionality_reduction.py +++ b/GPy/examples/dimensionality_reduction.py @@ -187,10 +187,10 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False): _np.random.seed(1234) x = _np.linspace(0, 4 * _np.pi, N)[:, None] - s1 = _np.vectorize(lambda x: -_np.sin(_np.exp(x))) + s1 = _np.vectorize(lambda x: _np.sin(x)) s2 = _np.vectorize(lambda x: _np.cos(x)**2) s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x))) - sS = _np.vectorize(lambda x: x*_np.sin(x)) + sS = _np.vectorize(lambda x: _np.cos(x)) s1 = s1(x) s2 = s2(x) @@ -202,7 +202,7 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False): s3 -= s3.mean(); s3 /= s3.std(0) sS -= sS.mean(); sS /= sS.std(0) - S1 = _np.hstack([s1, s2, sS]) + S1 = _np.hstack([s1, sS]) S2 = _np.hstack([s2, s3, sS]) S3 = _np.hstack([s3, sS]) @@ -270,7 +270,7 @@ def bgplvm_simulation(optimize=True, verbose=1, from GPy import kern from GPy.models import BayesianGPLVM - D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 5, 9 + D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) @@ -294,7 +294,7 @@ def bgplvm_simulation_missing_data(optimize=True, verbose=1, from GPy.models import BayesianGPLVM from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData - D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 5, 9 + D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 7, 9 _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim) Y = Ylist[0] k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py index fec61204..64707298 100644 --- a/GPy/inference/latent_function_inference/var_dtc.py +++ b/GPy/inference/latent_function_inference/var_dtc.py @@ -60,8 +60,7 @@ class VarDTC(object): _, output_dim = Y.shape #see whether we've got a different noise variance for each datum - beta = 1./np.squeeze(likelihood.variance) - + beta = 1./np.fmax(likelihood.variance, 1e-6) # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency! #self.YYTfactor = self.get_YYTfactor(Y) #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) @@ -214,7 +213,7 @@ class VarDTCMissingData(object): psi2_all = None Ys, traces = self._Y(Y) - beta_all = 1./likelihood.variance + beta_all = 1./np.fmax(likelihood.variance, 1e-6) het_noise = beta_all.size != 1 import itertools diff --git a/GPy/kern/_src/kern.py b/GPy/kern/_src/kern.py index eb3291e0..14e6ae49 100644 --- a/GPy/kern/_src/kern.py +++ b/GPy/kern/_src/kern.py @@ -112,10 +112,12 @@ class Kern(Parameterized): """ assert isinstance(other, Kern), "only kernels can be added to kernels..." from add import Add - return Add([self, other], tensor) - - def __call__(self, X, X2=None): - return self.K(X, X2) + kernels = [] + if not tensor and isinstance(self, Add): kernels.extend(self._parameters_) + else: kernels.append(self) + if not tensor and isinstance(other, Add): kernels.extend(other._parameters_) + else: kernels.append(other) + return Add(kernels, tensor) def __mul__(self, other): """ Here we overload the '*' operator. See self.prod for more information""" diff --git a/GPy/kern/_src/rbf.py b/GPy/kern/_src/rbf.py index 24b70671..baa5b932 100644 --- a/GPy/kern/_src/rbf.py +++ b/GPy/kern/_src/rbf.py @@ -20,7 +20,7 @@ class RBF(Stationary): """ - def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='RBF'): + def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='rbf'): super(RBF, self).__init__(input_dim, variance, lengthscale, ARD, name) self.weave_options = {} @@ -48,7 +48,7 @@ class RBF(Stationary): if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): psi2, _, _, _, _, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob) else: - _, _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) + _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) return psi2 def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): @@ -70,37 +70,39 @@ class RBF(Stationary): self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0) return - l2 = self.lengthscale **2 + elif isinstance(variational_posterior, variational.NormalPosterior): + + l2 = self.lengthscale **2 - #contributions from psi0: - self.variance.gradient = np.sum(dL_dpsi0) - self.lengthscale.gradient = 0. + #contributions from psi0: + self.variance.gradient = np.sum(dL_dpsi0) + self.lengthscale.gradient = 0. + + #from psi1 + denom, _, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + d_length = psi1[:,:,None] * ((dist_sq - 1.)/(self.lengthscale*denom) +1./self.lengthscale) + dpsi1_dlength = d_length * dL_dpsi1[:, :, None] + if self.ARD: + self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) + else: + self.lengthscale.gradient += dpsi1_dlength.sum() + self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance + + #from psi2 + S = variational_posterior.variance + _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * (2.*S[:,None,None,:]/l2 + 1.) + mudist_sq + S[:, None, None, :] / l2) / (2.*S[:,None,None,:] + l2)*self.lengthscale + + dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] + if not self.ARD: + self.lengthscale.gradient += dpsi2_dlength.sum() + else: + self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) + + self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance - #from psi1 - denom, _, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) - d_length = psi1[:,:,None] * ((dist_sq - 1.)/(self.lengthscale*denom) +1./self.lengthscale) - dpsi1_dlength = d_length * dL_dpsi1[:, :, None] - if not self.ARD: - self.lengthscale.gradient += dpsi1_dlength.sum() else: - self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) - self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance - - #from psi2 - S = variational_posterior.variance - denom, _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - d_length = 2.*psi2[:, :, :, None] * (Zdist_sq * denom + mudist_sq + S[:, None, None, :] / l2) / (self.lengthscale * denom) - #TODO: combine denom and l2 as denom_l2?? - #TODO: tidy the above! - #TODO: tensordot below? - - dpsi2_dlength = d_length * dL_dpsi2[:, :, :, None] - if not self.ARD: - self.lengthscale.gradient += dpsi2_dlength.sum() - else: - self.lengthscale.gradient += dpsi2_dlength.sum(0).sum(0).sum(0) - - self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance + raise ValueError, "unknown distriubtion received for psi-statistics" def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): # Spike-and-Slab GPLVM @@ -115,23 +117,26 @@ class RBF(Stationary): grad += (dL_dpsi2[:, :, :, None] * _dpsi2_dZ).sum(axis=0).sum(axis=1) return grad - - l2 = self.lengthscale **2 - #psi1 - denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) - denominator = l2 * denom - dpsi1_dZ = -psi1[:, :, None] * (dist / denominator) - grad = np.sum(dL_dpsi1[:, :, None] * dpsi1_dZ, 0) + elif isinstance(variational_posterior, variational.NormalPosterior): + + l2 = self.lengthscale **2 - #psi2 - denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - term1 = Zdist / l2 # M, M, Q - term2 = mudist / denom / l2 # N, M, M, Q - dZ = psi2[:, :, :, None] * (term1[None, :, :, :] + term2) #N,M,M,Q - grad += 2*(dL_dpsi2[:, :, :, None] * dZ).sum(0).sum(0) + #psi1 + denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + grad = np.einsum('ij,ij,ijk,ijk->jk', dL_dpsi1, psi1, dist, -1./(denom*l2)) - return grad + #psi2 + Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + term1 = Zdist / l2 # M, M, Q + S = variational_posterior.variance + term2 = mudist / (2.*S[:,None,None,:] + l2) # N, M, M, Q + + grad += 2.*np.einsum('ijk,ijk,ijkl->kl', dL_dpsi2, psi2, term1[None,:,:,:] + term2) + + return grad + else: + raise ValueError, "unknown distriubtion received for psi-statistics" def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): # Spike-and-Slab GPLVM @@ -151,18 +156,24 @@ class RBF(Stationary): grad_gamma += (dL_dpsi2[:,:,:, None] * _dpsi2_dgamma).reshape(ndata,-1,self.input_dim).sum(axis=1) return grad_mu, grad_S, grad_gamma + + elif isinstance(variational_posterior, variational.NormalPosterior): - l2 = self.lengthscale **2 - #psi1 - denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) - tmp = psi1[:, :, None] / l2 / denom - grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) - grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) - #psi2 - denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) - tmp = psi2[:, :, :, None] / l2 / denom - grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) - grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) + l2 = self.lengthscale **2 + #psi1 + denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) + tmp = psi1[:, :, None] / l2 / denom + grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) + grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) + #psi2 + _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) + S = variational_posterior.variance + tmp = psi2[:, :, :, None] / (2.*S[:,None,None,:] + l2) + grad_mu += -2.*np.einsum('ijk,ijkl,ijkl->il', dL_dpsi2, tmp , mudist) + grad_S += np.einsum('ijk,ijkl,ijkl->il', dL_dpsi2 , tmp , (2.*mudist_sq - 1)) + + else: + raise ValueError, "unknown distriubtion received for psi-statistics" return grad_mu, grad_S @@ -170,61 +181,6 @@ class RBF(Stationary): # Precomputations # #---------------------------------------# - #TODO: this function is unused, but it will be useful in the stationary class - def _dL_dlengthscales_via_K(self, dL_dK, X, X2): - """ - A helper function for update_gradients_* methods - - Computes the derivative of the objective L wrt the lengthscales via - - dL_dl = sum_{i,j}(dL_dK_{ij} dK_dl) - - assumes self._K_computations has just been called. - - This is only valid if self.ARD=True - """ - target = np.zeros(self.input_dim) - dvardLdK = self._K_dvar * dL_dK - var_len3 = self.variance / np.power(self.lengthscale, 3) - if X2 is None: - # save computation for the symmetrical case - dvardLdK = dvardLdK + dvardLdK.T - code = """ - int q,i,j; - double tmp; - for(q=0; qi', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)] ret /= self.lengthscale**2 return ret @@ -214,7 +245,7 @@ class Matern52(Stationary): .. math:: - k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) \ \ \ \ \ \\text{ where } r = \sqrt{\sum_{i=1}^input_dim \\frac{(x_i-y_i)^2}{\ell_i^2} } + k(r) = \sigma^2 (1 + \sqrt{5} r + \\frac53 r^2) \exp(- \sqrt{5} r) """ def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat52'): super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, name) @@ -225,7 +256,7 @@ class Matern52(Stationary): def dK_dr(self, r): return self.variance*(10./3*r -5.*r -5.*np.sqrt(5.)/3*r**2)*np.exp(-np.sqrt(5.)*r) - def Gram_matrix(self,F,F1,F2,F3,lower,upper): + def Gram_matrix(self, F, F1, F2, F3, lower, upper): """ Return the Gram matrix of the vector of functions F with respect to the RKHS norm. The use of this function is limited to input_dim=1. diff --git a/GPy/models/sparse_gp_regression.py b/GPy/models/sparse_gp_regression.py index 6a76df3f..4980d26a 100644 --- a/GPy/models/sparse_gp_regression.py +++ b/GPy/models/sparse_gp_regression.py @@ -8,7 +8,7 @@ from .. import likelihoods from .. import kern from ..inference.latent_function_inference import VarDTC from ..util.misc import param_to_array -from ..core.parameterization.variational import VariationalPosterior +from ..core.parameterization.variational import NormalPosterior class SparseGPRegression(SparseGP): """ @@ -47,7 +47,7 @@ class SparseGPRegression(SparseGP): likelihood = likelihoods.Gaussian() if not (X_variance is None): - X = VariationalPosterior(X,X_variance) + X = NormalPosterior(X,X_variance) SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC()) diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py index d72d2a3e..4ca4441e 100644 --- a/GPy/plotting/matplot_dep/models_plots.py +++ b/GPy/plotting/matplot_dep/models_plots.py @@ -56,10 +56,13 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', if ax is None: fig = pb.figure(num=fignum) ax = fig.add_subplot(111) - - X, Y = param_to_array(model.X, model.Y) - if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): X_variance = model.X_variance - + + if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): + X = model.X.mean + X_variance = param_to_array(model.X.variance) + else: + X = model.X + X, Y = param_to_array(X, model.Y) if hasattr(model, 'Z'): Z = param_to_array(model.Z) #work out what the inputs are for plotting (1D or 2D) @@ -98,10 +101,10 @@ def plot_fit(model, plot_limits=None, which_data_rows='all', #add error bars for uncertain (if input uncertainty is being modelled) - #if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): - # ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), - # xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()), - # ecolor='k', fmt=None, elinewidth=.5, alpha=.5) + if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): + ax.errorbar(X[which_data_rows, free_dims].flatten(), Y[which_data_rows, which_data_ycols].flatten(), + xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()), + ecolor='k', fmt=None, elinewidth=.5, alpha=.5) #set the limits of the plot to some sensible values diff --git a/GPy/testing/observable_tests.py b/GPy/testing/observable_tests.py new file mode 100644 index 00000000..6b4f1a87 --- /dev/null +++ b/GPy/testing/observable_tests.py @@ -0,0 +1,133 @@ +''' +Created on 27 Feb 2014 + +@author: maxz +''' +import unittest +from GPy.core.parameterization.parameterized import Parameterized +from GPy.core.parameterization.param import Param +import numpy + + +class ParamTestParent(Parameterized): + parent_changed_count = 0 + def parameters_changed(self): + self.parent_changed_count += 1 + +class ParameterizedTest(Parameterized): + params_changed_count = 0 + def parameters_changed(self): + self.params_changed_count += 1 + def _set_params(self, params, trigger_parent=True): + Parameterized._set_params(self, params, trigger_parent=trigger_parent) + +class Test(unittest.TestCase): + + def setUp(self): + self.parent = ParamTestParent('test parent') + self.par = ParameterizedTest('test model') + self.par2 = ParameterizedTest('test model 2') + self.p = Param('test parameter', numpy.random.normal(1,2,(10,3))) + + self.par.add_parameter(self.p) + self.par.add_parameter(Param('test1', numpy.random.normal(0,1,(1,)))) + self.par.add_parameter(Param('test2', numpy.random.normal(0,1,(1,)))) + + self.par2.add_parameter(Param('par2 test1', numpy.random.normal(0,1,(1,)))) + self.par2.add_parameter(Param('par2 test2', numpy.random.normal(0,1,(1,)))) + + self.parent.add_parameter(self.par) + self.parent.add_parameter(self.par2) + + self._observer_triggered = None + self._trigger_count = 0 + self._first = None + self._second = None + + def _trigger(self, which): + self._observer_triggered = float(which) + self._trigger_count += 1 + if self._first is not None: + self._second = self._trigger + else: + self._first = self._trigger + + def _trigger_priority(self, which): + if self._first is not None: + self._second = self._trigger_priority + else: + self._first = self._trigger_priority + + def test_observable(self): + self.par.add_observer(self, self._trigger, -1) + self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.p[0,1] = 3 # trigger observers + self.assertEqual(self._observer_triggered, 3, 'observer should have triggered') + self.assertEqual(self._trigger_count, 1, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 1, 'params changed once') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.par.remove_observer(self) + self.p[2,1] = 4 + self.assertEqual(self._observer_triggered, 3, 'observer should not have triggered') + self.assertEqual(self._trigger_count, 1, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 2, 'params changed second') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.par.add_observer(self, self._trigger, -1) + self.p[2,1] = 4 + self.assertEqual(self._observer_triggered, 4, 'observer should have triggered') + self.assertEqual(self._trigger_count, 2, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 3, 'params changed second') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + self.par.remove_observer(self, self._trigger) + self.p[0,1] = 3 + self.assertEqual(self._observer_triggered, 4, 'observer should not have triggered') + self.assertEqual(self._trigger_count, 2, 'observer should have triggered once') + self.assertEqual(self.par.params_changed_count, 4, 'params changed second') + self.assertEqual(self.par.params_changed_count, self.parent.parent_changed_count, 'parent should be triggered as often as param') + + def test_set_params(self): + self.assertEqual(self.par.params_changed_count, 0, 'no params changed yet') + self.par._set_params(numpy.ones(self.par.size)) + self.assertEqual(self.par.params_changed_count, 1, 'now params changed') + self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) + + self.parent._set_params(numpy.ones(self.parent.size) * 2) + self.assertEqual(self.par.params_changed_count, 2, 'now params changed') + self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) + + + def test_priority_notify(self): + self.assertEqual(self.par.params_changed_count, 0) + self.par._notify_observers(0, None) + self.assertEqual(self.par.params_changed_count, 1) + self.assertEqual(self.parent.parent_changed_count, self.par.params_changed_count) + + self.par._notify_observers(0, -numpy.inf) + self.assertEqual(self.par.params_changed_count, 2) + self.assertEqual(self.parent.parent_changed_count, 1) + + def test_priority(self): + self.par.add_observer(self, self._trigger, -1) + self.par.add_observer(self, self._trigger_priority, 0) + self.par._notify_observers(0) + self.assertEqual(self._first, self._trigger_priority, 'priority should be first') + self.assertEqual(self._second, self._trigger, 'priority should be first') + + self.par.remove_observer(self) + self._first = self._second = None + + self.par.add_observer(self, self._trigger, 1) + self.par.add_observer(self, self._trigger_priority, 0) + self.par._notify_observers(0) + self.assertEqual(self._first, self._trigger, 'priority should be second') + self.assertEqual(self._second, self._trigger_priority, 'priority should be second') + + +if __name__ == "__main__": + #import sys;sys.argv = ['', 'Test.testName'] + unittest.main() \ No newline at end of file diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py index 6f13d294..0da3f3ae 100644 --- a/GPy/testing/parameterized_tests.py +++ b/GPy/testing/parameterized_tests.py @@ -6,6 +6,7 @@ Created on Feb 13, 2014 import unittest import GPy import numpy as np +from GPy.core.parameterization.parameter_core import HierarchyError class Test(unittest.TestCase): @@ -65,7 +66,7 @@ class Test(unittest.TestCase): self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1]) def test_add_parameter_already_in_hirarchy(self): - self.test1.add_parameter(self.white._parameters_[0]) + 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) diff --git a/GPy/util/caching.py b/GPy/util/caching.py index 2899cb33..250efe11 100644 --- a/GPy/util/caching.py +++ b/GPy/util/caching.py @@ -1,4 +1,5 @@ from ..core.parameterization.parameter_core import Observable +import itertools class Cacher(object): """ @@ -38,8 +39,11 @@ class Cacher(object): if not all([isinstance(arg, Observable) for arg in observable_args]): return self.operation(*args) + # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING + # return self.operation(*args) + #if the result is cached, return the cached computation - state = [all(a is b for a, b in zip(args, cached_i)) for cached_i in self.cached_inputs] + 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]: