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

linK2_functions2 merged

Browse source
This commit is contained in:
Ricardo 2013-09-13 18:09:59 +01:00
parent f8c9e6b982
commit 1bc9374717
10 changed files with 113 additions and 75 deletions
Download patch file Download diff file Expand all files Collapse all files

44
GPy/core/gp_base.py
View file

@ -58,30 +58,30 @@ class GPBase(Model):
Model.setstate(self, state)
def plot_f(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, full_cov=False, fignum=None, ax=None,output=None):
"""
Plot the GP's view of the world, where the data is normalized and the
- In one dimension, the function is plotted with a shaded region identifying two standard deviations.
- In two dimsensions, a contour-plot shows the mean predicted function
- Not implemented in higher dimensions
"""
Plot the GP's view of the world, where the data is normalized and the
- In one dimension, the function is plotted with a shaded region identifying two standard deviations.
- In two dimsensions, a contour-plot shows the mean predicted function
- Not implemented in higher dimensions
:param samples: the number of a posteriori samples to plot
:param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
:param which_data: which if the training data to plot (default all)
:type which_data: 'all' or a slice object to slice self.X, self.Y
:param which_parts: which of the kernel functions to plot (additively)
:type which_parts: 'all', or list of bools
:param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
:type resolution: int
:param full_cov:
:type full_cov: bool
:param fignum: figure to plot on.
:type fignum: figure number
:param ax: axes to plot on.
:type ax: axes handle
:param samples: the number of a posteriori samples to plot
:param plot_limits: The limits of the plot. If 1D [xmin,xmax], if 2D [[xmin,ymin],[xmax,ymax]]. Defaluts to data limits
:param which_data: which if the training data to plot (default all)
:type which_data: 'all' or a slice object to slice self.X, self.Y
:param which_parts: which of the kernel functions to plot (additively)
:type which_parts: 'all', or list of bools
:param resolution: the number of intervals to sample the GP on. Defaults to 200 in 1D and 50 (a 50x50 grid) in 2D
:type resolution: int
:param full_cov:
:type full_cov: bool
:param fignum: figure to plot on.
:type fignum: figure number
:param ax: axes to plot on.
:type ax: axes handle
:param output: which output to plot (for multiple output models only)
:type output: integer (first output is 0)
"""
:param output: which output to plot (for multiple output models only)
:type output: integer (first output is 0)
"""
if which_data == 'all':
which_data = slice(None)

8
GPy/core/sparse_gp.py
View file

@ -165,7 +165,7 @@ class SparseGP(GPBase):
raise NotImplementedError, "heteroscedatic derivates with uncertain inputs not implemented"
else:
Lmi_psi1, nil = dtrtrs(self.Lm, np.asfortranarray(self.psi1.T), lower=1, trans=0)
Lmi_psi1, nil = dtrtrs(self._Lm, np.asfortranarray(self.psi1.T), lower=1, trans=0)
_LBi_Lmi_psi1, _ = dtrtrs(self.LB, np.asfortranarray(Lmi_psi1), lower=1, trans=0)
_Bi_Lmi_psi1, _ = dtrtrs(self.LB.T, np.asfortranarray(_LBi_Lmi_psi1), lower=1, trans=0)
@ -427,13 +427,13 @@ class SparseGP(GPBase):
"""
Bi, _ = dpotri(self.LB, lower=0) # WTH? this lower switch should be 1, but that doesn't work!
symmetrify(Bi)
Kmmi_LmiBLmi = backsub_both_sides(self.Lm, np.eye(self.num_inducing) - Bi)
Kmmi_LmiBLmi = backsub_both_sides(self._Lm, np.eye(self.num_inducing) - Bi)
if self.Cpsi1V is None:
psi1V = np.dot(self.psi1.T,self.likelihood.V)
tmp, _ = dtrtrs(self.Lm, np.asfortranarray(psi1V), lower=1, trans=0)
tmp, _ = dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
tmp, _ = dpotrs(self.LB, tmp, lower=1)
self.Cpsi1V, _ = dtrtrs(self.Lm, tmp, lower=1, trans=1)
self.Cpsi1V, _ = dtrtrs(self._Lm, tmp, lower=1, trans=1)
assert hasattr(self,'multioutput')
index = np.ones_like(_Xnew)*output

66
GPy/examples/regression.py
View file

@ -46,32 +46,23 @@ def coregionalisation_toy(max_iters=100):
"""
X1 = np.random.rand(50, 1) * 8
X2 = np.random.rand(30, 1) * 5
index = np.vstack((np.zeros_like(X1), np.ones_like(X2)))
X = np.hstack((np.vstack((X1, X2)), index))
X = np.vstack((X1, X2))
Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05
Y = np.vstack((Y1, Y2))
k1 = GPy.kern.rbf(1)
k2 = GPy.kern.coregionalise(2, 2)
k = k1**k2 #k1.prod(k2, tensor=True)
m = GPy.models.GPRegression(X, Y, kernel=k)
m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
m.constrain_fixed('.*rbf_var', 1.)
# m.constrain_positive('kappa')
m.optimize(max_iters=max_iters)
pb.figure()
Xtest1 = np.hstack((np.linspace(0, 9, 100)[:, None], np.zeros((100, 1))))
Xtest2 = np.hstack((np.linspace(0, 9, 100)[:, None], np.ones((100, 1))))
mean, var, low, up = m.predict(Xtest1)
GPy.util.plot.gpplot(Xtest1[:, 0], mean, low, up)
mean, var, low, up = m.predict(Xtest2)
GPy.util.plot.gpplot(Xtest2[:, 0], mean, low, up)
pb.plot(X1[:, 0], Y1[:, 0], 'rx', mew=2)
pb.plot(X2[:, 0], Y2[:, 0], 'gx', mew=2)
fig, axes = pb.subplots(2,1)
m.plot(output=0,ax=axes[0])
m.plot(output=1,ax=axes[1])
axes[0].set_title('Output 0')
axes[1].set_title('Output 1')
return m
def coregionalisation_sparse(max_iters=100):
"""
A simple demonstration of coregionalisation on two sinusoidal functions using sparse approximations.
@ -86,30 +77,39 @@ def coregionalisation_sparse(max_iters=100):
num_inducing = 40
Z = np.hstack((np.random.rand(num_inducing, 1) * 8, np.random.randint(0, 2, num_inducing)[:, None]))
Z = np.hstack((np.random.rand(num_inducing, 1) * 8, np.random.randint(0, 2, num_inducing)[:, None]))
k1 = GPy.kern.rbf(1)
k2 = GPy.kern.coregionalise(2, 2)
k = k1**k2 #.prod(k2, tensor=True) # + GPy.kern.white(2,0.001)
m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z)
m = GPy.models.SparseGPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1],num_inducing=20)
#k2 = GPy.kern.coregionalise(2, 2)
#k = k1**k2 #.prod(k2, tensor=True) # + GPy.kern.white(2,0.001)
#m = GPy.models.SparseGPRegression(X, Y, kernel=k, Z=Z)
m.constrain_fixed('.*rbf_var', 1.)
m.constrain_fixed('iip')
m.constrain_bounded('noise_variance', 1e-3, 1e-1)
#m.constrain_fixed('iip')
#m.constrain_bounded('noise_variance', 1e-3, 1e-1)
# m.optimize_restarts(5, robust=True, messages=1, max_iters=max_iters, optimizer='bfgs')
m.optimize(max_iters=max_iters)
fig, axes = pb.subplots(2,1)
m.plot(output=0,ax=axes[0])
m.plot(output=1,ax=axes[1])
axes[0].set_title('Output 0')
axes[1].set_title('Output 1')
# plotting:
pb.figure()
Xtest1 = np.hstack((np.linspace(0, 9, 100)[:, None], np.zeros((100, 1))))
Xtest2 = np.hstack((np.linspace(0, 9, 100)[:, None], np.ones((100, 1))))
mean, var, low, up = m.predict(Xtest1)
GPy.util.plot.gpplot(Xtest1[:, 0], mean, low, up)
mean, var, low, up = m.predict(Xtest2)
GPy.util.plot.gpplot(Xtest2[:, 0], mean, low, up)
pb.plot(X1[:, 0], Y1[:, 0], 'rx', mew=2)
pb.plot(X2[:, 0], Y2[:, 0], 'gx', mew=2)
y = pb.ylim()[0]
pb.plot(Z[:, 0][Z[:, 1] == 0], np.zeros(np.sum(Z[:, 1] == 0)) + y, 'r|', mew=2)
pb.plot(Z[:, 0][Z[:, 1] == 1], np.zeros(np.sum(Z[:, 1] == 1)) + y, 'g|', mew=2)
#pb.figure()
#Xtest1 = np.hstack((np.linspace(0, 9, 100)[:, None], np.zeros((100, 1))))
#Xtest2 = np.hstack((np.linspace(0, 9, 100)[:, None], np.ones((100, 1))))
#mean, var, low, up = m.predict(Xtest1)
#GPy.util.plot.gpplot(Xtest1[:, 0], mean, low, up)
#mean, var, low, up = m.predict(Xtest2)
#GPy.util.plot.gpplot(Xtest2[:, 0], mean, low, up)
#pb.plot(X1[:, 0], Y1[:, 0], 'rx', mew=2)
#pb.plot(X2[:, 0], Y2[:, 0], 'gx', mew=2)
#y = pb.ylim()[0]
#pb.plot(Z[:, 0][Z[:, 1] == 0], np.zeros(np.sum(Z[:, 1] == 0)) + y, 'r|', mew=2)
#pb.plot(Z[:, 0][Z[:, 1] == 1], np.zeros(np.sum(Z[:, 1] == 1)) + y, 'g|', mew=2)
return m
def epomeo_gpx(max_iters=100):

30
GPy/kern/constructors.py
View file

@ -340,7 +340,7 @@ def symmetric(k):
k_.parts = [symmetric.Symmetric(p) for p in k.parts]
return k_
def coregionalise(num_outpus,W_columns=1, W=None, kappa=None):
def coregionalise(num_outputs,W_columns=1, W=None, kappa=None):
"""
Coregionlization matrix B, of the form:
.. math::
@ -422,3 +422,31 @@ def hierarchical(k):
# assert (sl.start is None) and (sl.stop is None), "cannot adjust input slices! (TODO)"
_parts = [parts.hierarchical.Hierarchical(k.parts)]
return kern(k.input_dim+len(k.parts),_parts)
def build_lcm(input_dim, num_outputs, kernel_list = [], W_columns=1,W=None,kappa=None):
"""
Builds a kernel of a linear coregionalization model
:input_dim: Input dimensionality
:num_outputs: Number of outputs
:kernel_list: List of coregionalized kernels, each element in the list will be multiplied by a different corregionalization matrix
:type kernel_list: list of GPy kernels
:param W_columns: number tuples of the corregionalization parameters 'coregion_W'
:type W_columns: integer
..Note the kernels dimensionality is overwritten to fit input_dim
"""
for k in kernel_list:
if k.input_dim <> input_dim:
k.input_dim = input_dim
warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
k_coreg = coregionalise(num_outputs,W_columns,W,kappa)
kernel = kernel_list[0]**k_coreg.copy()
for k in kernel_list[1:]:
k_coreg = coregionalise(num_outputs,W_columns,W,kappa)
kernel += k**k_coreg.copy()
return kernel

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

@ -38,16 +38,16 @@ class Coregionalise(Kernpart):
self.num_outputs = num_outputs
self.W_columns = W_columns
if W is None:
self.W = 0.5*np.random.randn(self.output_dim,self.rank)/np.sqrt(self.rank)
self.W = 0.5*np.random.randn(self.num_outputs,self.W_columns)/np.sqrt(self.W_columns)
else:
assert W.shape==(self.output_dim,self.rank)
assert W.shape==(self.num_outputs,self.W_columns)
self.W = W
if kappa is None:
kappa = 0.5*np.ones(self.output_dim)
kappa = 0.5*np.ones(self.num_outputs)
else:
assert kappa.shape==(self.output_dim,)
assert kappa.shape==(self.num_outputs,)
self.kappa = kappa
self.num_params = self.output_dim*(self.rank + 1)
self.num_params = self.num_outputs*(self.W_columns + 1)
self._set_params(np.hstack([self.W.flatten(),self.kappa]))
def _get_params(self):

6
GPy/models/gp_multioutput_regression.py
View file

@ -6,7 +6,7 @@ import numpy as np
from ..core import GP
from .. import likelihoods
from .. import kern
from ..util import multioutput
#from ..util import multioutput
class GPMultioutputRegression(GP):
"""
@ -51,8 +51,8 @@ class GPMultioutputRegression(GP):
#Coregionalization kernel definition
if kernel_list is None:
kernel_list = [[kern.rbf(original_dim)],[]]
mkernel = multioutput.build_lcm(input_dim=original_dim, num_outputs=self.num_outputs, CK = kernel_list[0], NC = kernel_list[1], W_columns=W_columns)
kernel_list = [kern.rbf(original_dim)]
mkernel = kern.build_lcm(input_dim=original_dim, num_outputs=self.num_outputs, kernel_list = kernel_list, W_columns=W_columns)
self.multioutput = True
GP.__init__(self, X, likelihood, mkernel, normalize_X=normalize_X)

4
GPy/models/sparse_gp_multioutput_regression.py
View file

@ -71,8 +71,8 @@ class SparseGPMultioutputRegression(SparseGP):
#Coregionalization kernel definition
if kernel_list is None:
kernel_list = [[kern.rbf(original_dim)],[]]
mkernel = multioutput.build_lcm(input_dim=original_dim, num_outputs=self.num_outputs, CK = kernel_list[0], NC = kernel_list[1], W_columns=W_columns)
kernel_list = [kern.rbf(original_dim)]
mkernel = kern.build_lcm(input_dim=original_dim, num_outputs=self.num_outputs, kernel_list = kernel_list, W_columns=W_columns)
self.multioutput = True
SparseGP.__init__(self, X, likelihood, mkernel, Z=Z, normalize_X=normalize_X)

3
GPy/testing/kernel_tests.py
View file

@ -4,7 +4,6 @@
import unittest
import numpy as np
import GPy
class KernelTests(unittest.TestCase):
@ -12,7 +11,6 @@ class KernelTests(unittest.TestCase):
K = GPy.kern.rbf(5, ARD=True)
K.tie_params('.*[01]')
K.constrain_fixed('2')
X = np.random.rand(5,5)
Y = np.ones((5,1))
m = GPy.models.GPRegression(X,Y,K)
@ -68,7 +66,6 @@ class KernelTests(unittest.TestCase):
self.assertTrue(m.checkgrad())
if __name__ == "__main__":
print "Running unit tests, please be (very) patient..."
unittest.main()

15
GPy/testing/unit_tests.py
View file

@ -5,7 +5,6 @@
import unittest
import numpy as np
import GPy
from GPy.likelihoods.likelihood_functions import Binomial
class GradientTests(unittest.TestCase):
def setUp(self):
@ -226,6 +225,20 @@ class GradientTests(unittest.TestCase):
m.update_likelihood_approximation()
self.assertTrue(m.checkgrad())
def multioutput_regression_1D(self):
X1 = np.random.rand(50, 1) * 8
X2 = np.random.rand(30, 1) * 5
X = np.vstack((X1, X2))
Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
Y2 = -np.sin(X2) + np.random.randn(*X2.shape) * 0.05
Y = np.vstack((Y1, Y2))
k1 = GPy.kern.rbf(1)
m = GPy.models.GPMultioutputRegression(X_list=[X1,X2],Y_list=[Y1,Y2],kernel_list=[k1])
m.constrain_fixed('.*rbf_var', 1.)
self.assertTrue(m.checkgrad())
if __name__ == "__main__":
print "Running unit tests, please be (very) patient..."
unittest.main()

2
GPy/util/__init__.py
View file

@ -14,4 +14,4 @@ import visualize
import decorators
import classification
import latent_space_visualizations
import multioutput
#import multioutput