diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py index 29e845ad..4513ddac 100644 --- a/GPy/core/gp_base.py +++ b/GPy/core/gp_base.py @@ -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) diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py index 01ae1806..48a7d523 100644 --- a/GPy/core/sparse_gp.py +++ b/GPy/core/sparse_gp.py @@ -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 diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py index d12643f6..df7b92b8 100644 --- a/GPy/examples/regression.py +++ b/GPy/examples/regression.py @@ -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): diff --git a/GPy/kern/constructors.py b/GPy/kern/constructors.py index a90e6ef8..57c264ef 100644 --- a/GPy/kern/constructors.py +++ b/GPy/kern/constructors.py @@ -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 diff --git a/GPy/kern/parts/coregionalise.py b/GPy/kern/parts/coregionalise.py index 9a1c31c6..66e14052 100644 --- a/GPy/kern/parts/coregionalise.py +++ b/GPy/kern/parts/coregionalise.py @@ -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): diff --git a/GPy/models/gp_multioutput_regression.py b/GPy/models/gp_multioutput_regression.py index c0a5b557..d51f3bae 100644 --- a/GPy/models/gp_multioutput_regression.py +++ b/GPy/models/gp_multioutput_regression.py @@ -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) diff --git a/GPy/models/sparse_gp_multioutput_regression.py b/GPy/models/sparse_gp_multioutput_regression.py index 62712e15..041204b6 100644 --- a/GPy/models/sparse_gp_multioutput_regression.py +++ b/GPy/models/sparse_gp_multioutput_regression.py @@ -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) diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py index 65a8da77..49445276 100644 --- a/GPy/testing/kernel_tests.py +++ b/GPy/testing/kernel_tests.py @@ -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() diff --git a/GPy/testing/unit_tests.py b/GPy/testing/unit_tests.py index fd6e8911..6bb624df 100644 --- a/GPy/testing/unit_tests.py +++ b/GPy/testing/unit_tests.py @@ -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() diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py index 9b21f8f5..c205548b 100644 --- a/GPy/util/__init__.py +++ b/GPy/util/__init__.py @@ -14,4 +14,4 @@ import visualize import decorators import classification import latent_space_visualizations -import multioutput +#import multioutput