apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-05-24 14:15:14 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

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

Browse source
This commit is contained in:
Zhenwen Dai 2014-02-28 18:11:05 +00:00
commit c96d9ffc4c
20 changed files with 612 additions and 405 deletions
Download patch file Download diff file Expand all files Collapse all files

26
GPy/core/model.py
View file

@ -60,20 +60,6 @@ class Model(Parameterized):
self.priors = state.pop() self.priors = state.pop()
Parameterized._setstate(self, state) 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): 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 Perform random restarts of the model, and set the model to the best
@ -240,6 +226,11 @@ class Model(Parameterized):
TODO: valid args 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: if optimizer is None:
optimizer = self.preferred_optimizer optimizer = self.preferred_optimizer
@ -279,7 +270,7 @@ class Model(Parameterized):
and numerical gradients is within <tolerance> of unity. and numerical gradients is within <tolerance> of unity.
""" """
x = self._get_params_transformed().copy() x = self._get_params_transformed()
if not verbose: if not verbose:
# make sure only to test the selected parameters # make sure only to test the selected parameters
@ -297,7 +288,7 @@ class Model(Parameterized):
return return
# just check the global ratio # 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)) dx[transformed_index] = step * np.sign(np.random.uniform(-1, 1, transformed_index.size))
# evaulate around the point x # evaulate around the point x
@ -308,9 +299,8 @@ class Model(Parameterized):
dx = dx[transformed_index] dx = dx[transformed_index]
gradient = gradient[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))) 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: else:
# check the gradient of each parameter individually, and do some pretty printing # check the gradient of each parameter individually, and do some pretty printing
try: try:

16
GPy/core/parameterization/array_core.py
View file

@ -6,19 +6,6 @@ __updated__ = '2013-12-16'
import numpy as np import numpy as np
from parameter_core import Observable 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): class ObservableArray(np.ndarray, Observable):
""" """
An ndarray which reports changes to its observers. An ndarray which reports changes to its observers.
@ -62,10 +49,11 @@ class ObservableArray(np.ndarray, Observable):
def __setitem__(self, s, val): def __setitem__(self, s, val):
if self._s_not_empty(s): if self._s_not_empty(s):
super(ObservableArray, self).__setitem__(s, val) super(ObservableArray, self).__setitem__(s, val)
self._notify_observers() self._notify_observers(self[s])
def __getslice__(self, start, stop): def __getslice__(self, start, stop):
return self.__getitem__(slice(start, stop)) return self.__getitem__(slice(start, stop))
def __setslice__(self, start, stop, val): def __setslice__(self, start, stop, val):
return self.__setitem__(slice(start, stop), val) return self.__setitem__(slice(start, stop), val)

52
GPy/core/parameterization/index_operations.py
View file

@ -5,47 +5,7 @@ Created on Oct 2, 2013
''' '''
import numpy import numpy
from numpy.lib.function_base import vectorize from numpy.lib.function_base import vectorize
from param import Param from lists_and_dicts import IntArrayDict
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_([])
class ParameterIndexOperations(object): class ParameterIndexOperations(object):
''' '''
@ -194,8 +154,12 @@ class ParameterIndexOperationsView(object):
def shift_right(self, start, size): 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): def clear(self):
for i, ind in self.items(): for i, ind in self.items():
@ -232,9 +196,7 @@ class ParameterIndexOperationsView(object):
def __getitem__(self, prop): def __getitem__(self, prop):
ind = self._filter_index(self._param_index_ops[prop]) ind = self._filter_index(self._param_index_ops[prop])
if ind.size > 0: return ind
return ind
raise KeyError, prop
def __str__(self, *args, **kwargs): def __str__(self, *args, **kwargs):
import pprint import pprint

35
GPy/core/parameterization/lists_and_dicts.py Normal file
View file

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

65
GPy/core/parameterization/param.py
View file

@ -3,8 +3,8 @@
import itertools import itertools
import numpy import numpy
from parameter_core import Constrainable, Gradcheckable, Indexable, Parentable, adjust_name_for_printing from parameter_core import OptimizationHandlable, Gradcheckable, adjust_name_for_printing
from array_core import ObservableArray, ParamList from array_core import ObservableArray
###### printing ###### printing
__constraints_name__ = "Constraint" __constraints_name__ = "Constraint"
@ -15,7 +15,7 @@ __precision__ = numpy.get_printoptions()['precision'] # numpy printing precision
__print_threshold__ = 5 __print_threshold__ = 5
###### ######
class Param(Constrainable, ObservableArray, Gradcheckable): class Param(OptimizationHandlable, ObservableArray, Gradcheckable):
""" """
Parameter object for GPy models. Parameter object for GPy models.
@ -50,7 +50,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable):
obj._realsize_ = obj.size obj._realsize_ = obj.size
obj._realndim_ = obj.ndim obj._realndim_ = obj.ndim
obj._updated_ = False obj._updated_ = False
from index_operations import SetDict from lists_and_dicts import SetDict
obj._tied_to_me_ = SetDict() obj._tied_to_me_ = SetDict()
obj._tied_to_ = [] obj._tied_to_ = []
obj._original_ = True obj._original_ = True
@ -148,8 +148,11 @@ class Param(Constrainable, ObservableArray, Gradcheckable):
#=========================================================================== #===========================================================================
# get/set parameters # get/set parameters
#=========================================================================== #===========================================================================
def _set_params(self, param, update=True): def _set_params(self, param, trigger_parent=True):
self.flat = param self.flat = param
if trigger_parent: min_priority = None
else: min_priority = -numpy.inf
self._notify_observers(None, min_priority)
def _get_params(self): def _get_params(self):
return self.flat 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 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 except AttributeError: pass # returning 0d array or float, double etc
return new_arr return new_arr
def __setitem__(self, s, val): def __setitem__(self, s, val):
super(Param, self).__setitem__(s, val) super(Param, self).__setitem__(s, val)
if self.has_parent():
self._direct_parent_._notify_parameters_changed()
#self._notify_observers()
#=========================================================================== #===========================================================================
# Index Operations: # Index Operations:
@ -204,6 +205,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable):
ind = self._indices(slice_index) ind = self._indices(slice_index)
if ind.ndim < 2: ind = ind[:, None] 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) 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): def _expand_index(self, slice_index=None):
# this calculates the full indexing arrays from the slicing objects given by get_item for _real..._ attributes # 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 # it basically translates slices to their respective index arrays and turns negative indices around
@ -230,7 +232,8 @@ class Param(Constrainable, ObservableArray, Gradcheckable):
#=========================================================================== #===========================================================================
@property @property
def is_fixed(self): 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): #def round(self, decimals=0, out=None):
# view = super(Param, self).round(decimals, out).view(Param) # view = super(Param, self).round(decimals, out).view(Param)
# view.__array_finalize__(self) # view.__array_finalize__(self)
@ -267,7 +270,7 @@ class Param(Constrainable, ObservableArray, Gradcheckable):
return [t._short() for t in self._tied_to_] or [''] return [t._short() for t in self._tied_to_] or ['']
def __repr__(self, *args, **kwargs): def __repr__(self, *args, **kwargs):
name = "\033[1m{x:s}\033[0;0m:\n".format( 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) return name + super(Param, self).__repr__(*args, **kwargs)
def _ties_for(self, rav_index): def _ties_for(self, rav_index):
# size = sum(p.size for p in self._tied_to_) # 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) gen = map(lambda x: " ".join(map(str, x)), gen)
return reduce(lambda a, b:max(a, len(b)), gen, len(header)) return reduce(lambda a, b:max(a, len(b)), gen, len(header))
def _max_len_values(self): 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): def _max_len_index(self, ind):
return reduce(lambda a, b:max(a, len(str(b))), ind, len(__index_name__)) return reduce(lambda a, b:max(a, len(str(b))), ind, len(__index_name__))
def _short(self): def _short(self):
# short string to print # short string to print
name = self.hirarchy_name() name = self.hierarchy_name()
if self._realsize_ < 2: if self._realsize_ < 2:
return name return name
ind = self._indices() ind = self._indices()
@ -329,8 +332,8 @@ class Param(Constrainable, ObservableArray, Gradcheckable):
if lp is None: lp = self._max_len_names(prirs, __tie_name__) if lp is None: lp = self._max_len_names(prirs, __tie_name__)
sep = '-' 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}" 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 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.hirarchy_name(), c=__constraints_name__, i=__index_name__, t=__tie_name__, p=__priors_name__) # 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(['']) 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 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__() # 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. See :py:class:`GPy.core.parameter.Param` for more details on constraining.
""" """
# self.params = params # self.params = params
self.params = ParamList([]) from lists_and_dicts import ArrayList
self.params = ArrayList([])
for p in params: for p in params:
for p in p.flattened_parameters: for p in p.flattened_parameters:
if p not in self.params: if p not in self.params:
@ -353,6 +357,21 @@ class ParamConcatenation(object):
self._param_sizes = [p.size for p in self.params] self._param_sizes = [p.size for p in self.params]
startstops = numpy.cumsum([0] + self._param_sizes) startstops = numpy.cumsum([0] + self._param_sizes)
self._param_slices_ = [slice(start, stop) for start,stop in zip(startstops, startstops[1:])] 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 # Get/set items, enable broadcasting
#=========================================================================== #===========================================================================
@ -366,24 +385,26 @@ class ParamConcatenation(object):
val = val._vals() val = val._vals()
ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True; ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True;
vals = self._vals(); vals[s] = val; del val 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_)] for p, ps in zip(self.params, self._param_slices_)]
if update:
self.update_all_params()
def _vals(self): def _vals(self):
return numpy.hstack([p._get_params() for p in self.params]) return numpy.hstack([p._get_params() for p in self.params])
#=========================================================================== #===========================================================================
# parameter operations: # parameter operations:
#=========================================================================== #===========================================================================
def update_all_params(self): def update_all_params(self):
for p in self.params: for par in self.parents:
p._notify_observers() par._notify_observers(-numpy.inf)
def constrain(self, constraint, warning=True): 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() self.update_all_params()
constrain.__doc__ = Param.constrain.__doc__ constrain.__doc__ = Param.constrain.__doc__
def constrain_positive(self, warning=True): 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() self.update_all_params()
constrain_positive.__doc__ = Param.constrain_positive.__doc__ constrain_positive.__doc__ = Param.constrain_positive.__doc__
@ -393,12 +414,12 @@ class ParamConcatenation(object):
fix = constrain_fixed fix = constrain_fixed
def constrain_negative(self, warning=True): 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() self.update_all_params()
constrain_negative.__doc__ = Param.constrain_negative.__doc__ constrain_negative.__doc__ = Param.constrain_negative.__doc__
def constrain_bounded(self, lower, upper, warning=True): 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() self.update_all_params()
constrain_bounded.__doc__ = Param.constrain_bounded.__doc__ constrain_bounded.__doc__ = Param.constrain_bounded.__doc__

191
GPy/core/parameterization/parameter_core.py
View file

@ -2,34 +2,58 @@
# Licensed under the BSD 3-clause license (see LICENSE.txt) # Licensed under the BSD 3-clause license (see LICENSE.txt)
from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED
import heapq
__updated__ = '2013-12-16' __updated__ = '2013-12-16'
class HierarchyError(Exception):
"""
Gets thrown when something is wrong with the parameter hierarchy
"""
def adjust_name_for_printing(name): def adjust_name_for_printing(name):
if name is not None: if name is not None:
return name.replace(" ", "_").replace(".", "_").replace("-", "").replace("+", "").replace("!", "").replace("*", "").replace("/", "") return name.replace(" ", "_").replace(".", "_").replace("-", "").replace("+", "").replace("!", "").replace("*", "").replace("/", "")
return '' return ''
class Observable(object): class Observable(object):
_updated = True
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
from collections import defaultdict self._observer_callables_ = []
self._observer_callables_ = defaultdict(list)
def add_observer(self, observer, callble): def add_observer(self, observer, callble, priority=0):
self._observer_callables_[observer].append(callble) heapq.heappush(self._observer_callables_, (priority, observer, callble))
def remove_observer(self, observer, callble=None): def remove_observer(self, observer, callble=None):
if observer in self._observer_callables_: to_remove = []
if callble is None: for p, obs, clble in self._observer_callables_:
del self._observer_callables_[observer] if callble is not None:
elif callble in self._observer_callables_[observer]: if (obs == observer) and (callble == clble):
self._observer_callables_[observer].remove(callble) to_remove.append((p, obs, clble))
if len(self._observer_callables_[observer]) == 0: else:
self.remove_observer(observer) if obs is observer:
to_remove.append((p, obs, clble))
for r in to_remove:
self._observer_callables_.remove(r)
def _notify_observers(self): def _notify_observers(self, which=None, min_priority=None):
[[callble(self) for callble in callables] """
for callables in self._observer_callables_.itervalues()] 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): class Pickleable(object):
def _getstate(self): def _getstate(self):
@ -95,11 +119,11 @@ class Nameable(Parentable):
self._name = name self._name = name
if self.has_parent(): if self.has_parent():
self._direct_parent_._name_changed(self, from_name) 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) if adjust_for_printing: adjust = lambda x: adjust_name_for_printing(x)
else: adjust = lambda x: x else: adjust = lambda x: x
if self.has_parent(): 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) return adjust(self.name)
@ -156,7 +180,7 @@ class Constrainable(Nameable, Indexable):
#=========================================================================== #===========================================================================
# Fixing Parameters: # 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. 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: if value is not None:
self[:] = value 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) rav_i = self._highest_parent_._raveled_index_for(self)
self._highest_parent_._set_fixed(rav_i) self._highest_parent_._set_fixed(rav_i)
fix = constrain_fixed fix = constrain_fixed
@ -205,9 +229,9 @@ class Constrainable(Nameable, Indexable):
#=========================================================================== #===========================================================================
# Prior Operations # 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() 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): def unset_priors(self, *priors):
return self._remove_from_index_operations(self.priors, priors) return self._remove_from_index_operations(self.priors, priors)
@ -233,7 +257,7 @@ class Constrainable(Nameable, Indexable):
# Constrain operations -> done # 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` :param transform: the :py:class:`GPy.core.transformations.Transformation`
to constrain the this parameter to. to constrain the this parameter to.
@ -243,9 +267,9 @@ class Constrainable(Nameable, Indexable):
:py:class:`GPy.core.transformations.Transformation`. :py:class:`GPy.core.transformations.Transformation`.
""" """
if isinstance(transform, 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() 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): def unconstrain(self, *transforms):
""" """
@ -256,30 +280,30 @@ class Constrainable(Nameable, Indexable):
""" """
return self._remove_from_index_operations(self.constraints, transforms) 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. :param warning: print a warning if re-constraining parameters.
Constrain this parameter to the default positive constraint. 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. :param warning: print a warning if re-constraining parameters.
Constrain this parameter to the default negative constraint. 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 lower, upper: the limits to bound this parameter to
:param warning: print a warning if re-constraining parameters. :param warning: print a warning if re-constraining parameters.
Constrain this parameter to lie within the given range. 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): def unconstrain_positive(self):
""" """
@ -309,12 +333,11 @@ class Constrainable(Nameable, Indexable):
for p in self._parameters_: for p in self._parameters_:
p._parent_changed(parent) 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: 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) print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name)
which.add(transform, self._raveled_index()) which.add(transform, self._raveled_index())
if update:
self._notify_observers()
def _remove_from_index_operations(self, which, transforms): def _remove_from_index_operations(self, which, transforms):
if len(transforms) == 0: if len(transforms) == 0:
@ -329,12 +352,76 @@ class Constrainable(Nameable, Indexable):
return removed 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
class Parameterizable(Constrainable, Observable): 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])
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): def __init__(self, *args, **kwargs):
super(Parameterizable, self).__init__(*args, **kwargs) super(Parameterizable, self).__init__(*args, **kwargs)
from GPy.core.parameterization.array_core import ParamList from GPy.core.parameterization.lists_and_dicts import ArrayList
_parameters_ = ParamList() _parameters_ = ArrayList()
self._added_names_ = set() self._added_names_ = set()
def parameter_names(self, add_self=False, adjust_for_printing=False, recursive=True): 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_: if pname in self._added_names_:
del self.__dict__[pname] del self.__dict__[pname]
self._add_parameter_name(param) self._add_parameter_name(param)
else: elif pname not in dir(self):
self.__dict__[pname] = param self.__dict__[pname] = param
self._added_names_.add(pname) self._added_names_.add(pname)
@ -377,28 +464,26 @@ class Parameterizable(Constrainable, Observable):
import itertools import itertools
[p._collect_gradient(target[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] [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): def _set_gradient(self, g):
import itertools import itertools
[p._set_gradient(g[s]) for p, s in itertools.izip(self._parameters_, self._param_slices_)] [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): def copy(self):
"""Returns a (deep) copy of the current model""" """Returns a (deep) copy of the current model"""
import copy import copy
from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView
from .array_core import ParamList from .lists_and_dicts import ArrayList
dc = dict() dc = dict()
for k, v in self.__dict__.iteritems(): for k, v in self.__dict__.iteritems():
@ -412,7 +497,7 @@ class Parameterizable(Constrainable, Observable):
dc['_direct_parent_'] = None dc['_direct_parent_'] = None
dc['_parent_index_'] = None dc['_parent_index_'] = None
dc['_parameters_'] = ParamList() dc['_parameters_'] = ArrayList()
dc['constraints'].clear() dc['constraints'].clear()
dc['priors'].clear() dc['priors'].clear()
dc['size'] = 0 dc['size'] = 0
@ -425,12 +510,6 @@ class Parameterizable(Constrainable, Observable):
return s 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): def parameters_changed(self):
""" """
This method gets called when parameters have changed. This method gets called when parameters have changed.

105
GPy/core/parameterization/parameterized.py
View file

@ -7,9 +7,9 @@ import cPickle
import itertools import itertools
from re import compile, _pattern_type from re import compile, _pattern_type
from param import ParamConcatenation 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 transformations import __fixed__
from array_core import ParamList from lists_and_dicts import ArrayList
class Parameterized(Parameterizable, Pickleable, Gradcheckable): class Parameterized(Parameterizable, Pickleable, Gradcheckable):
""" """
@ -56,8 +56,9 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
def __init__(self, name=None, *a, **kw): def __init__(self, name=None, *a, **kw):
super(Parameterized, self).__init__(name=name, parent=None, parent_index=None, *a, **kw) super(Parameterized, self).__init__(name=name, parent=None, parent_index=None, *a, **kw)
self._in_init_ = True self._in_init_ = True
self._parameters_ = ParamList() self._parameters_ = ArrayList()
self.size = sum(p.size for p in self._parameters_) self.size = sum(p.size for p in self._parameters_)
self.add_observer(self, self._parameters_changed_notification, -100)
if not self._has_fixes(): if not self._has_fixes():
self._fixes_ = None self._fixes_ = None
self._param_slices_ = [] self._param_slices_ = []
@ -65,7 +66,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
del self._in_init_ del self._in_init_
def build_pydot(self, G=None): def build_pydot(self, G=None):
import pydot import pydot # @UnresolvedImport
iamroot = False iamroot = False
if G is None: if G is None:
G = pydot.Dot(graph_type='digraph') G = pydot.Dot(graph_type='digraph')
@ -104,6 +105,14 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
self.remove_parameter(param) self.remove_parameter(param)
self.add_parameter(param, index) self.add_parameter(param, index)
elif param not in self._parameters_: 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 # make sure the size is set
if index is None: if index is None:
self.constraints.update(param.constraints, self.size) self.constraints.update(param.constraints, self.size)
@ -116,12 +125,16 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
self.constraints.update(param.constraints, start) self.constraints.update(param.constraints, start)
self.priors.update(param.priors, start) self.priors.update(param.priors, start)
self._parameters_.insert(index, param) self._parameters_.insert(index, param)
param.add_observer(self, self._pass_through_notify_observers, -np.inf)
self.size += param.size self.size += param.size
self._connect_parameters()
self._notify_parent_change()
self._connect_fixes()
else: else:
raise RuntimeError, """Parameter exists already added and no copy made""" 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): def add_parameters(self, *parameters):
@ -144,12 +157,19 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
del self._parameters_[param._parent_index_] del self._parameters_[param._parent_index_]
param._disconnect_parent() 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.constraints.shift_left(start, param.size)
self._connect_fixes() self._connect_fixes()
self._connect_parameters() self._connect_parameters()
self._notify_parent_change() 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): def _connect_parameters(self):
# connect parameterlist to this parameterized object # connect parameterlist to this parameterized object
@ -170,6 +190,13 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
self._add_parameter_name(p) 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 # Pickling operations
#=========================================================================== #===========================================================================
def pickle(self, f, protocol=-1): 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)] g[self._offset_for(p) + numpy.array(list(i))] += g[self._raveled_index_for(t)]
if self._has_fixes(): return g[self._fixes_] if self._has_fixes(): return g[self._fixes_]
return g 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): def _offset_for(self, param):
# get the offset in the parameterized index array for param # get the offset in the parameterized index array for param
if param.has_parent(): if param.has_parent():
@ -297,34 +289,22 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
this is not in the global view of things! this is not in the global view of things!
""" """
return numpy.r_[:self.size] 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: # 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 @property
def is_fixed(self): def is_fixed(self):
for p in self._parameters_: for p in self._parameters_:
if not p.is_fixed: return False if not p.is_fixed: return False
return True return True
def _get_original(self, param): def _get_original(self, param):
# if advanced indexing is activated it happens that the array is a copy # if advanced indexing is activated it happens that the array is a copy
# you can retrieve the original param through this method, by passing # you can retrieve the original param through this method, by passing
# the copy here # the copy here
return self._parameters_[param._parent_index_] return self._parameters_[param._parent_index_]
#=========================================================================== #===========================================================================
# Get/set parameters: # Get/set parameters:
#=========================================================================== #===========================================================================
@ -365,7 +345,7 @@ class Parameterized(Parameterizable, Pickleable, Gradcheckable):
# Printing: # Printing:
#=========================================================================== #===========================================================================
def _short(self): def _short(self):
return self.hirarchy_name() return self.hierarchy_name()
@property @property
def flattened_parameters(self): def flattened_parameters(self):
return [xi for x in self._parameters_ for xi in x.flattened_parameters] 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): def _parameter_sizes_(self):
return [x.size for x in self._parameters_] return [x.size for x in self._parameters_]
@property @property
def size_transformed(self):
if self._has_fixes():
return sum(self._fixes_)
return self.size
@property
def parameter_shapes(self): def parameter_shapes(self):
return [xi for x in self._parameters_ for xi in x.parameter_shapes] return [xi for x in self._parameters_ for xi in x.parameter_shapes]
@property @property

30
GPy/core/parameterization/priors.py
View file

@ -64,6 +64,36 @@ class Gaussian(Prior):
return np.random.randn(n) * self.sigma + self.mu 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): class LogGaussian(Prior):
""" """
Implementation of the univariate *log*-Gaussian probability function, coupled with random variables. Implementation of the univariate *log*-Gaussian probability function, coupled with random variables.

8
GPy/core/parameterization/transformations.py
View file

@ -6,8 +6,11 @@ import numpy as np
from domains import _POSITIVE,_NEGATIVE, _BOUNDED from domains import _POSITIVE,_NEGATIVE, _BOUNDED
import weakref import weakref
import sys
#_lim_val = -np.log(sys.float_info.epsilon)
_exp_lim_val = np.finfo(np.float64).max _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 # Fixing constants
@ -35,7 +38,6 @@ class Transformation(object):
""" produce a sensible initial value for f(x)""" """ produce a sensible initial value for f(x)"""
raise NotImplementedError raise NotImplementedError
def plot(self, xlabel=r'transformed $\theta$', ylabel=r'$\theta$', axes=None, *args,**kw): 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." assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from ...plotting.matplot_dep import base_plots from ...plotting.matplot_dep import base_plots
@ -52,7 +54,7 @@ class Transformation(object):
class Logexp(Transformation): class Logexp(Transformation):
domain = _POSITIVE domain = _POSITIVE
def f(self, x): 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))) #raises overflow warning: return np.where(x>_lim_val, x, np.log(1. + np.exp(x)))
def finv(self, f): def finv(self, f):
return np.where(f>_lim_val, f, np.log(np.exp(f) - 1.)) return np.where(f>_lim_val, f, np.log(np.exp(f) - 1.))

12
GPy/core/sparse_gp.py
View file

@ -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_dKmm'], self.Z)
self.Z.gradient += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X) 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 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) Kx = self.kern.K(self.Z, Xnew)
mu = np.dot(Kx.T, self.posterior.woodbury_vector) mu = np.dot(Kx.T, self.posterior.woodbury_vector)
if full_cov: if full_cov:
@ -100,13 +100,13 @@ class SparseGP(GP):
Kxx = self.kern.Kdiag(Xnew) Kxx = self.kern.Kdiag(Xnew)
var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T var = (Kxx - np.sum(np.dot(np.atleast_3d(self.posterior.woodbury_inv).T, Kx) * Kx[None,:,:], 1)).T
else: else:
Kx = self.kern.psi1(self.Z, Xnew, X_variance_new) Kx = self.kern.psi1(self.Z, Xnew)
mu = np.dot(Kx, self.Cpsi1V) mu = np.dot(Kx, self.posterior.woodbury_vector)
if full_cov: if full_cov:
raise NotImplementedError, "TODO" raise NotImplementedError, "TODO"
else: else:
Kxx = self.kern.psi0(self.Z, Xnew, X_variance_new) Kxx = self.kern.psi0(self.Z, Xnew)
psi2 = self.kern.psi2(self.Z, Xnew, X_variance_new) psi2 = self.kern.psi2(self.Z, Xnew)
var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1) var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1)
return mu, var return mu, var

10
GPy/examples/dimensionality_reduction.py
View file

@ -187,10 +187,10 @@ def _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim=False):
_np.random.seed(1234) _np.random.seed(1234)
x = _np.linspace(0, 4 * _np.pi, N)[:, None] 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) s2 = _np.vectorize(lambda x: _np.cos(x)**2)
s3 = _np.vectorize(lambda x:-_np.exp(-_np.cos(2 * x))) 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) s1 = s1(x)
s2 = s2(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) s3 -= s3.mean(); s3 /= s3.std(0)
sS -= sS.mean(); sS /= sS.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]) S2 = _np.hstack([s2, s3, sS])
S3 = _np.hstack([s3, sS]) S3 = _np.hstack([s3, sS])
@ -270,7 +270,7 @@ def bgplvm_simulation(optimize=True, verbose=1,
from GPy import kern from GPy import kern
from GPy.models import BayesianGPLVM 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) _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
Y = Ylist[0] Y = Ylist[0]
k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) 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.models import BayesianGPLVM
from GPy.inference.latent_function_inference.var_dtc import VarDTCMissingData 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) _, _, Ylist = _simulate_sincos(D1, D2, D3, N, num_inducing, Q, plot_sim)
Y = Ylist[0] Y = Ylist[0]
k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q) k = kern.Linear(Q, ARD=True)# + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)

5
GPy/inference/latent_function_inference/var_dtc.py
View file

@ -60,8 +60,7 @@ class VarDTC(object):
_, output_dim = Y.shape _, output_dim = Y.shape
#see whether we've got a different noise variance for each datum #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! # VVT_factor is a matrix such that tdot(VVT_factor) = VVT...this is for efficiency!
#self.YYTfactor = self.get_YYTfactor(Y) #self.YYTfactor = self.get_YYTfactor(Y)
#VVT_factor = self.get_VVTfactor(self.YYTfactor, beta) #VVT_factor = self.get_VVTfactor(self.YYTfactor, beta)
@ -214,7 +213,7 @@ class VarDTCMissingData(object):
psi2_all = None psi2_all = None
Ys, traces = self._Y(Y) 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 het_noise = beta_all.size != 1
import itertools import itertools

10
GPy/kern/_src/kern.py
View file

@ -112,10 +112,12 @@ class Kern(Parameterized):
""" """
assert isinstance(other, Kern), "only kernels can be added to kernels..." assert isinstance(other, Kern), "only kernels can be added to kernels..."
from add import Add from add import Add
return Add([self, other], tensor) kernels = []
if not tensor and isinstance(self, Add): kernels.extend(self._parameters_)
def __call__(self, X, X2=None): else: kernels.append(self)
return self.K(X, X2) if not tensor and isinstance(other, Add): kernels.extend(other._parameters_)
else: kernels.append(other)
return Add(kernels, tensor)
def __mul__(self, other): def __mul__(self, other):
""" Here we overload the '*' operator. See self.prod for more information""" """ Here we overload the '*' operator. See self.prod for more information"""

182
GPy/kern/_src/rbf.py
View file

@ -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) super(RBF, self).__init__(input_dim, variance, lengthscale, ARD, name)
self.weave_options = {} self.weave_options = {}
@ -48,7 +48,7 @@ class RBF(Stationary):
if isinstance(variational_posterior, variational.SpikeAndSlabPosterior): 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) psi2, _, _, _, _, _, _ = ssrbf_psi_comp._psi2computations(self.variance, self.lengthscale, Z, variational_posterior.mean, variational_posterior.variance, variational_posterior.binary_prob)
else: else:
_, _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior) _, _, _, _, psi2 = self._psi2computations(Z, variational_posterior)
return psi2 return psi2
def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior): 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) self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).reshape(-1,self.input_dim).sum(axis=0)
return return
l2 = self.lengthscale **2 elif isinstance(variational_posterior, variational.NormalPosterior):
#contributions from psi0: l2 = self.lengthscale **2
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: else:
self.lengthscale.gradient += dpsi1_dlength.sum(0).sum(0) raise ValueError, "unknown distriubtion received for psi-statistics"
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
def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
# Spike-and-Slab GPLVM # Spike-and-Slab GPLVM
@ -116,22 +118,25 @@ class RBF(Stationary):
return grad return grad
l2 = self.lengthscale **2 elif isinstance(variational_posterior, variational.NormalPosterior):
#psi1 l2 = self.lengthscale **2
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)
#psi2 #psi1
denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior)
term1 = Zdist / l2 # M, M, Q grad = np.einsum('ij,ij,ijk,ijk->jk', dL_dpsi1, psi1, dist, -1./(denom*l2))
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)
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): def gradients_qX_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
# Spike-and-Slab GPLVM # Spike-and-Slab GPLVM
@ -152,17 +157,23 @@ class RBF(Stationary):
return grad_mu, grad_S, grad_gamma return grad_mu, grad_S, grad_gamma
l2 = self.lengthscale **2 elif isinstance(variational_posterior, variational.NormalPosterior):
#psi1
denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior) l2 = self.lengthscale **2
tmp = psi1[:, :, None] / l2 / denom #psi1
grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1) denom, dist, dist_sq, psi1 = self._psi1computations(Z, variational_posterior)
grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1) tmp = psi1[:, :, None] / l2 / denom
#psi2 grad_mu = np.sum(dL_dpsi1[:, :, None] * tmp * dist, 1)
denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) grad_S = np.sum(dL_dpsi1[:, :, None] * 0.5 * tmp * (dist_sq - 1), 1)
tmp = psi2[:, :, :, None] / l2 / denom #psi2
grad_mu += -2.*(dL_dpsi2[:, :, :, None] * tmp * mudist).sum(1).sum(1) _, _, mudist, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior)
grad_S += (dL_dpsi2[:, :, :, None] * tmp * (2.*mudist_sq - 1)).sum(1).sum(1) 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 return grad_mu, grad_S
@ -170,61 +181,6 @@ class RBF(Stationary):
# Precomputations # # 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; q<input_dim; q++){
tmp = 0;
for(i=0; i<num_data; i++){
for(j=0; j<i; j++){
tmp += (X(i,q)-X(j,q))*(X(i,q)-X(j,q))*dvardLdK(i,j);
}
}
target(q) += var_len3(q)*tmp;
}
"""
num_data, num_inducing, input_dim = X.shape[0], X.shape[0], self.input_dim
X, dvardLdK, var_len3 = param_to_array(X, dvardLdK, var_len3)
weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
else:
code = """
int q,i,j;
double tmp;
for(q=0; q<input_dim; q++){
tmp = 0;
for(i=0; i<num_data; i++){
for(j=0; j<num_inducing; j++){
tmp += (X(i,q)-X2(j,q))*(X(i,q)-X2(j,q))*dvardLdK(i,j);
}
}
target(q) += var_len3(q)*tmp;
}
"""
num_data, num_inducing, input_dim = X.shape[0], X2.shape[0], self.input_dim
X, X2, dvardLdK, var_len3 = param_to_array(X, X2, dvardLdK, var_len3)
weave.inline(code, arg_names=['num_data', 'num_inducing', 'input_dim', 'X', 'X2', 'target', 'dvardLdK', 'var_len3'], type_converters=weave.converters.blitz, **self.weave_options)
return target
@Cache_this(limit=1) @Cache_this(limit=1)
def _psi1computations(self, Z, vp): def _psi1computations(self, Z, vp):
mu, S = vp.mean, vp.variance mu, S = vp.mean, vp.variance
@ -237,7 +193,7 @@ class RBF(Stationary):
return denom, dist, dist_sq, psi1 return denom, dist, dist_sq, psi1
#@cache_this(ignore_args=(1,)) @Cache_this(limit=1, ignore_args=(0,))
def _Z_distances(self, Z): def _Z_distances(self, Z):
Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q Zhat = 0.5 * (Z[:, None, :] + Z[None, :, :]) # M,M,Q
Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q Zdist = 0.5 * (Z[:, None, :] - Z[None, :, :]) # M,M,Q
@ -309,8 +265,4 @@ class RBF(Stationary):
arg_names=['N', 'M', 'Q', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'denom_l2', 'Zdist_sq', 'half_log_denom', 'psi2', 'variance_sq'], arg_names=['N', 'M', 'Q', 'mu', 'Zhat', 'mudist_sq', 'mudist', 'denom_l2', 'Zdist_sq', 'half_log_denom', 'psi2', 'variance_sq'],
type_converters=weave.converters.blitz, **self.weave_options) type_converters=weave.converters.blitz, **self.weave_options)
return denom, Zdist, Zdist_sq, mudist, mudist_sq, psi2 return Zdist, Zdist_sq, mudist, mudist_sq, psi2
def input_sensitivity(self):
if self.ARD: return 1./self.lengthscale
else: return (1./self.lengthscale).repeat(self.input_dim)

93
GPy/kern/_src/stationary.py
View file

@ -12,6 +12,35 @@ from scipy import integrate
from ...util.caching import Cache_this from ...util.caching import Cache_this
class Stationary(Kern): class Stationary(Kern):
"""
Stationary kernels (covariance functions).
Stationary covariance fucntion depend only on r, where r is defined as
r = \sqrt{ \sum_{q=1}^Q (x_q - x'_q)^2 }
The covariance function k(x, x' can then be written k(r).
In this implementation, r is scaled by the lengthscales parameter(s):
r = \sqrt{ \sum_{q=1}^Q \frac{(x_q - x'_q)^2}{\ell_q^2} }.
By default, there's only one lengthscale: seaprate lengthscales for each
dimension can be enables by setting ARD=True.
To implement a stationary covariance function using this class, one need
only define the covariance function k(r), and it derivative.
...
def K_of_r(self, r):
return foo
def dK_dr(self, r):
return bar
The lengthscale(s) and variance parameters are added to the structure automatically.
"""
def __init__(self, input_dim, variance, lengthscale, ARD, name): def __init__(self, input_dim, variance, lengthscale, ARD, name):
super(Stationary, self).__init__(input_dim, name) super(Stationary, self).__init__(input_dim, name)
self.ARD = ARD self.ARD = ARD
@ -20,11 +49,11 @@ class Stationary(Kern):
lengthscale = np.ones(1) lengthscale = np.ones(1)
else: else:
lengthscale = np.asarray(lengthscale) lengthscale = np.asarray(lengthscale)
assert lengthscale.size == 1, "Only lengthscale needed for non-ARD kernel" assert lengthscale.size == 1, "Only 1 lengthscale needed for non-ARD kernel"
else: else:
if lengthscale is not None: if lengthscale is not None:
lengthscale = np.asarray(lengthscale) lengthscale = np.asarray(lengthscale)
assert lengthscale.size in [1, input_dim], "Bad lengthscales" assert lengthscale.size in [1, input_dim], "Bad number of lengthscales"
if lengthscale.size != input_dim: if lengthscale.size != input_dim:
lengthscale = np.ones(input_dim)*lengthscale lengthscale = np.ones(input_dim)*lengthscale
else: else:
@ -35,26 +64,25 @@ class Stationary(Kern):
self.add_parameters(self.variance, self.lengthscale) self.add_parameters(self.variance, self.lengthscale)
def K_of_r(self, r): def K_of_r(self, r):
raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" raise NotImplementedError, "implement the covariance function as a fn of r to use this class"
def dK_dr(self, r): def dK_dr(self, r):
raise NotImplementedError, "implement the covaraiance function as a fn of r to use this class" raise NotImplementedError, "implement derivative of the covariance function wrt r to use this class"
#@Cache_this(limit=5, ignore_args=()) @Cache_this(limit=5, ignore_args=())
def K(self, X, X2=None): def K(self, X, X2=None):
r = self._scaled_dist(X, X2) r = self._scaled_dist(X, X2)
return self.K_of_r(r) return self.K_of_r(r)
#@Cache_this(limit=5, ignore_args=(0,)) @Cache_this(limit=3, ignore_args=())
def _dist(self, X, X2): def dK_dr_via_X(self, X, X2):
if X2 is None: #a convenience function, so we can cache dK_dr
X2 = X return self.dK_dr(self._scaled_dist(X, X2))
return X[:, None, :] - X2[None, :, :]
#@Cache_this(limit=5, ignore_args=(0,)) @Cache_this(limit=5, ignore_args=(0,))
def _unscaled_dist(self, X, X2=None): def _unscaled_dist(self, X, X2=None):
""" """
Compute the square distance between each row of X and X2, or between Compute the Euclidean distance between each row of X and X2, or between
each pair of rows of X if X2 is None. each pair of rows of X if X2 is None.
""" """
if X2 is None: if X2 is None:
@ -65,12 +93,12 @@ class Stationary(Kern):
X2sq = np.sum(np.square(X2),1) X2sq = np.sum(np.square(X2),1)
return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:])) return np.sqrt(-2.*np.dot(X, X2.T) + (X1sq[:,None] + X2sq[None,:]))
#@Cache_this(limit=5, ignore_args=()) @Cache_this(limit=5, ignore_args=())
def _scaled_dist(self, X, X2=None): def _scaled_dist(self, X, X2=None):
""" """
Efficiently compute the scaled distance, r. Efficiently compute the scaled distance, r.
r = \sum_{q=1}^Q (x_q - x'q)^2/l_q^2 r = \sqrt( \sum_{q=1}^Q (x_q - x'q)^2/l_q^2 )
Note that if thre is only one lengthscale, l comes outside the sum. In Note that if thre is only one lengthscale, l comes outside the sum. In
this case we compute the unscaled distance first (in a separate this case we compute the unscaled distance first (in a separate
@ -84,7 +112,6 @@ class Stationary(Kern):
else: else:
return self._unscaled_dist(X, X2)/self.lengthscale return self._unscaled_dist(X, X2)/self.lengthscale
def Kdiag(self, X): def Kdiag(self, X):
ret = np.empty(X.shape[0]) ret = np.empty(X.shape[0])
ret[:] = self.variance ret[:] = self.variance
@ -95,20 +122,23 @@ class Stationary(Kern):
self.lengthscale.gradient = 0. self.lengthscale.gradient = 0.
def update_gradients_full(self, dL_dK, X, X2=None): def update_gradients_full(self, dL_dK, X, X2=None):
r = self._scaled_dist(X, X2)
K = self.K_of_r(r)
rinv = self._inv_dist(X, X2) self.variance.gradient = np.einsum('ij,ij,i', self.K(X, X2), dL_dK, 1./self.variance)
dL_dr = self.dK_dr(r) * dL_dK
#now the lengthscale gradient(s)
dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
if self.ARD: if self.ARD:
x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 #rinv = self._inv_dis# this is rather high memory? Should we loop instead?t(X, X2)
self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0) #d = X[:, None, :] - X2[None, :, :]
#x_xl3 = np.square(d)
#self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3
tmp = dL_dr*self._inv_dist(X, X2)
if X2 is None: X2 = X
self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(X[:,q:q+1] - X2[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)])
else: else:
x_xl3 = np.square(self._dist(X, X2)) / self.lengthscale**3 r = self._scaled_dist(X, X2)
self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum() self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale
self.variance.gradient = np.sum(K * dL_dK)/self.variance
def _inv_dist(self, X, X2=None): def _inv_dist(self, X, X2=None):
""" """
@ -116,7 +146,7 @@ class Stationary(Kern):
diagonal, where we return zero (the distance on the diagonal is zero). diagonal, where we return zero (the distance on the diagonal is zero).
This term appears in derviatives. This term appears in derviatives.
""" """
dist = self._scaled_dist(X, X2) dist = self._scaled_dist(X, X2).copy()
if X2 is None: if X2 is None:
nondiag = util.diag.offdiag_view(dist) nondiag = util.diag.offdiag_view(dist)
nondiag[:] = 1./nondiag nondiag[:] = 1./nondiag
@ -128,10 +158,11 @@ class Stationary(Kern):
""" """
Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X Given the derivative of the objective wrt K (dL_dK), compute the derivative wrt X
""" """
r = self._scaled_dist(X, X2)
invdist = self._inv_dist(X, X2) invdist = self._inv_dist(X, X2)
dL_dr = self.dK_dr(r) * dL_dK dL_dr = self.dK_dr_via_X(X, X2) * dL_dK
#The high-memory numpy way: ret = np.sum((invdist*dL_dr)[:,:,None]*self._dist(X, X2),1)/self.lengthscale**2 #The high-memory numpy way:
#d = X[:, None, :] - X2[None, :, :]
#ret = np.sum((invdist*dL_dr)[:,:,None]*d,1)/self.lengthscale**2
#if X2 is None: #if X2 is None:
#ret *= 2. #ret *= 2.
@ -141,7 +172,7 @@ class Stationary(Kern):
tmp *= 2. tmp *= 2.
X2 = X X2 = X
ret = np.empty(X.shape, dtype=np.float64) ret = np.empty(X.shape, dtype=np.float64)
[np.copyto(ret[:,q], np.sum(tmp*(X[:,q][:,None]-X2[:,q][None,:]), 1)) for q in xrange(self.input_dim)] [np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) for q in xrange(self.input_dim)]
ret /= self.lengthscale**2 ret /= self.lengthscale**2
return ret return ret
@ -214,7 +245,7 @@ class Matern52(Stationary):
.. math:: .. 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'): def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='Mat52'):
super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, name) super(Matern52, self).__init__(input_dim, variance, lengthscale, ARD, name)
@ -225,7 +256,7 @@ class Matern52(Stationary):
def dK_dr(self, r): 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) 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. 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.

4
GPy/models/sparse_gp_regression.py
View file

@ -8,7 +8,7 @@ from .. import likelihoods
from .. import kern from .. import kern
from ..inference.latent_function_inference import VarDTC from ..inference.latent_function_inference import VarDTC
from ..util.misc import param_to_array from ..util.misc import param_to_array
from ..core.parameterization.variational import VariationalPosterior from ..core.parameterization.variational import NormalPosterior
class SparseGPRegression(SparseGP): class SparseGPRegression(SparseGP):
""" """
@ -47,7 +47,7 @@ class SparseGPRegression(SparseGP):
likelihood = likelihoods.Gaussian() likelihood = likelihoods.Gaussian()
if not (X_variance is None): 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()) SparseGP.__init__(self, X, Y, Z, kernel, likelihood, inference_method=VarDTC())

17
GPy/plotting/matplot_dep/models_plots.py
View file

@ -57,9 +57,12 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
fig = pb.figure(num=fignum) fig = pb.figure(num=fignum)
ax = fig.add_subplot(111) 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():
if hasattr(model, 'has_uncertain_inputs') and model.has_uncertain_inputs(): X_variance = model.X_variance 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) if hasattr(model, 'Z'): Z = param_to_array(model.Z)
#work out what the inputs are for plotting (1D or 2D) #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) #add error bars for uncertain (if input uncertainty is being modelled)
#if hasattr(model,"has_uncertain_inputs") and model.has_uncertain_inputs(): 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(), 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()), xerr=2 * np.sqrt(X_variance[which_data_rows, free_dims].flatten()),
# ecolor='k', fmt=None, elinewidth=.5, alpha=.5) ecolor='k', fmt=None, elinewidth=.5, alpha=.5)
#set the limits of the plot to some sensible values #set the limits of the plot to some sensible values

133
GPy/testing/observable_tests.py Normal file
View file

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

3
GPy/testing/parameterized_tests.py
View file

@ -6,6 +6,7 @@ Created on Feb 13, 2014
import unittest import unittest
import GPy import GPy
import numpy as np import numpy as np
from GPy.core.parameterization.parameter_core import HierarchyError
class Test(unittest.TestCase): class Test(unittest.TestCase):
@ -65,7 +66,7 @@ class Test(unittest.TestCase):
self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1]) self.assertListEqual(self.test1.constraints[Logexp()].tolist(), [0,1])
def test_add_parameter_already_in_hirarchy(self): 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): def test_default_constraints(self):
self.assertIs(self.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops) self.assertIs(self.rbf.variance.constraints._param_index_ops, self.rbf.constraints._param_index_ops)

6
GPy/util/caching.py
View file

@ -1,4 +1,5 @@
from ..core.parameterization.parameter_core import Observable from ..core.parameterization.parameter_core import Observable
import itertools
class Cacher(object): class Cacher(object):
""" """
@ -38,8 +39,11 @@ class Cacher(object):
if not all([isinstance(arg, Observable) for arg in observable_args]): if not all([isinstance(arg, Observable) for arg in observable_args]):
return self.operation(*args) return self.operation(*args)
# TODO: WARNING !!! Cache OFFSWITCH !!! WARNING
# return self.operation(*args)
#if the result is cached, return the cached computation #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): if any(state):
i = state.index(True) i = state.index(True)
if self.inputs_changed[i]: if self.inputs_changed[i]: