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

changing all parameterized objects to be compatible with the new parameterization

Browse source
This commit is contained in:
Max Zwiessele 2013-10-25 15:29:04 +01:00
parent e1bee4536a
commit d3721b76a8
21 changed files with 645 additions and 529 deletions
Download patch file Download diff file Expand all files Collapse all files

34
GPy/core/gp.py
View file

@ -23,6 +23,7 @@ class GP(GPBase):
"""
def __init__(self, X, likelihood, kernel, normalize_X=False):
GPBase.__init__(self, X, likelihood, kernel, normalize_X=normalize_X)
#self._set_params(self._get_params())
def getstate(self):
@ -89,21 +90,26 @@ class GP(GPBase):
return (-0.5 * self.num_data * self.output_dim * np.log(2.*np.pi) -
0.5 * self.output_dim * self.K_logdet + self._model_fit_term() + self.likelihood.Z)
# def _log_likelihood_gradients(self):
# """
# The gradient of all parameters.
#
# Note, we use the chain rule: dL_dtheta = dL_dK * d_K_dtheta
# """
# #return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
#
# if not isinstance(self.likelihood,EP):
# tmp = np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
# else:
# tmp = np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
# return tmp
def dL_dtheta(self):
return self.kern.dK_dtheta(self.dL_dK, self.X)
def _log_likelihood_gradients(self):
"""
The gradient of all parameters.
Note, we use the chain rule: dL_dtheta = dL_dK * d_K_dtheta
"""
#return np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
if not isinstance(self.likelihood,EP):
tmp = np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
else:
tmp = np.hstack((self.kern.dK_dtheta(dL_dK=self.dL_dK, X=self.X), self.likelihood._gradients(partial=np.diag(self.dL_dK))))
return tmp
def dL_dlikelihood(self):
return self.likelihood._gradients(partial=np.diag(self.dL_dK))
def _raw_predict(self, _Xnew, which_parts='all', full_cov=False, stop=False):
"""
Internal helper function for making predictions, does not account

9
GPy/core/gp_base.py
View file

@ -31,9 +31,9 @@ class GPBase(Model):
self._Xoffset = np.zeros((1, self.input_dim))
self._Xscale = np.ones((1, self.input_dim))
self.set_as_parameters(*self.kern._parameters_)
self.set_as_parameters(*self.likelihood._parameters_)
self.add_parameter(self.kern, gradient=self.dL_dtheta)
self.add_parameter(self.likelihood, gradient=self.dL_dlikelihood)
# Model.__init__(self)
# All leaf nodes should call self._set_params(self._get_params()) at
# the end
@ -53,7 +53,8 @@ class GPBase(Model):
self.likelihood,
self.output_dim,
self._Xoffset,
self._Xscale]
self._Xscale,
]
def setstate(self, state):
self._Xscale = state.pop()

1
GPy/core/index_operations.py
View file

@ -34,7 +34,6 @@ class ParamDict(defaultdict):
for a in self.iterkeys():
if numpy.all(a==key) and a._parent_index_==key._parent_index_:
return super(ParamDict, self).__setitem__(a, value)
raise KeyError, key
defaultdict.__setitem__(self, key, value)

46
GPy/core/model.py
View file

@ -12,6 +12,7 @@ import numpy as np
from domains import _POSITIVE, _REAL
from numpy.linalg.linalg import LinAlgError
from index_operations import ParameterIndexOperations
import itertools
# import numdifftools as ndt
class Model(Parameterized):
@ -28,8 +29,15 @@ class Model(Parameterized):
def log_likelihood(self):
raise NotImplementedError, "this needs to be implemented to use the model class"
def _log_likelihood_gradients(self):
#def dK_d(self, param, dL_dK, X, X2)
g = np.zeros(self.size)
try:
[g.__setitem__(s, self.gradient_mapping[p]().flat) for p, s in itertools.izip(self._parameters_, self._param_slices_) if not p.is_fixed]
except KeyError:
raise KeyError, 'Gradient for {} not defined, please specify gradients for parameters to optimize'.format(p.name)
return g
raise NotImplementedError, "this needs to be implemented to use the model class"
def getstate(self):
"""
Get the current state of the class.
@ -153,19 +161,19 @@ class Model(Parameterized):
return ret
def _transform_gradients(self, g):
x = self._get_params()
for constraint, index in self._constraints_.iteritems():
if constraint != __fixed__:
g[index] = g[index] * constraint.gradfactor(x[index])
#[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
[np.put(g, i, v) for i, v in [[i, t.sum()] for p in self._parameters_ for t,i in p._tied_to_me_.iteritems()]]
# if len(self.tied_indices) or len(self.fixed_indices):
# to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
# return np.delete(g, to_remove)
# else:
if self._fixes_ is not None: return g[self._fixes_]
return g
# def _transform_gradients(self, g):
# x = self._get_params()
# for constraint, index in self.constraints.iteritems():
# if constraint != __fixed__:
# g[index] = g[index] * constraint.gradfactor(x[index])
# #[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
# [np.put(g, i, v) for i, v in [[i, t.sum()] for p in self.flattened_parameters for t,i in p._tied_to_me_.iteritems()]]
# # if len(self.tied_indices) or len(self.fixed_indices):
# # to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
# # return np.delete(g, to_remove)
# # else:
# if self._fixes_ is not None: return g[self._fixes_]
# return g
def randomize(self):
"""
@ -173,15 +181,15 @@ class Model(Parameterized):
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(x.size)
#x = self._get_params_transformed()
x = np.random.randn(self.size_transformed)
self._set_params_transformed(x)
# now draw from prior where possible
x = self._get_params()
if self.priors is not None:
[np.put(x, i, p.rvs(1)) for i, p in enumerate(self.priors) 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...)
#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):
@ -480,7 +488,7 @@ class Model(Parameterized):
names = self._get_param_names_transformed()
except NotImplementedError:
names = ['Variable %i' % i for i in range(len(x))]
import ipdb;ipdb.set_trace()
# Prepare for pretty-printing
header = ['Name', 'Ratio', 'Difference', 'Analytical', 'Numerical']
max_names = max([len(names[i]) for i in range(len(names))] + [len(header[0])])
@ -522,7 +530,7 @@ class Model(Parameterized):
d = '%.6f' % float(difference)
g = '%.6f' % gradient
ng = '%.6f' % float(numerical_gradient)
grad_string = "{0:^{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4])
grad_string = "{0:<{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name, r, d, g, ng, c0=cols[0] + 9, c1=cols[1], c2=cols[2], c3=cols[3], c4=cols[4])
print grad_string
def input_sensitivity(self):

212
GPy/core/parameter.py
View file

@ -6,7 +6,7 @@ Created on 4 Sep 2013
import itertools
import numpy
from transformations import Logexp, NegativeLogexp, Logistic
from GPy.util.misc import fast_array_equal
from parameterized import Parentable
###### printing
__constraints_name__ = "Constraint"
@ -14,9 +14,9 @@ __index_name__ = "Index"
__tie_name__ = "Tied to"
__precision__ = numpy.get_printoptions()['precision'] # numpy printing precision used, sublassing numpy ndarray after all
__print_threshold__ = 5
######
######
class Param(numpy.ndarray):
class Param(numpy.ndarray, Parentable):
"""
Parameter object for GPy models.
@ -40,13 +40,12 @@ class Param(numpy.ndarray):
See :py:class:`GPy.core.parameterized.Parameterized` for more details.
"""
__array_priority__ = -numpy.inf # Never give back Param
def __new__(cls, name, input_array, gradient, *args, **kwargs):
def __new__(cls, name, input_array, *args, **kwargs):
obj = numpy.atleast_1d(numpy.array(input_array)).view(cls)
obj._name_ = name
obj._parent_ = None
obj._direct_parent_ = None
obj._parent_index_ = None
obj._gradient_parent_ = None
obj._gradient_ = gradient
obj._highest_parent_ = None
obj._current_slice_ = (slice(obj.shape[0]),)
obj._realshape_ = obj.shape
obj._realsize_ = obj.size
@ -62,9 +61,9 @@ class Param(numpy.ndarray):
if obj is None: return
self._name_ = getattr(obj, '_name_', None)
self._current_slice_ = getattr(obj, '_current_slice_', None)
self._parent_ = getattr(obj, '_parent_', None)
self._direct_parent_ = getattr(obj, '_direct_parent_', None)
self._parent_index_ = getattr(obj, '_parent_index_', None)
self._gradient_ = getattr(obj, '_gradient_', None)
self._highest_parent_ = getattr(obj, '_highest_parent_', None)
self._tied_to_me_ = getattr(obj, '_tied_to_me_', None)
self._tied_to_ = getattr(obj, '_tied_to_', None)
self._realshape_ = getattr(obj, '_realshape_', None)
@ -72,7 +71,11 @@ class Param(numpy.ndarray):
self._realndim_ = getattr(obj, '_realndim_', None)
self._updated_ = getattr(obj, '_updated_', None)
self._original_ = getattr(obj, '_original_', None)
self._gradient_parent_ = getattr(obj, '_gradient_parent_', None)
def __eq__(self, other):
return other is self
if other is self:
return True
return super(Param, self).__eq__(other)
def __array_wrap__(self, out_arr, context=None):
return out_arr.view(numpy.ndarray)
@ -83,10 +86,9 @@ class Param(numpy.ndarray):
func, args, state = super(Param, self).__reduce__()
return func, args, (state,
(self._name_,
self._parent_,
self._direct_parent_,
self._parent_index_,
self._gradient_parent_,
self._gradient_,
self._highest_parent_,
self._current_slice_,
self._realshape_,
self._realsize_,
@ -106,16 +108,16 @@ class Param(numpy.ndarray):
self._realsize_ = state.pop()
self._realshape_ = state.pop()
self._current_slice_ = state.pop()
self._gradient_ = state.pop()
self._gradient_parent_ = state.pop()
self._highest_parent_ = state.pop()
self._parent_index_ = state.pop()
self._parent_ = state.pop()
self._direct_parent_ = state.pop()
self._name_ = state.pop()
#===========================================================================
# get/set parameters
#===========================================================================
def _set_params(self, param):
self.flat = param
self._fire_changed()
def _get_params(self):
return self.flat
@property
@ -132,7 +134,10 @@ class Param(numpy.ndarray):
def name(self, new_name):
from_name = self.name
self._name_ = new_name
self._parent_._name_changed(self, from_name)
self._direct_parent_._name_changed(self, from_name)
@property
def _parameters_(self):
return []
#===========================================================================
# Fixing Parameters:
#===========================================================================
@ -142,21 +147,15 @@ class Param(numpy.ndarray):
:param warning: print a warning for overwriting constraints.
"""
self._parent_._fix(self,warning)
self._highest_parent_._fix(self,warning)
fix = constrain_fixed
def unconstrain_fixed(self):
"""
This parameter will no longer be fixed.
"""
self._parent_._unfix(self)
self._highest_parent_._unfix(self)
unfix = unconstrain_fixed
#===========================================================================
# Convenience methods:
#===========================================================================
@property
def is_fixed(self):
return self._parent_._is_fixed(self)
#===========================================================================
# Constrain operations -> done
#===========================================================================
def constrain(self, transform, warning=True, update=True):
@ -171,10 +170,10 @@ class Param(numpy.ndarray):
if self._original_: # this happens when indexing created a copy of the array
self.__setitem__(slice(None), transform.initialize(self), update=False)
else:
self._parent_._get_original(self).__setitem__(self._current_slice_, transform.initialize(self), update=False)
self._parent_._add_constrain(self, transform, warning)
self._direct_parent_._get_original(self).__setitem__(self._current_slice_, transform.initialize(self), update=False)
self._highest_parent_._add_constrain(self, transform, warning)
if update:
self._parent_.parameters_changed()
self._highest_parent_.parameters_changed()
def constrain_positive(self, warning=True):
"""
:param warning: print a warning if re-constraining parameters.
@ -204,7 +203,7 @@ class Param(numpy.ndarray):
remove all :py:class:`GPy.core.transformations.Transformation`
transformats of this parameter object.
"""
self._parent_._remove_constrain(self, *transforms)
self._highest_parent_._remove_constrain(self, *transforms)
def unconstrain_positive(self):
"""
Remove positive constraint of this parameter.
@ -239,23 +238,22 @@ class Param(numpy.ndarray):
if self._original_: # this happens when indexing created a copy of the array
self[:] = param
else:
self._parent_._get_original(self)[self._current_slice_] = param
self._direct_parent_._get_original(self)[self._current_slice_] = param
except ValueError:
raise ValueError("Trying to tie {} with shape {} to {} with shape {}".format(self.name, self.shape, param.name, param.shape))
self._parent_._get_original(self)._tied_to_ += [param]
self._direct_parent_._get_original(self)._tied_to_ += [param]
param._add_tie_listener(self)
self._parent_._set_fixed(self)
# self._parent_._add_tie(self, param)
self._highest_parent_._set_fixed(self)
# self._direct_parent_._add_tie(self, param)
def untie(self, *ties):
"""
remove tie of this parameter to ties it was tied to.
"""
[self._parent_._get_original(t)._remove_tie_listener(self) for t in self._tied_to_]
self._tied_to_ = [tied_to for tied_to in self._tied_to_ for t in tied_to._tied_to_me_ if self._parent_index_==self._parent_._get_original(t)._parent_index_]
self._parent_._set_unfixed(self)
# self._parent_._remove_tie(self, *params)
[t._direct_parent_._get_original(t)._remove_tie_listener(self) for t in self._tied_to_]
self._tied_to_ = [tied_to for tied_to in self._tied_to_ for t in tied_to._tied_to_me_ if self._parent_index_==t._direct_parent_._get_original(t)._parent_index_]
self._highest_parent_._set_unfixed(self)
# self._direct_parent_._remove_tie(self, *params)
def _fire_changed(self):
for tied, ind in self._tied_to_me_.iteritems():
tied._on_change(self.base, list(ind))
@ -266,7 +264,7 @@ class Param(numpy.ndarray):
if t._parent_index_ == to_remove._parent_index_:
new_index = list(set(t._raveled_index()) - set(to_remove._raveled_index()))
if new_index:
tmp = self._parent_._get_original(t)[numpy.unravel_index(new_index,t._realshape_)]
tmp = t._direct_parent_._get_original(t)[numpy.unravel_index(new_index,t._realshape_)]
self._tied_to_me_[tmp] = self._tied_to_me_[t]
del self._tied_to_me_[t]
if len(self._tied_to_me_[tmp]) == 0:
@ -279,7 +277,7 @@ class Param(numpy.ndarray):
if self._original_:
self.__setitem__(slice(None), val[ind], update=False)
else: # this happens when indexing created a copy of the array
self._parent_._get_original(self).__setitem__(self._current_slice_, val[ind], update=False)
self._direct_parent_._get_original(self).__setitem__(self._current_slice_, val[ind], update=False)
self._fire_changed()
self._updated_ = False
#===========================================================================
@ -291,16 +289,16 @@ class Param(numpy.ndarray):
Set prior for this parameter.
"""
if not hasattr(self._parent_, '_set_prior'):
raise AttributeError("Parent of type {} does not support priors".format(self._parent_.__class__))
self._parent_._set_prior(self, prior)
if not hasattr(self._highest_parent_, '_set_prior'):
raise AttributeError("Parent of type {} does not support priors".format(self._highest_parent_.__class__))
self._highest_parent_._set_prior(self, prior)
def unset_prior(self, *priors):
"""
:param priors: priors to remove from this parameter
Remove all priors from this parameter
"""
self._parent_._remove_prior(self, *priors)
self._highest_parent_._remove_prior(self, *priors)
#===========================================================================
# Array operations -> done
#===========================================================================
@ -321,7 +319,7 @@ class Param(numpy.ndarray):
numpy.ndarray.__setitem__(self, s, val)
self._fire_changed()
if update:
self._parent_.parameters_changed()
self._highest_parent_.parameters_changed()
#===========================================================================
# Index Operations:
#===========================================================================
@ -368,22 +366,48 @@ class Param(numpy.ndarray):
return numpy.r_[:b]
return itertools.imap(f, itertools.izip_longest(slice_index[:self._realndim_], self._realshape_, fillvalue=slice(self.size)))
#===========================================================================
# Printing -> done
# Convienience
#===========================================================================
@property
def _desc(self):
if self.size <= 1: return "%f"%self
else: return str(self.shape)
@property
def _constr(self):
return ' '.join(map(lambda c: str(c[0]) if c[1].size==self._realsize_ else "{"+str(c[0])+"}", self._parent_._constraints_iter_items(self)))
def is_fixed(self):
return self._highest_parent_._is_fixed(self)
def round(self, decimals=0, out=None):
view = super(Param, self).round(decimals, out).view(Param)
view.__array_finalize__(self)
return view
round.__doc__ = numpy.round.__doc__
def _get_original(self, param):
return self
#===========================================================================
# Printing -> done
#===========================================================================
@property
def _description_str(self):
if self.size <= 1: return ["%f"%self]
else: return [str(self.shape)]
def _parameter_names(self, add_name):
return [self.name]
@property
def flattened_parameters(self):
return [self]
@property
def parameter_shapes(self):
return [self.shape]
@property
def _constraints_str(self):
return [' '.join(map(lambda c: str(c[0]) if c[1].size==self._realsize_ else "{"+str(c[0])+"}", self._highest_parent_._constraints_iter_items(self)))]
@property
def _ties_str(self):
return [t._short() for t in self._tied_to_] or ['']
@property
def name_hirarchical(self):
if self.has_parent():
return self._direct_parent_.hirarchy_name()+self.name
return self.name
def __repr__(self, *args, **kwargs):
return "\033[1m{x:s}\033[0;0m:\n".format(x=self.name)+super(Param, self).__repr__(*args,**kwargs)
name = "\033[1m{x:s}\033[0;0m:\n".format(
x=self.name_hirarchical)
return name + super(Param, self).__repr__(*args,**kwargs)
def _ties_for(self, rav_index):
ties = numpy.empty(shape=(len(self._tied_to_), numpy.size(rav_index)), dtype=Param)
for i, tied_to in enumerate(self._tied_to_):
@ -395,7 +419,7 @@ class Param(numpy.ndarray):
#[ties.__setitem__(i, ties[i] + [tied_to]) for i in t._raveled_index()]
return map(lambda a: sum(a,[]), zip(*[[[tie.flatten()] if tx!=None else [] for tx in t] for t,tie in zip(ties,self._tied_to_)]))
def _constraints_for(self, rav_index):
return self._parent_._constraints_for(self, rav_index)
return self._highest_parent_._constraints_for(self, rav_index)
def _indices(self, slice_index=None):
# get a int-array containing all indices in the first axis.
if slice_index is None:
@ -412,17 +436,18 @@ class Param(numpy.ndarray):
def _max_len_names(self, gen, header):
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.name))
return reduce(lambda a, b:max(a, len("{x:=.{0}g}".format(__precision__, x=b))), self.flat, len(self.name))
def _max_len_index(self, ind):
return reduce(lambda a, b:max(a, len(str(b))), ind, len(__index_name__))
def _short(self, slice_index=None):
def _short(self):
# short string to print
name = self._direct_parent_.hirarchy_name() + self.name
if self._realsize_ < 2:
return self.name
ind = self._indices(slice_index)
return name
ind = self._indices()
if ind.size > 4: indstr = ','.join(map(str,ind[:2])) + "..." + ','.join(map(str,ind[-2:]))
else: indstr = ','.join(map(str,ind))
return self.name+'['+indstr+']'
return name+'['+indstr+']'
def __str__(self, constr_matrix=None, indices=None, ties=None, lc=None, lx=None, li=None, lt=None):
filter_ = self._current_slice_
vals = self.flat
@ -437,7 +462,7 @@ class Param(numpy.ndarray):
if lt is None: lt = self._max_len_names(ties, __tie_name__)
header = " {i:^{2}s} | \033[1m{x:^{1}s}\033[0;0m | {c:^{0}s} | {t:^{3}s}".format(lc,lx,li,lt, x=self.name, c=__constraints_name__, i=__index_name__, t=__tie_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} | {t:^{4}s} ".format(lc,lx,__precision__,li,lt, x=x, c=" ".join(map(str,c)), t=(t or ''), i=i) for i,x,c,t in itertools.izip(indices,vals,constr_matrix,ties)]) # return all the constraints with right indices
return "\n".join([header]+[" {i!s:^{3}s} | {x: >{1}.{2}g} | {c:^{0}s} | {t:^{4}s} ".format(lc,lx,__precision__,li,lt, x=x, c=" ".join(map(str,c)), t=(t or ''), i=i) for i,x,c,t in itertools.izip(indices,vals,constr_matrix,ties)]) # return all the constraints with right indices
#except: return super(Param, self).__str__()
class ParamConcatenation(object):
@ -449,7 +474,12 @@ class ParamConcatenation(object):
See :py:class:`GPy.core.parameter.Param` for more details on constraining.
"""
self.params = params
#self.params = params
self.params = []
for p in params:
for p in p.flattened_parameters:
if p not in self.params:
self.params.append(p)
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:])]
@ -458,7 +488,7 @@ class ParamConcatenation(object):
#===========================================================================
def __getitem__(self, s):
ind = numpy.zeros(sum(self._param_sizes), dtype=bool); ind[s] = True;
params = [p.flatten()[ind[ps]] for p,ps in zip(self.params, self._param_slices_) if numpy.any(p.flat[ind[ps]])]
params = [p._get_params()[ind[ps]] for p,ps in zip(self.params, self._param_slices_) if numpy.any(p._get_params()[ind[ps]])]
if len(params)==1: return params[0]
return ParamConcatenation(params)
def __setitem__(self, s, val, update=True):
@ -467,7 +497,7 @@ class ParamConcatenation(object):
[numpy.place(p, ind[ps], vals[ps]) and p._fire_changed()
for p, ps in zip(self.params, self._param_slices_)]
if update:
self.params[0]._parent_.parameters_changed()
self.params[0]._highest_parent_.parameters_changed()
def _vals(self):
return numpy.hstack([p._get_params() for p in self.params])
#===========================================================================
@ -505,6 +535,8 @@ class ParamConcatenation(object):
def unconstrain_bounded(self, lower, upper):
[param.unconstrain_bounded(lower, upper) for param in self.params]
unconstrain_bounded.__doc__ = Param.unconstrain_bounded.__doc__
def untie(self, *ties):
[param.untie(*ties) for param in self.params]
__lt__ = lambda self, val: self._vals()<val
__le__ = lambda self, val: self._vals()<=val
__eq__ = lambda self, val: self._vals()==val
@ -514,39 +546,49 @@ class ParamConcatenation(object):
def __str__(self, *args, **kwargs):
def f(p):
ind = p._raveled_index()
return p._constraints_for(ind), p._ties_for(ind)
constr_matrices, ties_matrices = zip(*map(f, self.params))
indices = [p._indices() for p in self.params]
lc = max([p._max_len_names(cm, __constraints_name__) for p, cm in itertools.izip(self.params, constr_matrices)])
lx = max([p._max_len_values() for p in self.params])
li = max([p._max_len_index(i) for p, i in itertools.izip(self.params, indices)])
lt = max([p._max_len_names(tm, __tie_name__) for p, tm in itertools.izip(self.params, ties_matrices)])
strings = [p.__str__(cm, i, tm, lc, lx, li, lt) for p, cm, i, tm in itertools.izip(self.params,constr_matrices,indices,ties_matrices)]
return "\n{}\n".format(" -"+"- | -".join(['-'*l for l in [li,lx,lc,lt]])).join(strings)
return p._constraints_for(ind), p._ties_for(ind)
params = self.params
constr_matrices, ties_matrices = zip(*map(f, params))
indices = [p._indices() for p in params]
lc = max([p._max_len_names(cm, __constraints_name__) for p, cm in itertools.izip(params, constr_matrices)])
lx = max([p._max_len_values() for p in params])
li = max([p._max_len_index(i) for p, i in itertools.izip(params, indices)])
lt = max([p._max_len_names(tm, __tie_name__) for p, tm in itertools.izip(params, ties_matrices)])
strings = [p.__str__(cm, i, tm, lc, lx, li, lt) for p, cm, i, tm in itertools.izip(params,constr_matrices,indices,ties_matrices)]
return "\n".join(strings)
return "\n{}\n".format(" -"+"- | -".join(['-'*l for l in [li,lx,lc,lt]])).join(strings)
def __repr__(self):
return "\n".join(map(repr,self.params))
if __name__ == '__main__':
from GPy.core.parameterized import Parameterized
from GPy.core.parameter import Param
#X = numpy.random.randn(2,3,1,5,2,4,3)
X = numpy.random.randn(1000,20)
X = numpy.random.randn(3,2)
print "random done"
p = Param("q_mean", X, None)
p1 = Param("q_variance", numpy.random.rand(*p.shape), None)
p2 = Param("Y", numpy.random.randn(p.shape[0],1), None)
p3 = Param("rbf_variance", numpy.random.rand(), None)
p4 = Param("rbf_lengthscale", numpy.random.rand(2), None)
p = Param("q_mean", X)
p1 = Param("q_variance", numpy.random.rand(*p.shape))
p2 = Param("Y", numpy.random.randn(p.shape[0],1))
p3 = Param("variance", numpy.random.rand())
p4 = Param("lengthscale", numpy.random.rand(2))
m = Parameterized()
rbf = Parameterized(name='rbf')
rbf.add_parameter(p3,p4)
m.add_parameter(p,p1,rbf)
print "setting params"
m.set_as_parameters(p,p1,p2,p3,p4)
#print m.q_v[3:5,[1,4,5]]
print "constraining variance"
m[".*variance"].constrain_positive()
print "constraining rbf"
m.rbf_l.constrain_positive()
m.q_variance[1,[0,5,11,19,2]].tie_to(m.rbf_v)
m.rbf_v.tie_to(m.rbf_l[0])
m.rbf_l[0].tie_to(m.rbf_l[1])
#m[".*variance"].constrain_positive()
#print "constraining rbf"
#m.rbf_l.constrain_positive()
#m.q_variance[1,[0,5,11,19,2]].tie_to(m.rbf_v)
#m.rbf_v.tie_to(m.rbf_l[0])
#m.rbf_l[0].tie_to(m.rbf_l[1])
#m.q_v.tie_to(m.rbf_v)
# m.rbf_l.tie_to(m.rbf_va)
# pt = numpy.array(params._get_params_transformed())

282
GPy/core/parameterized.py
View file

@ -5,14 +5,36 @@
import numpy; np = numpy
import copy
import cPickle
from parameter import ParamConcatenation, Param
from index_operations import ParameterIndexOperations,\
index_empty
import transformations
import itertools
from re import compile, _pattern_type
import re
from GPy.util.misc import fast_array_equal
class Parentable(object):
_direct_parent_ = None
_parent_index_ = None
def has_parent(self):
return self._direct_parent_ is not None
class Nameable(Parentable):
_name = None
def __init__(self, name):
self._name = name or self.__class__.__name__
self.name = name
@property
def name(self):
return self._name
@name.setter
def name(self, name):
from_name = self.name
self._name = name
if self.has_parent():
self._direct_parent_._name_changed(self, from_name)
from parameter import ParamConcatenation
from index_operations import ParameterIndexOperations,\
index_empty
#===============================================================================
# Printing:
@ -25,7 +47,7 @@ FIXED = False
UNFIXED = True
#===============================================================================
class Parameterized(object):
class Parameterized(Nameable):
"""
Parameterized class
@ -64,65 +86,85 @@ class Parameterized(object):
- m.name[0].tie_to(m.name[1])
Fixing parameters will fix them to the value they are right now. If you change
the paramters value, the parameter will be fixed to the new value!
the parameters value, the parameter will be fixed to the new value!
If you want to operate on all parameters use m[''] to wildcard select all paramters
and concatenate them. Printing m[''] will result in printing of all parameters in detail.
"""
def __init__(self):
self._constraints_ = ParameterIndexOperations()
#self._fixes_ = TieIndexOperations(self)
#self._fixes_ = None
self._fixes_ = None
def __init__(self, name=None):
super(Parameterized, self).__init__(name)
self._in_init_ = True
self._constraints_ = None#ParameterIndexOperations()
self._fixes_ = None
if not hasattr(self, "_parameters_"):
self._parameters_ = []
#else:
# self._parameters_.extend(parameters)
self._connect_parameters()
self.gradient_mapping = {}
del self._in_init_
@property
def constraints(self):
if self._constraints_ is None:
self._constraints_ = ParameterIndexOperations()
return self._constraints_
#===========================================================================
# Parameter connection for model creation:
#===========================================================================
def set_as_parameter(self, name, array, gradient, index=None, gradient_parent=None):
# def set_as_parameter(self, name, array, gradient, index=None, gradient_parent=None):
# """
# :param name: name of the parameter (in print and plots), can be callable without parameters
# :type name: str, callable
# :param array: array which the parameter consists of
# :type array: array-like
# :param gradient: gradient method of the parameter
# :type gradient: callable
# :param index: (optional) index of the parameter when printing
#
# (:param gradient_parent: connect these parameters to this class, but tell
# updates to highest_parent, this is needed when parameterized classes
# contain parameterized classes, but want to access the parameters
# of their children)
#
#
# Set array (e.g. self.X) as parameter with name and gradient.
# I.e: self.set_as_parameter('curvature', self.lengthscale, self.dK_dlengthscale)
#
# Note: the order in which parameters are added can be adjusted by
# giving an index, of where to put this parameter in printing
# """
# if index is None:
# self._parameters_.append(Param(name, array, gradient))
# else:
# self._parameters_.insert(index, Param(name, array, gradient))
# self._connect_parameters(gradient_parent=gradient_parent)
def add_parameter(self, parameter, gradient=None, index=None):
"""
:param name: name of the parameter (in print and plots), can be callable without parameters
:type name: str, callable
:param array: array which the parameter consists of
:type array: array-like
:param gradient: gradient method of the parameter
:type gradient: callable
:param index: (optional) index of the parameter when printing
(:param gradient_parent: connect these parameters to this class, but tell
updates to highest_parent, this is needed when parameterized classes
contain parameterized classes, but want to access the parameters
of their children)
Set array (e.g. self.X) as parameter with name and gradient.
I.e: self.set_as_parameter('curvature', self.lengthscale, self.dK_dlengthscale)
Note: the order in which parameters are added can be adjusted by
giving an index, of where to put this parameter in printing
"""
if index is None:
self._parameters_.append(Param(name, array, gradient))
else:
self._parameters_.insert(index, Param(name, array, gradient))
self._connect_parameters(gradient_parent=gradient_parent)
def set_as_parameters(self, *parameters, **kwargs):
"""
:param parameters: the parameters to add
:param index: index of where to put parameters
:param parameters: the parameters to add
:type parameters: list of or one :py:class:`GPy.core.parameter.Param`
:param [gradients]: gradients for each parameter,
one gradient per parameter
:param [index]: index of where to put parameters
Add all parameters to this parameter class, you can insert parameters
at any given index using the :py:func:`list.insert` syntax
"""
if kwargs.get('index',None) is None:
self._parameters_.extend(parameters)
if index is None:
self._parameters_.append(parameter)
else:
self._parameters_.insert(kwargs['index'], parameters)
self._connect_parameters(gradient_parent=kwargs.get('gradient_parent', self))
self._parameters_.insert(index, parameter)
self._connect_parameters()
if gradient:
self.gradient_mapping[parameter] = gradient
def add_parameters(self, *parameters):
"""
convinience method for adding several
parameters without gradient specification
"""
[self.add_parameter(p) for p in parameters]
# def remove_parameter(self, *names_params_indices):
# """
# :param names_params_indices: mix of parameter_names, parameter objects, or indices
@ -138,23 +180,21 @@ class Parameterized(object):
def parameters_changed(self):
# will be called as soon as paramters have changed
pass
def _connect_parameters(self, gradient_parent=None):
def _connect_parameters(self):
# connect parameterlist to this parameterized object
# This just sets up the right connection for the params objects
# to be used as parameters
if not hasattr(self, "_parameters_") or len(self._parameters_) < 1:
# no parameters for this class
return
sizes = numpy.cumsum([0] + self._parameter_sizes_)
self._parameter_size_ = sizes[-1]
self._param_slices_ = [slice(start, stop) for start,stop in zip(sizes, sizes[1:])]
i = 0
for p in self._parameters_:
#if p._parent_ is None:
p._parent_ = self
p._direct_parent_ = self
p._parent_index_ = i
i += 1
p._gradient_parent_ = gradient_parent or p._gradient_parent_ or self
for pi in p.flattened_parameters:
pi._highest_parent_ = self
not_unique = []
# for k,v in self.__dict__.iteritems():
# try:
@ -163,13 +203,15 @@ class Parameterized(object):
# except: # parameter comparison, just for convenience
# pass
if p.name in self.__dict__:
if p.base is self.__dict__[p.name] or p is self.__dict__[p.name]:
self.__dict__[p.name] = p
else:
if not p is self.__dict__[p.name]:
not_unique.append(p.name)
del self.__dict__[p.name]
elif not (p.name in not_unique):
self.__dict__[p.name] = p
sizes = numpy.cumsum([0] + self._parameter_sizes_)
self.size = sizes[-1]
self._param_slices_ = [slice(start, stop) for start,stop in zip(sizes, sizes[1:])]
self.parameters_changed()
#===========================================================================
# Pickling operations
#===========================================================================
@ -214,9 +256,13 @@ class Parameterized(object):
self._constraints_,
self._priors_,
self._parameters_,
self._name,
self.gradient_mapping,
]
def setstate(self, state):
self.gradient_mapping = state.pop(),
self._name = state.pop()
self._parameters_ = state.pop()
self._connect_parameters()
self._priors = state.pop()
@ -224,10 +270,36 @@ class Parameterized(object):
self._fixes_ = state.pop()
self.parameters_changed()
#===========================================================================
# Gradient control
#===========================================================================
def _transform_gradients(self, g):
if self.has_parent():
return g
x = self._get_params()
#g = g.copy()
#for constraint, index in self.constraints.iteritems():
# if constraint != __fixed__:
# g[index] = g[index] * constraint.gradfactor(x[index])
[numpy.put(g, i, g[i]*c.gradfactor(x[i])) for c,i in self.constraints.iteritems() if c != __fixed__]
#[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
for p in self.flattened_parameters:
for t,i in p._tied_to_me_.iteritems():
g[self._offset_for(p) + numpy.array(list(i))] += g[self._raveled_index_for(t)]
#[g[self._offset_for(t) + numpy.array(list(i))].__iadd__(v) for i, v in [[i, g[self._raveled_index_for(p)].sum()] for p in self.flattened_parameters for t,i in p._tied_to_me_.iteritems()]]
# if len(self.tied_indices) or len(self.fixed_indices):
# to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
# return np.delete(g, to_remove)
# else:
if self._fixes_ is not None: return g[self._fixes_]
return g
#===========================================================================
# Optimization handles:
#===========================================================================
def _get_param_names_transformed(self):
return numpy.array([p.name+str(i) for p in self._parameters_ for i in p._indices()])[self._fixes_][0]
n = numpy.array([p.name_hirarchical+'['+str(i)+']' for p in self.flattened_parameters for i in p._indices()])
if self._fixes_ is not None:
return n[self._fixes_]
return n
def _get_params(self):
# don't overwrite this anymore!
return numpy.hstack([x._get_params() for x in self._parameters_])#numpy.fromiter(itertools.chain(*itertools.imap(lambda x: x._get_params(), self._parameters_)), dtype=numpy.float64, count=sum(self._parameter_sizes_))
@ -237,14 +309,14 @@ class Parameterized(object):
self.parameters_changed()
def _get_params_transformed(self):
p = self._get_params()
[numpy.put(p, ind, c.finv(p[ind])) for c,ind in self._constraints_.iteritems() if c != __fixed__]
[numpy.put(p, ind, c.finv(p[ind])) for c,ind in self.constraints.iteritems() if c != __fixed__]
if self._fixes_ is not None:
return p[self._fixes_]
return p
def _set_params_transformed(self, p):
p = p.copy()
if self._fixes_ is not None: 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__]
[numpy.put(p, ind, c.f(p[ind])) for c,ind in self.constraints.iteritems() if c != __fixed__]
self._set_params(p)
def _name_changed(self, param, old_name):
if hasattr(self, old_name):
@ -262,14 +334,19 @@ class Parameterized(object):
return ind
def _offset_for(self, param):
# get the offset in the parameterized index array for param
return self._param_slices_[param._parent_index_].start
if param._direct_parent_._get_original(param) in self._parameters_:
return self._param_slices_[param._direct_parent_._get_original(param)._parent_index_].start
if param.has_parent():
return self._offset_for(param._direct_parent_) + param._direct_parent_._offset_for(param)
return 0
def _raveled_index_for(self, param):
return param._raveled_index() + self._offset_for(param)
#===========================================================================
# Handle ties:
#===========================================================================
def _set_fixed(self, param_or_index):
if self._fixes_ is None: self._fixes_ = numpy.ones(self._parameter_size_, dtype=bool)
if self._fixes_ is None: self._fixes_ = numpy.ones(self.size, dtype=bool)
try:
param_or_index = self._raveled_index_for(param_or_index)
except AttributeError:
@ -277,7 +354,7 @@ class Parameterized(object):
self._fixes_[param_or_index] = FIXED
if numpy.all(self._fixes_): self._fixes_ = None # ==UNFIXED
def _set_unfixed(self, param_or_index):
if self._fixes_ is None: self._fixes_ = numpy.ones(self._parameter_size_, dtype=bool)
if self._fixes_ is None: self._fixes_ = numpy.ones(self.size, dtype=bool)
try:
param_or_index = self._raveled_index_for(param_or_index)
except AttributeError:
@ -288,7 +365,7 @@ class Parameterized(object):
# # tie param to tie_to, if the values match (with broadcasting)
# self._remove_tie(param) # delete if multiple ties should be allowed
# f, _ = self._fixes_.add(param, tied_to)
# if self._fixes_ is None: self._fixes_ = numpy.ones(self._parameter_size_, dtype=bool)
# if self._fixes_ is None: self._fixes_ = numpy.ones(self.size, dtype=bool)
# self._fixes_[f] = False
# def _remove_tie(self, param, *params):
# # remove the tie from param to all *params (can be None, so all ties get deleted for param)
@ -329,11 +406,20 @@ class Parameterized(object):
if self._fixes_ is None:
return False
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 parameter through this method, by passing
# the copy here
return self._parameters_[param._parent_index_]
def hirarchy_name(self):
if self.has_parent():
return self._direct_parent_.hirarchy_name() + self.name + "."
return ''
#===========================================================================
# Constraint Handling:
#===========================================================================
@ -341,7 +427,7 @@ class Parameterized(object):
rav_i = self._raveled_index_for(param)
reconstrained = self._remove_constrain(param, index=rav_i) # remove constraints before
# if removing constraints before adding new is not wanted, just delete the above line!
self._constraints_.add(transform, rav_i)
self.constraints.add(transform, rav_i)
if warning and any(reconstrained):
# if you want to print the whole params object, which was reconstrained use:
# m = str(param[self._backtranslate_index(param, reconstrained)])
@ -349,19 +435,19 @@ class Parameterized(object):
return rav_i
def _remove_constrain(self, param, *transforms, **kwargs):
if not transforms:
transforms = self._constraints_.properties()
transforms = self.constraints.properties()
removed_indices = numpy.array([]).astype(int)
if "index" in kwargs: index = kwargs['index']
else: index = self._raveled_index_for(param)
for constr in transforms:
removed = self._constraints_.remove(constr, index)
removed = self.constraints.remove(constr, index)
if constr is __fixed__:
self._set_unfixed(removed)
removed_indices = numpy.union1d(removed_indices, removed)
return removed_indices
# convienience for iterating over items
def _constraints_iter_items(self, param):
for constr, ind in self._constraints_.iteritems():
for constr, ind in self.constraints.iteritems():
ind = self._backtranslate_index(param, ind)
if not index_empty(ind):
yield constr, ind
@ -372,9 +458,9 @@ class Parameterized(object):
for _, ind in self._constraints_iter_items(param):
yield ind
def _constraint_indices(self, param, constraint):
return self._backtranslate_index(param, self._constraints_[constraint])
return self._backtranslate_index(param, self.constraints[constraint])
def _constraints_for(self, param, rav_index):
return self._constraints_.properties_for(rav_index+self._offset_for(param))
return self.constraints.properties_for(rav_index+self._offset_for(param))
#===========================================================================
# Get/set parameters:
#===========================================================================
@ -383,6 +469,13 @@ class Parameterized(object):
create a list of parameters, matching regular expression regexp
"""
if not isinstance(regexp, _pattern_type): regexp = compile(regexp)
found_params = []
for p in self._parameters_:
if regexp.match(p.name) is not None:
found_params.append(p)
if isinstance(p, Parameterized):
found_params.extend(p.grep_param_names(regexp))
return found_params
return [param for param in self._parameters_ if regexp.match(param.name) is not None]
def __getitem__(self, name, paramlist=None):
if paramlist is None:
@ -394,8 +487,8 @@ class Parameterized(object):
try: param = self.__getitem__(name, paramlist)
except AttributeError as a: raise a
param[:] = value
def __getattr__(self, name):
return self.__getitem__(name)
# def __getattr__(self, name):
# return self.__getitem__(name)
# def __getattribute__(self, name):
# #try:
# return object.__getattribute__(self, name)
@ -412,41 +505,52 @@ class Parameterized(object):
if len(paramlist) == 1: self.__setitem__(name, val, paramlist); return
object.__setattr__(self, name, val);
#===========================================================================
# Printing:
# Printing:
#===========================================================================
def _parameter_names(self, add_name=False):
if add_name:
return [self.name + "." + xi for x in self._parameters_ for xi in x._parameter_names(add_name=True)]
return [xi for x in self._parameters_ for xi in x._parameter_names(add_name=True)]
parameter_names = property(_parameter_names, doc="Names for all parameters handled by this parameterization object")
@property
def parameter_names(self):
return [x.name for x in self._parameters_]
def flattened_parameters(self):
return [xi for x in self._parameters_ for xi in x.flattened_parameters]
@property
def _parameter_sizes_(self):
return [x.size for x in self._parameters_]
@property
def _parameter_size_transformed_(self):
return sum(self._fixes_)
def size_transformed(self):
if self._fixes_ is not None:
return sum(self._fixes_)
return self.size
@property
def _parameter_shapes_(self):
return [x.shape for x in self._parameters_]
def parameter_shapes(self):
return [xi for x in self._parameters_ for xi in x.parameter_shapes]
@property
def _constrs(self):
return [p._constr for p in self._parameters_]
def _constraints_str(self):
return [cs for p in self._parameters_ for cs in p._constraints_str]
@property
def _descs(self):
return [x._desc for x in self._parameters_]
def _description_str(self):
return [xi for x in self._parameters_ for xi in x._description_str]
@property
def _ts(self):
return [' '.join([t._short() for t in x._tied_to_]) for x in self._parameters_]
def _ties_str(self):
return [xi for x in self._parameters_ for xi in x._ties_str]
def __str__(self, header=True):
nl = max([len(str(x)) for x in self.parameter_names + ["Name"]])
sl = max([len(str(x)) for x in self._descs + ["Value"]])
constrs = self._constrs; ts = self._ts
constrs = self._constraints_str; ts = self._ties_str
desc = self._description_str; names = self.parameter_names
nl = max([len(str(x)) for x in names + ["Name"]])
sl = max([len(str(x)) for x in desc + ["Value"]])
cl = max([len(str(x)) if x else 0 for x in constrs + ["Constraint"]])
tl = max([len(str(x)) if x else 0 for x in ts + ["Tied to"]])
format_spec = " \033[1m{{p.name:^{0}s}}\033[0;0m | {{p._desc:^{1}s}} | {{const:^{2}s}} | {{t:^{3}s}}".format(nl, sl, cl, tl)
to_print = [format_spec.format(p=p, const=c, t=t) for p, c, t in itertools.izip(self._parameters_, constrs, ts)]
format_spec = " \033[1m{{name:<{0}s}}\033[0;0m | {{desc:^{1}s}} | {{const:^{2}s}} | {{t:^{3}s}}".format(nl, sl, cl, tl)
to_print = []
for n, d, c, t in itertools.izip(names, desc, constrs, ts):
to_print.append(format_spec.format(name=n, desc=d, const=c, t=t))
#to_print = [format_spec.format(p=p, const=c, t=t) if isinstance(p, Param) else p.__str__(header=False) for p, c, t in itertools.izip(self._parameters_, constrs, ts)]
sep = '-'*(nl+sl+cl+tl+8*2+3)
if header:
header = " {{0:^{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}}".format(nl, sl, cl, tl).format("Name", "Value", "Constraint", "Tied to")
header += '\n' + sep
header = " {{0:<{0}s}} | {{1:^{1}s}} | {{2:^{2}s}} | {{3:^{3}s}}".format(nl, sl, cl, tl).format("Name", "Value", "Constraint", "Tied to")
#header += '\n' + sep
to_print.insert(0, header)
return '\n'.format(sep).join(to_print)
pass

9
GPy/examples/regression.py
View file

@ -255,7 +255,7 @@ def toy_rbf_1d(optimizer='tnc', max_nb_eval_optim=100):
print(m)
return m
def toy_rbf_1d_50(max_iters=100):
def toy_rbf_1d_50(max_iters=100, optimize=True):
"""Run a simple demonstration of a standard Gaussian process fitting it to data sampled from an RBF covariance."""
data = GPy.util.datasets.toy_rbf_1d_50()
@ -263,14 +263,15 @@ def toy_rbf_1d_50(max_iters=100):
m = GPy.models.GPRegression(data['X'], data['Y'])
# optimize
m.optimize(max_iters=max_iters)
if optimize:
m.optimize(max_iters=max_iters)
# plot
m.plot()
print(m)
return m
def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4):
def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4, optimize=True):
# Create an artificial dataset where the values in the targets (Y)
# only depend in dimensions 1 and 3 of the inputs (X). Run ARD to
# see if this dependency can be recovered
@ -300,7 +301,7 @@ def toy_ARD(max_iters=1000, kernel_type='linear', num_samples=300, D=4):
# len_prior = GPy.priors.inverse_gamma(1,18) # 1, 25
# m.set_prior('.*lengthscale',len_prior)
m.optimize(optimizer='scg', max_iters=max_iters, messages=1)
if optimize: m.optimize(optimizer='scg', max_iters=max_iters, messages=1)
m.kern.plot_ARD()
print(m)

281
GPy/kern/kern.py
View file

@ -32,57 +32,44 @@ class kern(Parameterized):
:type input_slices: list of slice objects, or list of bools
"""
self.parts = parts
self.num_parts = len(parts)
self.num_params = sum([p.num_params for p in self.parts])
super(kern, self).__init__('kern')
self.add_parameters(*parts)
self.input_dim = input_dim
part_names = [k.name for k in self.parts]
self.name=''
for name in part_names:
self.name += name + '+'
self.name = self.name[:-1]
# deal with input_slices
if input_slices is None:
self.input_slices = [slice(None) for p in self.parts]
self.input_slices = [slice(None) for p in self._parameters_]
else:
assert len(input_slices) == len(self.parts)
assert len(input_slices) == len(self._parameters_)
self.input_slices = [sl if type(sl) is slice else slice(None) for sl in input_slices]
for p in self.parts:
for p in self._parameters_:
assert isinstance(p, Kernpart), "bad kernel part"
self.compute_param_slices()
self._parameters_ = []
for p in self.parts:
self._parameters_.extend(p._parameters_)
super(kern, self).__init__()
#Parameterized_old.__init__(self)
def parameters_changed(self):
[p.parameters_changed() for p in self.parts]
[p.parameters_changed() for p in self._parameters_]
def getstate(self):
"""
Get the current state of the class,
here just all the indices, rest can get recomputed
"""
return Parameterized.getstate(self) + [self.parts,
self.num_parts,
self.num_params,
return Parameterized.getstate(self) + [self._parameters_,
#self.num_params,
self.input_dim,
self.input_slices,
self.param_slices
self._param_slices_
]
def setstate(self, state):
self.param_slices = state.pop()
self._param_slices_ = state.pop()
self.input_slices = state.pop()
self.input_dim = state.pop()
self.num_params = state.pop()
self.num_parts = state.pop()
self.parts = state.pop()
#self.num_params = state.pop()
self._parameters_ = state.pop()
Parameterized.setstate(self, state)
@ -107,7 +94,7 @@ class kern(Parameterized):
xticklabels = []
bars = []
x0 = 0
for p in self.parts:
for p in self._parameters_:
c = Tango.nextMedium()
if hasattr(p, 'ARD') and p.ARD:
if title is None:
@ -158,28 +145,28 @@ class kern(Parameterized):
ax.legend()
return ax
def _transform_gradients(self, g):
"""
Apply the transformations of the kernel so that the returned vector
represents the gradient in the transformed space (i.e. that given by
get_params_transformed())
:param g: the gradient vector for the current model, usually created by dK_dtheta
"""
x = self._get_params()
[np.place(g, index, g[index] * constraint.gradfactor(x[index]))
for constraint, index in self._constraints_.iteritems() if constraint is not __fixed__]
# for constraint, index in self._constraints_.iteritems():
# if constraint != __fixed__:
# g[index] = g[index] * constraint.gradfactor(x[index])
#[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
[np.put(g, i, v) for i, v in [[i, t.sum()] for p in self._parameters_ for t,i in p._tied_to_me_.iteritems()]]
# if len(self.tied_indices) or len(self.fixed_indices):
# to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
# return np.delete(g, to_remove)
# else:
if self._fixes_ is not None: return g[self._fixes_]
return g
# def _transform_gradients(self, g):
# """
# Apply the transformations of the kernel so that the returned vector
# represents the gradient in the transformed space (i.e. that given by
# get_params_transformed())
#
# :param g: the gradient vector for the current model, usually created by dK_dtheta
# """
# x = self._get_params()
# [np.place(g, index, g[index] * constraint.gradfactor(x[index]))
# for constraint, index in self.constraints.iteritems() if constraint is not __fixed__]
# # for constraint, index in self.constraints.iteritems():
# # if constraint != __fixed__:
# # g[index] = g[index] * constraint.gradfactor(x[index])
# #[np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
# [np.put(g, i, v) for i, v in [[i, t.sum()] for p in self._parameters_ for t,i in p._tied_to_me_.iteritems()]]
# # if len(self.tied_indices) or len(self.fixed_indices):
# # to_remove = np.hstack((self.fixed_indices + [t[1:] for t in self.tied_indices]))
# # return np.delete(g, to_remove)
# # else:
# if self._fixes_ is not None: return g[self._fixes_]
# return g
# x = self._get_params()
# [np.put(x, i, x * t.gradfactor(x[i])) for i, t in zip(self.constrained_indices, self.constraints)]
# [np.put(g, i, v) for i, v in [(t[0], np.sum(g[t])) for t in self.tied_indices]]
@ -189,16 +176,6 @@ class kern(Parameterized):
# else:
# return g
def compute_param_slices(self):
"""
Create a set of slices that can index the parameters of each part.
"""
self.param_slices = []
count = 0
for p in self.parts:
self.param_slices.append(slice(count, count + p.num_params))
count += p.num_params
def __add__(self, other):
""" Overloading of the '+' operator. for more control, see self.add """
return self.add(other)
@ -224,7 +201,7 @@ class kern(Parameterized):
other_input_indices = [sl.indices(other.input_dim) for sl in other.input_slices]
other_input_slices = [slice(i[0] + self.input_dim, i[1] + self.input_dim, i[2]) for i in other_input_indices]
newkern = kern(D, self.parts + other.parts, self_input_slices + other_input_slices)
newkern = kern(D, self._parameters_ + other._parameters_, self_input_slices + other_input_slices)
# transfer constraints:
# newkern.constrained_indices = self.constrained_indices + [x + self.num_params for x in other.constrained_indices]
@ -235,7 +212,7 @@ class kern(Parameterized):
# newkern.tied_indices = self.tied_indices + [self.num_params + x for x in other.tied_indices]
else:
assert self.input_dim == other.input_dim
newkern = kern(self.input_dim, self.parts + other.parts, self.input_slices + other.input_slices)
newkern = kern(self.input_dim, self._parameters_ + other._parameters_, self.input_slices + other.input_slices)
# transfer constraints:
# newkern.constrained_indices = self.constrained_indices + [i + self.num_params for i in other.constrained_indices]
# newkern.constraints = self.constraints + other.constraints
@ -244,8 +221,8 @@ class kern(Parameterized):
# newkern.tied_indices = self.tied_indices + [self.num_params + x for x in other.tied_indices]
[newkern._add_constrain(param, transform, warning=False)
for param, transform in itertools.izip(
*itertools.chain(self._constraints_.iteritems(),
other._constraints_.iteritems()))]
*itertools.chain(self.constraints.iteritems(),
other.constraints.iteritems()))]
newkern._fixes_ = ((self._fixes_ or 0) + (other._fixes_ or 0)) or None
return newkern
@ -279,73 +256,73 @@ class kern(Parameterized):
s1[sl1], s2[sl2] = [True], [True]
slices += [s1 + s2]
newkernparts = [prod(k1, k2, tensor) for k1, k2 in itertools.product(K1.parts, K2.parts)]
newkernparts = [prod(k1, k2, tensor) for k1, k2 in itertools.product(K1._parameters_, K2._parameters_)]
if tensor:
newkern = kern(K1.input_dim + K2.input_dim, newkernparts, slices)
else:
newkern = kern(K1.input_dim, newkernparts, slices)
newkern._follow_constrains(K1, K2)
#newkern._follow_constrains(K1, K2)
return newkern
def _follow_constrains(self, K1, K2):
# Build the array that allows to go from the initial indices of the param to the new ones
K1_param = []
n = 0
for k1 in K1.parts:
K1_param += [range(n, n + k1.num_params)]
n += k1.num_params
n = 0
K2_param = []
for k2 in K2.parts:
K2_param += [range(K1.num_params + n, K1.num_params + n + k2.num_params)]
n += k2.num_params
index_param = []
for p1 in K1_param:
for p2 in K2_param:
index_param += p1 + p2
index_param = np.array(index_param)
# Get the ties and constrains of the kernels before the multiplication
prev_ties = K1.tied_indices + [arr + K1.num_params for arr in K2.tied_indices]
prev_constr_ind = [K1.constrained_indices] + [K1.num_params + i for i in K2.constrained_indices]
prev_constr = K1.constraints + K2.constraints
# prev_constr_fix = K1.fixed_indices + [arr + K1.num_params for arr in K2.fixed_indices]
# prev_constr_fix_values = K1.fixed_values + K2.fixed_values
# follow the previous ties
for arr in prev_ties:
for j in arr:
index_param[np.where(index_param == j)[0]] = arr[0]
# ties and constrains
for i in range(K1.num_params + K2.num_params):
index = np.where(index_param == i)[0]
if index.size > 1:
self.tie_params(index)
for i, t in zip(prev_constr_ind, prev_constr):
self.constrain(np.where(index_param == i)[0], t)
def _get_params(self):
return np.hstack(self._parameters_)
return np.hstack([p._get_params() for p in self.parts])
# def _follow_constrains(self, K1, K2):
#
# # Build the array that allows to go from the initial indices of the param to the new ones
# K1_param = []
# n = 0
# for k1 in K1.parts:
# K1_param += [range(n, n + k1.num_params)]
# n += k1.num_params
# n = 0
# K2_param = []
# for k2 in K2.parts:
# K2_param += [range(K1.num_params + n, K1.num_params + n + k2.num_params)]
# n += k2.num_params
# index_param = []
# for p1 in K1_param:
# for p2 in K2_param:
# index_param += p1 + p2
# index_param = np.array(index_param)
#
# # Get the ties and constrains of the kernels before the multiplication
# prev_ties = K1.tied_indices + [arr + K1.num_params for arr in K2.tied_indices]
#
# prev_constr_ind = [K1.constrained_indices] + [K1.num_params + i for i in K2.constrained_indices]
# prev_constr = K1.constraints + K2.constraints
#
# # prev_constr_fix = K1.fixed_indices + [arr + K1.num_params for arr in K2.fixed_indices]
# # prev_constr_fix_values = K1.fixed_values + K2.fixed_values
#
# # follow the previous ties
# for arr in prev_ties:
# for j in arr:
# index_param[np.where(index_param == j)[0]] = arr[0]
#
# # ties and constrains
# for i in range(K1.num_params + K2.num_params):
# index = np.where(index_param == i)[0]
# if index.size > 1:
# self.tie_params(index)
# for i, t in zip(prev_constr_ind, prev_constr):
# self.constrain(np.where(index_param == i)[0], t)
#
# def _get_params(self):
# return np.hstack(self._parameters_)
# return np.hstack([p._get_params() for p in self._parameters_])
def _set_params(self, x):
import ipdb;ipdb.set_trace()
[p._set_params(x[s]) for p, s in zip(self.parts, self.param_slices)]
# def _set_params(self, x):
# import ipdb;ipdb.set_trace()
# [p._set_params(x[s]) for p, s in zip(self._parameters_, self._param_slices_)]
def _get_param_names(self):
# this is a bit nasty: we want to distinguish between parts with the same name by appending a count
part_names = np.array([k.name for k in self.parts], dtype=np.str)
counts = [np.sum(part_names == ni) for i, ni in enumerate(part_names)]
cum_counts = [np.sum(part_names[i:] == ni) for i, ni in enumerate(part_names)]
names = [name + '_' + str(cum_count) if count > 1 else name for name, count, cum_count in zip(part_names, counts, cum_counts)]
return sum([[name + '_' + n for n in k._get_param_names()] for name, k in zip(names, self.parts)], [])
# def _get_param_names(self):
# # this is a bit nasty: we want to distinguish between parts with the same name by appending a count
# part_names = np.array([k.name for k in self._parameters_], dtype=np.str)
# counts = [np.sum(part_names == ni) for i, ni in enumerate(part_names)]
# cum_counts = [np.sum(part_names[i:] == ni) for i, ni in enumerate(part_names)]
# names = [name + '_' + str(cum_count) if count > 1 else name for name, count, cum_count in zip(part_names, counts, cum_counts)]
#
# return sum([[name + '_' + n for n in k._get_param_names()] for name, k in zip(names, self._parameters_)], [])
def K(self, X, X2=None, which_parts='all'):
"""
@ -357,17 +334,17 @@ class kern(Parameterized):
handles this as X2 == X.
:param which_parts: a list of booleans detailing whether to include
each of the part functions. By default, 'all'
indicates [True]*self.num_parts
indicates all parts
"""
if which_parts == 'all':
which_parts = [True] * self.num_parts
which_parts = [True] * self.size
assert X.shape[1] == self.input_dim
if X2 is None:
target = np.zeros((X.shape[0], X.shape[0]))
[p.K(X[:, i_s], None, target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
[p.K(X[:, i_s], None, target=target) for p, i_s, part_i_used in zip(self._parameters_, self.input_slices, which_parts) if part_i_used]
else:
target = np.zeros((X.shape[0], X2.shape[0]))
[p.K(X[:, i_s], X2[:, i_s], target=target) for p, i_s, part_i_used in zip(self.parts, self.input_slices, which_parts) if part_i_used]
[p.K(X[:, i_s], X2[:, i_s], target=target) for p, i_s, part_i_used in zip(self._parameters_, self.input_slices, which_parts) if part_i_used]
return target
def dK_dtheta(self, dL_dK, X, X2=None):
@ -384,11 +361,11 @@ class kern(Parameterized):
returns: dL_dtheta
"""
assert X.shape[1] == self.input_dim
target = np.zeros(self.num_params)
target = np.zeros(self.size)
if X2 is None:
[p.dK_dtheta(dL_dK, X[:, i_s], None, target[ps]) for p, i_s, ps, in zip(self.parts, self.input_slices, self.param_slices)]
[p.dK_dtheta(dL_dK, X[:, i_s], None, target[ps]) for p, i_s, ps, in zip(self._parameters_, self.input_slices, self._param_slices_)]
else:
[p.dK_dtheta(dL_dK, X[:, i_s], X2[:, i_s], target[ps]) for p, i_s, ps, in zip(self.parts, self.input_slices, self.param_slices)]
[p.dK_dtheta(dL_dK, X[:, i_s], X2[:, i_s], target[ps]) for p, i_s, ps, in zip(self._parameters_, self.input_slices, self._param_slices_)]
return self._transform_gradients(target)
@ -404,18 +381,18 @@ class kern(Parameterized):
target = np.zeros_like(X)
if X2 is None:
[p.dK_dX(dL_dK, X[:, i_s], None, target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dK_dX(dL_dK, X[:, i_s], None, target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
else:
[p.dK_dX(dL_dK, X[:, i_s], X2[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dK_dX(dL_dK, X[:, i_s], X2[:, i_s], target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
return target
def Kdiag(self, X, which_parts='all'):
"""Compute the diagonal of the covariance function for inputs X."""
if which_parts == 'all':
which_parts = [True] * self.num_parts
which_parts = [True] * self.size
assert X.shape[1] == self.input_dim
target = np.zeros(X.shape[0])
[p.Kdiag(X[:, i_s], target=target) for p, i_s, part_on in zip(self.parts, self.input_slices, which_parts) if part_on]
[p.Kdiag(X[:, i_s], target=target) for p, i_s, part_on in zip(self._parameters_, self.input_slices, which_parts) if part_on]
return target
def dKdiag_dtheta(self, dL_dKdiag, X):
@ -423,49 +400,49 @@ class kern(Parameterized):
assert X.shape[1] == self.input_dim
assert dL_dKdiag.size == X.shape[0]
target = np.zeros(self.num_params)
[p.dKdiag_dtheta(dL_dKdiag, X[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)]
[p.dKdiag_dtheta(dL_dKdiag, X[:, i_s], target[ps]) for p, i_s, ps in zip(self._parameters_, self.input_slices, self._param_slices_)]
return self._transform_gradients(target)
def dKdiag_dX(self, dL_dKdiag, X):
assert X.shape[1] == self.input_dim
target = np.zeros_like(X)
[p.dKdiag_dX(dL_dKdiag, X[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dKdiag_dX(dL_dKdiag, X[:, i_s], target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
return target
def psi0(self, Z, mu, S):
target = np.zeros(mu.shape[0])
[p.psi0(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
[p.psi0(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self._parameters_, self.input_slices)]
return target
def dpsi0_dtheta(self, dL_dpsi0, Z, mu, S):
target = np.zeros(self.num_params)
[p.dpsi0_dtheta(dL_dpsi0, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)]
[p.dpsi0_dtheta(dL_dpsi0, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self._parameters_, self._param_slices_, self.input_slices)]
return self._transform_gradients(target)
def dpsi0_dmuS(self, dL_dpsi0, Z, mu, S):
target_mu, target_S = np.zeros_like(mu), np.zeros_like(S)
[p.dpsi0_dmuS(dL_dpsi0, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dpsi0_dmuS(dL_dpsi0, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
return target_mu, target_S
def psi1(self, Z, mu, S):
target = np.zeros((mu.shape[0], Z.shape[0]))
[p.psi1(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
[p.psi1(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self._parameters_, self.input_slices)]
return target
def dpsi1_dtheta(self, dL_dpsi1, Z, mu, S):
target = np.zeros((self.num_params))
[p.dpsi1_dtheta(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self.parts, self.param_slices, self.input_slices)]
[p.dpsi1_dtheta(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, ps, i_s in zip(self._parameters_, self._param_slices_, self.input_slices)]
return self._transform_gradients(target)
def dpsi1_dZ(self, dL_dpsi1, Z, mu, S):
target = np.zeros_like(Z)
[p.dpsi1_dZ(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dpsi1_dZ(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
return target
def dpsi1_dmuS(self, dL_dpsi1, Z, mu, S):
"""return shapes are num_samples,num_inducing,input_dim"""
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
[p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dpsi1_dmuS(dL_dpsi1, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
return target_mu, target_S
def psi2(self, Z, mu, S):
@ -478,13 +455,13 @@ class kern(Parameterized):
"""
target = np.zeros((mu.shape[0], Z.shape[0], Z.shape[0]))
[p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self.parts, self.input_slices)]
[p.psi2(Z[:, i_s], mu[:, i_s], S[:, i_s], target) for p, i_s in zip(self._parameters_, self.input_slices)]
# compute the "cross" terms
# TODO: input_slices needed
crossterms = 0
for [p1, i_s1], [p2, i_s2] in itertools.combinations(zip(self.parts, self.input_slices), 2):
for [p1, i_s1], [p2, i_s2] in itertools.combinations(zip(self._parameters_, self.input_slices), 2):
if i_s1 == i_s2:
# TODO psi1 this must be faster/better/precached/more nice
tmp1 = np.zeros((mu.shape[0], Z.shape[0]))
@ -501,14 +478,14 @@ class kern(Parameterized):
def dpsi2_dtheta(self, dL_dpsi2, Z, mu, S):
"""Gradient of the psi2 statistics with respect to the parameters."""
target = np.zeros(self.num_params)
[p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self.parts, self.input_slices, self.param_slices)]
[p.dpsi2_dtheta(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[ps]) for p, i_s, ps in zip(self._parameters_, self.input_slices, self._param_slices_)]
# compute the "cross" terms
# TODO: better looping, input_slices
for i1, i2 in itertools.permutations(range(len(self.parts)), 2):
p1, p2 = self.parts[i1], self.parts[i2]
for i1, i2 in itertools.permutations(range(len(self._parameters_)), 2):
p1, p2 = self._parameters_[i1], self._parameters_[i2]
# ipsl1, ipsl2 = self.input_slices[i1], self.input_slices[i2]
ps1, ps2 = self.param_slices[i1], self.param_slices[i2]
ps1, ps2 = self._param_slices_[i1], self._param_slices_[i2]
tmp = np.zeros((mu.shape[0], Z.shape[0]))
p1.psi1(Z, mu, S, tmp)
@ -518,12 +495,12 @@ class kern(Parameterized):
def dpsi2_dZ(self, dL_dpsi2, Z, mu, S):
target = np.zeros_like(Z)
[p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dpsi2_dZ(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
# target *= 2
# compute the "cross" terms
# TODO: we need input_slices here.
for p1, p2 in itertools.permutations(self.parts, 2):
for p1, p2 in itertools.permutations(self._parameters_, 2):
# if p1.name == 'linear' and p2.name == 'linear':
# raise NotImplementedError("We don't handle linear/linear cross-terms")
tmp = np.zeros((mu.shape[0], Z.shape[0]))
@ -534,11 +511,11 @@ class kern(Parameterized):
def dpsi2_dmuS(self, dL_dpsi2, Z, mu, S):
target_mu, target_S = np.zeros((2, mu.shape[0], mu.shape[1]))
[p.dpsi2_dmuS(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self.parts, self.input_slices)]
[p.dpsi2_dmuS(dL_dpsi2, Z[:, i_s], mu[:, i_s], S[:, i_s], target_mu[:, i_s], target_S[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
# compute the "cross" terms
# TODO: we need input_slices here.
for p1, p2 in itertools.permutations(self.parts, 2):
for p1, p2 in itertools.permutations(self._parameters_, 2):
# if p1.name == 'linear' and p2.name == 'linear':
# raise NotImplementedError("We don't handle linear/linear cross-terms")
tmp = np.zeros((mu.shape[0], Z.shape[0]))
@ -549,7 +526,7 @@ class kern(Parameterized):
def plot(self, x=None, plot_limits=None, which_parts='all', resolution=None, *args, **kwargs):
if which_parts == 'all':
which_parts = [True] * self.num_parts
which_parts = [True] * self.size
if self.input_dim == 1:
if x is None:
x = np.zeros((1, 1))

10
GPy/kern/parts/bias.py
View file

@ -15,13 +15,9 @@ class Bias(Kernpart):
:param variance: the variance of the kernel
:type variance: float
"""
super(Bias, self).__init__(input_dim)
self.input_dim = input_dim
self.num_params = 1
self.name = 'bias'
self.variance = Param(lambda: self.name+"_variance", variance, None)
self.set_as_parameters(self.variance)
super(Bias, self).__init__(input_dim, 'bias')
self.variance = Param("variance", variance, None)
self.add_parameter(self.variance)
#self._set_params(np.array([variance]).flatten())
# def _get_params(self):

15
GPy/kern/parts/kernpart.py
View file

@ -1,10 +1,10 @@
# Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt)
from ...core.parameter import Param
#from ...core.parameterized.Parameterized import set_as_parameter
from GPy.core.parameterized import Parameterized
class Kernpart(object):
def __init__(self,input_dim):
class Kernpart(Parameterized):
def __init__(self,input_dim,name):
"""
The base class for a kernpart: a positive definite function
which forms part of a covariance function (kernel).
@ -14,21 +14,14 @@ class Kernpart(object):
Do not instantiate.
"""
super(Kernpart, self).__init__(name)
# the input dimensionality for the covariance
self.input_dim = input_dim
# the number of optimisable parameters
self.num_params = 1
# the name of the covariance function.
self.name = 'unnamed'
# link to parameterized objects
self._parameters_ = []
def set_as_parameters(self, *params):
self._parameters_.extend(params)
def parameters_changed(self):
pass
# def set_as_parameter_named(self, name, gradient, index=None, *args, **kwargs):
# """
# :param names: name of parameter to set as parameter

38
GPy/kern/parts/linear.py
View file

@ -7,6 +7,7 @@ import numpy as np
from ...util.linalg import tdot
from ...util.misc import fast_array_equal
from scipy import weave
from GPy.core.parameter import Param
class Linear(Kernpart):
"""
@ -26,11 +27,9 @@ class Linear(Kernpart):
"""
def __init__(self, input_dim, variances=None, ARD=False):
self.input_dim = input_dim
super(Linear, self).__init__(input_dim, 'linear')
self.ARD = ARD
if ARD == False:
self.num_params = 1
self.name = 'linear'
if variances is not None:
variances = np.asarray(variances)
assert variances.size == 1, "Only one variance needed for non-ARD kernel"
@ -38,32 +37,33 @@ class Linear(Kernpart):
variances = np.ones(1)
self._Xcache, self._X2cache = np.empty(shape=(2,))
else:
self.num_params = self.input_dim
self.name = 'linear'
if variances is not None:
variances = np.asarray(variances)
assert variances.size == self.input_dim, "bad number of lengthscales"
else:
variances = np.ones(self.input_dim)
self._set_params(variances.flatten())
self.variances = Param('variances', variances)
self.add_parameters(self.variances)
# initialize cache
self._Z, self._mu, self._S = np.empty(shape=(3, 1))
self._X, self._X2, self._params = np.empty(shape=(3, 1))
self._X, self._X2 = np.empty(shape=(2, 1))
def _get_params(self):
return self.variances
def _set_params(self, x):
assert x.size == (self.num_params)
self.variances = x
# def _get_params(self):
# return self.variances
#
# def _set_params(self, x):
# assert x.size == (self.num_params)
# self.variances = x
def parameters_changed(self):
self.variances2 = np.square(self.variances)
def _get_param_names(self):
if self.num_params == 1:
return ['variance']
else:
return ['variance_%i' % i for i in range(self.variances.size)]
#
# def _get_param_names(self):
# if self.num_params == 1:
# return ['variance']
# else:
# return ['variance_%i' % i for i in range(self.variances.size)]
def K(self, X, X2, target):
if self.ARD:

50
GPy/kern/parts/prod.py
View file

@ -18,35 +18,37 @@ class Prod(Kernpart):
"""
def __init__(self,k1,k2,tensor=False):
self.num_params = k1.num_params + k2.num_params
self.name = '['+k1.name + '**' + k2.name +']'
self.k1 = k1
self.k2 = k2
if tensor:
self.input_dim = k1.input_dim + k2.input_dim
self.slice1 = slice(0,self.k1.input_dim)
self.slice2 = slice(self.k1.input_dim,self.k1.input_dim+self.k2.input_dim)
super(Prod, self).__init__(k1.input_dim + k2.input_dim, '['+k1.name + '**' + k2.name +']')
else:
assert k1.input_dim == k2.input_dim, "Error: The input spaces of the kernels to sum don't have the same dimension."
self.input_dim = k1.input_dim
self.slice1 = slice(0,self.input_dim)
self.slice2 = slice(0,self.input_dim)
super(Prod, self).__init__(k1.input_dim, '['+k1.name + '*' + k2.name +']')
#self.num_params = k1.num_params + k2.num_params
self.k1 = k1
self.k2 = k2
# if tensor:
# self.slice1 = slice(0,self.k1.input_dim)
# self.slice2 = slice(self.k1.input_dim,self.k1.input_dim+self.k2.input_dim)
# else:
# self.slice1 = slice(0,self.input_dim)
# self.slice2 = slice(0,self.input_dim)
self._X, self._X2, self._params = np.empty(shape=(3,1))
self._set_params(np.hstack((k1._get_params(),k2._get_params())))
self._X, self._X2 = np.empty(shape=(2,1))
self.add_parameters(self.k1, self.k2)
# self._set_params(np.hstack((k1._get_params(),k2._get_params())))
def _get_params(self):
"""return the value of the parameters."""
return np.hstack((self.k1._get_params(), self.k2._get_params()))
def _set_params(self,x):
"""set the value of the parameters."""
self.k1._set_params(x[:self.k1.num_params])
self.k2._set_params(x[self.k1.num_params:])
def _get_param_names(self):
"""return parameter names."""
return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
# def _get_params(self):
# """return the value of the parameters."""
# return np.hstack((self.k1._get_params(), self.k2._get_params()))
#
# def _set_params(self,x):
# """set the value of the parameters."""
# self.k1._set_params(x[:self.k1.num_params])
# self.k2._set_params(x[self.k1.num_params:])
#
# def _get_param_names(self):
# """return parameter names."""
# return [self.k1.name + '_' + param_name for param_name in self.k1._get_param_names()] + [self.k2.name + '_' + param_name for param_name in self.k2._get_param_names()]
def K(self,X,X2,target):
self._K_computations(X,X2)

62
GPy/kern/parts/rbf.py
View file

@ -33,20 +33,17 @@ class RBF(Kernpart):
"""
def __init__(self, input_dim, variance=1., lengthscale=None, ARD=False, name='rbf'):
super(RBF, self).__init__(input_dim)
super(RBF, self).__init__(input_dim, name)
self.input_dim = input_dim
self.name = name
self.ARD = ARD
if not ARD:
self.num_params = 2
if lengthscale is not None:
lengthscale = np.asarray(lengthscale)
assert lengthscale.size == 1, "Only one lengthscale needed for non-ARD kernel"
else:
lengthscale = np.ones(1)
else:
self.num_params = self.input_dim + 1
if lengthscale is not None:
lengthscale = np.asarray(lengthscale)
assert lengthscale.size == self.input_dim, "bad number of lengthscales"
@ -54,9 +51,10 @@ class RBF(Kernpart):
lengthscale = np.ones(self.input_dim)
#self._set_params(np.hstack((variance, lengthscale.flatten())))
self.variance = Param(lambda: self.name+'_variance', variance, None)
self.lengthscale = Param(lambda: self.name+'_lengthscale', lengthscale, None)
self.set_as_parameters(self.variance, self.lengthscale)
self.variance = Param('variance', variance, None)
self.lengthscale = Param('lengthscale', lengthscale, None)
self.add_parameters(self.variance, self.lengthscale)
# self.set_as_parameter('variance', self.variance, None)
# self.set_as_parameter('lengthscale', self.lengthscale, None)
@ -68,37 +66,35 @@ class RBF(Kernpart):
self.weave_options = {'headers' : ['<omp.h>'],
'extra_compile_args': ['-fopenmp -O3'], # -march=native'],
'extra_link_args' : ['-lgomp']}
def parameters_changed(self):
self.lengthscale2 = np.square(self.lengthscale)
# reset cached results
self._X, self._X2, self._params_save = np.empty(shape=(3, 1))
self._Z, self._mu, self._S = np.empty(shape=(3, 1)) # cached versions of Z,mu,S
def _get_params(self):
return np.hstack((self.variance, self.lengthscale))
#
def _set_params(self, x):
assert x.size == (self.num_params)
#self.variance = x[0]
#self.lengthscale = x[1:]
#self.lengthscale2 = np.square(self.lengthscale)
# reset cached results
#self._X, self._X2, self._params_save = np.empty(shape=(3, 1))
#self._Z, self._mu, self._S = np.empty(shape=(3, 1)) # cached versions of Z,mu,S
def _get_param_names(self):
if self.num_params == 2:
return ['variance', 'lengthscale']
else:
return ['variance'] + ['lengthscale_%i' % i for i in range(self.lengthscale.size)]
self._X, self._X2 = np.empty(shape=(2, 1))
self._Z, self._mu, self._S = np.empty(shape=(3, 1)) # cached versions of Z,mu,S
def _dK_dvariance(self, model, target):
pass
# def _get_params(self):
# return np.hstack((self.variance, self.lengthscale))
# #
# def _set_params(self, x):
# assert x.size == (self.num_params)
# #self.variance = x[0]
# #self.lengthscale = x[1:]
#
# #self.lengthscale2 = np.square(self.lengthscale)
#
# # reset cached results
# #self._X, self._X2, self._params_save = np.empty(shape=(3, 1))
# #self._Z, self._mu, self._S = np.empty(shape=(3, 1)) # cached versions of Z,mu,S
#
# def _get_param_names(self):
# if self.num_params == 2:
# return ['variance', 'lengthscale']
# else:
# return ['variance'] + ['lengthscale_%i' % i for i in range(self.lengthscale.size)]
def K(self, X, X2, target):
self._K_computations(X, X2)
@ -243,10 +239,10 @@ class RBF(Kernpart):
#---------------------------------------#
def _K_computations(self, X, X2):
params = self._get_params()
if not (fast_array_equal(X, self._X) and fast_array_equal(X2, self._X2) and fast_array_equal(self._params_save , params)):
#params = self._get_params()
if not (fast_array_equal(X, self._X) and fast_array_equal(X2, self._X2)):# and fast_array_equal(self._params_save , params)):
self._X = X.copy()
self._params_save = params.copy()
#self._params_save = params.copy()
if X2 is None:
self._X2 = None
X = X / self.lengthscale

53
GPy/kern/parts/rbfcos.py
View file

@ -5,6 +5,7 @@
from kernpart import Kernpart
import numpy as np
from ...core.parameter import Param
class RBFCos(Kernpart):
def __init__(self,input_dim,variance=1.,frequencies=None,bandwidths=None,ARD=False):
@ -41,29 +42,30 @@ class RBFCos(Kernpart):
else:
bandwidths = np.ones(1)
self.variance = Param('variance', variance)
self.frequencies = Param('frequencies', frequencies)
self.bandwidths = Param('bandwidths', bandwidths)
#initialise cache
self._X, self._X2, self._params = np.empty(shape=(3,1))
self._X, self._X2 = np.empty(shape=(3,1))
self._set_params(np.hstack((variance,frequencies.flatten(),bandwidths.flatten())))
# def _get_params(self):
# return np.hstack((self.variance,self.frequencies, self.bandwidths))
# def _set_params(self,x):
# assert x.size==(self.num_params)
# if self.ARD:
# self.variance = x[0]
# self.frequencies = x[1:1+self.input_dim]
# self.bandwidths = x[1+self.input_dim:]
# else:
# self.variance, self.frequencies, self.bandwidths = x
def _get_params(self):
return np.hstack((self.variance,self.frequencies, self.bandwidths))
def _set_params(self,x):
assert x.size==(self.num_params)
if self.ARD:
self.variance = x[0]
self.frequencies = x[1:1+self.input_dim]
self.bandwidths = x[1+self.input_dim:]
else:
self.variance, self.frequencies, self.bandwidths = x
def _get_param_names(self):
if self.num_params == 3:
return ['variance','frequency','bandwidth']
else:
return ['variance']+['frequency_%i'%i for i in range(self.input_dim)]+['bandwidth_%i'%i for i in range(self.input_dim)]
# def _get_param_names(self):
# if self.num_params == 3:
# return ['variance','frequency','bandwidth']
# else:
# return ['variance']+['frequency_%i'%i for i in range(self.input_dim)]+['bandwidth_%i'%i for i in range(self.input_dim)]
def K(self,X,X2,target):
self._K_computations(X,X2)
@ -94,6 +96,11 @@ class RBFCos(Kernpart):
def dKdiag_dX(self,dL_dKdiag,X,target):
pass
def parameters_changed(self):
self._rbf_part = np.exp(-2.*np.pi**2*np.sum(self._dist2*self.bandwidths,-1))
self._cos_part = np.prod(np.cos(2.*np.pi*self._dist*self.frequencies),-1)
self._dvar = self._rbf_part*self._cos_part
def _K_computations(self,X,X2):
if not (np.all(X==self._X) and np.all(X2==self._X2)):
if X2 is None: X2 = X
@ -107,11 +114,3 @@ class RBFCos(Kernpart):
#ensure the next section is computed:
self._params = np.empty(self.num_params)
if not np.all(self._params == self._get_params()):
self._params == self._get_params().copy()
self._rbf_part = np.exp(-2.*np.pi**2*np.sum(self._dist2*self.bandwidths,-1))
self._cos_part = np.prod(np.cos(2.*np.pi*self._dist*self.frequencies),-1)
self._dvar = self._rbf_part*self._cos_part

23
GPy/kern/parts/spline.py
View file

@ -5,6 +5,7 @@
from kernpart import Kernpart
import numpy as np
import hashlib
from ...core.parameter import Param
def theta(x):
"""Heaviside step function"""
return np.where(x>=0.,1.,0.)
@ -25,16 +26,18 @@ class Spline(Kernpart):
assert self.input_dim==1
self.num_params = 1
self.name = 'spline'
self._set_params(np.squeeze(variance))
def _get_params(self):
return self.variance
def _set_params(self,x):
self.variance = x
def _get_param_names(self):
return ['variance']
self.variance = Param('variance', variance)
self.lengthscale = Param('lengthscale', lengthscale)
self.add_parameters(self.variance, self.lengthscale)
# def _get_params(self):
# return self.variance
#
# def _set_params(self,x):
# self.variance = x
#
# def _get_param_names(self):
# return ['variance']
def K(self,X,X2,target):
assert np.all(X>0), "Spline covariance is for +ve domain only. TODO: symmetrise"

15
GPy/kern/parts/symmetric.py
View file

@ -24,19 +24,8 @@ class Symmetric(Kernpart):
self.num_params = k.num_params
self.name = k.name + '_symm'
self.k = k
self._set_params(k._get_params())
def _get_params(self):
"""return the value of the parameters."""
return self.k._get_params()
def _set_params(self,x):
"""set the value of the parameters."""
self.k._set_params(x)
def _get_param_names(self):
"""return parameter names."""
return self.k._get_param_names()
self.add_parameter(k)
#self._set_params(k._get_params())
def K(self,X,X2,target):
"""Compute the covariance matrix between X and X2."""

8
GPy/kern/parts/white.py
View file

@ -15,12 +15,10 @@ class White(Kernpart):
:type variance: float
"""
def __init__(self,input_dim,variance=1.):
super(White, self).__init__(input_dim)
super(White, self).__init__(input_dim, 'white')
self.input_dim = input_dim
self.num_params = 1
self.name = 'white'
self.variance = Param(lambda: self.name+'_variance', variance, None)
self.set_as_parameters(self.variance)
self.variance = Param('variance', variance, None)
self.add_parameters(self.variance)
# self._set_params(np.array([variance]).flatten())
self._psi1 = 0 # TODO: more elegance here

10
GPy/likelihoods/gaussian.py
View file

@ -15,7 +15,7 @@ class Gaussian(likelihood):
:type normalize: False|True
"""
def __init__(self, data, variance=1., normalize=False):
super(Gaussian, self).__init__()
super(Gaussian, self).__init__('gaussian')
self.is_heteroscedastic = False
self.num_params = 1
self.Z = 0. # a correction factor which accounts for the approximation made
@ -33,11 +33,12 @@ class Gaussian(likelihood):
self.set_data(data)
self.variance = Param('noise_variance', variance, None)
self.set_as_parameters(self.variance)
self.variance = Param('variance', variance)
self._variance = variance + 1
self.add_parameter(self.variance)
# self._set_params(np.asarray(variance))
@ -62,6 +63,7 @@ class Gaussian(likelihood):
# return ["noise_variance"]
#
# def _set_params(self, x):
# self.variance = x[0]
def parameters_changed(self):
if np.any(self._variance != self.variance):
if np.all(self.variance == 0.):#special case of zero noise

4
GPy/likelihoods/likelihood.py
View file

@ -21,8 +21,8 @@ class likelihood(Parameterized):
- self.V : self.precision * self.Y
"""
def __init__(self):
Parameterized.__init__(self)
def __init__(self, name=None):
Parameterized.__init__(self, name)
# def _get_params(self):
# raise NotImplementedError

2
GPy/testing/kernel_tests.py
View file

@ -78,7 +78,7 @@ class KernelTests(unittest.TestCase):
self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose))
def test_heterokernel(self):
kern = GPy.kern.hetero(5, mapping=GPy.mappings.Linear(5, 1), transform=GPy.core.transformations.logexp())
kern = GPy.kern.hetero(5, mapping=GPy.mappings.Linear(5, 1), transform=GPy.core.transformations.Logexp())
self.assertTrue(GPy.kern.kern_test(kern, verbose=verbose))
def test_mlpkernel(self):

10
GPy/testing/parameter_testing.py
View file

@ -31,10 +31,10 @@ class Test(unittest.TestCase):
self.Y = numpy.random.multivariate_normal(numpy.zeros(self.N), K + numpy.eye(self.N) * .2, self.D).T
self.bgplvm = BayesianGPLVM(Gaussian(self.Y, variance=self.noise_variance), self.Q, self.X, self.X_variance, kernel=self.kern)
self.bgplvm.ensure_default_constraints(warning=False)
self.bgplvm.tie_params("noise_variance|white_variance")
self.bgplvm.constrain_fixed("rbf_var", warning=False)
# self.bgplvm = BayesianGPLVM(Gaussian(self.Y, variance=self.noise_variance), self.Q, self.X, self.X_variance, kernel=self.kern)
# self.bgplvm.ensure_default_constraints(warning=False)
# self.bgplvm.tie_params("noise_variance|white_variance")
# self.bgplvm.constrain_fixed("rbf_var", warning=False)
self.parameter = Parameterized([
Parameterized([
Param('X', self.X),
@ -110,7 +110,7 @@ class Test(unittest.TestCase):
self.parameter[''].unconstrain()
self.parameter.X.constrain_positive()
self.parameter.X[:,numpy.s_[0::2]].unconstrain_positive()
assert(numpy.alltrue(self.parameter._constraints_.indices()[0] == numpy.r_[1:self.N*self.Q:2]))
assert(numpy.alltrue(self.parameter.constraints.indices()[0] == numpy.r_[1:self.N*self.Q:2]))
def testNdarrayFunc(self):
assert(numpy.alltrue(self.parameter.X * self.parameter.X == self.X * self.X))