apunkt/GPy
1
0
Fork 0
mirror of https://github.com/SheffieldML/GPy.git synced 2026-07-23 17:01:06 +02:00
Code Issues Projects Releases Packages Wiki Activity Actions

q

Browse source 
Merge branch 'params' of github.com:SheffieldML/GPy into params
This commit is contained in:
James Hensman 2014-03-12 13:08:29 +00:00
commit adc79b1027
25 changed files with 493 additions and 328 deletions
Download patch file Download diff file Expand all files Collapse all files

7
GPy/core/gp.py
View file

@ -48,7 +48,7 @@ class GP(Model):
self.Y_metadata = None self.Y_metadata = None
assert isinstance(kernel, kern.Kern) assert isinstance(kernel, kern.Kern)
assert self.input_dim == kernel.input_dim #assert self.input_dim == kernel.input_dim
self.kern = kernel self.kern = kernel
assert isinstance(likelihood, likelihoods.Likelihood) assert isinstance(likelihood, likelihoods.Likelihood)
@ -68,6 +68,7 @@ class GP(Model):
def parameters_changed(self): def parameters_changed(self):
self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata) self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata)
self.likelihood.update_gradients(np.diag(grad_dict['dL_dK']))
self.kern.update_gradients_full(grad_dict['dL_dK'], self.X) self.kern.update_gradients_full(grad_dict['dL_dK'], self.X)
def log_likelihood(self): def log_likelihood(self):
@ -185,7 +186,7 @@ class GP(Model):
""" """
assert "matplotlib" in sys.modules, "matplotlib package has not been imported." assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
from ..plotting.matplot_dep import models_plots from ..plotting.matplot_dep import models_plots
models_plots.plot_fit_f(self,*args,**kwargs) return models_plots.plot_fit_f(self,*args,**kwargs)
def plot(self, *args, **kwargs): def plot(self, *args, **kwargs):
""" """
@ -206,7 +207,7 @@ class GP(Model):
""" """
assert "matplotlib" in sys.modules, "matplotlib package has not been imported." assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
from ..plotting.matplot_dep import models_plots from ..plotting.matplot_dep import models_plots
models_plots.plot_fit(self,*args,**kwargs) return models_plots.plot_fit(self,*args,**kwargs)
def _getstate(self): def _getstate(self):
""" """

13
GPy/core/model.py
View file

@ -253,7 +253,7 @@ class Model(Parameterized):
sgd.run() sgd.run()
self.optimization_runs.append(sgd) self.optimization_runs.append(sgd)
def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3): def _checkgrad(self, target_param=None, verbose=False, step=1e-6, tolerance=1e-3, _debug=False):
""" """
Check the gradient of the ,odel by comparing to a numerical Check the gradient of the ,odel by comparing to a numerical
estimate. If the verbose flag is passed, invividual estimate. If the verbose flag is passed, invividual
@ -303,7 +303,7 @@ class Model(Parameterized):
global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator) global_ratio = (f1 - f2) / np.where(denominator==0., 1e-32, denominator)
gloabl_diff = (f1 - f2) - denominator gloabl_diff = (f1 - f2) - denominator
return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gloabl_diff) < tolerance) return (np.abs(1. - global_ratio) < tolerance) or (np.abs(gloabl_diff) == 0)
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:
@ -339,7 +339,7 @@ class Model(Parameterized):
print "No free parameters to check" print "No free parameters to check"
return return
gradient = self.objective_function_gradients(x) gradient = self.objective_function_gradients(x).copy()
np.where(gradient == 0, 1e-312, gradient) np.where(gradient == 0, 1e-312, gradient)
ret = True ret = True
for nind, xind in itertools.izip(param_index, transformed_index): for nind, xind in itertools.izip(param_index, transformed_index):
@ -349,6 +349,13 @@ class Model(Parameterized):
xx[xind] -= 2.*step xx[xind] -= 2.*step
f2 = self.objective_function(xx) f2 = self.objective_function(xx)
numerical_gradient = (f1 - f2) / (2 * step) numerical_gradient = (f1 - f2) / (2 * step)
if _debug:
for p in self.kern.flattened_parameters:
p._parent_._debug=True
self.gradient[xind] = numerical_gradient
self._set_params_transformed(x)
for p in self.kern.flattened_parameters:
p._parent_._debug=False
if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind] if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind]
else: ratio = (f1 - f2) / (2 * step * gradient[xind]) else: ratio = (f1 - f2) / (2 * step * gradient[xind])
difference = np.abs((f1 - f2) / 2 / step - gradient[xind]) difference = np.abs((f1 - f2) / 2 / step - gradient[xind])

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

@ -6,12 +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 _Array(np.ndarray):
def __init__(self, dtype=float, buffer=None, offset=0,
strides=None, order=None, *args, **kwargs):
super(_Array, self).__init__(dtype=dtype, buffer=buffer, offset=offset,
strides=strides, order=order, *args, **kwargs)
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.
@ -22,7 +16,7 @@ class ObservableArray(np.ndarray, Observable):
__array_priority__ = -1 # Never give back ObservableArray __array_priority__ = -1 # Never give back ObservableArray
def __new__(cls, input_array, *a, **kw): def __new__(cls, input_array, *a, **kw):
if not isinstance(input_array, ObservableArray): if not isinstance(input_array, ObservableArray):
obj = np.atleast_1d(np.require(input_array, dtype=np.float64, requirements=['C', 'W'])).view(cls) obj = np.atleast_1d(np.require(input_array, dtype=np.float64, requirements=['W', 'C'])).view(cls)
else: obj = input_array else: obj = input_array
cls.__name__ = "ObservableArray\n " cls.__name__ = "ObservableArray\n "
super(ObservableArray, obj).__init__(*a, **kw) super(ObservableArray, obj).__init__(*a, **kw)

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

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

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

@ -15,9 +15,8 @@ Observable Pattern for patameterization
from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED from transformations import Transformation, Logexp, NegativeLogexp, Logistic, __fixed__, FIXED, UNFIXED
import numpy as np import numpy as np
import itertools
__updated__ = '2013-12-16' __updated__ = '2014-03-12'
class HierarchyError(Exception): class HierarchyError(Exception):
""" """
@ -40,8 +39,9 @@ class Observable(object):
as an observer. Every time the observable changes, it sends a notification with as an observer. Every time the observable changes, it sends a notification with
self as only argument to all its observers. self as only argument to all its observers.
""" """
_updated = True
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(Observable, self).__init__() super(Observable, self).__init__(*args, **kwargs)
self._observer_callables_ = [] self._observer_callables_ = []
def add_observer(self, observer, callble, priority=0): def add_observer(self, observer, callble, priority=0):
@ -160,7 +160,7 @@ class Parentable(object):
_parent_ = None _parent_ = None
_parent_index_ = None _parent_index_ = None
def __init__(self, *args, **kwargs): def __init__(self, *args, **kwargs):
super(Parentable, self).__init__() super(Parentable, self).__init__(*args, **kwargs)
def has_parent(self): def has_parent(self):
""" """
@ -206,7 +206,7 @@ class Gradcheckable(Parentable):
def __init__(self, *a, **kw): def __init__(self, *a, **kw):
super(Gradcheckable, self).__init__(*a, **kw) super(Gradcheckable, self).__init__(*a, **kw)
def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3): def checkgrad(self, verbose=0, step=1e-6, tolerance=1e-3, _debug=False):
""" """
Check the gradient of this parameter with respect to the highest parent's Check the gradient of this parameter with respect to the highest parent's
objective function. objective function.
@ -220,14 +220,15 @@ class Gradcheckable(Parentable):
:param flaot tolerance: the tolerance for the gradient ratio or difference. :param flaot tolerance: the tolerance for the gradient ratio or difference.
""" """
if self.has_parent(): if self.has_parent():
return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance) return self._highest_parent_._checkgrad(self, verbose=verbose, step=step, tolerance=tolerance, _debug=_debug)
return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance) return self._checkgrad(self[''], verbose=verbose, step=step, tolerance=tolerance, _debug=_debug)
def _checkgrad(self, param):
def _checkgrad(self, param, verbose=0, step=1e-6, tolerance=1e-3, _debug=False):
""" """
Perform the checkgrad on the model. Perform the checkgrad on the model.
TODO: this can be done more efficiently, when doing it inside here TODO: this can be done more efficiently, when doing it inside here
""" """
raise NotImplementedError, "Need log likelihood to check gradient against" raise HierarchyError, "This parameter is not in a model with a likelihood, and, therefore, cannot be gradient checked!"
class Nameable(Gradcheckable): class Nameable(Gradcheckable):
@ -274,7 +275,7 @@ class Indexable(object):
The raveled index of an object is the index for its parameters in a flattened int array. The raveled index of an object is the index for its parameters in a flattened int array.
""" """
def __init__(self, *a, **kw): def __init__(self, *a, **kw):
super(Indexable, self).__init__() super(Indexable, self).__init__(*a, **kw)
def _raveled_index(self): def _raveled_index(self):
""" """
@ -318,7 +319,7 @@ class Constrainable(Nameable, Indexable):
:func:`constrain()` and :func:`unconstrain()` are main methods here :func:`constrain()` and :func:`unconstrain()` are main methods here
""" """
def __init__(self, name, default_constraint=None, *a, **kw): def __init__(self, name, default_constraint=None, *a, **kw):
super(Constrainable, self).__init__(name=name, default_constraint=default_constraint, *a, **kw) super(Constrainable, self).__init__(name=name, *a, **kw)
self._default_constraint_ = default_constraint self._default_constraint_ = default_constraint
from index_operations import ParameterIndexOperations from index_operations import ParameterIndexOperations
self.constraints = ParameterIndexOperations() self.constraints = ParameterIndexOperations()
@ -830,16 +831,13 @@ class Parameterizable(OptimizationHandlable):
self._gradient_array_ = np.empty(self.size, dtype=np.float64) self._gradient_array_ = np.empty(self.size, dtype=np.float64)
self._param_slices_ = [] self._param_slices_ = []
for i, p in enumerate(self._parameters_): for i, p in enumerate(self._parameters_):
p._parent_ = self p._parent_ = self
p._parent_index_ = i p._parent_index_ = i
pslice = slice(old_size, old_size+p.size) pslice = slice(old_size, old_size+p.size)
# first connect all children # first connect all children
p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice]) p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice])
# then connect children to self # then connect children to self
self._param_array_[pslice] = p._param_array_.ravel()#, requirements=['C', 'W']).ravel(order='C') self._param_array_[pslice] = p._param_array_.ravel()#, requirements=['C', 'W']).ravel(order='C')
self._gradient_array_[pslice] = p._gradient_array_.ravel()#, requirements=['C', 'W']).ravel(order='C') self._gradient_array_[pslice] = p._gradient_array_.ravel()#, requirements=['C', 'W']).ravel(order='C')
@ -867,6 +865,7 @@ class Parameterizable(OptimizationHandlable):
#=========================================================================== #===========================================================================
def copy(self): def copy(self):
"""Returns a (deep) copy of the current model""" """Returns a (deep) copy of the current model"""
raise NotImplementedError, "Copy is not yet implemented, TODO: Observable hierarchy"
import copy import copy
from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView from .index_operations import ParameterIndexOperations, ParameterIndexOperationsView
from .lists_and_dicts import ArrayList from .lists_and_dicts import ArrayList

1
GPy/core/parameterization/variational.py
View file

@ -71,6 +71,7 @@ class VariationalPosterior(Parameterized):
self.shape = self.mean.shape self.shape = self.mean.shape
self.num_data, self.input_dim = self.mean.shape self.num_data, self.input_dim = self.mean.shape
self.add_parameters(self.mean, self.variance) self.add_parameters(self.mean, self.variance)
self.num_data, self.input_dim = self.mean.shape
if self.has_uncertain_inputs(): if self.has_uncertain_inputs():
assert self.variance.shape == self.mean.shape, "need one variance per sample and dimenion" assert self.variance.shape == self.mean.shape, "need one variance per sample and dimenion"

8
GPy/core/sparse_gp.py
View file

@ -64,8 +64,8 @@ class SparseGP(GP):
self.kern.gradient += target self.kern.gradient += target
#gradients wrt Z #gradients wrt Z
self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z) self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(dL_dKmm, self.Z)
self.Z.gradient += self.kern.gradients_Z_expectations( self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_Z_expectations(
self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X) self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X)
else: else:
#gradients wrt kernel #gradients wrt kernel
@ -77,8 +77,8 @@ class SparseGP(GP):
self.kern.gradient += target self.kern.gradient += target
#gradients wrt Z #gradients wrt Z
self.Z.gradient = self.kern.gradients_X(self.grad_dict['dL_dKmm'], self.Z) self.Z.gradient[:,self.kern.active_dims] = 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.active_dims] += self.kern.gradients_X(self.grad_dict['dL_dKnm'].T, self.Z, self.X)
def _raw_predict(self, Xnew, full_cov=False): def _raw_predict(self, Xnew, full_cov=False):
""" """

2
GPy/examples/regression.py
View file

@ -468,7 +468,7 @@ def sparse_GP_regression_2D(num_samples=400, num_inducing=50, max_iters=100, opt
def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True): def uncertain_inputs_sparse_regression(max_iters=200, optimize=True, plot=True):
"""Run a 1D example of a sparse GP regression with uncertain inputs.""" """Run a 1D example of a sparse GP regression with uncertain inputs."""
fig, axes = pb.subplots(1, 2, figsize=(12, 5)) fig, axes = pb.subplots(1, 2, figsize=(12, 5), sharex=True, sharey=True)
# sample inputs and outputs # sample inputs and outputs
S = np.ones((20, 1)) S = np.ones((20, 1))

3
GPy/inference/latent_function_inference/exact_gaussian_inference.py
View file

@ -49,9 +49,6 @@ class ExactGaussianInference(object):
dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi) dL_dK = 0.5 * (tdot(alpha) - Y.shape[1] * Wi)
#TODO: does this really live here?
likelihood.update_gradients(np.diag(dL_dK))
return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK} return Posterior(woodbury_chol=LW, woodbury_vector=alpha, K=K), log_marginal, {'dL_dK':dL_dK}

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

@ -3,6 +3,7 @@
from posterior import Posterior from posterior import Posterior
from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify from ...util.linalg import jitchol, backsub_both_sides, tdot, dtrtrs, dtrtri, dpotri, dpotrs, symmetrify
from ...util import diag
from ...core.parameterization.variational import VariationalPosterior from ...core.parameterization.variational import VariationalPosterior
import numpy as np import numpy as np
from ...util.misc import param_to_array from ...util.misc import param_to_array
@ -78,9 +79,9 @@ class VarDTC(object):
num_data = Y.shape[0] num_data = Y.shape[0]
# kernel computations, using BGPLVM notation # kernel computations, using BGPLVM notation
Kmm = kern.K(Z) +np.eye(Z.shape[0]) * self.const_jitter Kmm = kern.K(Z).copy()
diag.add(Kmm, self.const_jitter)
Lm = jitchol(Kmm+np.eye(Z.shape[0])*self.const_jitter) Lm = jitchol(Kmm)
# The rather complex computations of A # The rather complex computations of A
if uncertain_inputs: if uncertain_inputs:
@ -169,7 +170,6 @@ class VarDTC(object):
Bi, _ = dpotri(LB, lower=1) Bi, _ = dpotri(LB, lower=1)
symmetrify(Bi) symmetrify(Bi)
Bi = -dpotri(LB, lower=1)[0] Bi = -dpotri(LB, lower=1)[0]
from ...util import diag
diag.add(Bi, 1) diag.add(Bi, 1)
woodbury_inv = backsub_both_sides(Lm, Bi) woodbury_inv = backsub_both_sides(Lm, Bi)
@ -238,7 +238,8 @@ class VarDTCMissingData(object):
dL_dKmm = 0 dL_dKmm = 0
log_marginal = 0 log_marginal = 0
Kmm = kern.K(Z) Kmm = kern.K(Z).copy()
diag.add(Kmm, self.const_jitter)
#factor Kmm #factor Kmm
Lm = jitchol(Kmm) Lm = jitchol(Kmm)
if uncertain_inputs: LmInv = dtrtri(Lm) if uncertain_inputs: LmInv = dtrtri(Lm)
@ -324,7 +325,6 @@ class VarDTCMissingData(object):
Bi, _ = dpotri(LB, lower=1) Bi, _ = dpotri(LB, lower=1)
symmetrify(Bi) symmetrify(Bi)
Bi = -dpotri(LB, lower=1)[0] Bi = -dpotri(LB, lower=1)[0]
from ...util import diag
diag.add(Bi, 1) diag.add(Bi, 1)
woodbury_inv_all[:, :, ind] = backsub_both_sides(Lm, Bi)[:,:,None] woodbury_inv_all[:, :, ind] = backsub_both_sides(Lm, Bi)[:,:,None]

160
GPy/kern/_src/add.py
View file

@ -4,46 +4,47 @@
import numpy as np import numpy as np
import itertools import itertools
from ...core.parameterization import Parameterized from ...core.parameterization import Parameterized
from kern import Kern from ...util.caching import Cache_this
from kern import CombinationKernel
class Add(Kern): class Add(CombinationKernel):
def __init__(self, subkerns, tensor): """
assert all([isinstance(k, Kern) for k in subkerns]) Add given list of kernels together.
if tensor: propagates gradients thorugh.
input_dim = sum([k.input_dim for k in subkerns]) """
self.input_slices = [] def __init__(self, subkerns, name='add'):
n = 0 super(Add, self).__init__(subkerns, name)
for k in subkerns:
self.input_slices.append(slice(n, n+k.input_dim))
n += k.input_dim
else:
assert all([k.input_dim == subkerns[0].input_dim for k in subkerns])
input_dim = subkerns[0].input_dim
self.input_slices = [slice(None) for k in subkerns]
super(Add, self).__init__(input_dim, 'add')
self.add_parameters(*subkerns)
@Cache_this(limit=2, force_kwargs=['which_parts'])
def K(self, X, X2=None): def K(self, X, X2=None, which_parts=None):
""" """
Compute the kernel function. Add all kernels together.
If a list of parts (of this kernel!) `which_parts` is given, only
:param X: the first set of inputs to the kernel the parts of the list are taken to compute the covariance.
:param X2: (optional) the second set of arguments to the kernel. If X2
is None, this is passed throgh to the 'part' object, which
handLes this as X2 == X.
""" """
assert X.shape[1] == self.input_dim assert X.shape[1] == self.input_dim
if X2 is None: if which_parts is None:
return sum([p.K(X[:, i_s], None) for p, i_s in zip(self._parameters_, self.input_slices)]) which_parts = self.parts
else: elif not isinstance(which_parts, (list, tuple)):
return sum([p.K(X[:, i_s], X2[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]) # if only one part is given
which_parts = [which_parts]
return reduce(np.add, (p.K(X, X2) for p in which_parts))
@Cache_this(limit=2, force_kwargs=['which_parts'])
def Kdiag(self, X, which_parts=None):
assert X.shape[1] == self.input_dim
if which_parts is None:
which_parts = self.parts
elif not isinstance(which_parts, (list, tuple)):
# if only one part is given
which_parts = [which_parts]
return reduce(np.add, (p.Kdiag(X) for p in which_parts))
def update_gradients_full(self, dL_dK, X, X2=None): def update_gradients_full(self, dL_dK, X, X2=None):
if X2 is None: [p.update_gradients_full(dL_dK, X, X2) for p in self.parts]
[p.update_gradients_full(dL_dK, X[:,i_s], X2) for p, i_s in zip(self._parameters_, self.input_slices)]
else: def update_gradients_diag(self, dL_dK, X):
[p.update_gradients_full(dL_dK, X[:,i_s], X2[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)] [p.update_gradients_diag(dL_dK, X) for p in self.parts]
def update_gradients_diag(self, dL_dKdiag, X): def update_gradients_diag(self, dL_dKdiag, X):
[p.update_gradients_diag(dL_dKdiag, X[:,i_s]) for p, i_s in zip(self._parameters_, self.input_slices)] [p.update_gradients_diag(dL_dKdiag, X[:,i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
@ -58,27 +59,19 @@ class Add(Kern):
:param X2: Observed data inputs (optional, defaults to X) :param X2: Observed data inputs (optional, defaults to X)
:type X2: np.ndarray (num_inducing x input_dim)""" :type X2: np.ndarray (num_inducing x input_dim)"""
target = np.zeros_like(X) target = np.zeros(X.shape)
if X2 is None: [target.__setitem__([Ellipsis, p.active_dims], target[:, p.active_dims]+p.gradients_X(dL_dK, X, X2)) for p in self.parts]
[np.add(target[:,i_s], p.gradients_X(dL_dK, X[:, i_s], None), target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
else:
[np.add(target[:,i_s], p.gradients_X(dL_dK, X[:, i_s], X2[:,i_s]), target[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)]
return target return target
def Kdiag(self, X):
assert X.shape[1] == self.input_dim
return sum([p.Kdiag(X[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)])
def psi0(self, Z, variational_posterior): def psi0(self, Z, variational_posterior):
return np.sum([p.psi0(Z[:, i_s], variational_posterior[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)],0) return reduce(np.add, (p.psi0(Z, variational_posterior) for p in self.parts))
def psi1(self, Z, variational_posterior): def psi1(self, Z, variational_posterior):
return np.sum([p.psi1(Z[:, i_s], variational_posterior[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)], 0) return reduce(np.add, (p.psi1(Z, variational_posterior) for p in self.parts))
def psi2(self, Z, variational_posterior): def psi2(self, Z, variational_posterior):
psi2 = np.sum([p.psi2(Z[:, i_s], variational_posterior[:, i_s]) for p, i_s in zip(self._parameters_, self.input_slices)], 0) psi2 = reduce(np.add, (p.psi2(Z, variational_posterior) for p in self.parts))
#return psi2
# compute the "cross" terms # compute the "cross" terms
from static import White, Bias from static import White, Bias
from rbf import RBF from rbf import RBF
@ -86,54 +79,52 @@ class Add(Kern):
from linear import Linear from linear import Linear
#ffrom fixed import Fixed #ffrom fixed import Fixed
for (p1, i1), (p2, i2) in itertools.combinations(itertools.izip(self._parameters_, self.input_slices), 2): for p1, p2 in itertools.combinations(self.parts, 2):
# i1, i2 = p1.active_dims, p2.active_dims
# white doesn;t combine with anything # white doesn;t combine with anything
if isinstance(p1, White) or isinstance(p2, White): if isinstance(p1, White) or isinstance(p2, White):
pass pass
# rbf X bias # rbf X bias
#elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, (RBF, RBFInv)): #elif isinstance(p1, (Bias, Fixed)) and isinstance(p2, (RBF, RBFInv)):
elif isinstance(p1, Bias) and isinstance(p2, (RBF, Linear)): elif isinstance(p1, Bias) and isinstance(p2, (RBF, Linear)):
tmp = p2.psi1(Z[:,i2], variational_posterior[:, i_s]) tmp = p2.psi1(Z, variational_posterior)
psi2 += p1.variance * (tmp[:, :, None] + tmp[:, None, :]) psi2 += p1.variance * (tmp[:, :, None] + tmp[:, None, :])
#elif isinstance(p2, (Bias, Fixed)) and isinstance(p1, (RBF, RBFInv)): #elif isinstance(p2, (Bias, Fixed)) and isinstance(p1, (RBF, RBFInv)):
elif isinstance(p2, Bias) and isinstance(p1, (RBF, Linear)): elif isinstance(p2, Bias) and isinstance(p1, (RBF, Linear)):
tmp = p1.psi1(Z[:,i1], variational_posterior[:, i_s]) tmp = p1.psi1(Z, variational_posterior)
psi2 += p2.variance * (tmp[:, :, None] + tmp[:, None, :]) psi2 += p2.variance * (tmp[:, :, None] + tmp[:, None, :])
elif isinstance(p2, (RBF, Linear)) and isinstance(p1, (RBF, Linear)):
assert np.intersect1d(p1.active_dims, p2.active_dims).size == 0, "only non overlapping kernel dimensions allowed so far"
tmp1 = p1.psi1(Z, variational_posterior)
tmp2 = p2.psi1(Z, variational_posterior)
psi2 += (tmp1[:, :, None] * tmp2[:, None, :]) + (tmp2[:, :, None] * tmp1[:, None, :])
else: else:
raise NotImplementedError, "psi2 cannot be computed for this kernel" raise NotImplementedError, "psi2 cannot be computed for this kernel"
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):
from static import White, Bias from static import White, Bias
mu, S = variational_posterior.mean, variational_posterior.variance for p1 in self.parts:
for p1, is1 in zip(self._parameters_, self.input_slices):
#compute the effective dL_dpsi1. Extra terms appear becaue of the cross terms in psi2! #compute the effective dL_dpsi1. Extra terms appear becaue of the cross terms in psi2!
eff_dL_dpsi1 = dL_dpsi1.copy() eff_dL_dpsi1 = dL_dpsi1.copy()
for p2, is2 in zip(self._parameters_, self.input_slices): for p2 in self.parts:
if p2 is p1: if p2 is p1:
continue continue
if isinstance(p2, White): if isinstance(p2, White):
continue continue
elif isinstance(p2, Bias): elif isinstance(p2, Bias):
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2.
else: else:# np.setdiff1d(p1.active_dims, ar2, assume_unique): # TODO: Careful, not correct for overlapping active_dims
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], variational_posterior[:, is1]) * 2. eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2.
p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior)
p1.update_gradients_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z[:,is1], variational_posterior[:, is1])
def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior): def gradients_Z_expectations(self, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
from static import White, Bias from static import White, Bias
target = np.zeros(Z.shape) target = np.zeros(Z.shape)
for p1, is1 in zip(self._parameters_, self.input_slices): for p1 in self.parts:
#compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2! #compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2!
eff_dL_dpsi1 = dL_dpsi1.copy() eff_dL_dpsi1 = dL_dpsi1.copy()
for p2, is2 in zip(self._parameters_, self.input_slices): for p2 in self.parts:
if p2 is p1: if p2 is p1:
continue continue
if isinstance(p2, White): if isinstance(p2, White):
@ -141,22 +132,18 @@ class Add(Kern):
elif isinstance(p2, Bias): elif isinstance(p2, Bias):
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2.
else: else:
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], variational_posterior[:, is2]) * 2. eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2.
target[:, p1.active_dims] += p1.gradients_Z_expectations(eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior)
target += p1.gradients_Z_expectations(eff_dL_dpsi1, dL_dpsi2, Z[:,is1], variational_posterior[:, is1])
return target return target
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):
from static import White, Bias from static import White, Bias
target_mu = np.zeros(variational_posterior.shape) target_mu = np.zeros(variational_posterior.shape)
target_S = np.zeros(variational_posterior.shape) target_S = np.zeros(variational_posterior.shape)
for p1, is1 in zip(self._parameters_, self.input_slices): for p1 in self._parameters_:
#compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2! #compute the effective dL_dpsi1. extra terms appear becaue of the cross terms in psi2!
eff_dL_dpsi1 = dL_dpsi1.copy() eff_dL_dpsi1 = dL_dpsi1.copy()
for p2, is2 in zip(self._parameters_, self.input_slices): for p2 in self._parameters_:
if p2 is p1: if p2 is p1:
continue continue
if isinstance(p2, White): if isinstance(p2, White):
@ -164,35 +151,20 @@ class Add(Kern):
elif isinstance(p2, Bias): elif isinstance(p2, Bias):
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2. eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.variance * 2.
else: else:
eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z[:,is2], variational_posterior[:, is2]) * 2. eff_dL_dpsi1 += dL_dpsi2.sum(1) * p2.psi1(Z, variational_posterior) * 2.
a, b = p1.gradients_qX_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z, variational_posterior)
target_mu[:, p1.active_dims] += a
a, b = p1.gradients_qX_expectations(dL_dpsi0, eff_dL_dpsi1, dL_dpsi2, Z[:,is1], variational_posterior[:, is1]) target_S[:, p1.active_dims] += b
target_mu += a
target_S += b
return target_mu, target_S return target_mu, target_S
def input_sensitivity(self):
in_sen = np.zeros((self.num_params, self.input_dim))
for i, [p, i_s] in enumerate(zip(self._parameters_, self.input_slices)):
in_sen[i, i_s] = p.input_sensitivity()
return in_sen
def _getstate(self): def _getstate(self):
""" """
Get the current state of the class, Get the current state of the class,
here just all the indices, rest can get recomputed here just all the indices, rest can get recomputed
""" """
return Parameterized._getstate(self) + [#self._parameters_, return super(Add, self)._getstate()
self.input_dim,
self.input_slices,
self._param_slices_
]
def _setstate(self, state): def _setstate(self, state):
self._param_slices_ = state.pop() super(Add, self)._setstate(state)
self.input_slices = state.pop()
self.input_dim = state.pop()
Parameterized._setstate(self, state)

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

@ -3,12 +3,19 @@
import sys import sys
import numpy as np import numpy as np
import itertools from ...core.parameterization.parameterized import Parameterized
from ...core.parameterization import Parameterized from kernel_slice_operations import KernCallsViaSlicerMeta
from ...core.parameterization.param import Param from ...util.caching import Cache_this
class Kern(Parameterized): class Kern(Parameterized):
#===========================================================================
# This adds input slice support. The rather ugly code for slicing can be
# found in kernel_slice_operations
__metaclass__ = KernCallsViaSlicerMeta
#===========================================================================
_debug=False
def __init__(self, input_dim, name, *a, **kw): def __init__(self, input_dim, name, *a, **kw):
""" """
The base class for a kernel: a positive definite function The base class for a kernel: a positive definite function
@ -20,11 +27,29 @@ class Kern(Parameterized):
Do not instantiate. Do not instantiate.
""" """
super(Kern, self).__init__(name=name, *a, **kw) super(Kern, self).__init__(name=name, *a, **kw)
self.input_dim = input_dim if isinstance(input_dim, int):
self.active_dims = np.r_[0:input_dim]
self.input_dim = input_dim
else:
self.active_dims = np.r_[input_dim]
self.input_dim = len(self.active_dims)
self._sliced_X = 0
@Cache_this(limit=10)#, ignore_args = (0,))
def _slice_X(self, X):
return X[:, self.active_dims]
def K(self, X, X2): def K(self, X, X2):
"""
Compute the kernel function.
:param X: the first set of inputs to the kernel
:param X2: (optional) the second set of arguments to the kernel. If X2
is None, this is passed throgh to the 'part' object, which
handLes this as X2 == X.
"""
raise NotImplementedError raise NotImplementedError
def Kdiag(self, Xa): def Kdiag(self, X):
raise NotImplementedError raise NotImplementedError
def psi0(self, Z, variational_posterior): def psi0(self, Z, variational_posterior):
raise NotImplementedError raise NotImplementedError
@ -34,7 +59,7 @@ class Kern(Parameterized):
raise NotImplementedError raise NotImplementedError
def gradients_X(self, dL_dK, X, X2): def gradients_X(self, dL_dK, X, X2):
raise NotImplementedError raise NotImplementedError
def gradients_X_diag(self, dL_dK, X): def gradients_X_diag(self, dL_dKdiag, X):
raise NotImplementedError raise NotImplementedError
def update_gradients_diag(self, dL_dKdiag, X): def update_gradients_diag(self, dL_dKdiag, X):
@ -44,7 +69,9 @@ class Kern(Parameterized):
def update_gradients_full(self, dL_dK, X, X2): def update_gradients_full(self, dL_dK, X, X2):
"""Set the gradients of all parameters when doing full (N) inference.""" """Set the gradients of all parameters when doing full (N) inference."""
raise NotImplementedError raise NotImplementedError
def update_gradients_diag(self, dL_dKdiag, X):
"""Set the gradients for all parameters for the derivative of the diagonal of the covariance w.r.t the kernel parameters."""
raise NotImplementedError
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):
""" """
Set the gradients of all parameters when doing inference with Set the gradients of all parameters when doing inference with
@ -99,17 +126,10 @@ class Kern(Parameterized):
""" Overloading of the '+' operator. for more control, see self.add """ """ Overloading of the '+' operator. for more control, see self.add """
return self.add(other) return self.add(other)
def add(self, other, tensor=False): def add(self, other, name='add'):
""" """
Add another kernel to this one. Add another kernel to this one.
If Tensor is False, both kernels are defined on the same _space_. then
the created kernel will have the same number of inputs as self and
other (which must be the same).
If Tensor is True, then the dimensions are stacked 'horizontally', so
that the resulting kernel has self.input_dim + other.input_dim
:param other: the other kernel to be added :param other: the other kernel to be added
:type other: GPy.kern :type other: GPy.kern
@ -117,23 +137,23 @@ 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
kernels = [] kernels = []
if not tensor and isinstance(self, Add): kernels.extend(self._parameters_) if isinstance(self, Add): kernels.extend(self._parameters_)
else: kernels.append(self) else: kernels.append(self)
if not tensor and isinstance(other, Add): kernels.extend(other._parameters_) if isinstance(other, Add): kernels.extend(other._parameters_)
else: kernels.append(other) else: kernels.append(other)
return Add(kernels, tensor) return Add(kernels, name=name)
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"""
return self.prod(other) return self.prod(other)
def __pow__(self, other): #def __pow__(self, other):
""" # """
Shortcut for tensor `prod`. # Shortcut for tensor `prod`.
""" # """
return self.prod(other, tensor=True) # return self.prod(other, tensor=True)
def prod(self, other, tensor=False, name=None): def prod(self, other, name=None):
""" """
Multiply two kernels (either on the same space, or on the tensor Multiply two kernels (either on the same space, or on the tensor
product of the input space). product of the input space).
@ -146,4 +166,42 @@ 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 prod import Prod from prod import Prod
return Prod(self, other, tensor, name) kernels = []
if isinstance(self, Prod): kernels.extend(self._parameters_)
else: kernels.append(self)
if isinstance(other, Prod): kernels.extend(other._parameters_)
else: kernels.append(other)
return Prod(self, other, name)
def _getstate(self):
"""
Get the current state of the class,
here just all the indices, rest can get recomputed
"""
return super(Kern, self)._getstate() + [
self.active_dims,
self.input_dim,
self._sliced_X]
def _setstate(self, state):
self._sliced_X = state.pop()
self.input_dim = state.pop()
self.active_dims = state.pop()
super(Kern, self)._setstate(state)
class CombinationKernel(Kern):
def __init__(self, kernels, name):
assert all([isinstance(k, Kern) for k in kernels])
input_dim = reduce(np.union1d, (x.active_dims for x in kernels))
super(CombinationKernel, self).__init__(input_dim, name)
self.add_parameters(*kernels)
@property
def parts(self):
return self._parameters_
def input_sensitivity(self):
in_sen = np.zeros((self.num_params, self.input_dim))
for i, p in enumerate(self.parts):
in_sen[i, p.active_dims] = p.input_sensitivity()
return in_sen

108
GPy/kern/_src/kernel_slice_operations.py Normal file
View file

@ -0,0 +1,108 @@
'''
Created on 11 Mar 2014
@author: maxz
'''
from ...core.parameterization.parameterized import ParametersChangedMeta
class KernCallsViaSlicerMeta(ParametersChangedMeta):
def __call__(self, *args, **kw):
instance = super(ParametersChangedMeta, self).__call__(*args, **kw)
instance.K = _slice_wrapper(instance, instance.K)
instance.Kdiag = _slice_wrapper(instance, instance.Kdiag, True)
instance.update_gradients_full = _slice_wrapper(instance, instance.update_gradients_full, False, True)
instance.update_gradients_diag = _slice_wrapper(instance, instance.update_gradients_diag, True, True)
instance.gradients_X = _slice_wrapper(instance, instance.gradients_X, False, True)
instance.gradients_X_diag = _slice_wrapper(instance, instance.gradients_X_diag, True, True)
instance.psi0 = _slice_wrapper(instance, instance.psi0, False, False)
instance.psi1 = _slice_wrapper(instance, instance.psi1, False, False)
instance.psi2 = _slice_wrapper(instance, instance.psi2, False, False)
instance.update_gradients_expectations = _slice_wrapper(instance, instance.update_gradients_expectations, psi_stat=True)
instance.gradients_Z_expectations = _slice_wrapper(instance, instance.gradients_Z_expectations, psi_stat_Z=True)
instance.gradients_qX_expectations = _slice_wrapper(instance, instance.gradients_qX_expectations, psi_stat=True)
instance.parameters_changed()
return instance
def _slice_wrapper(kern, operation, diag=False, derivative=False, psi_stat=False, psi_stat_Z=False):
"""
This method wraps the functions in kernel to make sure all kernels allways see their respective input dimension.
The different switches are:
diag: if X2 exists
derivative: if first arg is dL_dK
psi_stat: if first 3 args are dL_dpsi0..2
psi_stat_Z: if first 2 args are dL_dpsi1..2
"""
if derivative:
if diag:
def x_slice_wrapper(dL_dK, X):
X = kern._slice_X(X) if not kern._sliced_X else X
kern._sliced_X += 1
try:
ret = operation(dL_dK, X)
except:
raise
finally:
kern._sliced_X -= 1
return ret
else:
def x_slice_wrapper(dL_dK, X, X2=None):
X, X2 = kern._slice_X(X) if not kern._sliced_X else X, kern._slice_X(X2) if X2 is not None and not kern._sliced_X else X2
kern._sliced_X += 1
try:
ret = operation(dL_dK, X, X2)
except:
raise
finally:
kern._sliced_X -= 1
return ret
elif psi_stat:
def x_slice_wrapper(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior
kern._sliced_X += 1
try:
ret = operation(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior)
except:
raise
finally:
kern._sliced_X -= 1
return ret
elif psi_stat_Z:
def x_slice_wrapper(dL_dpsi1, dL_dpsi2, Z, variational_posterior):
Z, variational_posterior = kern._slice_X(Z) if not kern._sliced_X else Z, kern._slice_X(variational_posterior) if not kern._sliced_X else variational_posterior
kern._sliced_X += 1
try:
ret = operation(dL_dpsi1, dL_dpsi2, Z, variational_posterior)
except:
raise
finally:
kern._sliced_X -= 1
return ret
else:
if diag:
def x_slice_wrapper(X, *args, **kw):
X = kern._slice_X(X) if not kern._sliced_X else X
kern._sliced_X += 1
try:
ret = operation(X, *args, **kw)
except:
raise
finally:
kern._sliced_X -= 1
return ret
else:
def x_slice_wrapper(X, X2=None, *args, **kw):
X, X2 = kern._slice_X(X) if not kern._sliced_X else X, kern._slice_X(X2) if X2 is not None and not kern._sliced_X else X2
kern._sliced_X += 1
try:
ret = operation(X, X2, *args, **kw)
except: raise
finally:
kern._sliced_X -= 1
return ret
x_slice_wrapper._operation = operation
x_slice_wrapper.__name__ = ("slicer("+operation.__name__
+(","+str(bool(diag)) if diag else'')
+(','+str(bool(derivative)) if derivative else '')
+')')
x_slice_wrapper.__doc__ = "**sliced**\n" + (operation.__doc__ or "")
return x_slice_wrapper

1
GPy/kern/_src/linear.py
View file

@ -147,7 +147,6 @@ class Linear(Kern):
mu = variational_posterior.mean mu = variational_posterior.mean
S = variational_posterior.variance S = variational_posterior.variance
mu2S = np.square(mu)+S mu2S = np.square(mu)+S
_dpsi2_dvariance, _, _, _, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma) _dpsi2_dvariance, _, _, _, _ = linear_psi_comp._psi2computations(self.variances, Z, mu, S, gamma)
grad = np.einsum('n,nq,nq->q',dL_dpsi0,gamma,mu2S) + np.einsum('nm,nq,mq,nq->q',dL_dpsi1,gamma,Z,mu) +\ grad = np.einsum('n,nq,nq->q',dL_dpsi0,gamma,mu2S) + np.einsum('nm,nq,mq,nq->q',dL_dpsi1,gamma,Z,mu) +\
np.einsum('nmo,nmoq->q',dL_dpsi2,_dpsi2_dvariance) np.einsum('nmo,nmoq->q',dL_dpsi2,_dpsi2_dvariance)

51
GPy/kern/_src/prod.py
View file

@ -1,10 +1,12 @@
# Copyright (c) 2012, GPy authors (see AUTHORS.txt). # Copyright (c) 2012, GPy authors (see AUTHORS.txt).
# Licensed under the BSD 3-clause license (see LICENSE.txt) # Licensed under the BSD 3-clause license (see LICENSE.txt)
from kern import Kern
import numpy as np import numpy as np
from kern import CombinationKernel
from ...util.caching import Cache_this
import itertools
class Prod(Kern): class Prod(CombinationKernel):
""" """
Computes the product of 2 kernels Computes the product of 2 kernels
@ -15,34 +17,31 @@ class Prod(Kern):
:rtype: kernel object :rtype: kernel object
""" """
def __init__(self, k1, k2, tensor=False,name=None): def __init__(self, kernels, name='prod'):
if tensor: super(Prod, self).__init__(kernels, name)
name = k1.name + '_xx_' + k2.name if name is None else name
super(Prod, self).__init__(k1.input_dim + k2.input_dim, name)
self.slice1 = slice(0,k1.input_dim)
self.slice2 = slice(k1.input_dim,k1.input_dim+k2.input_dim)
else:
assert k1.input_dim == k2.input_dim, "Error: The input spaces of the kernels to multiply don't have the same dimension."
name = k1.name + '_x_' + k2.name if name is None else name
super(Prod, self).__init__(k1.input_dim, name)
self.slice1 = slice(0, self.input_dim)
self.slice2 = slice(0, self.input_dim)
self.k1 = k1
self.k2 = k2
self.add_parameters(self.k1, self.k2)
def K(self, X, X2=None): @Cache_this(limit=2, force_kwargs=['which_parts'])
if X2 is None: def K(self, X, X2=None, which_parts=None):
return self.k1.K(X[:,self.slice1], None) * self.k2.K(X[:,self.slice2], None) assert X.shape[1] == self.input_dim
else: if which_parts is None:
return self.k1.K(X[:,self.slice1], X2[:,self.slice1]) * self.k2.K(X[:,self.slice2], X2[:,self.slice2]) which_parts = self.parts
elif not isinstance(which_parts, (list, tuple)):
# if only one part is given
which_parts = [which_parts]
return reduce(np.multiply, (p.K(X, X2) for p in which_parts))
def Kdiag(self, X): @Cache_this(limit=2, force_kwargs=['which_parts'])
return self.k1.Kdiag(X[:,self.slice1]) * self.k2.Kdiag(X[:,self.slice2]) def Kdiag(self, X, which_parts=None):
assert X.shape[1] == self.input_dim
if which_parts is None:
which_parts = self.parts
return reduce(np.multiply, (p.Kdiag(X) for p in which_parts))
def update_gradients_full(self, dL_dK, X): def update_gradients_full(self, dL_dK, X):
self.k1.update_gradients_full(dL_dK*self.k2.K(X[:,self.slice2]), X[:,self.slice1]) for k1,k2 in itertools.combinations(self.parts, 2):
self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2]) k1._sliced_X = k1._sliced_X2 = k2._sliced_X = k2._sliced_X2 = True
k1.update_gradients_full(dL_dK*k2.K(X, X)
self.k2.update_gradients_full(dL_dK*self.k1.K(X[:,self.slice1]), X[:,self.slice2])
def gradients_X(self, dL_dK, X, X2=None): def gradients_X(self, dL_dK, X, X2=None):
target = np.zeros(X.shape) target = np.zeros(X.shape)

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

@ -19,7 +19,6 @@ class RBF(Stationary):
k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r^2 \\bigg) k(r) = \sigma^2 \exp \\bigg(- \\frac{1}{2} r^2 \\bigg)
""" """
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 = {}
@ -67,7 +66,6 @@ class RBF(Stationary):
else: else:
self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).sum() self.lengthscale.gradient = (dL_dpsi1[:,:,None]*_dpsi1_dlengthscale).sum()
#from psi2 #from psi2
self.variance.gradient += (dL_dpsi2 * _dpsi2_dvariance).sum() self.variance.gradient += (dL_dpsi2 * _dpsi2_dvariance).sum()
if self.ARD: if self.ARD:
@ -76,11 +74,14 @@ class RBF(Stationary):
self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).sum() self.lengthscale.gradient += (dL_dpsi2[:,:,:,None] * _dpsi2_dlengthscale).sum()
elif isinstance(variational_posterior, variational.NormalPosterior): elif isinstance(variational_posterior, variational.NormalPosterior):
l2 = self.lengthscale**2
l2 = self.lengthscale **2 if l2.size != self.input_dim:
l2 = l2*np.ones(self.input_dim)
#contributions from psi0: #contributions from psi0:
self.variance.gradient = np.sum(dL_dpsi0) self.variance.gradient = np.sum(dL_dpsi0)
if self._debug:
num_grad = self.lengthscale.gradient.copy()
self.lengthscale.gradient = 0. self.lengthscale.gradient = 0.
#from psi1 #from psi1
@ -92,16 +93,16 @@ class RBF(Stationary):
else: else:
self.lengthscale.gradient += dpsi1_dlength.sum() self.lengthscale.gradient += dpsi1_dlength.sum()
self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance self.variance.gradient += np.sum(dL_dpsi1 * psi1) / self.variance
#from psi2 #from psi2
S = variational_posterior.variance S = variational_posterior.variance
_, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior) _, Zdist_sq, _, mudist_sq, psi2 = self._psi2computations(Z, variational_posterior)
if not self.ARD: if not self.ARD:
self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2).sum() self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2).sum()
else: else:
self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2) self.lengthscale.gradient += self._weave_psi2_lengthscale_grads(dL_dpsi2, psi2, Zdist_sq, S, mudist_sq, l2)
if self._debug:
import ipdb;ipdb.set_trace()
self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance self.variance.gradient += 2.*np.sum(dL_dpsi2 * psi2)/self.variance
else: else:
@ -122,7 +123,6 @@ class RBF(Stationary):
return grad return grad
elif isinstance(variational_posterior, variational.NormalPosterior): elif isinstance(variational_posterior, variational.NormalPosterior):
l2 = self.lengthscale **2 l2 = self.lengthscale **2
#psi1 #psi1
@ -153,6 +153,7 @@ class RBF(Stationary):
grad_mu = (dL_dpsi1[:, :, None] * _dpsi1_dmu).sum(axis=1) grad_mu = (dL_dpsi1[:, :, None] * _dpsi1_dmu).sum(axis=1)
grad_S = (dL_dpsi1[:, :, None] * _dpsi1_dS).sum(axis=1) grad_S = (dL_dpsi1[:, :, None] * _dpsi1_dS).sum(axis=1)
grad_gamma = (dL_dpsi1[:,:,None] * _dpsi1_dgamma).sum(axis=1) grad_gamma = (dL_dpsi1[:,:,None] * _dpsi1_dgamma).sum(axis=1)
#psi2 #psi2
grad_mu += (dL_dpsi2[:, :, :, None] * _dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1) grad_mu += (dL_dpsi2[:, :, :, None] * _dpsi2_dmu).reshape(ndata,-1,self.input_dim).sum(axis=1)
grad_S += (dL_dpsi2[:, :, :, None] * _dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1) grad_S += (dL_dpsi2[:, :, :, None] * _dpsi2_dS).reshape(ndata,-1,self.input_dim).sum(axis=1)

28
GPy/kern/_src/static.py
View file

@ -89,3 +89,31 @@ class Bias(Static):
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):
self.variance.gradient = dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum() self.variance.gradient = dL_dpsi0.sum() + dL_dpsi1.sum() + 2.*self.variance*dL_dpsi2.sum()
class Fixed(Static):
def __init__(self, input_dim, covariance_matrix, variance=1., name='fixed'):
"""
:param input_dim: the number of input dimensions
:type input_dim: int
:param variance: the variance of the kernel
:type variance: float
"""
super(Bias, self).__init__(input_dim, variance, name)
self.fixed_K = covariance_matrix
def K(self, X, X2):
return self.variance * self.fixed_K
def Kdiag(self, X):
return self.variance * self.fixed_K.diag()
def update_gradients_full(self, dL_dK, X, X2=None):
self.variance.gradient = np.einsum('ij,ij', dL_dK, self.fixed_K)
def update_gradients_diag(self, dL_dKdiag, X):
self.variance.gradient = np.einsum('i,i', dL_dKdiag, self.fixed_K)
def psi2(self, Z, variational_posterior):
return np.zeros((variational_posterior.shape[0], Z.shape[0], Z.shape[0]), dtype=np.float64)
def update_gradients_expectations(self, dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, variational_posterior):
self.variance.gradient = dL_dpsi0.sum()

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

@ -85,12 +85,14 @@ class Stationary(Kern):
Compute the Euclidean 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.
""" """
#X, = self._slice_X(X)
if X2 is None: if X2 is None:
Xsq = np.sum(np.square(X),1) Xsq = np.sum(np.square(X),1)
r2 = -2.*tdot(X) + (Xsq[:,None] + Xsq[None,:]) r2 = -2.*tdot(X) + (Xsq[:,None] + Xsq[None,:])
util.diag.view(r2)[:,]= 0. # force diagnoal to be zero: sometime numerically a little negative util.diag.view(r2)[:,]= 0. # force diagnoal to be zero: sometime numerically a little negative
return np.sqrt(r2) return np.sqrt(r2)
else: else:
#X2, = self._slice_X(X2)
X1sq = np.sum(np.square(X),1) X1sq = np.sum(np.square(X),1)
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,:]))
@ -124,7 +126,6 @@ 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):
self.variance.gradient = np.einsum('ij,ij,i', self.K(X, X2), dL_dK, 1./self.variance) self.variance.gradient = np.einsum('ij,ij,i', self.K(X, X2), dL_dK, 1./self.variance)
#now the lengthscale gradient(s) #now the lengthscale gradient(s)
@ -136,7 +137,7 @@ class Stationary(Kern):
#self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3 #self.lengthscale.gradient = -((dL_dr*rinv)[:,:,None]*x_xl3).sum(0).sum(0)/self.lengthscale**3
tmp = dL_dr*self._inv_dist(X, X2) tmp = dL_dr*self._inv_dist(X, X2)
if X2 is None: X2 = X 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)]) self.lengthscale.gradient = np.array([np.einsum('ij,ij,...', tmp, np.square(self._slice_X(X)[:,q:q+1] - self._slice_X(X2)[:,q:q+1].T), -1./self.lengthscale[q]**3) for q in xrange(self.input_dim)])
else: else:
r = self._scaled_dist(X, X2) r = self._scaled_dist(X, X2)
self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale self.lengthscale.gradient = -np.sum(dL_dr*r)/self.lengthscale
@ -176,7 +177,6 @@ class Stationary(Kern):
ret = np.empty(X.shape, dtype=np.float64) ret = np.empty(X.shape, dtype=np.float64)
[np.einsum('ij,ij->i', tmp, X[:,q][:,None]-X2[:,q][None,:], out=ret[:,q]) 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
def gradients_X_diag(self, dL_dKdiag, X): def gradients_X_diag(self, dL_dKdiag, X):

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

@ -68,7 +68,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
#work out what the inputs are for plotting (1D or 2D) #work out what the inputs are for plotting (1D or 2D)
fixed_dims = np.array([i for i,v in fixed_inputs]) fixed_dims = np.array([i for i,v in fixed_inputs])
free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims) free_dims = np.setdiff1d(np.arange(model.input_dim),fixed_dims)
plots = {}
#one dimensional plotting #one dimensional plotting
if len(free_dims) == 1: if len(free_dims) == 1:
@ -89,20 +89,20 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
m, v, lower, upper = model.predict(Xgrid) m, v, lower, upper = model.predict(Xgrid)
Y = Y Y = Y
for d in which_data_ycols: for d in which_data_ycols:
gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], ax=ax, edgecol=linecol, fillcol=fillcol) plots['gpplot'] = gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], ax=ax, edgecol=linecol, fillcol=fillcol)
ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5) plots['dataplot'] = ax.plot(X[which_data_rows,free_dims], Y[which_data_rows, d], 'kx', mew=1.5)
#optionally plot some samples #optionally plot some samples
if samples: #NOTE not tested with fixed_inputs if samples: #NOTE not tested with fixed_inputs
Ysim = model.posterior_samples(Xgrid, samples) Ysim = model.posterior_samples(Xgrid, samples)
for yi in Ysim.T: for yi in Ysim.T:
ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25) plots['posterior_samples'] = ax.plot(Xnew, yi[:,None], Tango.colorsHex['darkBlue'], linewidth=0.25)
#ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs. #ax.plot(Xnew, yi[:,None], marker='x', linestyle='--',color=Tango.colorsHex['darkBlue']) #TODO apply this line for discrete outputs.
#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(), plots['xerrorbar'] = 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)
@ -118,7 +118,7 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
#Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims]
Zu = Z[:,free_dims] Zu = Z[:,free_dims]
z_height = ax.get_ylim()[0] z_height = ax.get_ylim()[0]
ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12) plots['inducing_inputs'] = ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12)
@ -143,8 +143,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
Y = Y Y = Y
for d in which_data_ycols: for d in which_data_ycols:
m_d = m[:,d].reshape(resolution, resolution).T m_d = m[:,d].reshape(resolution, resolution).T
ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet) plots['contour'] = ax.contour(x, y, m_d, levels, vmin=m.min(), vmax=m.max(), cmap=pb.cm.jet)
ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.) plots['dataplot'] = ax.scatter(X[which_data_rows, free_dims[0]], X[which_data_rows, free_dims[1]], 40, Y[which_data_rows, d], cmap=pb.cm.jet, vmin=m.min(), vmax=m.max(), linewidth=0.)
#set the limits of the plot to some sensible values #set the limits of the plot to some sensible values
ax.set_xlim(xmin[0], xmax[0]) ax.set_xlim(xmin[0], xmax[0])
@ -157,11 +157,11 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
if hasattr(model,"Z"): if hasattr(model,"Z"):
#Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims] #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims]
Zu = Z[:,free_dims] Zu = Z[:,free_dims]
ax.plot(Zu[:,free_dims[0]], Zu[:,free_dims[1]], 'wo') plots['inducing_inputs'] = ax.plot(Zu[:,free_dims[0]], Zu[:,free_dims[1]], 'wo')
else: else:
raise NotImplementedError, "Cannot define a frame with more than two input dimensions" raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
return plots
def plot_fit_f(model, *args, **kwargs): def plot_fit_f(model, *args, **kwargs):
""" """

30
GPy/testing/kernel_tests.py
View file

@ -6,7 +6,9 @@ import numpy as np
import GPy import GPy
import sys import sys
verbose = True verbose = 0
class Kern_check_model(GPy.core.Model): class Kern_check_model(GPy.core.Model):
""" """
@ -91,7 +93,7 @@ class Kern_check_dKdiag_dX(Kern_check_dK_dX):
def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False): def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verbose=False):
""" """
This function runs on kernels to check the correctness of their This function runs on kernels to check the correctness of their
implementation. It checks that the covariance function is positive definite implementation. It checks that the covariance function is positive definite
@ -210,7 +212,7 @@ def kern_test(kern, X=None, X2=None, output_ind=None, verbose=False):
class KernelTestsContinuous(unittest.TestCase): class KernelGradientTestsContinuous(unittest.TestCase):
def setUp(self): def setUp(self):
self.X = np.random.randn(100,2) self.X = np.random.randn(100,2)
self.X2 = np.random.randn(110,2) self.X2 = np.random.randn(110,2)
@ -220,16 +222,34 @@ class KernelTestsContinuous(unittest.TestCase):
def test_Matern32(self): def test_Matern32(self):
k = GPy.kern.Matern32(2) k = GPy.kern.Matern32(2)
self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
def test_Matern52(self): def test_Matern52(self):
k = GPy.kern.Matern52(2) k = GPy.kern.Matern52(2)
self.assertTrue(kern_test(k, X=self.X, X2=self.X2, verbose=verbose)) self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
#TODO: turn off grad checkingwrt X for indexed kernels liek coregionalize #TODO: turn off grad checkingwrt X for indexed kernels liek coregionalize
class KernelTestsMiscellaneous(unittest.TestCase):
def setUp(self):
N, D = 100, 10
self.X = np.linspace(-np.pi, +np.pi, N)[:,None] * np.ones(D)
self.rbf = GPy.kern.RBF(range(2))
self.linear = GPy.kern.Linear((3,5,6))
self.matern = GPy.kern.Matern32(np.array([2,4,7]))
self.sumkern = self.rbf + self.linear
self.sumkern += self.matern
self.sumkern.randomize()
def test_active_dims(self):
self.assertListEqual(self.sumkern.active_dims.tolist(), range(8))
def test_which_parts(self):
self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=[self.linear, self.matern]), self.linear.K(self.X)+self.matern.K(self.X)))
self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=[self.linear, self.rbf]), self.linear.K(self.X)+self.rbf.K(self.X)))
self.assertTrue(np.allclose(self.sumkern.K(self.X, which_parts=self.sumkern.parts[0]), self.rbf.K(self.X)))
if __name__ == "__main__": if __name__ == "__main__":
print "Running unit tests, please be (very) patient..." print "Running unit tests, please be (very) patient..."

58
GPy/testing/psi_stat_expectation_tests.py
View file

@ -12,6 +12,7 @@ import numpy
from GPy.kern import RBF from GPy.kern import RBF
from GPy.kern import Linear from GPy.kern import Linear
from copy import deepcopy from copy import deepcopy
from GPy.core.parameterization.variational import NormalPosterior
__test__ = lambda: 'deep' in sys.argv __test__ = lambda: 'deep' in sys.argv
# np.random.seed(0) # np.random.seed(0)
@ -28,53 +29,21 @@ def ard(p):
class Test(unittest.TestCase): class Test(unittest.TestCase):
input_dim = 9 input_dim = 9
num_inducing = 13 num_inducing = 13
N = 300 N = 1000
Nsamples = 1e6 Nsamples = 1e6
def setUp(self): def setUp(self):
i_s_dim_list = [2,4,3]
indices = numpy.cumsum(i_s_dim_list).tolist()
input_slices = [slice(a,b) for a,b in zip([None]+indices, indices)]
#input_slices[2] = deepcopy(input_slices[1])
input_slice_kern = GPy.kern.kern(9,
[
RBF(i_s_dim_list[0], np.random.rand(), np.random.rand(i_s_dim_list[0]), ARD=True),
RBF(i_s_dim_list[1], np.random.rand(), np.random.rand(i_s_dim_list[1]), ARD=True),
Linear(i_s_dim_list[2], np.random.rand(i_s_dim_list[2]), ARD=True)
],
input_slices = input_slices
)
self.kerns = ( self.kerns = (
# input_slice_kern, #GPy.kern.RBF([0,1,2], ARD=True)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim),
# (GPy.kern.rbf(self.input_dim, ARD=True) + #GPy.kern.RBF(self.input_dim)+GPy.kern.Bias(self.input_dim)+GPy.kern.White(self.input_dim),
# GPy.kern.linear(self.input_dim, ARD=True) + #GPy.kern.Linear(self.input_dim) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim),
# GPy.kern.bias(self.input_dim) + #GPy.kern.Linear(self.input_dim, ARD=True) + GPy.kern.Bias(self.input_dim) + GPy.kern.White(self.input_dim),
# GPy.kern.white(self.input_dim)), GPy.kern.Linear([1,3,6,7], ARD=True) + GPy.kern.RBF([0,5,8], ARD=True) + GPy.kern.White(self.input_dim),
(#GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
GPy.kern.Linear(self.input_dim, np.random.rand(self.input_dim), ARD=True)
+GPy.kern.RBF(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
# +GPy.kern.bias(self.input_dim)
# +GPy.kern.white(self.input_dim)),
),
# (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True) +
# GPy.kern.bias(self.input_dim, np.random.rand())),
# (GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
# +GPy.kern.rbf(self.input_dim, np.random.rand(), np.random.rand(self.input_dim), ARD=True)
# #+GPy.kern.bias(self.input_dim, np.random.rand())
# #+GPy.kern.white(self.input_dim, np.random.rand())),
# ),
# GPy.kern.white(self.input_dim, np.random.rand())),
# GPy.kern.rbf(self.input_dim), GPy.kern.rbf(self.input_dim, ARD=True),
# GPy.kern.linear(self.input_dim, ARD=False), GPy.kern.linear(self.input_dim, ARD=True),
# GPy.kern.linear(self.input_dim) + GPy.kern.bias(self.input_dim),
# GPy.kern.rbf(self.input_dim) + GPy.kern.bias(self.input_dim),
# GPy.kern.linear(self.input_dim) + GPy.kern.bias(self.input_dim) + GPy.kern.white(self.input_dim),
# GPy.kern.rbf(self.input_dim) + GPy.kern.bias(self.input_dim) + GPy.kern.white(self.input_dim),
# GPy.kern.bias(self.input_dim), GPy.kern.white(self.input_dim),
) )
self.q_x_mean = np.random.randn(self.input_dim) self.q_x_mean = np.random.randn(self.input_dim)[None]
self.q_x_variance = np.exp(np.random.randn(self.input_dim)) self.q_x_variance = np.exp(.5*np.random.randn(self.input_dim))[None]
self.q_x_samples = np.random.randn(self.Nsamples, self.input_dim) * np.sqrt(self.q_x_variance) + self.q_x_mean self.q_x_samples = np.random.randn(self.Nsamples, self.input_dim) * np.sqrt(self.q_x_variance) + self.q_x_mean
self.q_x = NormalPosterior(self.q_x_mean, self.q_x_variance)
self.Z = np.random.randn(self.num_inducing, self.input_dim) self.Z = np.random.randn(self.num_inducing, self.input_dim)
self.q_x_mean.shape = (1, self.input_dim) self.q_x_mean.shape = (1, self.input_dim)
self.q_x_variance.shape = (1, self.input_dim) self.q_x_variance.shape = (1, self.input_dim)
@ -114,8 +83,9 @@ class Test(unittest.TestCase):
def test_psi2(self): def test_psi2(self):
for kern in self.kerns: for kern in self.kerns:
kern.randomize()
Nsamples = int(np.floor(self.Nsamples/self.N)) Nsamples = int(np.floor(self.Nsamples/self.N))
psi2 = kern.psi2(self.Z, self.q_x_mean, self.q_x_variance) psi2 = kern.psi2(self.Z, self.q_x)
K_ = np.zeros((self.num_inducing, self.num_inducing)) K_ = np.zeros((self.num_inducing, self.num_inducing))
diffs = [] diffs = []
for i, q_x_sample_stripe in enumerate(np.array_split(self.q_x_samples, self.Nsamples / Nsamples)): for i, q_x_sample_stripe in enumerate(np.array_split(self.q_x_samples, self.Nsamples / Nsamples)):
@ -130,8 +100,8 @@ class Test(unittest.TestCase):
pylab.figure(msg) pylab.figure(msg)
pylab.plot(diffs, marker='x', mew=.2) pylab.plot(diffs, marker='x', mew=.2)
# print msg, np.allclose(psi2.squeeze(), K_, rtol=1e-1, atol=.1) # print msg, np.allclose(psi2.squeeze(), K_, rtol=1e-1, atol=.1)
self.assertTrue(np.allclose(psi2.squeeze(), K_), self.assertTrue(np.allclose(psi2.squeeze(), K_,
#rtol=1e-1, atol=.1), atol=.1, rtol=1),
msg=msg + ": not matching") msg=msg + ": not matching")
# sys.stdout.write(".") # sys.stdout.write(".")
except: except:

12
GPy/testing/psi_stat_gradient_tests.py
View file

@ -11,6 +11,7 @@ import itertools
from GPy.core import Model from GPy.core import Model
from GPy.core.parameterization.param import Param from GPy.core.parameterization.param import Param
from GPy.core.parameterization.transformations import Logexp from GPy.core.parameterization.transformations import Logexp
from GPy.core.parameterization.variational import NormalPosterior
class PsiStatModel(Model): class PsiStatModel(Model):
def __init__(self, which, X, X_variance, Z, num_inducing, kernel): def __init__(self, which, X, X_variance, Z, num_inducing, kernel):
@ -18,23 +19,24 @@ class PsiStatModel(Model):
self.which = which self.which = which
self.X = Param("X", X) self.X = Param("X", X)
self.X_variance = Param('X_variance', X_variance, Logexp()) self.X_variance = Param('X_variance', X_variance, Logexp())
self.q = NormalPosterior(self.X, self.X_variance)
self.Z = Param("Z", Z) self.Z = Param("Z", Z)
self.N, self.input_dim = X.shape self.N, self.input_dim = X.shape
self.num_inducing, input_dim = Z.shape self.num_inducing, input_dim = Z.shape
assert self.input_dim == input_dim, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape) assert self.input_dim == input_dim, "shape missmatch: Z:{!s} X:{!s}".format(Z.shape, X.shape)
self.kern = kernel self.kern = kernel
self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance) self.psi_ = self.kern.__getattribute__(self.which)(self.Z, self.q)
self.add_parameters(self.X, self.X_variance, self.Z, self.kern) self.add_parameters(self.q, self.Z, self.kern)
def log_likelihood(self): def log_likelihood(self):
return self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance).sum() return self.kern.__getattribute__(self.which)(self.Z, self.X, self.X_variance).sum()
def parameters_changed(self): def parameters_changed(self):
psimu, psiS = self.kern.__getattribute__("d" + self.which + "_dmuS")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance) psimu, psiS = self.kern.__getattribute__("d" + self.which + "_dmuS")(numpy.ones_like(self.psi_), self.Z, self.q)
self.X.gradient = psimu self.X.gradient = psimu
self.X_variance.gradient = psiS self.X_variance.gradient = psiS
#psimu, psiS = numpy.ones(self.N * self.input_dim), numpy.ones(self.N * self.input_dim) #psimu, psiS = numpy.ones(self.N * self.input_dim), numpy.ones(self.N * self.input_dim)
try: psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.X, self.X_variance) try: psiZ = self.kern.__getattribute__("d" + self.which + "_dZ")(numpy.ones_like(self.psi_), self.Z, self.q)
except AttributeError: psiZ = numpy.zeros_like(self.Z) except AttributeError: psiZ = numpy.zeros_like(self.Z)
self.Z.gradient = psiZ self.Z.gradient = psiZ
#psiZ = numpy.ones(self.num_inducing * self.input_dim) #psiZ = numpy.ones(self.num_inducing * self.input_dim)
@ -176,6 +178,6 @@ if __name__ == "__main__":
+GPy.kern.White(input_dim) +GPy.kern.White(input_dim)
) )
) )
m2.ensure_default_constraints() #m2.ensure_default_constraints()
else: else:
unittest.main() unittest.main()

37
GPy/util/caching.py
View file

@ -9,24 +9,27 @@ class Cacher(object):
""" """
def __init__(self, operation, limit=5, ignore_args=()): def __init__(self, operation, limit=5, ignore_args=(), force_kwargs=()):
self.limit = int(limit) self.limit = int(limit)
self.ignore_args = ignore_args self.ignore_args = ignore_args
self.force_kwargs = force_kwargs
self.operation=operation self.operation=operation
self.cached_inputs = [] self.cached_inputs = []
self.cached_outputs = [] self.cached_outputs = []
self.inputs_changed = [] self.inputs_changed = []
def __call__(self, *args): def __call__(self, *args, **kw):
""" """
A wrapper function for self.operation, A wrapper function for self.operation,
""" """
#ensure that specified arguments are ignored #ensure that specified arguments are ignored
items = sorted(kw.items(), key=lambda x: x[0])
oa_all = args + tuple(a for _,a in items)
if len(self.ignore_args) != 0: if len(self.ignore_args) != 0:
oa = [a for i,a in enumerate(args) if i not in self.ignore_args] oa = [a for i,a in itertools.chain(enumerate(args), items) if i not in self.ignore_args and i not in self.force_kwargs]
else: else:
oa = args oa = oa_all
# this makes sure we only add an observer once, and that None can be in args # this makes sure we only add an observer once, and that None can be in args
observable_args = [] observable_args = []
@ -37,8 +40,13 @@ class Cacher(object):
#make sure that all the found argument really are observable: #make sure that all the found argument really are observable:
#otherswise don't cache anything, pass args straight though #otherswise don't cache anything, pass args straight though
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, **kw)
if len(self.force_kwargs) != 0:
# check if there are force args, which force reloading
for k in self.force_kwargs:
if k in kw and kw[k] is not None:
return self.operation(*args, **kw)
# TODO: WARNING !!! Cache OFFSWITCH !!! WARNING # TODO: WARNING !!! Cache OFFSWITCH !!! WARNING
# return self.operation(*args) # return self.operation(*args)
@ -48,7 +56,7 @@ class Cacher(object):
i = state.index(True) i = state.index(True)
if self.inputs_changed[i]: if self.inputs_changed[i]:
#(elements of) the args have changed since we last computed: update #(elements of) the args have changed since we last computed: update
self.cached_outputs[i] = self.operation(*args) self.cached_outputs[i] = self.operation(*args, **kw)
self.inputs_changed[i] = False self.inputs_changed[i] = False
return self.cached_outputs[i] return self.cached_outputs[i]
else: else:
@ -62,11 +70,11 @@ class Cacher(object):
self.cached_outputs.pop(0) self.cached_outputs.pop(0)
#compute #compute
self.cached_inputs.append(args) self.cached_inputs.append(oa_all)
self.cached_outputs.append(self.operation(*args)) self.cached_outputs.append(self.operation(*args, **kw))
self.inputs_changed.append(False) self.inputs_changed.append(False)
[a.add_observer(self, self.on_cache_changed) for a in observable_args] [a.add_observer(self, self.on_cache_changed) for a in observable_args]
return self.cached_outputs[-1]#Max says return. return self.cached_outputs[-1]#return
def on_cache_changed(self, arg): def on_cache_changed(self, arg):
""" """
@ -90,15 +98,16 @@ class Cache_this(object):
""" """
A decorator which can be applied to bound methods in order to cache them A decorator which can be applied to bound methods in order to cache them
""" """
def __init__(self, limit=5, ignore_args=()): def __init__(self, limit=5, ignore_args=(), force_kwargs=()):
self.limit = limit self.limit = limit
self.ignore_args = ignore_args self.ignore_args = ignore_args
self.force_args = force_kwargs
self.c = None self.c = None
def __call__(self, f): def __call__(self, f):
def f_wrap(*args): def f_wrap(*args, **kw):
if self.c is None: if self.c is None:
self.c = Cacher(f, self.limit, ignore_args=self.ignore_args) self.c = Cacher(f, self.limit, ignore_args=self.ignore_args, force_kwargs=self.force_args)
return self.c(*args) return self.c(*args, **kw)
f_wrap._cacher = self f_wrap._cacher = self
f_wrap.__doc__ = "**cached**\n\n" + (f.__doc__ or "") f_wrap.__doc__ = "**cached**" + (f.__doc__ or "")
return f_wrap return f_wrap