diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index bd31edd8..fb5fe789 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -173,7 +173,7 @@ class GP(GPBase):
            This is to allow for different normalizations of the output dimensions.
 
         """
-        assert isinstance(self.likelihood,EP_Mixed_Noise)
+        assert hasattr(self,'multioutput')
         index = np.ones_like(Xnew)*output
         Xnew = np.hstack((Xnew,index))
 
@@ -182,7 +182,10 @@ class GP(GPBase):
         mu, var = self._raw_predict(Xnew, full_cov=full_cov, which_parts=which_parts)
 
         # now push through likelihood
-        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
+        if isinstance(self.likelihood,EP_Mixed_Noise):
+            mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
+        else:
+            mean, var, _025pm, _975pm = self.likelihood_list[output].predictive_values(mu, var, full_cov)
         return mean, var, _025pm, _975pm
 
     def _raw_predict_single_output(self, _Xnew, output=0, which_parts='all', full_cov=False,stop=False):
diff --git a/GPy/core/gp_base.py b/GPy/core/gp_base.py
index 8f63e256..2820a447 100644
--- a/GPy/core/gp_base.py
+++ b/GPy/core/gp_base.py
@@ -106,13 +106,16 @@ class GPBase(Model):
                 gpplot(Xnew, m, m - 2 * np.sqrt(np.diag(v)[:, None]), m + 2 * np.sqrt(np.diag(v))[:, None, ], axes=ax)
                 for i in range(samples):
                     ax.plot(Xnew, Ysim[i, :], Tango.colorsHex['darkBlue'], linewidth=0.25)
-            #ax.plot(self.X[which_data], self.likelihood.Y[which_data], 'kx', mew=1.5)
-            #ax.plot(Xu[which_data], self.likelihood.Y[self.likelihood.index==output][:,None], 'kx', mew=1.5)
             ax.set_xlim(xmin, xmax)
             ymin, ymax = min(np.append(self.likelihood.Y, m - 2 * np.sqrt(np.diag(v)[:, None]))), max(np.append(self.likelihood.Y, m + 2 * np.sqrt(np.diag(v)[:, None])))
             ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin)
             ax.set_ylim(ymin, ymax)
 
+            if hasattr(self,'Z'):
+                Zu = self.Z[self.Z[:,-1]==output,:]
+                Zu = self.Z * self._Xscale + self._Xoffset
+                Zu = self.Z[self.Z[:,-1]==output ,0:1] #??
+                ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12)
 
         else:
             raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
@@ -120,7 +123,7 @@ class GPBase(Model):
     def plot(self, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, samples=0, fignum=None, ax=None, output=None):
         """
         TODO: Docstrings!
-        
+
         :param levels: for 2D plotting, the number of contour levels to use
         is ax is None, create a new figure
         """
@@ -132,7 +135,7 @@ class GPBase(Model):
             fig = pb.figure(num=fignum)
             ax = fig.add_subplot(111)
 
-        if self.X.shape[1] == 1 and not isinstance(self.likelihood,EP_Mixed_Noise):
+        if self.X.shape[1] == 1 and not hasattr(self,'multioutput'):
 
             Xu = self.X * self._Xscale + self._Xoffset # NOTE self.X are the normalized values now
 
@@ -146,7 +149,7 @@ class GPBase(Model):
             ax.set_xlim(xmin, xmax)
             ax.set_ylim(ymin, ymax)
 
-        elif self.X.shape[1] == 2 and not isinstance(self.likelihood,EP_Mixed_Noise): # FIXME
+        elif self.X.shape[1] == 2 and not hasattr(self,'multioutput'):
             resolution = resolution or 50
             Xnew, _, _, xmin, xmax = x_frame2D(self.X, plot_limits, resolution)
             x, y = np.linspace(xmin[0], xmax[0], resolution), np.linspace(xmin[1], xmax[1], resolution)
@@ -158,17 +161,23 @@ class GPBase(Model):
             ax.set_xlim(xmin[0], xmax[0])
             ax.set_ylim(xmin[1], xmax[1])
 
-        elif self.X.shape[1] == 2 and isinstance(self.likelihood,EP_Mixed_Noise):
-            Xu = self.X[self.X[:,-1]==output,:]
+        elif self.X.shape[1] == 2 and hasattr(self,'multioutput'):
+            Xu = self.X[self.X[:,-1]==output,:] #keep the output of interest
             Xu = self.X * self._Xscale + self._Xoffset
-            Xu = self.X[self.X[:,-1]==output ,0:1]
+            Xu = self.X[self.X[:,-1]==output ,0:1] #get rid of the index column
 
             Xnew, xmin, xmax = x_frame1D(Xu, plot_limits=plot_limits)
+
             m, _, lower, upper = self.predict_single_output(Xnew, which_parts=which_parts,output=output)
+            #if not isinstance(self.likelihood,EP_Mixed_Noise):
+            #    m, _, lower, upper = self.predict(np.hstack([Xnew,np.repeat(output,Xnew.size)[:,None]]), which_parts=which_parts)
+            #else:
+            #    m, _, lower, upper = self.predict_single_output(Xnew, which_parts=which_parts,output=output)
+
             for d in range(m.shape[1]):
                 gpplot(Xnew, m[:, d], lower[:, d], upper[:, d], axes=ax)
-                #ax.plot(Xu[which_data], self.likelihood.data[which_data, d], 'kx', mew=1.5)
-                ax.plot(Xu[which_data], self.likelihood.data[self.likelihood.index==output][:,None], 'kx', mew=1.5)
+                #ax.plot(Xu[which_data], self.likelihood.data[self.likelihood.index==output][:,None], 'kx', mew=1.5)
+                ax.plot(Xu[which_data], self.likelihood_list[output].data, 'kx', mew=1.5)
             ymin, ymax = min(np.append(self.likelihood.data, lower)), max(np.append(self.likelihood.data, upper))
             ymin, ymax = ymin - 0.1 * (ymax - ymin), ymax + 0.1 * (ymax - ymin)
             ax.set_xlim(xmin, xmax)
diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py
index 04762522..6d9761c4 100644
--- a/GPy/core/sparse_gp.py
+++ b/GPy/core/sparse_gp.py
@@ -293,7 +293,7 @@ class SparseGP(GPBase):
 
         return mean, var, _025pm, _975pm
 
-    def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None):
+    def plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, fignum=None, ax=None, output=None):
         if ax is None:
             fig = pb.figure(num=fignum)
             ax = fig.add_subplot(111)
@@ -301,8 +301,8 @@ class SparseGP(GPBase):
         if which_data is 'all':
             which_data = slice(None)
 
-        GPBase.plot(self, samples=0, plot_limits=None, which_data='all', which_parts='all', resolution=None, levels=20, ax=ax)
-        if self.X.shape[1] == 1:
+        GPBase.plot(self, samples=0, plot_limits=plot_limits, which_data='all', which_parts='all', resolution=None, levels=20, ax=ax, output=output)
+        if self.X.shape[1] == 1 and not hasattr(self,'multioutput'):
             if self.has_uncertain_inputs:
                 Xu = self.X * self._Xscale + self._Xoffset  # NOTE self.X are the normalized values now
                 ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0],
@@ -311,10 +311,31 @@ class SparseGP(GPBase):
             Zu = self.Z * self._Xscale + self._Xoffset
             ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12)
 
-        elif self.X.shape[1] == 2:
+        elif self.X.shape[1] == 2 and not hasattr(self,'multioutput'):
             Zu = self.Z * self._Xscale + self._Xoffset
             ax.plot(Zu[:, 0], Zu[:, 1], 'wo')
 
+        elif self.X.shape[1] == 2 and hasattr(self,'multioutput'):
+            Xu = self.X[self.X[:,-1]==output,:]
+            if self.has_uncertain_inputs:
+                Xu = self.X * self._Xscale + self._Xoffset  # NOTE self.X are the normalized values now
+
+                Xu = self.X[self.X[:,-1]==output ,0:1] #??
+
+                ax.errorbar(Xu[which_data, 0], self.likelihood.data[which_data, 0],
+                            xerr=2 * np.sqrt(self.X_variance[which_data, 0]),
+                            ecolor='k', fmt=None, elinewidth=.5, alpha=.5)
+
+            Zu = self.Z[self.Z[:,-1]==output,:]
+            Zu = self.Z * self._Xscale + self._Xoffset
+            Zu = self.Z[self.Z[:,-1]==output ,0:1] #??
+            ax.plot(Zu, np.zeros_like(Zu) + ax.get_ylim()[0], 'r|', mew=1.5, markersize=12)
+            #ax.set_ylim(ax.get_ylim()[0],)
+
+        else:
+            raise NotImplementedError, "Cannot define a frame with more than two input dimensions"
+
+
 
     def predict_single_output(self, Xnew, output=0, which_parts='all', full_cov=False):
         """
@@ -336,7 +357,7 @@ class SparseGP(GPBase):
            This is to allow for different normalizations of the output dimensions.
 
         """
-        assert isinstance(self.likelihood,EP_Mixed_Noise)
+        assert hasattr(self,'multioutput')
         index = np.ones_like(Xnew)*output
         Xnew = np.hstack((Xnew,index))
 
@@ -345,6 +366,51 @@ class SparseGP(GPBase):
         mu, var = self._raw_predict(Xnew, full_cov=full_cov, which_parts=which_parts)
 
         # now push through likelihood
-        mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
+        if isinstance(self.likelihood,EP_Mixed_Noise):
+            mean, var, _025pm, _975pm = self.likelihood.predictive_values(mu, var, full_cov, noise_model = output)
+        else:
+            mean, var, _025pm, _975pm = self.likelihood_list[output].predictive_values(mu, var, full_cov)
         return mean, var, _025pm, _975pm
 
+
+
+    def _raw_predict_single_output(self, _Xnew, output=0, X_variance_new=None, which_parts='all', full_cov=False,stop=False):
+        """
+        Internal helper function for making predictions, does not account
+        for normalization or likelihood
+        """
+        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)
+
+        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, _ = dpotrs(self.LB, tmp, lower=1)
+            self.Cpsi1V, _ = dtrtrs(self.Lm, tmp, lower=1, trans=1)
+
+        assert hasattr(self,'multioutput')
+        index = np.ones_like(_Xnew)*output
+        _Xnew = np.hstack((_Xnew,index))
+
+        if X_variance_new is None:
+            Kx = self.kern.K(self.Z, _Xnew, which_parts=which_parts)
+            mu = np.dot(Kx.T, self.Cpsi1V)
+            if full_cov:
+                Kxx = self.kern.K(_Xnew, which_parts=which_parts)
+                var = Kxx - mdot(Kx.T, Kmmi_LmiBLmi, Kx) # NOTE this won't work for plotting
+            else:
+                Kxx = self.kern.Kdiag(_Xnew, which_parts=which_parts)
+                var = Kxx - np.sum(Kx * np.dot(Kmmi_LmiBLmi, Kx), 0)
+        else:
+            # assert which_p.Tarts=='all', "swithching out parts of variational kernels is not implemented"
+            Kx = self.kern.psi1(self.Z, _Xnew, X_variance_new) # , which_parts=which_parts) TODO: which_parts
+            mu = np.dot(Kx, self.Cpsi1V)
+            if full_cov:
+                raise NotImplementedError, "TODO"
+            else:
+                Kxx = self.kern.psi0(self.Z, _Xnew, X_variance_new)
+                psi2 = self.kern.psi2(self.Z, _Xnew, X_variance_new)
+                var = Kxx - np.sum(np.sum(psi2 * Kmmi_LmiBLmi[None, :, :], 1), 1)
+
+        return mu, var[:, None]
diff --git a/GPy/kern/parts/prod.py b/GPy/kern/parts/prod.py
index db31c626..8b4e832b 100644
--- a/GPy/kern/parts/prod.py
+++ b/GPy/kern/parts/prod.py
@@ -18,7 +18,7 @@ class Prod(Kernpart):
     """
     def __init__(self,k1,k2,tensor=False):
         self.num_params = k1.num_params + k2.num_params
-        self.name = k1.name + '<times>' + k2.name
+        self.name = '['+k1.name + '(x)' + k2.name +']'
         self.k1 = k1
         self.k2 = k2
         if tensor:
diff --git a/GPy/likelihoods/ep_mixed_noise.py b/GPy/likelihoods/ep_mixed_noise.py
index 150de09c..a00b0643 100644
--- a/GPy/likelihoods/ep_mixed_noise.py
+++ b/GPy/likelihoods/ep_mixed_noise.py
@@ -249,7 +249,7 @@ class EP_Mixed_Noise(likelihood):
                 self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
                 self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
                 #Marginal moments
-                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
+                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model_list[self.index[i]].moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
                 #Site parameters update
                 Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
                 Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
@@ -344,7 +344,7 @@ class EP_Mixed_Noise(likelihood):
                 self.tau_[i] = 1./Sigma_diag[i] - self.eta*self.tau_tilde[i]
                 self.v_[i] = mu[i]/Sigma_diag[i] - self.eta*self.v_tilde[i]
                 #Marginal moments
-                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model.moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
+                self.Z_hat[i], mu_hat[i], sigma2_hat[i] = self.noise_model_list[self.index[i]].moments_match(self._transf_data[i],self.tau_[i],self.v_[i])
                 #Site parameters update
                 Delta_tau = self.delta/self.eta*(1./sigma2_hat[i] - 1./Sigma_diag[i])
                 Delta_v = self.delta/self.eta*(mu_hat[i]/sigma2_hat[i] - mu[i]/Sigma_diag[i])
diff --git a/GPy/models/gp_multioutput.py b/GPy/models/gp_multioutput.py
index 0fdad786..96138bbf 100644
--- a/GPy/models/gp_multioutput.py
+++ b/GPy/models/gp_multioutput.py
@@ -6,6 +6,7 @@ import numpy as np
 from ..core import GP
 from .. import likelihoods
 from .. import kern
+from ..util import multioutput
 
 
 import pylab as pb
@@ -29,7 +30,7 @@ class GPMultioutput(GP):
 
     """
 
-    def __init__(self,X_list,Y_list,noise_list=[],kernel_list=None,normalize_X=False,normalize_Y=False,W=1): #TODO W
+    def __init__(self,X_list,Y_list,kernel_list=None,normalize_X=False,normalize_Y=False,W=1,mixed_noise_list=[]): #TODO W
 
         assert len(X_list) == len(Y_list)
         index = []
@@ -40,53 +41,30 @@ class GPMultioutput(GP):
             i += 1
         index = np.vstack(index)
 
-        if noise_list == []:
-            likelihood_list = []
+        self.likelihood_list = []
+        if mixed_noise_list == []:
             for Y in Y_list:
-                likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
+                self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
+
+            Y = np.vstack([l_.Y for l_ in self.likelihood_list])
+            likelihood = likelihoods.Gaussian(Y,normalize=False)
+            likelihood.index = index
+
+        else:
+            assert len(Y_list) == len(mixed_noise_list)
+            for noise,Y in zip(mixed_noise_list,Y_list):
+                self.likelihood_list.append(likelihoods.EP(Y,noise))
+            likelihood = likelihoods.EP_Mixed_Noise(Y_list, mixed_noise_list)
 
-        Y = np.vstack([l_.Y for l_ in likelihood_list])
-        likelihood = likelihoods.Gaussian(Y,normalize=False)
-        likelihood.index = index
 
         X = np.hstack([np.vstack(X_list),index])
+        original_dim = X.shape[1] - 1
 
         if kernel_list is None:
-            original_dim = X.shape[1]-1
-            kernel_list = [kern.rbf(original_dim) + kern.white(original_dim)]
+            kernel_list = [[kern.rbf(original_dim)],[kern.white(original_dim+1)]]
+
+        mkernel = multioutput.build_cor_kernel(input_dim=original_dim, Nout=len(X_list), CK = kernel_list[0], NC = kernel_list[1], W=1)
 
-        mkernel = kernel_list[0].prod(kern.coregionalise(len(X_list),W),tensor=True)
-        for k in kernel_list[1:]:
-            mkernel += k.prod(kern.coregionalise(len(X_list),W),tensor=True)
         self.multioutput = True
         GP.__init__(self, X, likelihood, mkernel, normalize_X=normalize_X)
         self.ensure_default_constraints()
-
-
-"""
-if likelihood is None:
-noise_model_list = []
-for Y in Y_list:
-noise_model_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
-#noise_model_list = [likelihoods.gaussian(variance=1.) for Y in Y_list]
-#likelihood = likelihoods.EP_Mixed_Noise(Y_list, noise_model_list)
-
-elif Y_list is not None:
-if not all(np.vstack(Y_list).flatten() == likelihood.data.flatten()):
-raise Warning, 'likelihood.data and Y_list values are different.'
-
-X = np.hstack([np.vstack(X_list),likelihood.index])
-
-if kernel_list is None:
-original_dim = X.shape[1]-1
-kernel_list = [kern.rbf(original_dim) + kern.white(original_dim)]
-
-mkernel = kernel_list[0].prod(kern.coregionalise(len(X_list),W),tensor=True) #TODO W
-for k in kernel_list[1:]:
-mkernel += k.prod(kern.coregionalise(len(X_list),W),tensor=True) #TODO W
-
-#kern1 = kern.rbf(1) + kern.white(1)
-#kern2 = kern.coregionalise(2,1)
-#kern3 = kern1.prod(kern2,tensor=True)
-"""
-
diff --git a/GPy/models/sparse_gp_multioutput.py b/GPy/models/sparse_gp_multioutput.py
new file mode 100644
index 00000000..ae99421f
--- /dev/null
+++ b/GPy/models/sparse_gp_multioutput.py
@@ -0,0 +1,97 @@
+# Copyright (c) 2013, Ricardo Andrade
+# Licensed under the BSD 3-clause license (see LICENSE.txt)
+
+
+import numpy as np
+from ..core import SparseGP
+from .. import likelihoods
+from .. import kern
+from ..util import multioutput
+
+
+import pylab as pb
+
+class SparseGPMultioutput(SparseGP):
+    """
+    Multiple output Gaussian process
+
+    This is a thin wrapper around the models.GP class, with a set of sensible defaults
+
+    :param X_list: input observations
+    :param Y_list: observed values
+    :param L_list: a GPy likelihood, defaults to Binomial with probit link_function
+    :param kernel_list: a GPy kernel, defaults to rbf
+    :param normalize_X:  whether to normalize the input data before computing (predictions will be in original scales)
+    :type normalize_X: False|True
+    :param normalize_Y:  whether to normalize the input data before computing (predictions will be in original scales)
+    :type normalize_Y: False|True
+
+    .. Note:: Multiple independent outputs are allowed using columns of Y
+
+    """
+
+    def __init__(self,X_list,Y_list,kernel_list=None,normalize_X=False,normalize_Y=False,Z_list=None,num_inducing_list=10,X_variance=None,W=1,mixed_noise_list=[]): #TODO W
+
+        assert len(X_list) == len(Y_list)
+        index = []
+        for x,y,j in zip(X_list,Y_list,range(len(X_list))):
+            assert x.shape[0] == y.shape[0]
+            index.append(np.repeat(j,y.size)[:,None])
+        index = np.vstack(index)
+
+
+        self.likelihood_list = []
+        if mixed_noise_list == []:
+            for Y in Y_list:
+                self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
+
+            Y = np.vstack([l_.Y for l_ in self.likelihood_list])
+            likelihood = likelihoods.Gaussian(Y,normalize=False)
+            likelihood.index = index
+
+        else:
+            assert len(Y_list) == len(mixed_noise_list)
+            for noise,Y in zip(mixed_noise_list,Y_list):
+                self.likelihood_list.append(likelihoods.EP(Y,noise))
+            likelihood = likelihoods.EP_Mixed_Noise(Y_list, mixed_noise_list)
+
+        """
+        if noise_list == []:
+            self.likelihood_list = []
+            for Y in Y_list:
+                self.likelihood_list.append(likelihoods.Gaussian(Y,normalize = normalize_Y))
+
+        Y = np.vstack([l_.Y for l_ in self.likelihood_list])
+        likelihood = likelihoods.Gaussian(Y,normalize=False)
+        likelihood.index = index
+        """
+        X = np.hstack([np.vstack(X_list),index])
+        original_dim = X.shape[1] - 1
+
+        if kernel_list is None:
+            kernel_list = [[kern.rbf(original_dim)],[kern.white(original_dim+1)]]
+
+        mkernel = multioutput.build_cor_kernel(input_dim=original_dim, Nout=len(X_list), CK = kernel_list[0], NC = kernel_list[1], W=1)
+
+        z_index = []
+        if Z_list is None:
+            if isinstance(num_inducing_list,int):
+                num_inducing_list = [num_inducing_list for Xj in X_list]
+            Z_list = []
+            for Xj,nj,j in zip(X_list,num_inducing_list,range(len(X_list))):
+                i = np.random.permutation(Xj.shape[0])[:nj]
+                z_index.append(np.repeat(j,nj)[:,None])
+                Z_list.append(Xj[i].copy())
+        else:
+            assert len(Z_list) == len(X_list)
+            for Zj,Xj,j in zip(Z_list,X_list,range(len(Z_list))):
+                assert Zj.shape[1] == Xj.shape[1]
+                z_index.append(np.repeat(j,Zj.shape[0])[:,None])
+
+        Z = np.hstack([np.vstack(Z_list),np.vstack(z_index)])
+
+
+        self.multioutput = True
+        SparseGP.__init__(self, X, likelihood, mkernel, Z=Z, normalize_X=normalize_X, X_variance=X_variance)
+        self.constrain_fixed('.*iip_\d+_1')
+        self.ensure_default_constraints()
diff --git a/GPy/util/__init__.py b/GPy/util/__init__.py
index 27d25518..f5384356 100644
--- a/GPy/util/__init__.py
+++ b/GPy/util/__init__.py
@@ -14,3 +14,4 @@ import mocap
 import visualize
 import decorators
 import classification
+import multioutput
diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py
new file mode 100644
index 00000000..44b70b6f
--- /dev/null
+++ b/GPy/util/multioutput.py
@@ -0,0 +1,35 @@
+import numpy as np
+import warnings
+from .. import kern
+
+def build_cor_kernel(input_dim, Nout, CK = [], NC = [], W=1):
+    """
+    Builds an appropiate coregionalized kernel
+
+    :input_dim: Input dimensionality
+    :Nout: Number of outputs
+    :param CK: List of coregionalized kernels (i.e., this will be multiplied by a coregionalise kernel).
+    :param K: List of kernels that won't be multiplied by a coregionalise kernel
+    :W:
+    """
+
+    for k in CK:
+        if k.input_dim <> input_dim:
+            k.input_dim = input_dim
+            #raise Warning("kernel's input dimension overwritten to fit input_dim parameter.")
+            warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
+
+    for k in NC:
+        if k.input_dim <> input_dim + 1:
+            k.input_dim = input_dim + 1
+            #raise Warning("kernel's input dimension overwritten to fit input_dim parameter.")
+            warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
+
+    kernel = CK[0].prod(kern.coregionalise(Nout,W),tensor=True)
+    for k in CK[1:]:
+        kernel += k.prod(kern.coregionalise(Nout,W),tensor=True)
+
+    for k in NC:
+        kernel += k
+
+    return kernel