diff --git a/GPy/core/parameterization/parameter_core.py b/GPy/core/parameterization/parameter_core.py
index 6a8f1b1d..48fe69c2 100644
--- a/GPy/core/parameterization/parameter_core.py
+++ b/GPy/core/parameterization/parameter_core.py
@@ -541,12 +541,12 @@ class Constrainable(Nameable, Indexable):
             print "WARNING: reconstraining parameters {}".format(self.parameter_names() or self.name)
         which.add(what, self._raveled_index())
 
-    def _remove_from_index_operations(self, which, what):
+    def _remove_from_index_operations(self, which, transforms):
         """
         Helper preventing copy code.
         Remove given what (transform prior etc) from which param index ops.
         """
-        if len(what) == 0:
+        if len(transforms) == 0:
             transforms = which.properties()
         removed = np.empty((0,), dtype=int)
         for t in transforms:
@@ -567,10 +567,10 @@ class OptimizationHandlable(Constrainable):
         super(OptimizationHandlable, self).__init__(name, default_constraint=default_constraint, *a, **kw)
 
     def transform(self):
-        [np.put(self._param_array_, ind, c.finv(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
+        [np.put(self._param_array_, ind, c.finv(self._param_array_.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
 
     def untransform(self):
-        [np.put(self._param_array_, ind, c.f(self._param_array_[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
+        [np.put(self._param_array_, ind, c.f(self._param_array_.flat[ind])) for c, ind in self.constraints.iteritems() if c != __fixed__]
 
     def _get_params_transformed(self):
         # transformed parameters (apply transformation rules)
@@ -590,9 +590,9 @@ class OptimizationHandlable(Constrainable):
         if self.has_parent() and self.constraints[__fixed__].size != 0:
             fixes = np.ones(self.size).astype(bool)
             fixes[self.constraints[__fixed__]] = FIXED
-            self._param_array_[fixes] = p
-        elif self._has_fixes(): self._param_array_[self._fixes_] = p
-        else: self._param_array_[:] = p
+            self._param_array_.flat[fixes] = p
+        elif self._has_fixes(): self._param_array_.flat[self._fixes_] = p
+        else: self._param_array_.flat = p
         self.untransform()
         self._trigger_params_changed()
 
@@ -670,8 +670,8 @@ class OptimizationHandlable(Constrainable):
         for pi in self._parameters_:
             pislice = slice(pi_old_size, pi_old_size+pi.size)
 
-            self._param_array_[pislice] = pi._param_array_.ravel()#, requirements=['C', 'W']).flat
-            self._gradient_array_[pislice] = pi._gradient_array_.ravel()#, requirements=['C', 'W']).flat
+            self._param_array_[pislice] = pi._param_array_.flat#, requirements=['C', 'W']).flat
+            self._gradient_array_[pislice] = pi._gradient_array_.flat#, requirements=['C', 'W']).flat
 
             pi._param_array_.data = parray[pislice].data
             pi._gradient_array_.data = garray[pislice].data
@@ -878,8 +878,8 @@ class Parameterizable(OptimizationHandlable):
             # first connect all children
             p._propagate_param_grad(self._param_array_[pslice], self._gradient_array_[pslice])
             # then connect children to self
-            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._param_array_[pslice] = p._param_array_.flat#, requirements=['C', 'W']).ravel(order='C')
+            self._gradient_array_[pslice] = p._gradient_array_.flat#, requirements=['C', 'W']).ravel(order='C')
 
             if not p._param_array_.flags['C_CONTIGUOUS']:
                 import ipdb;ipdb.set_trace()
diff --git a/GPy/core/sparse_gp.py b/GPy/core/sparse_gp.py
index a0b09564..0b796171 100644
--- a/GPy/core/sparse_gp.py
+++ b/GPy/core/sparse_gp.py
@@ -64,8 +64,8 @@ class SparseGP(GP):
             self.kern.gradient += target
 
             #gradients wrt Z
-            self.Z.gradient[:,self.kern.active_dims] = self.kern.gradients_X(dL_dKmm, self.Z)
-            self.Z.gradient[:,self.kern.active_dims] += self.kern.gradients_Z_expectations(
+            self.Z.gradient = self.kern.gradients_X(dL_dKmm, self.Z)
+            self.Z.gradient += self.kern.gradients_Z_expectations(
                                self.grad_dict['dL_dpsi1'], self.grad_dict['dL_dpsi2'], Z=self.Z, variational_posterior=self.X)
         else:
             #gradients wrt kernel
@@ -76,8 +76,8 @@ class SparseGP(GP):
             self.kern.update_gradients_full(self.grad_dict['dL_dKmm'], self.Z, None)
             self.kern.gradient += target
             #gradients wrt Z
-            self.Z.gradient[:,self.kern.active_dims] = 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_dKnm'].T, self.Z, self.X)
+            self.Z.gradient = 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)
 
     def _raw_predict(self, Xnew, full_cov=False):
         """
diff --git a/GPy/examples/coreg_example.py b/GPy/examples/coreg_example.py
index 8f4cfcc6..66ba143d 100644
--- a/GPy/examples/coreg_example.py
+++ b/GPy/examples/coreg_example.py
@@ -23,13 +23,10 @@ K = Bias.prod(Coreg,name='X')
 
 #K.coregion.W = 0
 #print K.coregion.W
-
 #print Bias.K(_X,_X)
 #print K.K(X,X)
-
 #pb.matshow(K.K(X,X))
 
-
 Mlist = [GPy.kern.Matern32(1,lengthscale=20.,name="Mat")]
 kern = GPy.util.multioutput.LCM(input_dim=1,num_outputs=2,kernels_list=Mlist,name='H')
 kern.B.W = 0
@@ -37,16 +34,22 @@ kern.B.kappa = 1.
 #kern.B.W.fix()
 #kern.B.kappa.fix()
 #m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2], kernel=kern)
-m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], kernel=kern)
+
+
+Z1 = np.array([1.5,2.5])[:,None]
+
+m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1], Y_list=[Y1], Z_list = [Z1], kernel=kern)
 #m.optimize()
 m.checkgrad(verbose=1)
 
+"""
 fig = pb.figure()
 ax0 = fig.add_subplot(211)
 ax1 = fig.add_subplot(212)
 slices = GPy.util.multioutput.get_slices([Y1,Y2])
 m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],ax=ax0)
 #m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],ax=ax1)
+"""
 
 
 
diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index dfd922f0..ea997d63 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -160,6 +160,7 @@ def swiss_roll(optimize=True, verbose=1, plot=True, N=1000, num_inducing=15, Q=4
 def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40, max_iters=1000, **k):
     import GPy
     from matplotlib import pyplot as plt
+    from ..util.misc import param_to_array
 
     _np.random.seed(0)
     data = GPy.util.datasets.oil()
@@ -173,11 +174,11 @@ def bgplvm_oil(optimize=True, verbose=1, plot=True, N=200, Q=7, num_inducing=40,
         m.optimize('scg', messages=verbose, max_iters=max_iters, gtol=.05)
 
     if plot:
-        y = m.Y[0, :]
+        y = m.Y
         fig, (latent_axes, sense_axes) = plt.subplots(1, 2)
         m.plot_latent(ax=latent_axes)
         data_show = GPy.plotting.matplot_dep.visualize.vector_show(y)
-        lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(m.X[0, :], # @UnusedVariable
+        lvm_visualizer = GPy.plotting.matplot_dep.visualize.lvm_dimselect(param_to_array(m.X.mean), # @UnusedVariable
             m, data_show, latent_axes=latent_axes, sense_axes=sense_axes)
         raw_input('Press enter to finish')
         plt.close(fig)
diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index 190af93b..2a4b91b3 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -25,80 +25,51 @@ def olympic_marathon_men(optimize=True, plot=True):
 
     return m
 
-def coregionalization_toy2(optimize=True, plot=True):
+def coregionalization_toy(optimize=True, plot=True):
     """
     A simple demonstration of coregionalization on two sinusoidal functions.
     """
     #build a design matrix with a column of integers indicating the output
     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))
 
     #build a suitable set of observed variables
     Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
     Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
-    Y = np.vstack((Y1, Y2))
 
-    #build the kernel
-    k1 = GPy.kern.RBF(1) + GPy.kern.Bias(1)
-    k2 = GPy.kern.Coregionalize(2,1)
-    k = k1**k2
-    m = GPy.models.GPRegression(X, Y, kernel=k)
+    m = GPy.models.GPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
 
     if optimize:
         m.optimize('bfgs', max_iters=100)
 
     if plot:
-        m.plot(fixed_inputs=[(1,0)])
-        m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
-
+        slices = GPy.util.multioutput.get_slices([X1,X2])
+        m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
+        m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
     return m
 
-#FIXME: Needs recovering once likelihoods are consolidated
-#def coregionalization_toy(optimize=True, plot=True):
-#    """
-#    A simple demonstration of coregionalization on two sinusoidal functions.
-#    """
-#    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.)
-#    m.optimize(max_iters=100)
-#
-#    fig, axes = pb.subplots(2,1)
-#    m.plot(fixed_inputs=[(1,0)],ax=axes[0])
-#    m.plot(fixed_inputs=[(1,1)],ax=axes[1])
-#    axes[0].set_title('Output 0')
-#    axes[1].set_title('Output 1')
-#    return m
-
 def coregionalization_sparse(optimize=True, plot=True):
     """
     A simple demonstration of coregionalization on two sinusoidal functions using sparse approximations.
     """
-    #fetch the data from the non sparse examples
-    m = coregionalization_toy2(optimize=False, plot=False)
-    X, Y = m.X, m.Y
+    #build a design matrix with a column of integers indicating the output
+    X1 = np.random.rand(50, 1) * 8
+    X2 = np.random.rand(30, 1) * 5
 
-    k = GPy.kern.RBF(1)**GPy.kern.Coregionalize(2)
+    #build a suitable set of observed variables
+    Y1 = np.sin(X1) + np.random.randn(*X1.shape) * 0.05
+    Y2 = np.sin(X2) + np.random.randn(*X2.shape) * 0.05 + 2.
 
-    #construct a model
-    m = GPy.models.SparseGPRegression(X,Y, num_inducing=25, kernel=k)
-    m.Z[:,1].fix() # don't optimize the inducing input indexes
+    m = GPy.models.SparseGPCoregionalizedRegression(X_list=[X1,X2], Y_list=[Y1,Y2])
 
     if optimize:
-        m.optimize('bfgs', max_iters=100, messages=1)
+        m.optimize('bfgs', max_iters=100)
 
     if plot:
-        m.plot(fixed_inputs=[(1,0)])
-        m.plot(fixed_inputs=[(1,1)], ax=pb.gca())
+        slices = GPy.util.multioutput.get_slices([X1,X2])
+        m.plot(fixed_inputs=[(1,0)],which_data_rows=slices[0],Y_metadata={'output_index':0})
+        m.plot(fixed_inputs=[(1,1)],which_data_rows=slices[1],Y_metadata={'output_index':1},ax=pb.gca())
+        pb.ylim(-3,)
 
     return m
 
diff --git a/GPy/inference/latent_function_inference/expectation_propagation.py b/GPy/inference/latent_function_inference/expectation_propagation.py
index 514a6dc7..ff60d2e3 100644
--- a/GPy/inference/latent_function_inference/expectation_propagation.py
+++ b/GPy/inference/latent_function_inference/expectation_propagation.py
@@ -11,9 +11,9 @@ class EP(object):
 
         :param epsilon: Convergence criterion, maximum squared difference allowed between mean updates to stop iterations (float)
         :type epsilon: float
-        :param eta: Power EP thing TODO: Ricardo: what, exactly?
+        :param eta: parameter for fractional EP updates.
         :type eta: float64
-        :param delta: Power EP thing TODO: Ricardo: what, exactly?
+        :param delta: damping EP updates factor.
         :type delta: float64
         """
         self.epsilon, self.eta, self.delta = epsilon, eta, delta
diff --git a/GPy/inference/latent_function_inference/var_dtc.py b/GPy/inference/latent_function_inference/var_dtc.py
index e2aa95f5..669ccd08 100644
--- a/GPy/inference/latent_function_inference/var_dtc.py
+++ b/GPy/inference/latent_function_inference/var_dtc.py
@@ -176,7 +176,6 @@ class VarDTC(object):
 
         #construct a posterior object
         post = Posterior(woodbury_inv=woodbury_inv, woodbury_vector=woodbury_vector, K=Kmm, mean=None, cov=None, K_chol=Lm)
-
         return post, log_marginal, grad_dict
 
 class VarDTCMissingData(object):
@@ -365,7 +364,7 @@ class VarDTCMissingData(object):
         return post, log_marginal, grad_dict
 
 def _compute_dL_dpsi(num_inducing, num_data, output_dim, beta, Lm, VVT_factor, Cpsi1Vf, DBi_plus_BiPBi, psi1, het_noise, uncertain_inputs):
-    dL_dpsi0 = -0.5 * output_dim * (beta * np.ones([num_data, 1])).flatten()
+    dL_dpsi0 = -0.5 * output_dim * (beta[:,None] * np.ones([num_data, 1])).flatten()
     dL_dpsi1 = np.dot(VVT_factor, Cpsi1Vf.T)
     dL_dpsi2_beta = 0.5 * backsub_both_sides(Lm, output_dim * np.eye(num_inducing) - DBi_plus_BiPBi)
     if het_noise:
diff --git a/GPy/kern/_src/ODE_UY.py b/GPy/kern/_src/ODE_UY.py
index cc68416b..510b4f7c 100644
--- a/GPy/kern/_src/ODE_UY.py
+++ b/GPy/kern/_src/ODE_UY.py
@@ -139,11 +139,11 @@ class ODE_UY(Kern):
         dVdly = 0.5/np.sqrt(ly)*np.sqrt(2*Vy)
         dVdVy = 0.5/np.sqrt(Vy)*np.sqrt(2*ly)
 
-        rd=rdist.shape[0]
-        dktheta1 = np.zeros([rd,rd])
-        dktheta2 = np.zeros([rd,rd])
-        dkUdvar = np.zeros([rd,rd])
-        dkYdvar = np.zeros([rd,rd])
+        rd=rdist.shape
+        dktheta1 = np.zeros(rd)
+        dktheta2 = np.zeros(rd)
+        dkUdvar = np.zeros(rd)
+        dkYdvar = np.zeros(rd)
 
         # dk dtheta for UU
         UUdtheta1 = lambda dist: np.exp(-lu* dist)*dist + (-dist)*np.exp(-lu* dist)*(1+lu*dist)
diff --git a/GPy/kern/_src/add.py b/GPy/kern/_src/add.py
index cb73087e..57e611ed 100644
--- a/GPy/kern/_src/add.py
+++ b/GPy/kern/_src/add.py
@@ -23,7 +23,7 @@ class Add(CombinationKernel):
         If a list of parts (of this kernel!) `which_parts` is given, only
         the parts of the list are taken to compute the covariance.
         """
-        assert X.shape[1] == self.input_dim
+        assert X.shape[1] > max(np.r_[self.active_dims])
         if which_parts is None:
             which_parts = self.parts
         elif not isinstance(which_parts, (list, tuple)):
@@ -33,7 +33,7 @@ class Add(CombinationKernel):
 
     @Cache_this(limit=2, force_kwargs=['which_parts'])
     def Kdiag(self, X, which_parts=None):
-        assert X.shape[1] == self.input_dim
+        assert X.shape[1] > max(np.r_[self.active_dims])
         if which_parts is None:
             which_parts = self.parts
         elif not isinstance(which_parts, (list, tuple)):
@@ -172,7 +172,7 @@ class Add(CombinationKernel):
 
     def add(self, other, name='sum'):
         if isinstance(other, Add):
-            other_params = other._parameters_.copy()
+            other_params = other._parameters_[:]
             for p in other_params:
                 other.remove_parameter(p)
             self.add_parameters(*other_params)
diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py
index 15e23d5c..7d9eeac2 100644
--- a/GPy/kern/_src/linear.py
+++ b/GPy/kern/_src/linear.py
@@ -312,5 +312,4 @@ class Linear(Kern):
         return np.dot(ZA, inner).swapaxes(0, 1)  # NOTE: self.ZAinner \in [num_inducing x num_data x input_dim]!
 
     def input_sensitivity(self):
-        if self.ARD: return self.variances
-        else: return self.variances.repeat(self.input_dim)
+        return np.ones(self.input_dim) * self.variances
diff --git a/GPy/likelihoods/gaussian.py b/GPy/likelihoods/gaussian.py
index c7001278..6a5030bc 100644
--- a/GPy/likelihoods/gaussian.py
+++ b/GPy/likelihoods/gaussian.py
@@ -92,7 +92,7 @@ class Gaussian(Likelihood):
     def predictive_variance(self, mu, sigma, predictive_mean=None):
         return self.variance + sigma**2
 
-    def predictive_quantiles(self, mu, var, quantiles, Y_metadata):
+    def predictive_quantiles(self, mu, var, quantiles, Y_metadata=None):
         return  [stats.norm.ppf(q/100.)*np.sqrt(var) + mu for q in quantiles]
 
     def pdf_link(self, link_f, y, Y_metadata=None):
diff --git a/GPy/likelihoods/likelihood.py b/GPy/likelihoods/likelihood.py
index 6df573df..aabe93ef 100644
--- a/GPy/likelihoods/likelihood.py
+++ b/GPy/likelihoods/likelihood.py
@@ -411,10 +411,7 @@ class Likelihood(Parameterized):
         #compute the quantiles by sampling!!!
         N_samp = 1000
         s = np.random.randn(mu.shape[0], N_samp)*np.sqrt(var) + mu
-        #ss_f = s.flatten()
-        #ss_y = self.samples(ss_f, Y_metadata)
         ss_y = self.samples(s, Y_metadata)
-        #ss_y = ss_y.reshape(mu.shape[0], N_samp)
 
         return [np.percentile(ss_y ,q, axis=1)[:,None] for q in quantiles]
 
diff --git a/GPy/likelihoods/mixed_noise.py b/GPy/likelihoods/mixed_noise.py
index 5f4d0705..c2435508 100644
--- a/GPy/likelihoods/mixed_noise.py
+++ b/GPy/likelihoods/mixed_noise.py
@@ -11,7 +11,7 @@ import itertools
 
 class MixedNoise(Likelihood):
     def __init__(self, likelihoods_list, name='mixed_noise'):
-
+        #NOTE at the moment this likelihood only works for using a list of gaussians
         super(Likelihood, self).__init__(name=name)
 
         self.add_parameters(*likelihoods_list)
@@ -24,11 +24,11 @@ class MixedNoise(Likelihood):
         variance = np.zeros(ind.size)
         for lik, j in zip(self.likelihoods_list, range(len(self.likelihoods_list))):
             variance[ind==j] = lik.variance
-        return variance[:,None]
+        return variance
 
     def betaY(self,Y,Y_metadata):
         #TODO not here.
-        return Y/self.gaussian_variance(Y_metadata=Y_metadata)
+        return Y/self.gaussian_variance(Y_metadata=Y_metadata)[:,None]
 
     def update_gradients(self, gradients):
         self.gradient = gradients
@@ -39,24 +39,27 @@ class MixedNoise(Likelihood):
         return np.array([dL_dKdiag[ind==i].sum() for i in range(len(self.likelihoods_list))])
 
     def predictive_values(self, mu, var, full_cov=False, Y_metadata=None):
-        if all([isinstance(l, Gaussian) for l in self.likelihoods_list]):
-            ind = Y_metadata['output_index'].flatten()
-            _variance = np.array([self.likelihoods_list[j].variance for j in ind ])
-            if full_cov:
-                var += np.eye(var.shape[0])*_variance
-            else:
-                var += _variance
-            return mu, var
+        ind = Y_metadata['output_index'].flatten()
+        _variance = np.array([self.likelihoods_list[j].variance for j in ind ])
+        if full_cov:
+            var += np.eye(var.shape[0])*_variance
         else:
-            raise NotImplementedError
+            var += _variance
+        return mu, var
 
-    def predictive_variance(self, mu, sigma, **other_shit):
-        if isinstance(noise_index,int):
-            _variance = self.variance[noise_index]
-        else:
-            _variance = np.array([ self.variance[j] for j in noise_index ])[:,None]
+    def predictive_variance(self, mu, sigma, Y_metadata):
+        _variance = self.gaussian_variance(Y_metadata)
         return _variance + sigma**2
 
+    def predictive_quantiles(self, mu, var, quantiles, Y_metadata):
+        ind = Y_metadata['output_index'].flatten()
+        outputs = np.unique(ind)
+        Q = np.zeros( (mu.size,len(quantiles)) )
+        for j in outputs:
+            q = self.likelihoods_list[j].predictive_quantiles(mu[ind==j,:],
+                var[ind==j,:],quantiles,Y_metadata=None)
+            Q[ind==j,:] = np.hstack(q)
+        return [q[:,None] for q in Q.T]
 
     def samples(self, gp, Y_metadata):
         """
@@ -75,4 +78,3 @@ class MixedNoise(Likelihood):
             _ysim = np.array([np.random.normal(lik.gp_link.transf(gpj), scale=np.sqrt(lik.variance), size=1) for gpj in gp_filtered.flatten()])
             Ysim[flt,:] = _ysim.reshape(n1,N2)
         return Ysim
-
diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py
index e0818a2f..9a6a8f4c 100644
--- a/GPy/models/bayesian_gplvm.py
+++ b/GPy/models/bayesian_gplvm.py
@@ -72,15 +72,19 @@ class BayesianGPLVM(SparseGP):
         self.variational_prior.update_gradients_KL(self.X)
 
 
-    def plot_latent(self, plot_inducing=True, *args, **kwargs):
-        """
-        See GPy.plotting.matplot_dep.dim_reduction_plots.plot_latent
-        """
+    def plot_latent(self, labels=None, which_indices=None,
+                resolution=50, ax=None, marker='o', s=40,
+                fignum=None, plot_inducing=True, legend=True,
+                plot_limits=None, 
+                aspect='auto', updates=False, **kwargs):
         import sys
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
         from ..plotting.matplot_dep import dim_reduction_plots
 
-        return dim_reduction_plots.plot_latent(self, plot_inducing=plot_inducing, *args, **kwargs)
+        return dim_reduction_plots.plot_latent(self, labels, which_indices,
+                resolution, ax, marker, s,
+                fignum, plot_inducing, legend,
+                plot_limits, aspect, updates, **kwargs)
 
     def do_test_latents(self, Y):
         """
diff --git a/GPy/models/gp_coregionalized_regression.py b/GPy/models/gp_coregionalized_regression.py
index 6d478fd9..5854fe63 100644
--- a/GPy/models/gp_coregionalized_regression.py
+++ b/GPy/models/gp_coregionalized_regression.py
@@ -36,7 +36,7 @@ class GPCoregionalizedRegression(GP):
 
         #Kernel
         if kernel is None:
-            kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name)
+            kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name)
 
         #Likelihood
         likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list)
diff --git a/GPy/models/gplvm.py b/GPy/models/gplvm.py
index 5f7e3265..b85540ce 100644
--- a/GPy/models/gplvm.py
+++ b/GPy/models/gplvm.py
@@ -67,12 +67,22 @@ class GPLVM(GP):
         assert self.likelihood.Y.shape[1] == 2
         pb.scatter(self.likelihood.Y[:, 0], self.likelihood.Y[:, 1], 40, self.X[:, 0].copy(), linewidth=0, cmap=pb.cm.jet)  # @UndefinedVariable
         Xnew = np.linspace(self.X.min(), self.X.max(), 200)[:, None]
-        mu, var, upper, lower = self.predict(Xnew)
+        mu, _ = self.predict(Xnew)
         pb.plot(mu[:, 0], mu[:, 1], 'k', linewidth=1.5)
 
-    def plot_latent(self, *args, **kwargs):
+    def plot_latent(self, labels=None, which_indices=None,
+                resolution=50, ax=None, marker='o', s=40,
+                fignum=None, legend=True,
+                plot_limits=None, 
+                aspect='auto', updates=False, **kwargs):
+        import sys
+        assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
         from ..plotting.matplot_dep import dim_reduction_plots
 
-        return dim_reduction_plots.plot_latent(self, *args, **kwargs)
+        return dim_reduction_plots.plot_latent(self, labels, which_indices,
+                resolution, ax, marker, s,
+                fignum, False, legend,
+                plot_limits, aspect, updates, **kwargs)
+
     def plot_magnification(self, *args, **kwargs):
         return util.plot_latent.plot_magnification(self, *args, **kwargs)
diff --git a/GPy/models/sparse_gp_coregionalized_regression.py b/GPy/models/sparse_gp_coregionalized_regression.py
index a97696d2..ea6ec9f1 100644
--- a/GPy/models/sparse_gp_coregionalized_regression.py
+++ b/GPy/models/sparse_gp_coregionalized_regression.py
@@ -43,14 +43,14 @@ class SparseGPCoregionalizedRegression(SparseGP):
 
         #Kernel
         if kernel is None:
-            kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=GPy.kern.rbf(X.shape[1]-1), W_rank=1,name=kernel_name)
+            kernel = util.multioutput.ICM(input_dim=X.shape[1]-1, num_outputs=Ny, kernel=kern.RBF(X.shape[1]-1), W_rank=1,name=kernel_name)
 
         #Likelihood
         likelihood = util.multioutput.build_likelihood(Y_list,self.output_index,likelihoods_list)
 
         #Inducing inputs list
         if len(Z_list):
-            assert len(Z_list) == self.output_dim, 'Number of outputs do not match length of inducing inputs list.'
+            assert len(Z_list) == Ny, 'Number of outputs do not match length of inducing inputs list.'
         else:
             if isinstance(num_inducing,np.int):
                 num_inducing = [num_inducing] * Ny
diff --git a/GPy/plotting/matplot_dep/dim_reduction_plots.py b/GPy/plotting/matplot_dep/dim_reduction_plots.py
index bf9297b9..57d932cc 100644
--- a/GPy/plotting/matplot_dep/dim_reduction_plots.py
+++ b/GPy/plotting/matplot_dep/dim_reduction_plots.py
@@ -30,7 +30,8 @@ def most_significant_input_dimensions(model, which_indices):
 def plot_latent(model, labels=None, which_indices=None,
                 resolution=50, ax=None, marker='o', s=40,
                 fignum=None, plot_inducing=False, legend=True,
-                aspect='auto', updates=False):
+                plot_limits=None, 
+                aspect='auto', updates=False, **kwargs):
     """
     :param labels: a np.array of size model.num_data containing labels for the points (can be number, strings, etc)
     :param resolution: the resolution of the grid on which to evaluate the predictive variance
@@ -38,6 +39,8 @@ def plot_latent(model, labels=None, which_indices=None,
     if ax is None:
         fig = pb.figure(num=fignum)
         ax = fig.add_subplot(111)
+    else:
+        fig = ax.figure
     Tango.reset()
 
     if labels is None:
@@ -57,15 +60,28 @@ def plot_latent(model, labels=None, which_indices=None,
     def plot_function(x):
         Xtest_full = np.zeros((x.shape[0], model.X.shape[1]))
         Xtest_full[:, [input_1, input_2]] = x
-        mu, var, low, up = model.predict(Xtest_full)
+        _, var = model.predict(Xtest_full)
         var = var[:, :1]
         return np.log(var)
 
     #Create an IMshow controller that can re-plot the latent space shading at a good resolution
+    if plot_limits is None:
+        xmin, ymin = X[:, [input_1, input_2]].min(0)
+        xmax, ymax = X[:, [input_1, input_2]].max(0)
+        x_r, y_r = xmax-xmin, ymax-ymin
+        xmin -= .1*x_r
+        xmax += .1*x_r
+        ymin -= .1*y_r
+        ymax += .1*y_r
+    else:
+        try:
+            xmin, xmax, ymin, ymax = plot_limits
+        except (TypeError, ValueError) as e:
+            raise e.__class__, "Wrong plot limits: {} given -> need (xmin, xmax, ymin, ymax)".format(plot_limits)
     view = ImshowController(ax, plot_function,
-                            tuple(X[:, [input_1, input_2]].min(0)) + tuple(X[:, [input_1, input_2]].max(0)),
+                            (xmin, ymin, xmax, ymax),
                             resolution, aspect=aspect, interpolation='bilinear',
-                            cmap=pb.cm.binary)
+                            cmap=pb.cm.binary, **kwargs)
 
     # make sure labels are in order of input:
     ulabels = []
@@ -99,18 +115,31 @@ def plot_latent(model, labels=None, which_indices=None,
     if not np.all(labels == 1.) and legend:
         ax.legend(loc=0, numpoints=1)
 
-    #ax.set_xlim(xmin[0], xmax[0])
-    #ax.set_ylim(xmin[1], xmax[1])
     ax.grid(b=False) # remove the grid if present, it doesn't look good
     ax.set_aspect('auto') # set a nice aspect ratio
 
     if plot_inducing:
         Z = param_to_array(model.Z)
         ax.plot(Z[:, input_1], Z[:, input_2], '^w')
+    
+    ax.set_xlim((xmin, xmax))
+    ax.set_ylim((ymin, ymax))
 
+    try:
+        fig.canvas.draw()
+        fig.tight_layout()
+        fig.canvas.draw()
+    except Exception as e:
+        print "Could not invoke tight layout: {}".format(e)
+        pass
+    
     if updates:
-        ax.figure.canvas.show()
+        try:
+            ax.figure.canvas.show()
+        except Exception as e:
+            print "Could not invoke show: {}".format(e)
         raw_input('Enter to continue')
+        view.deactivate()
     return ax
 
 def plot_magnification(model, labels=None, which_indices=None,
@@ -186,7 +215,7 @@ def plot_magnification(model, labels=None, which_indices=None,
         ax.plot(model.Z[:, input_1], model.Z[:, input_2], '^w')
 
     if updates:
-        ax.figure.canvas.show()
+        fig.canvas.show()
         raw_input('Enter to continue')
 
     pb.title('Magnification Factor')
diff --git a/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/axis_event_controller.py b/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/axis_event_controller.py
index d5aaefd2..62b622c5 100644
--- a/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/axis_event_controller.py
+++ b/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/axis_event_controller.py
@@ -33,7 +33,7 @@ class AxisChangedController(AxisEventController):
         Constructor
         '''
         super(AxisChangedController, self).__init__(ax)
-        self._lim_ratio_threshold = update_lim or .8
+        self._lim_ratio_threshold = update_lim or .95
         self._x_lim = self.ax.get_xlim()
         self._y_lim = self.ax.get_ylim()
 
@@ -80,6 +80,10 @@ class AxisChangedController(AxisEventController):
 class BufferedAxisChangedController(AxisChangedController):
     def __init__(self, ax, plot_function, plot_limits, resolution=50, update_lim=None, **kwargs):
         """
+        Buffered axis changed controller. Controls the buffer and handles update events for when the axes changed.
+        
+        Updated plotting will be after first reload (first time will be within plot limits, after that the limits will be buffered)
+        
         :param plot_function:
             function to use for creating image for plotting (return ndarray-like)
             plot_function gets called with (2D!) Xtest grid if replotting required
@@ -91,11 +95,13 @@ class BufferedAxisChangedController(AxisChangedController):
         """
         super(BufferedAxisChangedController, self).__init__(ax, update_lim=update_lim)
         self.plot_function = plot_function
-        xmin, xmax = self._x_lim # self._compute_buffered(*self._x_lim)
-        ymin, ymax = self._y_lim # self._compute_buffered(*self._y_lim)
+        xmin, ymin, xmax, ymax = plot_limits#self._x_lim # self._compute_buffered(*self._x_lim)
+        # imshow acts on the limits of the plot, this is why we need to override the limits here, to make sure the right plot limits are used:
+        self._x_lim = xmin, xmax
+        self._y_lim = ymin, ymax
         self.resolution = resolution
         self._not_init = False
-        self.view = self._init_view(self.ax, self.recompute_X(), xmin, xmax, ymin, ymax, **kwargs)
+        self.view = self._init_view(self.ax, self.recompute_X(buffered=False), xmin, xmax, ymin, ymax, **kwargs)
         self._not_init = True
 
     def update(self, ax):
@@ -111,14 +117,16 @@ class BufferedAxisChangedController(AxisChangedController):
     def update_view(self, view, X, xmin, xmax, ymin, ymax):
         raise NotImplementedError('update view given in here')
 
-    def get_grid(self):
-        xmin, xmax = self._compute_buffered(*self._x_lim)
-        ymin, ymax = self._compute_buffered(*self._y_lim)
+    def get_grid(self, buffered=True):
+        if buffered: comp = self._compute_buffered
+        else: comp = lambda a,b: (a,b)
+        xmin, xmax = comp(*self._x_lim)
+        ymin, ymax = comp(*self._y_lim)
         x, y = numpy.mgrid[xmin:xmax:1j * self.resolution, ymin:ymax:1j * self.resolution]
         return numpy.hstack((x.flatten()[:, None], y.flatten()[:, None]))
 
-    def recompute_X(self):
-        X = self.plot_function(self.get_grid())
+    def recompute_X(self, buffered=True):
+        X = self.plot_function(self.get_grid(buffered))
         if isinstance(X, (tuple, list)):
             for x in X:
                 x.shape = [self.resolution, self.resolution]
diff --git a/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/imshow_controller.py b/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/imshow_controller.py
index b473dd96..de1114a2 100644
--- a/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/imshow_controller.py
+++ b/GPy/plotting/matplot_dep/latent_space_visualizations/controllers/imshow_controller.py
@@ -9,7 +9,7 @@ import numpy
 
 
 class ImshowController(BufferedAxisChangedController):
-    def __init__(self, ax, plot_function, plot_limits, resolution=50, update_lim=.5, **kwargs):
+    def __init__(self, ax, plot_function, plot_limits, resolution=50, update_lim=.8, **kwargs):
         """
         :param plot_function:
             function to use for creating image for plotting (return ndarray-like)
diff --git a/GPy/plotting/matplot_dep/models_plots.py b/GPy/plotting/matplot_dep/models_plots.py
index cbb213b1..b626758f 100644
--- a/GPy/plotting/matplot_dep/models_plots.py
+++ b/GPy/plotting/matplot_dep/models_plots.py
@@ -123,6 +123,8 @@ def plot_fit(model, plot_limits=None, which_data_rows='all',
         #add inducing inputs (if a sparse model is used)
         if hasattr(model,"Z"):
             #Zu = model.Z[:,free_dims] * model._Xscale[:,free_dims] + model._Xoffset[:,free_dims]
+            if isinstance(model,SparseGPCoregionalizedRegression):
+                Z = Z[Z[:,-1] == Y_metadata['output_index'],:]
             Zu = Z[:,free_dims]
             z_height = ax.get_ylim()[0]
             plots['inducing_inputs'] = ax.plot(Zu, np.zeros_like(Zu) + z_height, 'r|', mew=1.5, markersize=12)
diff --git a/GPy/plotting/matplot_dep/visualize.py b/GPy/plotting/matplot_dep/visualize.py
index fb085f39..f8bcc9f9 100644
--- a/GPy/plotting/matplot_dep/visualize.py
+++ b/GPy/plotting/matplot_dep/visualize.py
@@ -4,6 +4,8 @@ import GPy
 import numpy as np
 import matplotlib as mpl
 import time
+from ...util.misc import param_to_array
+from GPy.core.parameterization.variational import VariationalPosterior
 try:
     import visual
     visual_available = True
@@ -72,12 +74,13 @@ class vector_show(matplotlib_show):
     """
     def __init__(self, vals, axes=None):
         matplotlib_show.__init__(self, vals, axes)
-        self.handle = self.axes.plot(np.arange(0, len(vals))[:, None], self.vals.T)[0]
+        self.handle = self.axes.plot(np.arange(0, len(vals))[:, None], self.vals)
 
     def modify(self, vals):
         self.vals = vals.copy()
-        xdata, ydata = self.handle.get_data()
-        self.handle.set_data(xdata, self.vals.T)
+        for handle, vals in zip(self.handle, self.vals.T):
+            xdata, ydata = handle.get_data()
+            handle.set_data(xdata, vals)
         self.axes.figure.canvas.draw()
 
 
@@ -91,8 +94,12 @@ class lvm(matplotlib_show):
         :param latent_axes: the axes where the latent visualization should be plotted.
         """
         if vals == None:
-            vals = model.X[0]
-
+            if isinstance(model.X, VariationalPosterior):
+                vals = param_to_array(model.X.mean)
+            else:
+                vals = param_to_array(model.X)
+        
+        vals = param_to_array(vals)
         matplotlib_show.__init__(self, vals, axes=latent_axes)
 
         if isinstance(latent_axes,mpl.axes.Axes):
diff --git a/GPy/testing/kernel_tests.py b/GPy/testing/kernel_tests.py
index 3eef6768..9ed218d8 100644
--- a/GPy/testing/kernel_tests.py
+++ b/GPy/testing/kernel_tests.py
@@ -152,7 +152,12 @@ def check_kernel_gradient_functions(kern, X=None, X2=None, output_ind=None, verb
 
     if verbose:
         print("Checking gradients of Kdiag(X) wrt theta.")
-    result = Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=verbose)
+    try:
+        result = Kern_check_dKdiag_dtheta(kern, X=X).checkgrad(verbose=verbose)
+    except NotImplementedError:
+        result=True
+        if verbose:
+            print("update_gradients_diag not implemented for " + kern.name)
     if result and verbose:
         print("Check passed.")
     if not result:
@@ -240,9 +245,22 @@ class KernelGradientTestsContinuous(unittest.TestCase):
 
     def test_Add(self):
         k = GPy.kern.Matern32(2, active_dims=[2,3]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
+        k += GPy.kern.Matern32(2, active_dims=[2,3]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
         k.randomize()
         self.assertTrue(check_kernel_gradient_functions(k, X=self.X, X2=self.X2, verbose=verbose))
 
+    def test_Add_dims(self):
+        k = GPy.kern.Matern32(2, active_dims=[2,self.D]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
+        k.randomize()
+        self.assertRaises(AssertionError, k.K, self.X)
+        k = GPy.kern.Matern32(2, active_dims=[2,self.D-1]) + GPy.kern.RBF(2, active_dims=[0,4]) + GPy.kern.Linear(self.D)
+        k.randomize()
+        # assert it runs:
+        try:
+            k.K(self.X)
+        except AssertionError:
+            raise AssertionError, "k.K(X) should run on self.D-1 dimension"
+
     def test_Matern52(self):
         k = GPy.kern.Matern52(self.D)
         k.randomize()
@@ -329,17 +347,11 @@ class KernelTestsNonContinuous(unittest.TestCase):
         kern = GPy.kern.IndependentOutputs(k, -1, name='ind_split')
         self.assertTrue(check_kernel_gradient_functions(kern, X=self.X, X2=self.X2, verbose=verbose, fixed_X_dims=-1))
 
-class test_ODE_UY(unittest.TestCase):
-    def setUp(self):
-        self.k = GPy.kern.ODE_UY(2)
-        self.X = np.random.randn(50,2)
-        self.X[:,1] = np.random.randint(0,2,50)
-        i = np.argsort(X[:,1])
-        self.X = self.X[i]
-        self.Y = np.random.randn(50, 1)
-    def checkgrad(self):
-        m = GPy.models.GPRegression(X,Y,kernel=k)
-        self.assertTrue(m.checkgrad())
+    def test_ODE_UY(self):
+        kern = GPy.kern.ODE_UY(2, active_dims=[0, self.D])
+        X = self.X[self.X[:,-1]!=2]
+        X2 = self.X2[self.X2[:,-1]!=2]
+        self.assertTrue(check_kernel_gradient_functions(kern, X=X, X2=X2, verbose=verbose, fixed_X_dims=-1))
 
 
 if __name__ == "__main__":
diff --git a/GPy/testing/parameterized_tests.py b/GPy/testing/parameterized_tests.py
index cd5127c8..dc59449f 100644
--- a/GPy/testing/parameterized_tests.py
+++ b/GPy/testing/parameterized_tests.py
@@ -121,6 +121,11 @@ class ParameterizedTest(unittest.TestCase):
         self.test1.randomize()
         self.assertEqual(val, self.white.variance)
 
+    def test_randomize(self):
+        ps = self.test1.param.view(np.ndarray).copy()
+        self.test1.param.randomize()
+        self.assertFalse(np.all(ps==self.test1.param))
+
     def test_fixing_randomize_parameter_handling(self):
         self.rbf.fix(warning=True)
         val = float(self.rbf.variance)
diff --git a/GPy/util/multioutput.py b/GPy/util/multioutput.py
index d9e8b704..b769d6a5 100644
--- a/GPy/util/multioutput.py
+++ b/GPy/util/multioutput.py
@@ -56,8 +56,6 @@ def ICM(input_dim, num_outputs, kernel, W_rank=1,W=None,kappa=None,name='X'):
         warnings.warn("kernel's input dimension overwritten to fit input_dim parameter.")
 
     K = kernel.prod(GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B'),name=name)
-    #K = kernel * GPy.kern.Coregionalize(1, num_outputs, active_dims=[input_dim], rank=W_rank,W=W,kappa=kappa,name='B')
-    #K = kernel ** GPy.kern.Coregionalize(input_dim, num_outputs,W_rank,W,kappa, name= 'B')
     K['.*variance'] = 1.
     K['.*variance'].fix()
     return K