diff --git a/GPy/examples/dimensionality_reduction.py b/GPy/examples/dimensionality_reduction.py
index 9d6686f9..0d72949a 100644
--- a/GPy/examples/dimensionality_reduction.py
+++ b/GPy/examples/dimensionality_reduction.py
@@ -343,6 +343,29 @@ def bgplvm_simulation(optimize=True, verbose=1,
         m.kern.plot_ARD('BGPLVM Simulation ARD Parameters')
     return m
 
+def gplvm_simulation(optimize=True, verbose=1,
+                      plot=True, plot_sim=False,
+                      max_iters=2e4,
+                      ):
+    from GPy import kern
+    from GPy.models import GPLVM
+
+    D1, D2, D3, N, num_inducing, Q = 13, 5, 8, 45, 3, 9
+    _, _, Ylist = _simulate_matern(D1, D2, D3, N, num_inducing, plot_sim)
+    Y = Ylist[0]
+    k = kern.Linear(Q, ARD=True)  # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
+    # k = kern.RBF(Q, ARD=True, lengthscale=10.)
+    m = GPLVM(Y, Q, init="PCA", kernel=k)
+    m.likelihood.variance = .1
+
+    if optimize:
+        print("Optimizing model:")
+        m.optimize('bfgs', messages=verbose, max_iters=max_iters,
+                   gtol=.05)
+    if plot:
+        m.X.plot("BGPLVM Latent Space 1D")
+        m.kern.plot_ARD('BGPLVM Simulation ARD Parameters')
+    return m
 def ssgplvm_simulation(optimize=True, verbose=1,
                       plot=True, plot_sim=False,
                       max_iters=2e4, useGPU=False
diff --git a/GPy/kern/_src/linear.py b/GPy/kern/_src/linear.py
index 0a582ac8..b743fdf0 100644
--- a/GPy/kern/_src/linear.py
+++ b/GPy/kern/_src/linear.py
@@ -109,6 +109,9 @@ class Linear(Kern):
     def gradients_X_diag(self, dL_dKdiag, X):
         return 2.*self.variances*dL_dKdiag[:,None]*X
 
+    def gradients_XX_diag(self, dL_dKdiag, X):
+        return 2*np.ones(X.shape)*self.variances
+
     def input_sensitivity(self, summarize=True):
         return np.ones(self.input_dim) * self.variances
 
diff --git a/GPy/models/bayesian_gplvm.py b/GPy/models/bayesian_gplvm.py
index 88123227..fd02cb3e 100644
--- a/GPy/models/bayesian_gplvm.py
+++ b/GPy/models/bayesian_gplvm.py
@@ -99,137 +99,3 @@ class BayesianGPLVM(SparseGP_MPI):
         self.variational_prior.update_gradients_KL(self.X)
         self._Xgrad = self.X.gradient.copy()
 
-        #super(BayesianGPLVM, self).parameters_changed()
-        #self._log_marginal_likelihood -= self.variational_prior.KL_divergence(self.X)
-
-        #self.X.mean.gradient, self.X.variance.gradient = self.kern.gradients_qX_expectations(variational_posterior=self.X, Z=self.Z, dL_dpsi0=self.grad_dict['dL_dpsi0'], dL_dpsi1=self.grad_dict['dL_dpsi1'], dL_dpsi2=self.grad_dict['dL_dpsi2'])
-
-        # This is testing code -------------------------
-#         i = np.random.randint(self.X.shape[0])
-#         X_ = self.X.mean
-#         which = np.sqrt(((X_ - X_[i:i+1])**2).sum(1)).argsort()>(max(0, self.X.shape[0]-51))
-#         _, _, grad_dict = self.inference_method.inference(self.kern, self.X[which], self.Z, self.likelihood, self.Y[which], self.Y_metadata)
-#         grad = self.kern.gradients_qX_expectations(variational_posterior=self.X[which], Z=self.Z, dL_dpsi0=grad_dict['dL_dpsi0'], dL_dpsi1=grad_dict['dL_dpsi1'], dL_dpsi2=grad_dict['dL_dpsi2'])
-#
-#         self.X.mean.gradient[:] = 0
-#         self.X.variance.gradient[:] = 0
-#         self.X.mean.gradient[which] = grad[0]
-#         self.X.variance.gradient[which] = grad[1]
-
-        # update for the KL divergence
-#         self.variational_prior.update_gradients_KL(self.X, which)
-        # -----------------------------------------------
-
-        # update for the KL divergence
-        #self.variational_prior.update_gradients_KL(self.X)
-
-    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, predict_kwargs={}, imshow_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, labels, which_indices,
-                resolution, ax, marker, s,
-                fignum, plot_inducing, legend,
-                plot_limits, aspect, updates, predict_kwargs, imshow_kwargs)
-
-    def do_test_latents(self, Y):
-        """
-        Compute the latent representation for a set of new points Y
-
-        Notes:
-        This will only work with a univariate Gaussian likelihood (for now)
-        """
-        N_test = Y.shape[0]
-        input_dim = self.Z.shape[1]
-
-        means = np.zeros((N_test, input_dim))
-        covars = np.zeros((N_test, input_dim))
-
-        dpsi0 = -0.5 * self.input_dim / self.likelihood.variance
-        dpsi2 = self.grad_dict['dL_dpsi2'][0][None, :, :] # TODO: this may change if we ignore het. likelihoods
-        V = Y/self.likelihood.variance
-
-        #compute CPsi1V
-        #if self.Cpsi1V is None:
-        #    psi1V = np.dot(self.psi1.T, self.likelihood.V)
-        #    tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
-        #    tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1)
-        #    self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1)
-
-        dpsi1 = np.dot(self.posterior.woodbury_vector, V.T)
-
-        #start = np.zeros(self.input_dim * 2)
-
-
-        from scipy.optimize import minimize
-
-        for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]):
-            args = (input_dim, self.kern.copy(), self.Z, dpsi0, dpsi1_n.T, dpsi2)
-            res = minimize(latent_cost_and_grad, jac=True, x0=np.hstack((means[n], covars[n])), args=args, method='BFGS')
-            xopt = res.x
-            mu, log_S = xopt.reshape(2, 1, -1)
-            means[n] = mu[0].copy()
-            covars[n] = np.exp(log_S[0]).copy()
-
-        X = NormalPosterior(means, covars)
-
-        return X
-
-    def dmu_dX(self, Xnew):
-        """
-        Calculate the gradient of the prediction at Xnew w.r.t Xnew.
-        """
-        dmu_dX = np.zeros_like(Xnew)
-        for i in range(self.Z.shape[0]):
-            dmu_dX += self.kern.gradients_X(self.grad_dict['dL_dpsi1'][i:i + 1, :], Xnew, self.Z[i:i + 1, :])
-        return dmu_dX
-
-    def dmu_dXnew(self, Xnew):
-        """
-        Individual gradient of prediction at Xnew w.r.t. each sample in Xnew
-        """
-        gradients_X = np.zeros((Xnew.shape[0], self.num_inducing))
-        ones = np.ones((1, 1))
-        for i in range(self.Z.shape[0]):
-            gradients_X[:, i] = self.kern.gradients_X(ones, Xnew, self.Z[i:i + 1, :]).sum(-1)
-        return np.dot(gradients_X, self.grad_dict['dL_dpsi1'])
-
-    def plot_steepest_gradient_map(self, *args, ** kwargs):
-        """
-        See GPy.plotting.matplot_dep.dim_reduction_plots.plot_steepest_gradient_map
-        """
-        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_steepest_gradient_map(self,*args,**kwargs)
-
-
-def latent_cost_and_grad(mu_S, input_dim, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
-    """
-    objective function for fitting the latent variables for test points
-    (negative log-likelihood: should be minimised!)
-    """
-    mu = mu_S[:input_dim][None]
-    log_S = mu_S[input_dim:][None]
-    S = np.exp(log_S)
-
-    X = NormalPosterior(mu, S)
-
-    psi0 = kern.psi0(Z, X)
-    psi1 = kern.psi1(Z, X)
-    psi2 = kern.psi2(Z, X)
-
-    lik = dL_dpsi0 * psi0.sum() + np.einsum('ij,kj->...', dL_dpsi1, psi1) + np.einsum('ijk,lkj->...', dL_dpsi2, psi2) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
-
-    dLdmu, dLdS = kern.gradients_qX_expectations(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, X)
-    dmu = dLdmu - mu
-    # dS = S0 + S1 + S2 -0.5 + .5/S
-    dlnS = S * (dLdS - 0.5) + .5
-
-    return -lik, -np.hstack((dmu.flatten(), dlnS.flatten()))
diff --git a/GPy/models/bayesian_gplvm_minibatch.py b/GPy/models/bayesian_gplvm_minibatch.py
index 175c3e56..3ef43753 100644
--- a/GPy/models/bayesian_gplvm_minibatch.py
+++ b/GPy/models/bayesian_gplvm_minibatch.py
@@ -128,115 +128,4 @@ class BayesianGPLVMMiniBatch(SparseGPMiniBatch):
                 d = self.output_dim
                 self._log_marginal_likelihood -= kl_fctr*self.variational_prior.KL_divergence(self.X)*self.stochastics.batchsize/d
 
-        self._Xgrad = self.X.gradient.copy()
-
-    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, predict_kwargs={}, imshow_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, labels, which_indices,
-                resolution, ax, marker, s,
-                fignum, plot_inducing, legend,
-                plot_limits, aspect, updates, predict_kwargs, imshow_kwargs)
-
-    def do_test_latents(self, Y):
-        """
-        Compute the latent representation for a set of new points Y
-
-        Notes:
-        This will only work with a univariate Gaussian likelihood (for now)
-        """
-        N_test = Y.shape[0]
-        input_dim = self.Z.shape[1]
-
-        means = np.zeros((N_test, input_dim))
-        covars = np.zeros((N_test, input_dim))
-
-        dpsi0 = -0.5 * self.input_dim / self.likelihood.variance
-        dpsi2 = self.grad_dict['dL_dpsi2'][0][None, :, :] # TODO: this may change if we ignore het. likelihoods
-        V = Y/self.likelihood.variance
-
-        #compute CPsi1V
-        #if self.Cpsi1V is None:
-        #    psi1V = np.dot(self.psi1.T, self.likelihood.V)
-        #    tmp, _ = linalg.dtrtrs(self._Lm, np.asfortranarray(psi1V), lower=1, trans=0)
-        #    tmp, _ = linalg.dpotrs(self.LB, tmp, lower=1)
-        #    self.Cpsi1V, _ = linalg.dtrtrs(self._Lm, tmp, lower=1, trans=1)
-
-        dpsi1 = np.dot(self.posterior.woodbury_vector, V.T)
-
-        #start = np.zeros(self.input_dim * 2)
-
-
-        from scipy.optimize import minimize
-
-        for n, dpsi1_n in enumerate(dpsi1.T[:, :, None]):
-            args = (input_dim, self.kern.copy(), self.Z, dpsi0, dpsi1_n.T, dpsi2)
-            res = minimize(latent_cost_and_grad, jac=True, x0=np.hstack((means[n], covars[n])), args=args, method='BFGS')
-            xopt = res.x
-            mu, log_S = xopt.reshape(2, 1, -1)
-            means[n] = mu[0].copy()
-            covars[n] = np.exp(log_S[0]).copy()
-
-        X = NormalPosterior(means, covars)
-
-        return X
-
-    def dmu_dX(self, Xnew):
-        """
-        Calculate the gradient of the prediction at Xnew w.r.t Xnew.
-        """
-        dmu_dX = np.zeros_like(Xnew)
-        for i in range(self.Z.shape[0]):
-            dmu_dX += self.kern.gradients_X(self.grad_dict['dL_dpsi1'][i:i + 1, :], Xnew, self.Z[i:i + 1, :])
-        return dmu_dX
-
-    def dmu_dXnew(self, Xnew):
-        """
-        Individual gradient of prediction at Xnew w.r.t. each sample in Xnew
-        """
-        gradients_X = np.zeros((Xnew.shape[0], self.num_inducing))
-        ones = np.ones((1, 1))
-        for i in range(self.Z.shape[0]):
-            gradients_X[:, i] = self.kern.gradients_X(ones, Xnew, self.Z[i:i + 1, :]).sum(-1)
-        return np.dot(gradients_X, self.grad_dict['dL_dpsi1'])
-
-    def plot_steepest_gradient_map(self, *args, ** kwargs):
-        """
-        See GPy.plotting.matplot_dep.dim_reduction_plots.plot_steepest_gradient_map
-        """
-        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_steepest_gradient_map(self,*args,**kwargs)
-
-
-def latent_cost_and_grad(mu_S, input_dim, kern, Z, dL_dpsi0, dL_dpsi1, dL_dpsi2):
-    """
-    objective function for fitting the latent variables for test points
-    (negative log-likelihood: should be minimised!)
-    """
-    mu = mu_S[:input_dim][None]
-    log_S = mu_S[input_dim:][None]
-    S = np.exp(log_S)
-
-    X = NormalPosterior(mu, S)
-
-    psi0 = kern.psi0(Z, X)
-    psi1 = kern.psi1(Z, X)
-    psi2 = kern.psi2(Z, X)
-
-    lik = dL_dpsi0 * psi0.sum() + np.einsum('ij,kj->...', dL_dpsi1, psi1) + np.einsum('ijk,lkj->...', dL_dpsi2, psi2) - 0.5 * np.sum(np.square(mu) + S) + 0.5 * np.sum(log_S)
-
-    dLdmu, dLdS = kern.gradients_qX_expectations(dL_dpsi0, dL_dpsi1, dL_dpsi2, Z, X)
-    dmu = dLdmu - mu
-    # dS = S0 + S1 + S2 -0.5 + .5/S
-    dlnS = S * (dLdS - 0.5) + .5
-
-    return -lik, -np.hstack((dmu.flatten(), dlnS.flatten()))
+        self._Xgrad = self.X.gradient.copy()
\ No newline at end of file
diff --git a/GPy/plotting/__init__.py b/GPy/plotting/__init__.py
index 0b6c9c89..b90bf938 100644
--- a/GPy/plotting/__init__.py
+++ b/GPy/plotting/__init__.py
@@ -36,14 +36,20 @@ if config.get('plotting', 'library') is not 'none':
     GP.plot_samples = gpy_plot.gp_plots.plot_samples
     GP.plot = gpy_plot.gp_plots.plot
     GP.plot_f = gpy_plot.gp_plots.plot_f
-    GP.plot_magnificaion = gpy_plot.latent_plots.plot_magnification
+    GP.plot_magnification = gpy_plot.latent_plots.plot_magnification
     
     from ..core import SparseGP
     SparseGP.plot_inducing = gpy_plot.data_plots.plot_inducing
     
     from ..models import GPLVM
+    from ..models import BayesianGPLVM
+    from ..models import bayesian_gplvm_minibatch
     GPLVM.plot_prediction_fit = gpy_plot.latent_plots.plot_prediction_fit
-    #GPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
+    GPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
+    BayesianGPLVM.plot_latent = gpy_plot.latent_plots.plot_latent
+    BayesianGPLVM.plot_prediction_fit = gpy_plot.latent_plots.plot_prediction_fit
+    bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_latent = gpy_plot.latent_plots.plot_latent
+    bayesian_gplvm_minibatch.BayesianGPLVMMiniBatch.plot_prediction_fit = gpy_plot.latent_plots.plot_prediction_fit
     
     from ..kern import Kern
     #Kern.plot_covariance = gpy_plot.kern_plots.plot_kern
diff --git a/GPy/plotting/gpy_plot/__init__.py b/GPy/plotting/gpy_plot/__init__.py
index cec60bcf..6a4b7391 100644
--- a/GPy/plotting/gpy_plot/__init__.py
+++ b/GPy/plotting/gpy_plot/__init__.py
@@ -1,3 +1,3 @@
 from .. import plotting_library as pl
-from . import data_plots, gp_plots, latent_plots
+from . import data_plots, gp_plots, latent_plots, kernel_plots, plot_util
 
diff --git a/GPy/plotting/gpy_plot/latent_plots.py b/GPy/plotting/gpy_plot/latent_plots.py
index 23235c44..bd19fb29 100644
--- a/GPy/plotting/gpy_plot/latent_plots.py
+++ b/GPy/plotting/gpy_plot/latent_plots.py
@@ -33,6 +33,15 @@ from .plot_util import get_x_y_var, get_free_dims, get_which_data_ycols,\
     get_which_data_rows, update_not_existing_kwargs, helper_predict_with_model,\
     helper_for_plot_data
 import itertools
+from GPy.plotting.gpy_plot.plot_util import scatter_label_generator, subsample_X
+
+def _wait_for_updates(view, updates):
+    if updates:
+        clear = raw_input('yes or enter to deactivate updates - otherwise still do updates - use plots[imshow].deactivate() to clear')
+        if clear.lower() in 'yes' or clear == '':
+            view.deactivate()
+    else:
+        view.deactivate()
 
 def plot_prediction_fit(self, plot_limits=None,
         which_data_rows='all', which_data_ycols='all', 
@@ -88,7 +97,7 @@ def _plot_prediction_fit(self, canvas, plot_limits=None,
                 scatter_kwargs = {}
             update_not_existing_kwargs(scatter_kwargs, pl.defaults.data_y_1d)  # @UndefinedVariable
             plots['output'] = pl.scatter(canvas, Y[rows, ycols[0]], Y[rows, ycols[1]],
-                                      c=X[rows, free_dims[0]],
+                                      color=X[rows, free_dims[0]],
                                       **scatter_kwargs)
             if predict_kw is None:
                 predict_kw = {}
@@ -108,7 +117,9 @@ def plot_magnification(self, labels=None, which_indices=None,
                 plot_limits=None,
                 updates=False, 
                 mean=True, covariance=True, 
-                kern=None, marker='<>^vsd', imshow_kwargs=None, **kwargs):
+                kern=None, marker='<>^vsd', 
+                num_samples=1000,
+                imshow_kwargs=None, **kwargs):
     """
     Plot the magnification factor of the GP on the inputs. This is the 
     density of the GP as a gray scale.
@@ -124,95 +135,23 @@ def plot_magnification(self, labels=None, which_indices=None,
     :param bool covariance: use the covariance of the Wishart embedding for the magnification factor
     :param :py:class:`~GPy.kern.Kern` kern: the kernel to use for prediction
     :param str marker: markers to use - cycle if more labels then markers are given
+    :param int num_samples: the number of samples to plot maximally. We do a stratified subsample from the labels, if the number of samples (in X) is higher then num_samples. 
     :param imshow_kwargs: the kwargs for the imshow (magnification factor)
     :param kwargs: the kwargs for the scatter plots
     """
     input_1, input_2 = self.get_most_significant_input_dimensions(which_indices)
 
-    #fethch the data points X that we'd like to plot
-    X, _, _ = get_x_y_var(self)
-
-    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:
-            try:
-                xmin, xmax = plot_limits
-                ymin, ymax = xmin, xmax
-            except (TypeError, ValueError) as e:
-                raise e.__class__("Wrong plot limits: {} given -> need (xmin, xmax, ymin, ymax)".format(plot_limits))
-    xlim = (xmin, xmax)
-    ylim = (ymin, ymax)
-
     from .. import Tango
     Tango.reset()
     
     if labels is None:
         labels = np.ones(self.num_data)
-
-    if X.shape[0] > 1000:
-        print("Warning: subsampling X, as it has more samples then 1000. X.shape={!s}".format(X.shape))
-        subsample = np.random.choice(X.shape[0], size=1000, replace=False)
-        X = X[subsample]
-        labels = labels[subsample]
-        #=======================================================================
-        #     <<<WORK IN PROGRESS>>>
-        #     <<<DO NOT DELETE>>>
-        #     plt.close('all')
-        #     fig, ax = plt.subplots(1,1)
-        #     from GPy.plotting.matplot_dep.dim_reduction_plots import most_significant_input_dimensions
-        #     import matplotlib.patches as mpatches
-        #     i1, i2 = most_significant_input_dimensions(m, None)
-        #     xmin, xmax = 100, -100
-        #     ymin, ymax = 100, -100
-        #     legend_handles = []
-        #
-        #     X = m.X.mean[:, [i1, i2]]
-        #     X = m.X.variance[:, [i1, i2]]
-        #
-        #     xmin = X[:,0].min(); xmax = X[:,0].max()
-        #     ymin = X[:,1].min(); ymax = X[:,1].max()
-        #     range_ = [[xmin, xmax], [ymin, ymax]]
-        #     ul = np.unique(labels)
-        #
-        #     for i, l in enumerate(ul):
-        #         #cdict = dict(red  =[(0., colors[i][0], colors[i][0]), (1., colors[i][0], colors[i][0])],
-        #         #             green=[(0., colors[i][0], colors[i][1]), (1., colors[i][1], colors[i][1])],
-        #         #             blue =[(0., colors[i][0], colors[i][2]), (1., colors[i][2], colors[i][2])],
-        #         #             alpha=[(0., 0., .0), (.5, .5, .5), (1., .5, .5)])
-        #         #cmap = LinearSegmentedColormap('{}'.format(l), cdict)
-        #         cmap = LinearSegmentedColormap.from_list('cmap_{}'.format(str(l)), [colors[i], colors[i]], 255)
-        #         cmap._init()
-        #         #alphas = .5*(1+scipy.special.erf(np.linspace(-2,2, cmap.N+3)))#np.log(np.linspace(np.exp(0), np.exp(1.), cmap.N+3))
-        #         alphas = (scipy.special.erf(np.linspace(0,2.4, cmap.N+3)))#np.log(np.linspace(np.exp(0), np.exp(1.), cmap.N+3))
-        #         cmap._lut[:, -1] = alphas
-        #         print l
-        #         x, y = X[labels==l].T
-        #
-        #         heatmap, xedges, yedges = np.histogram2d(x, y, bins=300, range=range_)
-        #         #heatmap, xedges, yedges = np.histogram2d(x, y, bins=100)
-        #
-        #         im = ax.imshow(heatmap, extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap=cmap, aspect='auto', interpolation='nearest', label=str(l))
-        #         legend_handles.append(mpatches.Patch(color=colors[i], label=l))
-        #     ax.set_xlim(xmin, xmax)
-        #     ax.set_ylim(ymin, ymax)
-        #     plt.legend(legend_handles, [l.get_label() for l in legend_handles])
-        #     plt.draw()
-        #     plt.show()
-        #=======================================================================
-    
+        legend = False # No legend if there is no labels given
     
     canvas, kwargs = pl.get_new_canvas(xlabel='latent dimension %i' % input_1, ylabel='latent dimension %i' % input_2, **kwargs)
 
-    _, _, _, _, _, Xgrid, _, _, _, _, resolution = helper_for_plot_data(self, ((xmin, ymin), (xmax, ymax)), (input_1, input_2), None, resolution)
+    X, _, _, _, _, Xgrid, _, _, xmin, xmax, resolution = helper_for_plot_data(self, plot_limits, (input_1, input_2), None, resolution)
+    X, labels = subsample_X(X, labels)
     
     def plot_function(x):
         Xtest_full = np.zeros((x.shape[0], X.shape[1]))
@@ -223,44 +162,79 @@ def plot_magnification(self, labels=None, which_indices=None,
     imshow_kwargs = update_not_existing_kwargs(imshow_kwargs, pl.defaults.magnification)
     Y = plot_function(Xgrid[:, [input_1, input_2]]).reshape(resolution, resolution).T[::-1, :]
     view = pl.imshow(canvas, Y, 
-                     (xmin, ymin, xmax, ymax), 
+                     (xmin[0], xmin[1], xmax[1], xmax[1]), 
                      None, plot_function, resolution,
                      vmin=Y.min(), vmax=Y.max(), 
                      **imshow_kwargs)
-    
-    # make sure labels are in order of input:
-    ulabels = []
-    for lab in labels:
-        if not lab in ulabels:
-            ulabels.append(lab)
 
-    marker = itertools.cycle(list(marker))
-    scatters = []
-
-    for ul in ulabels:
-        if type(ul) is np.string_:
-            this_label = ul
-        elif type(ul) is np.int64:
-            this_label = 'class %i' % ul
-        else:
-            this_label = unicode(ul)
-        m = marker.next()
-
-        index = np.nonzero(labels == ul)[0]
-        if self.input_dim == 1:
-            x = X[index, input_1]
-            y = np.zeros(index.size)
-        else:
-            x = X[index, input_1]
-            y = X[index, input_2]
+    scatters = []    
+    for x, y, this_label, _, m in scatter_label_generator(labels, X, input_1, input_2, marker):
         update_not_existing_kwargs(kwargs, pl.defaults.latent_scatter)
         scatters.append(pl.scatter(canvas, x, y, marker=m, color=Tango.nextMedium(), label=this_label, **kwargs))
     
-    plots = pl.show_canvas(canvas, dict(scatter=scatters, imshow=view), legend=legend, xlim=xlim, ylim=ylim)
-    if updates:
-        clear = raw_input('yes or enter to deactivate updates - otherwise still do updates - use plots[imshow].deactivate() to clear')
-        if clear.lower() in 'yes' or clear == '':
-            view.deactivate()
-    else:
-        view.deactivate()
+    plots = pl.show_canvas(canvas, dict(scatter=scatters, imshow=view), legend=legend, xlim=(xmin[0], xmax[0]), ylim=(xmin[1], xmax[1]))
+    _wait_for_updates(view, updates)
+    return plots
+
+
+def plot_latent(self, labels=None, which_indices=None,
+                resolution=60, legend=True,
+                plot_limits=None,
+                updates=False, 
+                kern=None, marker='<>^vsd', 
+                num_samples=1000,
+                imshow_kwargs=None, **kwargs):
+    """
+    Plot the latent space of the GP on the inputs. This is the 
+    density of the GP posterior as a grey scale and the 
+    scatter plot of the input dimemsions selected by which_indices.
+    
+    :param array-like labels: a label for each data point (row) of the inputs
+    :param (int, int) which_indices: which input dimensions to plot against each other
+    :param int resolution: the resolution at which we predict the magnification factor
+    :param bool legend: whether to plot the legend on the figure
+    :param plot_limits: the plot limits for the plot
+    :type plot_limits: (xmin, xmax, ymin, ymax) or ((xmin, xmax), (ymin, ymax))
+    :param bool updates: if possible, make interactive updates using the specific library you are using
+    :param :py:class:`~GPy.kern.Kern` kern: the kernel to use for prediction
+    :param str marker: markers to use - cycle if more labels then markers are given
+    :param int num_samples: the number of samples to plot maximally. We do a stratified subsample from the labels, if the number of samples (in X) is higher then num_samples. 
+    :param imshow_kwargs: the kwargs for the imshow (magnification factor)
+    :param kwargs: the kwargs for the scatter plots
+    """
+    input_1, input_2 = self.get_most_significant_input_dimensions(which_indices)
+
+    from .. import Tango
+    Tango.reset()
+    
+    if labels is None:
+        labels = np.ones(self.num_data)
+        legend = False # No legend if there is no labels given
+    
+    canvas, kwargs = pl.get_new_canvas(xlabel='latent dimension %i' % input_1, ylabel='latent dimension %i' % input_2, **kwargs)
+
+    X, _, _, _, _, Xgrid, _, _, xmin, xmax, resolution = helper_for_plot_data(self, plot_limits, (input_1, input_2), None, resolution)
+    X, labels = subsample_X(X, labels)
+    
+    def plot_function(x):
+        Xtest_full = np.zeros((x.shape[0], X.shape[1]))
+        Xtest_full[:, [input_1, input_2]] = x
+        mf = np.log(self.predict(Xtest_full, kern=kern)[1])
+        return mf
+
+    imshow_kwargs = update_not_existing_kwargs(imshow_kwargs, pl.defaults.latent)
+    Y = plot_function(Xgrid[:, [input_1, input_2]]).reshape(resolution, resolution).T[::-1, :]
+    view = pl.imshow(canvas, Y, 
+                     (xmin[0], xmin[1], xmax[1], xmax[1]), 
+                     None, plot_function, resolution,
+                     vmin=Y.min(), vmax=Y.max(), 
+                     **imshow_kwargs)
+
+    scatters = []    
+    for x, y, this_label, _, m in scatter_label_generator(labels, X, input_1, input_2, marker):
+        update_not_existing_kwargs(kwargs, pl.defaults.latent_scatter)
+        scatters.append(pl.scatter(canvas, x, y, marker=m, color=Tango.nextMedium(), label=this_label, **kwargs))
+    
+    plots = pl.show_canvas(canvas, dict(scatter=scatters, imshow=view), legend=legend, xlim=(xmin[0], xmax[0]), ylim=(xmin[1], xmax[1]))
+    _wait_for_updates(view, updates)
     return plots
diff --git a/GPy/plotting/gpy_plot/plot_util.py b/GPy/plotting/gpy_plot/plot_util.py
index 0defc9df..4f3fdaf2 100644
--- a/GPy/plotting/gpy_plot/plot_util.py
+++ b/GPy/plotting/gpy_plot/plot_util.py
@@ -30,6 +30,7 @@
 
 import numpy as np
 from scipy import sparse
+import itertools
 
 def helper_predict_with_model(self, Xgrid, plot_raw, apply_link, percentiles, which_data_ycols, predict_kw, samples=0):
     """
@@ -117,6 +118,102 @@ def helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resoluti
             Xgrid[:,i] = v    
     return X, Xvar, Y, fixed_dims, free_dims, Xgrid, x, y, xmin, xmax, resolution
 
+def scatter_label_generator(labels, X, input_1, input_2=None, marker=None):
+    ulabels = []
+    for lab in labels:
+        if not lab in ulabels:
+            ulabels.append(lab)
+
+    if marker is not None:
+        marker = itertools.cycle(list(marker))    
+    else:
+        m = None
+    
+    for ul in ulabels:
+        if type(ul) is np.string_:
+            this_label = ul
+        elif type(ul) is np.int64:
+            this_label = 'class %i' % ul
+        else:
+            this_label = unicode(ul)
+        
+        if marker is not None:
+            m = marker.next()
+
+        index = np.nonzero(labels == ul)[0]
+        
+        if input_2 is None:
+            x = X[index, input_1]
+            y = np.zeros(index.size)
+        else:
+            x = X[index, input_1]
+            y = X[index, input_2]
+        yield x, y, this_label, index, m
+
+def subsample_X(X, labels, num_samples=1000):
+    """
+    Stratified subsampling if labels are given. 
+    This means due to rounding errors you might get a little differences between the 
+    num_samples and the returned subsampled X.
+    """
+    if X.shape[0] > num_samples:
+        print("Warning: subsampling X, as it has more samples then 1000. X.shape={!s}".format(X.shape))
+        if labels is not None:
+            subsample = []
+            for _, _, _, index, _ in scatter_label_generator(labels, X, 0):
+                subsample.append(np.random.choice(index, size=max(2, int(index.size*(float(num_samples)/X.shape[0]))), replace=False))
+            subsample = np.hstack(subsample)
+        else:
+            subsample = np.random.choice(X.shape[0], size=1000, replace=False)
+        X = X[subsample]
+        labels = labels[subsample]
+        #=======================================================================
+        #     <<<WORK IN PROGRESS>>>
+        #     <<<DO NOT DELETE>>>
+        #     plt.close('all')
+        #     fig, ax = plt.subplots(1,1)
+        #     from GPy.plotting.matplot_dep.dim_reduction_plots import most_significant_input_dimensions
+        #     import matplotlib.patches as mpatches
+        #     i1, i2 = most_significant_input_dimensions(m, None)
+        #     xmin, xmax = 100, -100
+        #     ymin, ymax = 100, -100
+        #     legend_handles = []
+        #
+        #     X = m.X.mean[:, [i1, i2]]
+        #     X = m.X.variance[:, [i1, i2]]
+        #
+        #     xmin = X[:,0].min(); xmax = X[:,0].max()
+        #     ymin = X[:,1].min(); ymax = X[:,1].max()
+        #     range_ = [[xmin, xmax], [ymin, ymax]]
+        #     ul = np.unique(labels)
+        #
+        #     for i, l in enumerate(ul):
+        #         #cdict = dict(red  =[(0., colors[i][0], colors[i][0]), (1., colors[i][0], colors[i][0])],
+        #         #             green=[(0., colors[i][0], colors[i][1]), (1., colors[i][1], colors[i][1])],
+        #         #             blue =[(0., colors[i][0], colors[i][2]), (1., colors[i][2], colors[i][2])],
+        #         #             alpha=[(0., 0., .0), (.5, .5, .5), (1., .5, .5)])
+        #         #cmap = LinearSegmentedColormap('{}'.format(l), cdict)
+        #         cmap = LinearSegmentedColormap.from_list('cmap_{}'.format(str(l)), [colors[i], colors[i]], 255)
+        #         cmap._init()
+        #         #alphas = .5*(1+scipy.special.erf(np.linspace(-2,2, cmap.N+3)))#np.log(np.linspace(np.exp(0), np.exp(1.), cmap.N+3))
+        #         alphas = (scipy.special.erf(np.linspace(0,2.4, cmap.N+3)))#np.log(np.linspace(np.exp(0), np.exp(1.), cmap.N+3))
+        #         cmap._lut[:, -1] = alphas
+        #         print l
+        #         x, y = X[labels==l].T
+        #
+        #         heatmap, xedges, yedges = np.histogram2d(x, y, bins=300, range=range_)
+        #         #heatmap, xedges, yedges = np.histogram2d(x, y, bins=100)
+        #
+        #         im = ax.imshow(heatmap, extent=[xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap=cmap, aspect='auto', interpolation='nearest', label=str(l))
+        #         legend_handles.append(mpatches.Patch(color=colors[i], label=l))
+        #     ax.set_xlim(xmin, xmax)
+        #     ax.set_ylim(ymin, ymax)
+        #     plt.legend(legend_handles, [l.get_label() for l in legend_handles])
+        #     plt.draw()
+        #     plt.show()
+        #=======================================================================
+    return X, labels
+    
 
 def update_not_existing_kwargs(to_update, update_from):
     """
diff --git a/GPy/plotting/matplot_dep/plot_definitions.py b/GPy/plotting/matplot_dep/plot_definitions.py
index f7d49fcb..46b67283 100644
--- a/GPy/plotting/matplot_dep/plot_definitions.py
+++ b/GPy/plotting/matplot_dep/plot_definitions.py
@@ -79,7 +79,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
     def scatter(self, ax, X, Y, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, marker='o', **kwargs):
         if Z is not None:
             return ax.scatter(X, Y, c=color, zs=Z, label=label, marker=marker, **kwargs)
-        return ax.scatter(X, Y, c=color, label=label, **kwargs)
+        return ax.scatter(X, Y, c=color, label=label, marker=marker, **kwargs)
     
     def plot(self, ax, X, Y, color=None, label=None, **kwargs):
         return ax.plot(X, Y, color=color, label=label, **kwargs)
diff --git a/GPy/testing/plotting_tests.py b/GPy/testing/plotting_tests.py
index 3f8da666..4d2cb39d 100644
--- a/GPy/testing/plotting_tests.py
+++ b/GPy/testing/plotting_tests.py
@@ -31,6 +31,7 @@ import numpy as np
 import GPy, os
 from nose import SkipTest
 from matplotlib.testing.compare import compare_images
+from matplotlib.testing.noseclasses import ImageComparisonFailure
 
 try:
     from matplotlib import cbook, pyplot as plt
@@ -41,14 +42,15 @@ except:
 
 extensions = ['png']
 
-def _image_directories(func):
+def _image_directories():
     """
     Compute the baseline and result image directories for testing *func*.
     Create the result directory if it doesn't exist.
     """
-    module_name = func.__module__
-    mods = module_name.split('.')
-    basedir = os.path.join(*mods)
+    basedir = os.path.splitext(os.path.relpath(os.path.abspath(__file__)))[0]
+    #module_name = __init__.__module__
+    #mods = module_name.split('.')
+    #basedir = os.path.join(*mods)
     result_dir = os.path.join(basedir, 'testresult')
     baseline_dir = os.path.join(basedir, 'baseline')
     if not os.path.exists(result_dir):
@@ -56,45 +58,44 @@ def _image_directories(func):
     return baseline_dir, result_dir
 
 
-def sequenceEqual(a, b):
+def _sequenceEqual(a, b):
     assert len(a) == len(b), "Sequences not same length"
     for i, [x, y], in enumerate(zip(a, b)):
         assert x == y, "element not matching {}".format(i) 
 
-def notFound(path):
+def _notFound(path):
     raise IOError('File {} not in baseline')
 
-class test_image_comparison(object):
-    def __init__(self, baseline_images=[], extensions=['pdf','svg','ong']):
-        self.baseline_images = baseline_images
-        self.extensions = extensions
-        self.f = None
-    
-    def __call__(self, func):
-        self.baseline_dir, self.result_dir = _image_directories(func)
-        def test_wrap():
-            func()
-            for num, base in zip(plt.get_fignums(), self.baseline_images):
-                for ext in self.extensions:
-                    fig = plt.figure(num)
-                    fig.axes[0].set_axis_off()
-                    fig.set_frameon(False)
-                    fig.savefig(os.path.join(self.result_dir, "{}.{}".format(base, ext)), frameon=False)
-                    actual = os.path.join(self.result_dir, "{}.{}".format(base, ext))
-                    expected = os.path.join(self.baseline_dir, "{}.{}".format(base, ext))
-                    yield compare_images, actual, expected, 1e-3
-            plt.close('all')
-                    #with open(os.path.join(self.result_dir, "{}.{}".format(base, ext)), 'r') as f:
-                    #    try:
-                    #        with open(os.path.join(self.baseline_dir, "{}.{}".format(base, ext)), 'r') as b:
-                    #    except:
-                    #        yield notFound, os.path.join(self.baseline_dir, "{}.{}".format(base, ext))
-            #plt.close(num)
+def _image_comparison(baseline_images, extensions=['pdf','svg','ong'], tol=1e-3):
+    baseline_dir, result_dir = _image_directories()
+    for num, base in zip(plt.get_fignums(), baseline_images):
+        for ext in extensions:
+            fig = plt.figure(num)
+            fig.axes[0].set_axis_off()
+            fig.set_frameon(False)
+            fig.canvas.draw()
+            fig.savefig(os.path.join(result_dir, "{}.{}".format(base, ext)))
+        for num, base in zip(plt.get_fignums(), baseline_images):
+            for ext in extensions:
+                #plt.close(num)
+                actual = os.path.join(result_dir, "{}.{}".format(base, ext))
+                expected = os.path.join(baseline_dir, "{}.{}".format(base, ext))
+                def do_test():
+                    err = compare_images(actual, expected, tol)
+                    try:
+                        if not os.path.exists(expected):
+                            raise ImageComparisonFailure(
+                                'image does not exist: %s' % expected)
+                        if err:
+                            raise ImageComparisonFailure(
+                                'images not close: %(actual)s vs. %(expected)s '
+                                '(RMS %(rms).3f)'%err)
+                    except ImageComparisonFailure:
+                        pass
+                yield do_test
+    plt.close('all')
 
-        return test_wrap
-        
-@test_image_comparison(baseline_images=['gp_{}'.format(sub) for sub in ["data", "mean", 'conf', 'density', 'error']], extensions=extensions)
-def Plot(self=None):
+def test_plot(self=None):
     np.random.seed(11111)
     X = np.random.uniform(0, 1, (40, 1))
     f = .2 * np.sin(1.3*X) + 1.3*np.cos(2*X)
@@ -106,24 +107,22 @@ def Plot(self=None):
     m.plot_confidence()
     m.plot_density()
     m.plot_errorbars_trainset()
+    m.plot_samples()
+    for do_test in _image_comparison(baseline_images=['gp_{}'.format(sub) for sub in ["data", "mean", 'conf', 'density', 'error', 'samples']], extensions=extensions):
+        yield (do_test, )
 
-@test_image_comparison(baseline_images=['sparse_gp_{}'.format(sub) for sub in ["data", "mean", 'conf', 'density', 'error', 'inducing']], extensions=extensions)
-def PlotSparse(self=None):
+def test_plot_sparse(self=None):
     np.random.seed(11111)
-    X = np.random.uniform(0, 1, (40, 1))
+    X = np.random.uniform(-1, 1, (40, 1))
     f = .2 * np.sin(1.3*X) + 1.3*np.cos(2*X)
     Y = f+np.random.normal(0, .1, f.shape)
     m = GPy.models.SparseGPRegression(X, Y)
     m.optimize()
-    m.plot_data()
-    m.plot_mean()
-    m.plot_confidence()
-    m.plot_density()
-    m.plot_errorbars_trainset()
     m.plot_inducing()
+    for do_test in _image_comparison(baseline_images=['sparse_gp_{}'.format(sub) for sub in ['inducing']], extensions=extensions):
+        yield (do_test, )
 
-@test_image_comparison(baseline_images=['gp_class_{}'.format(sub) for sub in ["", "raw", 'link', 'raw_link']], extensions=extensions)
-def PlotClassification(self=None):
+def test_plot_classification(self=None):
     np.random.seed(11111)
     X = np.random.uniform(0, 1, (40, 1))
     f = .2 * np.sin(1.3*X) + 1.3*np.cos(2*X)
@@ -134,9 +133,11 @@ def PlotClassification(self=None):
     m.plot(plot_raw=True)
     m.plot(plot_raw=False, apply_link=True)
     m.plot(plot_raw=True, apply_link=True)
+    for do_test in _image_comparison(baseline_images=['gp_class_{}'.format(sub) for sub in ["", "raw", 'link', 'raw_link']], extensions=extensions):
+        yield (do_test, )
 
-@test_image_comparison(baseline_images=['sparse_gp_class_{}'.format(sub) for sub in ["", "raw", 'link', 'raw_link']], extensions=extensions)
-def PlotSparseClassification(self=None):
+ 
+def test_plot_sparse_classification(self=None):
     np.random.seed(11111)
     X = np.random.uniform(0, 1, (40, 1))
     f = .2 * np.sin(1.3*X) + 1.3*np.cos(2*X)
@@ -147,3 +148,29 @@ def PlotSparseClassification(self=None):
     m.plot(plot_raw=True)
     m.plot(plot_raw=False, apply_link=True)
     m.plot(plot_raw=True, apply_link=True)
+    for do_test in _image_comparison(baseline_images=['sparse_gp_class_{}'.format(sub) for sub in ["", "raw", 'link', 'raw_link']], extensions=extensions):
+        yield (do_test, )
+
+def test_gplvm_plot(self=None):
+    from ..examples.dimensionality_reduction import _simulate_matern
+    from ..kern import RBF
+    from ..models import GPLVM
+    Q = 3
+    _, _, Ylist = _simulate_matern(5, 1, 1, 100, num_inducing=5, plot_sim=False)
+    Y = Ylist[0]
+    k = RBF(Q, ARD=True)  # + kern.white(Q, _np.exp(-2)) # + kern.bias(Q)
+    # k = kern.RBF(Q, ARD=True, lengthscale=10.)
+    m = GPLVM(Y, Q, init="PCA", kernel=k)
+    m.likelihood.variance = .1
+    m.optimize(messages=0)
+    labels = np.random.multinomial(1, np.random.dirichlet([.3333333, .3333333, .3333333]), size=(m.Y.shape[0])).nonzero()[1]
+    m.plot_prediction_fit(which_data_ycols=(0,1)) # ignore this test, as plotting is not consistent!!
+    plt.close('all')
+    m.plot_latent()
+    m.plot_magnification(labels=labels)
+    for do_test in _image_comparison(baseline_images=['gplvm_{}'.format(sub) for sub in ["latent", "magnification"]], extensions=extensions):
+        yield (do_test, )
+        
+if __name__ == '__main__':
+    import nose
+    nose.main()
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_-failed-diff.png
deleted file mode 100644
index e6f3f308..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_.png b/GPy/testing/plotting_tests/baseline/gp_class_.png
index 0d9b84a1..8a1ba4d2 100644
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_.png and b/GPy/testing/plotting_tests/baseline/gp_class_.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class__pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class__pdf-failed-diff.png
deleted file mode 100644
index 7254c426..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class__pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class__pdf.png b/GPy/testing/plotting_tests/baseline/gp_class__pdf.png
deleted file mode 100644
index bf7681cd..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class__pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_link-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_link-failed-diff.png
deleted file mode 100644
index 26a94894..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_link-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_link.png b/GPy/testing/plotting_tests/baseline/gp_class_link.png
index ecb14e1c..8a1ba4d2 100644
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_link.png and b/GPy/testing/plotting_tests/baseline/gp_class_link.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_link_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_link_pdf-failed-diff.png
deleted file mode 100644
index 14bf6173..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_link_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_link_pdf.png b/GPy/testing/plotting_tests/baseline/gp_class_link_pdf.png
deleted file mode 100644
index 9ab36f73..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_link_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_raw-failed-diff.png
deleted file mode 100644
index 504cb8b4..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw.png b/GPy/testing/plotting_tests/baseline/gp_class_raw.png
index e65873b4..b168bc91 100644
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw.png and b/GPy/testing/plotting_tests/baseline/gp_class_raw.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw_link-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_raw_link-failed-diff.png
deleted file mode 100644
index 6222d49f..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw_link-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw_link.png b/GPy/testing/plotting_tests/baseline/gp_class_raw_link.png
index c6deb9b1..85e3158b 100644
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw_link.png and b/GPy/testing/plotting_tests/baseline/gp_class_raw_link.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw_link_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_raw_link_pdf-failed-diff.png
deleted file mode 100644
index c89d8eda..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw_link_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw_link_pdf.png b/GPy/testing/plotting_tests/baseline/gp_class_raw_link_pdf.png
deleted file mode 100644
index b186fe9a..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw_link_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_class_raw_pdf-failed-diff.png
deleted file mode 100644
index 2057b7a9..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_class_raw_pdf.png b/GPy/testing/plotting_tests/baseline/gp_class_raw_pdf.png
deleted file mode 100644
index b4db6744..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_class_raw_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_conf_pdf.png b/GPy/testing/plotting_tests/baseline/gp_conf_pdf.png
deleted file mode 100644
index 9ab36f73..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_conf_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_data_pdf.png b/GPy/testing/plotting_tests/baseline/gp_data_pdf.png
deleted file mode 100644
index bf7681cd..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_data_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_density-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_density-failed-diff.png
deleted file mode 100644
index d9c7d872..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_density-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_density.png b/GPy/testing/plotting_tests/baseline/gp_density.png
index c6deb9b1..cd248420 100644
Binary files a/GPy/testing/plotting_tests/baseline/gp_density.png and b/GPy/testing/plotting_tests/baseline/gp_density.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_density_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/gp_density_pdf-failed-diff.png
deleted file mode 100644
index e935e772..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_density_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_density_pdf.png b/GPy/testing/plotting_tests/baseline/gp_density_pdf.png
deleted file mode 100644
index b186fe9a..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_density_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_error_pdf.png b/GPy/testing/plotting_tests/baseline/gp_error_pdf.png
deleted file mode 100644
index ec839872..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_error_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_mean_pdf.png b/GPy/testing/plotting_tests/baseline/gp_mean_pdf.png
deleted file mode 100644
index b4db6744..00000000
Binary files a/GPy/testing/plotting_tests/baseline/gp_mean_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/gp_samples.png b/GPy/testing/plotting_tests/baseline/gp_samples.png
new file mode 100644
index 00000000..22f92d05
Binary files /dev/null and b/GPy/testing/plotting_tests/baseline/gp_samples.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gplvm_latent.png b/GPy/testing/plotting_tests/baseline/gplvm_latent.png
new file mode 100644
index 00000000..f8167b13
Binary files /dev/null and b/GPy/testing/plotting_tests/baseline/gplvm_latent.png differ
diff --git a/GPy/testing/plotting_tests/baseline/gplvm_magnification.png b/GPy/testing/plotting_tests/baseline/gplvm_magnification.png
new file mode 100644
index 00000000..584578a0
Binary files /dev/null and b/GPy/testing/plotting_tests/baseline/gplvm_magnification.png differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_-failed-diff.png
deleted file mode 100644
index b43d0363..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_.png
index 0d9b84a1..aed79690 100644
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_.png and b/GPy/testing/plotting_tests/baseline/sparse_gp_class_.png differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class__pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class__pdf-failed-diff.png
deleted file mode 100644
index 421c9a65..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class__pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class__pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class__pdf.png
deleted file mode 100644
index bf7681cd..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class__pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_link-failed-diff.png
deleted file mode 100644
index 0e6e5adc..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_link.png
index ecb14e1c..aed79690 100644
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link.png and b/GPy/testing/plotting_tests/baseline/sparse_gp_class_link.png differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_link_pdf-failed-diff.png
deleted file mode 100644
index e9c5eea2..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_link_pdf.png
deleted file mode 100644
index 9ab36f73..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_link_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw-failed-diff.png
deleted file mode 100644
index 0733541b..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw.png
index e65873b4..c94d34e1 100644
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw.png and b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw.png differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link-failed-diff.png
deleted file mode 100644
index 05f6e64f..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link.png
index c6deb9b1..94c017e5 100644
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link.png and b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link.png differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link_pdf-failed-diff.png
deleted file mode 100644
index 12ff09e5..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link_pdf.png
deleted file mode 100644
index b186fe9a..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_pdf-failed-diff.png
deleted file mode 100644
index 8c9061b9..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_pdf.png
deleted file mode 100644
index b4db6744..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_conf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_conf-failed-diff.png
deleted file mode 100644
index b89e020e..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_conf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_conf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_conf.png
deleted file mode 100644
index ecb14e1c..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_conf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_conf_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_conf_pdf.png
deleted file mode 100644
index 9ab36f73..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_conf_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_data.png b/GPy/testing/plotting_tests/baseline/sparse_gp_data.png
deleted file mode 100644
index 0d9b84a1..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_data.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_data_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_data_pdf.png
deleted file mode 100644
index bf7681cd..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_data_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_density-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_density-failed-diff.png
deleted file mode 100644
index d9c7d872..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_density-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_density.png b/GPy/testing/plotting_tests/baseline/sparse_gp_density.png
deleted file mode 100644
index c6deb9b1..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_density.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_density_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_density_pdf-failed-diff.png
deleted file mode 100644
index 8666ff25..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_density_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_density_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_density_pdf.png
deleted file mode 100644
index b186fe9a..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_density_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_error.png b/GPy/testing/plotting_tests/baseline/sparse_gp_error.png
deleted file mode 100644
index 3dd5cddb..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_error.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_error_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_error_pdf.png
deleted file mode 100644
index ec839872..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_error_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_inducing-failed-diff.png
deleted file mode 100644
index 776f2fe0..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing.png b/GPy/testing/plotting_tests/baseline/sparse_gp_inducing.png
index 8a1ba4d2..4c3becae 100644
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing.png and b/GPy/testing/plotting_tests/baseline/sparse_gp_inducing.png differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing_pdf-failed-diff.png b/GPy/testing/plotting_tests/baseline/sparse_gp_inducing_pdf-failed-diff.png
deleted file mode 100644
index f548c535..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing_pdf-failed-diff.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_inducing_pdf.png
deleted file mode 100644
index 64589594..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_inducing_pdf.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_mean.png b/GPy/testing/plotting_tests/baseline/sparse_gp_mean.png
deleted file mode 100644
index e65873b4..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_mean.png and /dev/null differ
diff --git a/GPy/testing/plotting_tests/baseline/sparse_gp_mean_pdf.png b/GPy/testing/plotting_tests/baseline/sparse_gp_mean_pdf.png
deleted file mode 100644
index b4db6744..00000000
Binary files a/GPy/testing/plotting_tests/baseline/sparse_gp_mean_pdf.png and /dev/null differ