diff --git a/GPy/installation.cfg b/GPy/installation.cfg
index 901a7ef5..2056c0b9 100644
--- a/GPy/installation.cfg
+++ b/GPy/installation.cfg
@@ -1,2 +1,5 @@
 # This is the local installation configuration file for GPy
 
+[plotting]
+#library = plotly
+library = matplotlib
\ No newline at end of file
diff --git a/GPy/plotting/abstract_plotting_library.py b/GPy/plotting/abstract_plotting_library.py
index 4743f1bc..f46843cc 100644
--- a/GPy/plotting/abstract_plotting_library.py
+++ b/GPy/plotting/abstract_plotting_library.py
@@ -81,16 +81,6 @@ class AbstractPlottingLibrary(object):
 
         E.g. in matplotlib this means it deletes references to ax, as
         plotting is done on the axis itself and is not a kwarg. 
-        """
-        raise NotImplementedError("Implement all plot functions in AbstractPlottingLibrary in order to use your own plotting library")
-
-    def show_canvas(self, canvas, plots, legend=True, title=None, **kwargs):
-        """
-        Show the canvas given. 
-        plots is a dictionary with the plots
-        as the items.
-        
-        the kwargs are plotting library specific kwargs!
 
         :param xlabel: the label to put on the xaxis
         :param ylabel: the label to put on the yaxis
@@ -100,11 +90,32 @@ class AbstractPlottingLibrary(object):
         :param (float, float) xlim: the limits for the xaxis
         :param (float, float) ylim: the limits for the yaxis
         :param (float, float) zlim: the limits for the zaxis (if plotting in 3d)
-        
-        E.g. in matplotlib this does not have to do anything, we make the tight plot, though.
         """
         raise NotImplementedError("Implement all plot functions in AbstractPlottingLibrary in order to use your own plotting library")
 
+    def add_to_canvas(self, canvas, plots, legend=True, title=None, **kwargs):
+        """
+        Add plots is a dictionary with the plots as the 
+        items or a list of plots as items to canvas.
+        
+        The kwargs are plotting library specific kwargs!
+        
+        E.g. in matplotlib this does not have to do anything to add stuff, but
+        we set the legend and title.
+
+        !This function returns the updated canvas!
+
+        :param title: the title of the plot
+        :param legend: whether to plot a legend or not
+        """
+        raise NotImplementedError("Implement all plot functions in AbstractPlottingLibrary in order to use your own plotting library")
+
+    def show_canvas(self, canvas, **kwargs):
+        """
+        Draw/Plot the canvas given.
+        """
+        raise NotImplementedError
+
     def plot(self, cavas, X, Y, Z=None, color=None, label=None, **kwargs):
         """
         Make a line plot from for Y on X (Y = f(X)) on the canvas.
@@ -177,8 +188,8 @@ class AbstractPlottingLibrary(object):
     def yerrorbar(self, canvas, X, Y, error, color=None, label=None, **kwargs):
         """
         Make errorbars along the yaxis on the canvas given.
-        if error is two dimensional, the lower error is error[:,0] and
-        the upper error is error[:,1]
+        if error is two dimensional, the lower error is error[0, :] and
+        the upper error is error[1, :]
                 
         the kwargs are plotting library specific kwargs!
         """
diff --git a/GPy/plotting/gpy_plot/data_plots.py b/GPy/plotting/gpy_plot/data_plots.py
index f5743aaa..e2fb1fe7 100644
--- a/GPy/plotting/gpy_plot/data_plots.py
+++ b/GPy/plotting/gpy_plot/data_plots.py
@@ -56,7 +56,7 @@ def plot_data(self, which_data_rows='all',
     """
     canvas, plot_kwargs = pl.new_canvas(projection=projection, **plot_kwargs)
     plots = _plot_data(self, canvas, which_data_rows, which_data_ycols, visible_dims, projection, label, **plot_kwargs)
-    return pl.show_canvas(canvas, plots)
+    return pl.add_to_canvas(canvas, plots)
 
 def _plot_data(self, canvas, which_data_rows='all',
         which_data_ycols='all', visible_dims=None,
@@ -81,12 +81,12 @@ def _plot_data(self, canvas, which_data_rows='all',
             for d in ycols:
                 update_not_existing_kwargs(plot_kwargs, pl.defaults.data_2d)  # @UndefinedVariable
                 plots['dataplot'].append(pl.scatter(canvas, X[rows, free_dims[0]], X[rows, free_dims[1]], 
-                                               color=Y[rows, d], vmin=Y.min(), vmax=Y.max(), label=label, **plot_kwargs))
+                                               color=Y[rows, d], label=label, **plot_kwargs))
         else:
             for d in ycols:
                 update_not_existing_kwargs(plot_kwargs, pl.defaults.data_2d)  # @UndefinedVariable
                 plots['dataplot'].append(pl.scatter(canvas, X[rows, free_dims[0]], X[rows, free_dims[1]], 
-                                                    Z=Y[rows, d], vmin=Y.min(), color=Y[rows, d], vmax=Y.max(), label=label, **plot_kwargs))
+                                                    Z=Y[rows, d], color=Y[rows, d], label=label, **plot_kwargs))
     elif len(free_dims) == 0:
         pass #Nothing to plot!
     else:
@@ -117,9 +117,9 @@ def plot_data_error(self, which_data_rows='all',
     
     :returns list: of plots created.
     """
-    canvas, error_kwargs = pl.new_canvas(projection=='3d', **error_kwargs)
+    canvas, error_kwargs = pl.new_canvas(projection=projection, **error_kwargs)
     plots = _plot_data_error(self, canvas, which_data_rows, which_data_ycols, visible_dims, projection, label, **error_kwargs)
-    return pl.show_canvas(canvas, plots)
+    return pl.add_to_canvas(canvas, plots)
 
 def _plot_data_error(self, canvas, which_data_rows='all',
         which_data_ycols='all', visible_dims=None,
@@ -167,7 +167,7 @@ def plot_inducing(self, visible_dims=None, projection='2d', label=None, **plot_k
     """
     canvas, kwargs = pl.new_canvas(projection=projection, **plot_kwargs)
     plots = _plot_inducing(self, canvas, visible_dims, projection, label, **kwargs)
-    return pl.show_canvas(canvas, plots)
+    return pl.add_to_canvas(canvas, plots)
 
 def _plot_inducing(self, canvas, visible_dims, projection, label, **plot_kwargs):
     if visible_dims is None:
@@ -220,7 +220,7 @@ def plot_errorbars_trainset(self, which_data_rows='all',
     canvas, kwargs = pl.new_canvas(projection=projection, **plot_kwargs)
     plots = _plot_errorbars_trainset(self, canvas, which_data_rows, which_data_ycols, 
                                      fixed_inputs, plot_raw, apply_link, label, projection, predict_kw, **kwargs)
-    return pl.show_canvas(canvas, plots)
+    return pl.add_to_canvas(canvas, plots)
 
 def _plot_errorbars_trainset(self, canvas, 
         which_data_rows='all', which_data_ycols='all', 
@@ -245,25 +245,31 @@ def _plot_errorbars_trainset(self, canvas,
     
     if len(free_dims)<=2:
         update_not_existing_kwargs(plot_kwargs, pl.defaults.yerrorbar)
-        if len(free_dims)==1:
-            if predict_kw is None:
+        if predict_kw is None:
                 predict_kw = {}
-            if 'Y_metadata' not in predict_kw:
-                predict_kw['Y_metadata'] = self.Y_metadata or {}
-            _, percs, _ = helper_predict_with_model(self, Xgrid, plot_raw, 
-                                              apply_link, (2.5, 97.5), 
-                                              ycols, predict_kw)
+        if 'Y_metadata' not in predict_kw:
+            predict_kw['Y_metadata'] = self.Y_metadata or {}
+        mu, percs, _ = helper_predict_with_model(self, Xgrid, plot_raw, 
+                                          apply_link, (2.5, 97.5), 
+                                          ycols, predict_kw)
+        if len(free_dims)==1:
             for d in ycols:
-                plots.append(pl.yerrorbar(canvas, X[rows,free_dims[0]], Y[rows,d], 
-                                          np.vstack([Y[rows,d]-percs[0][rows,d], percs[1][rows,d]-Y[rows,d]]),
+                plots.append(pl.yerrorbar(canvas, X[rows,free_dims[0]], mu[rows,d], 
+                                          np.vstack([mu[rows, d] - percs[0][rows, d], percs[1][rows, d] - mu[rows,d]]),
                                           label=label, 
                                           **plot_kwargs))
         elif len(free_dims) == 2:
-            plots.append(pl.yerrorbar(canvas, X[rows,free_dims[0]], X[rows,free_dims[1]], 
-                          np.vstack([Y[rows,d]-percs[0][rows,d], percs[1][rows,d]-Y[rows,d]]),
-                          Y[rows,d],
-                          label=label, 
-                          **plot_kwargs))
+            for d in ycols:
+                plots.append(pl.yerrorbar(canvas, X[rows,free_dims[0]], X[rows,free_dims[1]], 
+                              np.vstack([mu[rows, d] - percs[0][rows, d], percs[1][rows, d] - mu[rows,d]]),
+                              color=Y[rows,d],
+                              label=label, 
+                              **plot_kwargs))
+                plots.append(pl.xerrorbar(canvas, X[rows,free_dims[0]], X[rows,free_dims[1]], 
+                              np.vstack([mu[rows, d] - percs[0][rows, d], percs[1][rows, d] - mu[rows,d]]),
+                              color=Y[rows,d],
+                              label=label, 
+                              **plot_kwargs))
             pass #Nothing to plot!
     else:
         raise NotImplementedError("Cannot plot in more then one dimension.")
diff --git a/GPy/plotting/gpy_plot/gp_plots.py b/GPy/plotting/gpy_plot/gp_plots.py
index 222a4318..d10e75ff 100644
--- a/GPy/plotting/gpy_plot/gp_plots.py
+++ b/GPy/plotting/gpy_plot/gp_plots.py
@@ -63,26 +63,23 @@ def plot_mean(self, plot_limits=None, fixed_inputs=None,
     :param dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=<specific kernel>) in here
     """
     canvas, kwargs = pl.new_canvas(projection=projection, **kwargs)
-    plots = _plot_mean(self, canvas, plot_limits, fixed_inputs, 
-                       resolution, plot_raw, 
-                       apply_link, visible_dims, which_data_ycols, 
-                       levels, projection, label, predict_kw, **kwargs)
-    return pl.show_canvas(canvas, plots)
-
-def _plot_mean(self, canvas, plot_limits=None, fixed_inputs=None,
-              resolution=None, plot_raw=False,
-              apply_link=False, visible_dims=None,
-              which_data_ycols='all',
-              levels=20, projection='2d', label=None,
-              predict_kw=None,
-              **kwargs):
-    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
-
-    if len(free_dims)<=2:
-        mu, _, _ = helper_predict_with_model(self, Xgrid, plot_raw, 
+    helper_data = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+    helper_prediction = helper_predict_with_model(self, helper_data[5], plot_raw, 
                                           apply_link, None, 
                                           get_which_data_ycols(self, which_data_ycols), 
                                           predict_kw)
+    plots = _plot_mean(self, canvas, helper_data, helper_prediction, 
+                       levels, projection, label, **kwargs)
+    pl.add_to_canvas(canvas, plots)
+    return pl.show_canvas(canvas)
+
+def _plot_mean(self, canvas, helper_data, helper_prediction, 
+              levels=20, projection='2d', label=None,
+              **kwargs):
+
+    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_data
+    if len(free_dims)<=2:
+        mu, _, _ = helper_prediction
         if len(free_dims)==1:
             # 1D plotting:
             update_not_existing_kwargs(kwargs, pl.defaults.meanplot_1d)  # @UndefinedVariable
@@ -108,7 +105,7 @@ def _plot_mean(self, canvas, plot_limits=None, fixed_inputs=None,
 def plot_confidence(self, lower=2.5, upper=97.5, plot_limits=None, fixed_inputs=None,
               resolution=None, plot_raw=False,
               apply_link=False, visible_dims=None,
-              which_data_ycols='all',
+              which_data_ycols='all', label=None,
               predict_kw=None, 
               **kwargs):
     """
@@ -132,33 +129,23 @@ def plot_confidence(self, lower=2.5, upper=97.5, plot_limits=None, fixed_inputs=
     :param dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=<specific kernel>) in here
     """
     canvas, kwargs = pl.new_canvas(**kwargs)
-    plots = _plot_confidence(self, canvas, lower, upper, plot_limits, 
-                             fixed_inputs, resolution, plot_raw, 
-                             apply_link, visible_dims, which_data_ycols, 
-                             predict_kw, **kwargs)
-    return pl.show_canvas(canvas, plots)
-
-def _plot_confidence(self, canvas, lower, upper, plot_limits=None, fixed_inputs=None,
-              resolution=None, plot_raw=False,
-              apply_link=False, visible_dims=None,
-              which_data_ycols=None, 
-              predict_kw=None, 
-              **kwargs):
-    _, _, _, _, free_dims, Xgrid, _, _, _, _, _ = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
-
     ycols = get_which_data_ycols(self, which_data_ycols)
-    
-    update_not_existing_kwargs(kwargs, pl.defaults.confidence_interval)  # @UndefinedVariable
-    
-    if len(free_dims)<=1:
-        if len(free_dims)==1:
-            _, percs, _ = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, 
+    helper_data = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+    helper_prediction = helper_predict_with_model(self, helper_data[5], plot_raw, apply_link, 
                                                  (lower, upper), 
                                                  ycols, predict_kw)
-    
+    plots = _plot_confidence(self, canvas, helper_data, helper_prediction, label, **kwargs)
+    return pl.add_to_canvas(canvas, plots)
+
+def _plot_confidence(self, canvas, helper_data, helper_prediction, label, **kwargs):
+    _, _, _, _, free_dims, Xgrid, _, _, _, _, _ = helper_data
+    update_not_existing_kwargs(kwargs, pl.defaults.confidence_interval)  # @UndefinedVariable
+    if len(free_dims)<=1:
+        if len(free_dims)==1:
+            percs = helper_prediction[1]
             fills = []
-            for d in ycols:
-                fills.append(pl.fill_between(canvas, Xgrid[:,free_dims[0]], percs[0][:,d], percs[1][:,d], **kwargs))
+            for d in range(helper_prediction[0].shape[1]):
+                fills.append(pl.fill_between(canvas, Xgrid[:,free_dims[0]], percs[0][:,d], percs[1][:,d], label=label, **kwargs))
             return dict(gpconfidence=fills)
         else:
             pass #Nothing to plot!
@@ -170,7 +157,8 @@ def plot_samples(self, plot_limits=None, fixed_inputs=None,
               resolution=None, plot_raw=True,
               apply_link=False, visible_dims=None,
               which_data_ycols='all',
-              samples=3, projection='2d', predict_kw=None,
+              samples=3, projection='2d', label=None,
+              predict_kw=None,
               **kwargs):
     """
     Plot the mean of the GP.
@@ -191,36 +179,31 @@ def plot_samples(self, plot_limits=None, fixed_inputs=None,
     :param int levels: for 2D plotting, the number of contour levels to use is 
     """
     canvas, kwargs = pl.new_canvas(projection=projection, **kwargs)
-    plots = _plot_samples(self, canvas, plot_limits, fixed_inputs, 
-                       resolution, plot_raw, 
-                       apply_link, visible_dims, which_data_ycols, samples, projection, 
-                       predict_kw, **kwargs)
-    return pl.show_canvas(canvas, plots)
+    ycols = get_which_data_ycols(self, which_data_ycols)
+    helper_data = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+    helper_prediction = helper_predict_with_model(self, helper_data[5], plot_raw, apply_link, 
+                                                 None, 
+                                                 ycols, predict_kw, samples)
+    plots = _plot_samples(self, canvas, helper_data, helper_prediction,
+                          projection, label, **kwargs)
+    return pl.add_to_canvas(canvas, plots)
 
-def _plot_samples(self, canvas, plot_limits=None, fixed_inputs=None,
-              resolution=None, plot_raw=True,
-              apply_link=False, visible_dims=None,
-              which_data_ycols=None,
-              samples=3, projection='2d',
-              label=None,
-              predict_kw=None, **kwargs):
-    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+def _plot_samples(self, canvas, helper_data, helper_prediction, projection,
+              label, **kwargs):
+    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_data
+    samples = helper_prediction[2]
 
     if len(free_dims)<=2:
         if len(free_dims)==1:
             # 1D plotting:
-            _, _, samples = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, 
-                                     None, get_which_data_ycols(self, which_data_ycols), predict_kw, samples)
             update_not_existing_kwargs(kwargs, pl.defaults.samples_1d)  # @UndefinedVariable
-            return dict(gpmean=[pl.plot(canvas, Xgrid[:, free_dims], samples, **kwargs)])
+            return dict(gpmean=[pl.plot(canvas, Xgrid[:, free_dims], samples, label=label, **kwargs)])
         elif len(free_dims)==2 and projection=='3d':
-            _, _, samples = helper_predict_with_model(self, Xgrid, plot_raw, apply_link, 
-                                     None, get_which_data_ycols(self, which_data_ycols), predict_kw, samples)
             update_not_existing_kwargs(kwargs, pl.defaults.samples_3d)  # @UndefinedVariable
             for s in range(samples.shape[-1]):
                 return dict(gpmean=[pl.surface(canvas, x, 
-                                            y, samples[:, s].reshape(resolution, resolution), 
-                                            **kwargs)])            
+                                    y, samples[:, s].reshape(resolution, resolution), 
+                                    **kwargs)])            
         else:
             pass # Nothing to plot!
     else:
@@ -231,7 +214,7 @@ def plot_density(self, plot_limits=None, fixed_inputs=None,
               resolution=None, plot_raw=False,
               apply_link=False, visible_dims=None, 
               which_data_ycols='all',
-              levels=35,
+              levels=35, label=None, 
               predict_kw=None, 
               **kwargs):
     """
@@ -254,35 +237,26 @@ def plot_density(self, plot_limits=None, fixed_inputs=None,
     :param dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=<specific kernel>) in here
     """
     canvas, kwargs = pl.new_canvas(**kwargs)
-    plots = _plot_density(self, canvas, plot_limits, 
-                             fixed_inputs, resolution, plot_raw,  
-                             apply_link, visible_dims, which_data_ycols, 
-                             levels,
-                             predict_kw, **kwargs)
-    return pl.show_canvas(canvas, plots)
+    helper_data = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+    helper_prediction = helper_predict_with_model(self, helper_data[5], plot_raw, 
+                                          apply_link, np.linspace(2.5, 97.5, levels*2), 
+                                          get_which_data_ycols(self, which_data_ycols), 
+                                          predict_kw)
+    plots = _plot_density(self, canvas, helper_data, helper_prediction, label, **kwargs)
+    return pl.add_to_canvas(canvas, plots)
 
-def _plot_density(self, canvas, plot_limits=None, fixed_inputs=None,
-              resolution=None, plot_raw=False,
-              apply_link=False, visible_dims=None,
-              which_data_ycols=None,
-              levels=35, 
-              predict_kw=None, **kwargs):
-    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
-
-    ycols = get_which_data_ycols(self, which_data_ycols)
+def _plot_density(self, canvas, helper_data, helper_prediction, label, **kwargs):
+    _, _, _, _, free_dims, Xgrid, _, _, _, _, _ = helper_data
+    mu, percs, _ = helper_prediction
 
     update_not_existing_kwargs(kwargs, pl.defaults.density)  # @UndefinedVariable
 
     if len(free_dims)<=1:
         if len(free_dims)==1:
-            _, percs, _ = helper_predict_with_model(self, Xgrid, plot_raw, 
-                                          apply_link, np.linspace(2.5, 97.5, levels*2), 
-                                          get_which_data_ycols(self, which_data_ycols), 
-                                          predict_kw)
             # 1D plotting:
             fills = []
-            for d in ycols:
-                fills.append(pl.fill_gradient(canvas, Xgrid[:, free_dims[0]], [p[:,d] for p in percs], **kwargs))
+            for d in range(mu.shape[1]):
+                fills.append(pl.fill_gradient(canvas, Xgrid[:, free_dims[0]], [p[:,d] for p in percs], label=label, **kwargs))
             return dict(gpdensity=fills)
         else:
             pass # Nothing to plot!
@@ -327,11 +301,28 @@ def plot(self, plot_limits=None, fixed_inputs=None,
     :param dict predict_kw: the keyword arguments for the prediction. If you want to plot a specific kernel give dict(kern=<specific kernel>) in here
     """
     canvas, _ = pl.new_canvas(projection=projection, **kwargs)
-    plots = _plot(self, canvas, plot_limits, fixed_inputs, resolution, plot_raw, 
-                  apply_link, which_data_ycols, which_data_rows, visible_dims, 
-                  levels, samples, samples_likelihood, lower, upper, plot_data, 
-                  plot_inducing, plot_density, projection, predict_kw)
-    return pl.show_canvas(canvas, plots)
+    helper_data = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+    helper_prediction = helper_predict_with_model(self, helper_data[5], plot_raw, 
+                                          apply_link, np.linspace(2.5, 97.5, levels*2) if plot_density else (lower,upper), 
+                                          get_which_data_ycols(self, which_data_ycols), 
+                                          predict_kw, samples)
+    if plot_raw and not apply_link:
+        # It does not make sense to plot the data (which lives not in the latent function space) into latent function space. 
+        plot_data = False
+    plots = {}
+    if plot_data:
+        plots.update(_plot_data(self, canvas, which_data_rows, which_data_ycols, visible_dims, projection))
+        plots.update(_plot_data_error(self, canvas, which_data_rows, which_data_ycols, visible_dims, projection))
+    plots.update(_plot(self, canvas, plots, helper_data, helper_prediction, levels, plot_inducing, plot_density, projection))
+    if plot_raw and (samples_likelihood > 0):
+        helper_prediction = helper_predict_with_model(self, helper_data[5], False, 
+                                      apply_link, None, 
+                                      get_which_data_ycols(self, which_data_ycols), 
+                                      predict_kw, samples_likelihood)
+        plots.update(_plot_samples(canvas, helper_data, helper_prediction, projection))
+    if hasattr(self, 'Z') and plot_inducing:
+        plots.update(_plot_inducing(self, canvas, visible_dims, projection, None))
+    return pl.add_to_canvas(canvas, plots)
 
 
 def plot_f(self, plot_limits=None, fixed_inputs=None,
@@ -375,50 +366,35 @@ def plot_f(self, plot_limits=None, fixed_inputs=None,
     :param dict error_kwargs: kwargs for the error plot for the plotting library you are using
     :param kwargs plot_kwargs: kwargs for the data plot for the plotting library you are using
     """
-    canvas, _ = pl.new_canvas(projection=='3d', **kwargs)
-    plots = _plot(self, canvas, plot_limits, fixed_inputs, resolution, 
-                  True, apply_link, which_data_ycols, which_data_rows, 
-                  visible_dims, levels, samples, 0, lower, upper, 
-                  plot_data, plot_inducing, plot_density, projection,
-                  predict_kw)
-    return pl.show_canvas(canvas, plots)
+    canvas, _ = pl.new_canvas(projection=projection, **kwargs)
+    helper_data = helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resolution)
+    helper_prediction = helper_predict_with_model(self, helper_data[5], True, 
+                                          apply_link, np.linspace(2.5, 97.5, levels*2) if plot_density else (lower,upper), 
+                                          get_which_data_ycols(self, which_data_ycols), 
+                                          predict_kw, samples)
+    if not apply_link:
+        # It does not make sense to plot the data (which lives not in the latent function space) into latent function space. 
+        plot_data = False
+    plots = {}
+    if plot_data:
+        plots.update(_plot_data(self, canvas, which_data_rows, which_data_ycols, visible_dims, projection))
+        plots.update(_plot_data_error(self, canvas, which_data_rows, which_data_ycols, visible_dims, projection))
+    plots.update(_plot(self, canvas, plots, helper_data, helper_prediction, levels, plot_inducing, plot_density, projection))
+    if hasattr(self, 'Z') and plot_inducing:
+        plots.update(_plot_inducing(self, canvas, visible_dims, projection, None))
+    return pl.add_to_canvas(canvas, plots)
 
 
 
-def _plot(self, canvas, plot_limits=None, fixed_inputs=None,
-              resolution=None, 
-              plot_raw=False, apply_link=False, 
-              which_data_ycols='all', which_data_rows='all',
-              visible_dims=None, 
-              levels=20, samples=0, samples_likelihood=0, lower=2.5, upper=97.5, 
-              plot_data=True, plot_inducing=True, plot_density=False, projection='2d',
-              predict_kw=None):  
-        plots = {}
-        if plot_raw and not apply_link:
-            # It does not make sense to plot the data (which lives not in the latent function space) into latent function space. 
-            plot_data = False
-            
-        if plot_data:
-            plots.update(_plot_data(self, canvas, which_data_rows, which_data_ycols, visible_dims,
-                                    projection, label=None))
-            plots.update(_plot_data_error(self, canvas, which_data_rows, which_data_ycols, visible_dims, 
-                                          projection, label=None))
-            
-        plots.update(_plot_mean(self, canvas, plot_limits, fixed_inputs, resolution, plot_raw, apply_link, visible_dims, which_data_ycols, levels, projection, label=None,
-              predict_kw=None))
+def _plot(self, canvas, plots, helper_data, helper_prediction, levels, plot_inducing=True, plot_density=False, projection='2d'):    
+        plots.update(_plot_mean(self, canvas, helper_data, helper_prediction, levels, projection, None))
         
         if projection=='2d':
             if not plot_density:
-                plots.update(_plot_confidence(self, canvas, lower, upper, plot_limits, fixed_inputs, resolution, plot_raw, apply_link, visible_dims, which_data_ycols, predict_kw))
+                plots.update(_plot_confidence(self, canvas, helper_data, helper_prediction, None))
             else:
-                plots.update(_plot_density(self, canvas, plot_limits, fixed_inputs, resolution, plot_raw, apply_link, visible_dims, which_data_ycols, levels, predict_kw))
-        
-        if samples > 0:
-            plots.update(_plot_samples(self, canvas, plot_limits, fixed_inputs, resolution, True, apply_link, visible_dims, which_data_ycols, samples, predict_kw))
-        if samples_likelihood > 0:
-            plots.update(_plot_samples(self, canvas, plot_limits, fixed_inputs, resolution, False, apply_link, visible_dims, which_data_ycols, samples, predict_kw))
-
-        if hasattr(self, 'Z') and plot_inducing:
-            plots.update(_plot_inducing(self, canvas, visible_dims, projection, None))
+                plots.update(_plot_density(self, canvas, helper_data, helper_prediction, None))
         
+        if helper_prediction[2] is not None:
+            plots.update(_plot_samples(self, canvas, helper_data, helper_prediction, projection, None))        
         return plots
\ No newline at end of file
diff --git a/GPy/plotting/gpy_plot/inference_plots.py b/GPy/plotting/gpy_plot/inference_plots.py
index bbb75c3b..8a0d5b04 100644
--- a/GPy/plotting/gpy_plot/inference_plots.py
+++ b/GPy/plotting/gpy_plot/inference_plots.py
@@ -13,7 +13,7 @@ def plot_optimizer(optimizer, **kwargs):
     else:
         canvas, kwargs = pl.new_canvas(**kwargs)
         plots = dict(trace=pl.plot(range(len(optimizer.trace)), optimizer.trace))
-        return pl.show_canvas(canvas, plots, xlabel='Iteration', ylabel='f(x)')
+        return pl.add_to_canvas(canvas, plots, xlabel='Iteration', ylabel='f(x)')
 
 def plot_sgd_traces(optimizer):
     figure = pl.figure(2,1)
@@ -21,8 +21,8 @@ def plot_sgd_traces(optimizer):
     plots = dict(lines=[])
     for k in optimizer.param_traces.keys():
         plots['lines'].append(pl.plot(canvas, range(len(optimizer.param_traces[k])), optimizer.param_traces[k], label=k))
-    pl.show_canvas(canvas, legend=True)
+    pl.add_to_canvas(canvas, legend=True)
     canvas, _ = pl.new_canvas(figure, 1, 2, title="Objective function")
     pl.plot(canvas, range(len(optimizer.fopt_trace)), optimizer.fopt_trace)
-    return pl.show_canvas(canvas, plots, legend=True)
+    return pl.add_to_canvas(canvas, plots, legend=True)
     
\ No newline at end of file
diff --git a/GPy/plotting/gpy_plot/kernel_plots.py b/GPy/plotting/gpy_plot/kernel_plots.py
index e9b20196..532faf14 100644
--- a/GPy/plotting/gpy_plot/kernel_plots.py
+++ b/GPy/plotting/gpy_plot/kernel_plots.py
@@ -30,6 +30,10 @@
 import numpy as np
 from . import pl
 from .. import Tango
+from .plot_util import get_x_y_var,\
+    update_not_existing_kwargs, \
+    helper_for_plot_data, scatter_label_generator, subsample_X,\
+    find_best_layout_for_subplots
 
 def plot_ARD(kernel, filtering=None, **kwargs):
     """
@@ -64,11 +68,95 @@ def plot_ARD(kernel, filtering=None, **kwargs):
         else:
             print("filtering out {}".format(kernel.parameters[i].name))
 
-    plt.show_canvas()
+    plt.add_to_canvas()
     ax.set_xlim(-.5, kernel.input_dim - .5)
     add_bar_labels(fig, ax, [bars[-1]], bottom=bottom-last_bottom)
 
     return dict(barplots=bars)
 
-def plot_covariance():
+def plot_covariance(kernel, x=None, label=None, plot_limits=None, visible_dims=None, resolution=None, projection=None, levels=20, **mpl_kwargs):
+    """
+    plot a kernel.
+    :param x: the value to use for the other kernel argument (kernels are a function of two variables!)
+    :param fignum: figure number of the plot
+    :param ax: matplotlib axis to plot on
+    :param title: the matplotlib title
+    :param plot_limits: the range over which to plot the kernel
+    :resolution: the resolution of the lines used in plotting
+    :mpl_kwargs avalid keyword arguments to pass through to matplotlib (e.g. lw=7)
+    """
+    canvas, error_kwargs = pl.new_canvas(projection=projection, **error_kwargs)
+    _, _, _, _, free_dims, Xgrid, x, y, _, _, resolution = helper_for_plot_data(kernel, plot_limits, visible_dims, None, resolution)
+
+    if len(free_dims)<=2:
+        if len(free_dims)==1:
+            if x is None: x = np.zeros((1, 1))
+            else:
+                x = np.asarray(x)
+                assert x.size == 1, "The size of the fixed variable x is not 1"
+                x = x.reshape((1, 1))
+            # 1D plotting:
+            update_not_existing_kwargs(kwargs, pl.defaults.meanplot_1d)  # @UndefinedVariable
+            plots = dict(covariance=[pl.plot(canvas, Xgrid[:, free_dims], mu, label=label, **kwargs)])
+        else:
+            if projection == '2d':
+                update_not_existing_kwargs(kwargs, pl.defaults.meanplot_2d)  # @UndefinedVariable
+                plots = dict(covariance=[pl.contour(canvas, x, y, 
+                                               mu.reshape(resolution, resolution).T, 
+                                               levels=levels, label=label, **kwargs)])
+            elif projection == '3d':
+                update_not_existing_kwargs(kwargs, pl.defaults.meanplot_3d)  # @UndefinedVariable
+                plots = dict(covariance=[pl.surface(canvas, x, y, 
+                                               mu.reshape(resolution, resolution).T, 
+                                               label=label, 
+                                               **kwargs)])
+                
+        return pl.add_to_canvas(canvas, plots)
+
+    if kernel.input_dim == 1:
+
+        if plot_limits == None:
+            xmin, xmax = (x - 5).flatten(), (x + 5).flatten()
+        elif len(plot_limits) == 2:
+            xmin, xmax = plot_limits
+        else:
+            raise ValueError("Bad limits for plotting")
+
+        Xnew = np.linspace(xmin, xmax, resolution or 201)[:, None]
+        Kx = kernel.K(Xnew, x)
+        ax.plot(Xnew, Kx, **mpl_kwargs)
+        ax.set_xlim(xmin, xmax)
+        ax.set_xlabel("x")
+        ax.set_ylabel("k(x,%0.1f)" % x)
+
+    elif kernel.input_dim == 2:
+        if x is None:
+            x = np.zeros((1, 2))
+        else:
+            x = np.asarray(x)
+            assert x.size == 2, "The size of the fixed variable x is not 2"
+            x = x.reshape((1, 2))
+
+        if plot_limits is None:
+            xmin, xmax = (x - 5).flatten(), (x + 5).flatten()
+        elif len(plot_limits) == 2:
+            xmin, xmax = plot_limits
+        else:
+            raise ValueError("Bad limits for plotting")
+
+
+        resolution = resolution or 51
+        xx, yy = np.mgrid[xmin[0]:xmax[0]:1j * resolution, xmin[1]:xmax[1]:1j * resolution]
+        Xnew = np.vstack((xx.flatten(), yy.flatten())).T
+        Kx = kernel.K(Xnew, x)
+        Kx = Kx.reshape(resolution, resolution).T
+        ax.contour(xx, yy, Kx, vmin=Kx.min(), vmax=Kx.max(), cmap=pb.cm.jet, **mpl_kwargs) # @UndefinedVariable
+        ax.set_xlim(xmin[0], xmax[0])
+        ax.set_ylim(xmin[1], xmax[1])
+        ax.set_xlabel("x1")
+        ax.set_ylabel("x2")
+        ax.set_title("k(x1,x2 ; %0.1f,%0.1f)" % (x[0, 0], x[0, 1]))
+    else:
+        raise NotImplementedError("Cannot plot a kernel with more than two input dimensions")
+    
     pass
\ No newline at end of file
diff --git a/GPy/plotting/gpy_plot/latent_plots.py b/GPy/plotting/gpy_plot/latent_plots.py
index 43e3f650..0f5021ea 100644
--- a/GPy/plotting/gpy_plot/latent_plots.py
+++ b/GPy/plotting/gpy_plot/latent_plots.py
@@ -94,7 +94,11 @@ def plot_latent_scatter(self, labels=None,
     else:
         legend = find_best_layout_for_subplots(len(np.unique(labels)))[1]
     scatters = _plot_latent_scatter(canvas, X, sig_dims, labels, marker, num_samples, projection=projection, **kwargs)
-    return pl.show_canvas(canvas, dict(scatter=scatters), legend=legend)
+    return pl.add_to_canvas(canvas, dict(scatter=scatters), legend=legend)
+
+
+
+
 def plot_latent_inducing(self,  
                         which_indices=None,
                         legend=False,
@@ -125,7 +129,12 @@ def plot_latent_inducing(self,
     Z = self.Z.values
     labels = np.array(['inducing'] * Z.shape[0])
     scatters = _plot_latent_scatter(canvas, Z, sig_dims, labels, marker, num_samples, projection=projection, **kwargs)
-    return pl.show_canvas(canvas, dict(scatter=scatters), legend=legend)
+    return pl.add_to_canvas(canvas, dict(scatter=scatters), legend=legend)
+    
+
+
+
+
 
 def _plot_magnification(self, canvas, which_indices, Xgrid, 
                         xmin, xmax, resolution, updates,
@@ -182,12 +191,12 @@ def plot_magnification(self, labels=None, which_indices=None,
         legend = False
     scatters = _plot_latent_scatter(canvas, X, which_indices, labels, marker, num_samples, projection='2d', **scatter_kwargs or {})
     view = _plot_magnification(self, canvas, which_indices[:2], Xgrid, xmin, xmax, resolution, updates, mean, covariance, kern, **imshow_kwargs)
-    plots = pl.show_canvas(canvas, dict(scatter=scatters, imshow=view), 
+    retval = pl.add_to_canvas(canvas, dict(scatter=scatters, imshow=view), 
                            legend=legend, 
                            )
     _wait_for_updates(view, updates)
-    return plots
-
+    return retval
+    
 
 
 
@@ -245,10 +254,10 @@ def plot_latent(self, labels=None, which_indices=None,
         legend = False
     scatters = _plot_latent_scatter(canvas, X, which_indices, labels, marker, num_samples, projection='2d', **scatter_kwargs or {})
     view = _plot_latent(self, canvas, which_indices, Xgrid, xmin, xmax, resolution, updates, kern, **imshow_kwargs)
-    plots = pl.show_canvas(canvas, dict(scatter=scatters, imshow=view), legend=legend)
+    retval = pl.add_to_canvas(canvas, dict(scatter=scatters, imshow=view), legend=legend)
     _wait_for_updates(view, updates)
-    return plots
-
+    return retval
+    
 def _plot_steepest_gradient_map(self, canvas, which_indices, Xgrid, 
                         xmin, xmax, resolution, output_labels, updates,
                         kern=None, annotation_kwargs=None, 
@@ -310,10 +319,9 @@ def plot_steepest_gradient_map(self, output_labels=None, data_labels=None, which
         legend = False
     plots = dict(scatter=_plot_latent_scatter(canvas, X, which_indices, data_labels, marker, num_samples, **scatter_kwargs or {}))
     plots.update(_plot_steepest_gradient_map(self, canvas, which_indices, Xgrid, xmin, xmax, resolution, output_labels, updates, kern, annotation_kwargs=annotation_kwargs, **imshow_kwargs))
-    show = pl.show_canvas(canvas, plots, legend=legend)
+    retval = pl.add_to_canvas(canvas, plots, legend=legend)
     _wait_for_updates(plots['annotation'], updates)
-    return show
-
+    return retval
 
 
 
diff --git a/GPy/plotting/gpy_plot/plot_util.py b/GPy/plotting/gpy_plot/plot_util.py
index 954593d6..749b4a15 100644
--- a/GPy/plotting/gpy_plot/plot_util.py
+++ b/GPy/plotting/gpy_plot/plot_util.py
@@ -32,6 +32,15 @@ import numpy as np
 from scipy import sparse
 import itertools
 
+def in_ipynb():
+    try:
+        cfg = get_ipython().config 
+        if cfg['IPKernelApp']['parent_appname'] == 'ipython-notebook':
+            return True
+        else:
+            return False
+    except NameError:
+        return False
 
 def find_best_layout_for_subplots(num_subplots):
     r, c = 1, 1
@@ -109,7 +118,7 @@ def helper_for_plot_data(self, plot_limits, visible_dims, fixed_inputs, resoluti
     if fixed_inputs is None:
         fixed_inputs = []
     fixed_dims = get_fixed_dims(self, fixed_inputs)
-    free_dims = get_free_dims(self, visible_dims, fixed_dims)[:2]
+    free_dims = get_free_dims(self, visible_dims, fixed_dims)
     
     if len(free_dims) == 1:
         #define the frame on which to plot
diff --git a/GPy/plotting/matplot_dep/plot_definitions.py b/GPy/plotting/matplot_dep/plot_definitions.py
index 4678f1a8..15578446 100644
--- a/GPy/plotting/matplot_dep/plot_definitions.py
+++ b/GPy/plotting/matplot_dep/plot_definitions.py
@@ -80,7 +80,7 @@ class MatplotlibPlots(AbstractPlottingLibrary):
             if zlabel is not None: ax.set_zlabel(zlabel)
         return ax, kwargs
     
-    def show_canvas(self, ax, plots, legend=False, title=None, **kwargs):
+    def add_to_canvas(self, ax, plots, legend=False, title=None, **kwargs):
         ax.autoscale_view()
         fontdict=dict(family='sans-serif', weight='light', size=9)
         if legend is True:
@@ -89,8 +89,13 @@ class MatplotlibPlots(AbstractPlottingLibrary):
             #ax.legend(prop=fontdict)
             legend_ontop(ax, ncol=legend, fontdict=fontdict)
         if title is not None: ax.figure.suptitle(title)
+        return ax
+    
+    def show_canvas(self, ax, tight_layout=False, **kwargs):
+        if tight_layout:
+            ax.figure.tight_layout()
         ax.figure.canvas.draw()
-        return plots
+        return ax.figure
     
     def scatter(self, ax, X, Y, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, marker='o', **kwargs):
         if Z is not None:
diff --git a/GPy/plotting/plotly_dep/defaults.py b/GPy/plotting/plotly_dep/defaults.py
index f67b6e14..d10201d6 100644
--- a/GPy/plotting/plotly_dep/defaults.py
+++ b/GPy/plotting/plotly_dep/defaults.py
@@ -28,7 +28,6 @@
 # OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 #===============================================================================
 
-from matplotlib import cm
 from .. import Tango
 
 '''
@@ -43,33 +42,33 @@ it gives back an empty default, when defaults are not defined.
 '''
 
 # Data plots:
-data_1d = dict(lw=1.5, marker='x', edgecolor='k')
-data_2d = dict(s=35, edgecolors='none', linewidth=0., cmap=cm.get_cmap('hot'), alpha=.5)
-inducing_1d = dict(lw=0, s=500, facecolors=Tango.colorsHex['darkRed'])
-inducing_2d = dict(s=14, edgecolors='k', linewidth=.4, facecolors='white', alpha=.5)
-inducing_3d = dict(lw=.3, s=500, facecolors='white', edgecolors='k')
-xerrorbar = dict(color='k', fmt='none', elinewidth=.5, alpha=.5)
-yerrorbar = dict(color=Tango.colorsHex['darkRed'], fmt='none', elinewidth=.5, alpha=.5)
-
-# GP plots:
-meanplot_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=2)
-meanplot_2d = dict(cmap='hot', linewidth=.5)
-meanplot_3d = dict(linewidth=0, antialiased=True, cstride=1, rstride=1, cmap='hot', alpha=.3)
-samples_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=.3)
-samples_3d = dict(cmap='hot', alpha=.1, antialiased=True, cstride=1, rstride=1, linewidth=0)
-confidence_interval = dict(edgecolor=Tango.colorsHex['darkBlue'], linewidth=.5, color=Tango.colorsHex['lightBlue'],alpha=.2)
-density = dict(alpha=.5, color=Tango.colorsHex['lightBlue'])
-
-# GPLVM plots:
-data_y_1d = dict(linewidth=0, cmap='RdBu', s=40)
-data_y_1d_plot = dict(color='k', linewidth=1.5)
-
-# Kernel plots:
-ard = dict(edgecolor='k', linewidth=1.2)
-
-# Input plots:
-latent = dict(aspect='auto', cmap='Greys', interpolation='bicubic')
-gradient = dict(aspect='auto', cmap='RdBu', interpolation='nearest', alpha=.7)
-magnification = dict(aspect='auto', cmap='Greys', interpolation='bicubic')
-latent_scatter = dict(s=40, linewidth=.2, edgecolor='k', alpha=.9)
-annotation = dict(fontdict=dict(family='sans-serif', weight='light', fontsize=9), zorder=.3, alpha=.7)
\ No newline at end of file
+data_1d = dict(marker_kwargs=dict(linewidth=.7, ), marker='x', color='black')
+data_2d = dict(marker='o', cmap='Hot', marker_kwargs=dict(opacity=.5))
+# inducing_1d = dict(lw=0, s=500, facecolors=Tango.colorsHex['darkRed'])
+# inducing_2d = dict(s=14, edgecolors='k', linewidth=.4, facecolors='white', alpha=.5)
+# inducing_3d = dict(lw=.3, s=500, facecolors='white', edgecolors='k')
+# xerrorbar = dict(color='k', fmt='none', elinewidth=.5, alpha=.5)
+yerrorbar = dict(color=Tango.colorsHex['darkRed'], error_kwargs=dict(thickness=.5), opacity=.5)
+# 
+# # GP plots:
+# meanplot_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=2)
+# meanplot_2d = dict(cmap='hot', linewidth=.5)
+# meanplot_3d = dict(linewidth=0, antialiased=True, cstride=1, rstride=1, cmap='hot', alpha=.3)
+# samples_1d = dict(color=Tango.colorsHex['mediumBlue'], linewidth=.3)
+# samples_3d = dict(cmap='hot', alpha=.1, antialiased=True, cstride=1, rstride=1, linewidth=0)
+# confidence_interval = dict(edgecolor=Tango.colorsHex['darkBlue'], linewidth=.5, color=Tango.colorsHex['lightBlue'],alpha=.2)
+# density = dict(alpha=.5, color=Tango.colorsHex['lightBlue'])
+# 
+# # GPLVM plots:
+# data_y_1d = dict(linewidth=0, cmap='RdBu', s=40)
+# data_y_1d_plot = dict(color='k', linewidth=1.5)
+# 
+# # Kernel plots:
+# ard = dict(edgecolor='k', linewidth=1.2)
+# 
+# # Input plots:
+# latent = dict(aspect='auto', cmap='Greys', interpolation='bicubic')
+# gradient = dict(aspect='auto', cmap='RdBu', interpolation='nearest', alpha=.7)
+# magnification = dict(aspect='auto', cmap='Greys', interpolation='bicubic')
+# latent_scatter = dict(s=40, linewidth=.2, edgecolor='k', alpha=.9)
+# annotation = dict(fontdict=dict(family='sans-serif', weight='light', fontsize=9), zorder=.3, alpha=.7)
\ No newline at end of file
diff --git a/GPy/plotting/plotly_dep/plot_definitions.py b/GPy/plotting/plotly_dep/plot_definitions.py
index 0de1c7c4..41a56d3b 100644
--- a/GPy/plotting/plotly_dep/plot_definitions.py
+++ b/GPy/plotting/plotly_dep/plot_definitions.py
@@ -34,7 +34,22 @@ from . import defaults
 import itertools
 from plotly import tools
 from plotly import plotly as py
-from plotly.graph_objs import Scatter, Line, Data
+from plotly import matplotlylib
+from plotly.graph_objs import Scatter, Scatter3d, Line, Marker, ErrorX, ErrorY, Bar
+
+SYMBOL_MAP = {
+    'o': 'dot',
+    'v': 'triangle-down',
+    '^': 'triangle-up',
+    '<': 'triangle-left',
+    '>': 'triangle-right',
+    's': 'square',
+    '+': 'cross',
+    'x': 'x',
+    '*': 'x',  # no star yet in plotly!!
+    'D': 'diamond',
+    'd': 'diamond',
+}
 
 class PlotlyPlots(AbstractPlottingLibrary):
     def __init__(self):
@@ -42,87 +57,61 @@ class PlotlyPlots(AbstractPlottingLibrary):
         self._defaults = defaults.__dict__
         self.current_states = dict()
     
-    def figure(self, rows=1, cols=1, **kwargs):
-        if 'filename' not in kwargs:
-            print('PlotlyWarning: filename was not given, this may clutter your plotly workspace')
-            filename = None
-        else: filename = kwargs.pop('filename')
-        figure = tools.make_subplots(rows, cols, **kwargs)
-        self.current_states[hex(id(figure))] = dict(filename=filename)
-        index = 1
-        for i in rows:
-            for j in cols:
-                self.current_states[hex(id(figure))][(i,j)] = ('xaxis{}'.format(index), 'yaxis{}'.format(index))
-                index += 1  
+    def figure(self, rows=1, cols=1, specs=None, is_3d=False):
+        if specs is None:
+            specs = [[{'is_3d': is_3d}]*cols]*rows
+        figure = tools.make_subplots(rows, cols, specs=specs)
         return figure
     
     def new_canvas(self, figure=None, row=1, col=1, projection='2d', xlabel=None, ylabel=None, zlabel=None, title=None, xlim=None, ylim=None, zlim=None, **kwargs):
+        if 'filename' not in kwargs:
+            print('PlotlyWarning: filename was not given, this may clutter your plotly workspace')
+            filename = None
+        else: 
+            filename = kwargs.pop('filename')
         if figure is None:
-            figure = self.figure(**kwargs)
-        return (figure, self.current_states[hex(id(figure))][(row,col)]), kwargs
+            figure = self.figure(is_3d=projection=='3d')
+        self.current_states[hex(id(figure))] = dict(filename=filename)
+        return (figure, row, col), kwargs
     
-    def show_canvas(self, canvas, traces, legend=False, **kwargs):
-        fig = canvas[0] 
-        axis = fig['layout'][canvas[1]]
-        axis.update(title=title, label)
-        figure.add_traces(traces, row, col, **kwargs)
-        
-        # If shared_axes is False (default) use list_of_domains
-        # This is used for insets and irregular layouts
-        if not shared_xaxes and not shared_yaxes:
-            x_dom = list_of_domains[::2]
-            y_dom = list_of_domains[1::2]
-            subtitle_pos_x = []
-            subtitle_pos_y = []
-            for x_domains in x_dom:
-                subtitle_pos_x.append(sum(x_domains) / 2)
-            for y_domains in y_dom:
-                subtitle_pos_y.append(y_domains[1])
-        # If shared_axes is True the domin of each subplot is not returned so the
-        # title position must be calculated for each subplot
+    def add_to_canvas(self, canvas, traces, legend=False, **kwargs):
+        figure, row, col = canvas
+        def recursive_append(traces):
+            for _, trace in traces.items():
+                if isinstance(trace, (tuple, list)):
+                    for t in trace:
+                        figure.append_trace(t, row, col)
+                elif isinstance(trace, dict):
+                    recursive_append(trace)
+                else:
+                    figure.append_trace(trace, row, col)
+        recursive_append(traces)
+        figure.layout['showlegend'] = legend
+        return canvas
+    
+    def show_canvas(self, canvas, **kwargs):
+        figure, _, _ = canvas
+        from ..gpy_plot.plot_util import in_ipynb
+        if in_ipynb():
+            py.iplot(figure, filename=self.current_states[hex(id(figure))]['filename'])
         else:
-            subtitle_pos_x = [None] * cols
-            subtitle_pos_y = [None] * rows
-            delt_x = (x_e - x_s)
-            for index in range(cols):
-                subtitle_pos_x[index] = ((delt_x / 2) +
-                                         ((delt_x + horizontal_spacing) * index))
-            subtitle_pos_x *= rows
-            for index in range(rows):
-                subtitle_pos_y[index] = (1 - ((y_e + vertical_spacing) * index))
-            subtitle_pos_y *= cols
-            subtitle_pos_y = sorted(subtitle_pos_y, reverse=True)
+            py.plot(figure, filename=self.current_states[hex(id(figure))]['filename'])
+        return figure
     
-        plot_titles = []
-        for index in range(len(subplot_titles)):
-            if not subplot_titles[index]:
-                pass
-            else:
-                plot_titles.append({'y': subtitle_pos_y[index],
-                                    'xref': 'paper',
-                                    'x': subtitle_pos_x[index],
-                                    'yref': 'paper',
-                                    'text': subplot_titles[index],
-                                    'showarrow': False,
-                                    'font': graph_objs.Font(size=16),
-                                    'xanchor': 'center',
-                                    'yanchor': 'bottom'
-                                    })
-    
-                layout['annotations'] = plot_titles
+    def scatter(self, ax, X, Y, Z=None, color=Tango.colorsHex['mediumBlue'], cmap=None, label=None, marker='o', marker_kwargs=None, **kwargs):
         try:
-            url = py.iplot(figure, self.current_filename)
+            marker = SYMBOL_MAP[marker]
         except:
-            url = py.plot(figure, self.current_filename)
-        return url
-    
-    def scatter(self, ax, X, Y, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, marker='o', **kwargs):
-        
-    
-    def plot(self, ax, X, Y, Z=None, color=None, label=None, **kwargs):
+            #not matplotlib marker
+            pass
         if Z is not None:
-            return ax.plot(X, Y, color=color, zs=Z, label=label, **kwargs)
-        return ax.plot(X, Y, color=color, label=label, **kwargs)
+            return Scatter3d(x=X, y=Y, z=Z, mode='markers', marker=Marker(color=color, symbol=marker, colorscale=cmap, **marker_kwargs or {}), name=label, **kwargs)
+        return Scatter(x=X, y=Y, mode='markers', marker=Marker(color=color, symbol=marker, colorscale=cmap, **marker_kwargs or {}), name=label, **kwargs)
+    
+    def plot(self, ax, X, Y, Z=None, color=None, label=None, line_kwargs=None, **kwargs):
+        if Z is not None:
+            return Scatter3d(x=X, y=Y, z=Z, mode='lines', line=Line(color=color, **line_kwargs or {}), name=label, **kwargs)
+        return Scatter(x=X, y=Y, mode='lines', line=Line(color=color, **line_kwargs or {}), name=label, **kwargs)
 
     def plot_axis_lines(self, ax, X, color=Tango.colorsHex['mediumBlue'], label=None, **kwargs):
         from matplotlib import transforms
@@ -137,26 +126,31 @@ class PlotlyPlots(AbstractPlottingLibrary):
             return ax.scatter(X[:,0], X[:,1], ax.get_zlim()[0], c=color, label=label, **kwargs)
         return ax.scatter(X, np.zeros_like(X), c=color, label=label, **kwargs)
 
-    def barplot(self, ax, x, height, width=0.8, bottom=0, color=Tango.colorsHex['mediumBlue'], label=None, **kwargs):
-        if 'align' not in kwargs:
-            kwargs['align'] = 'center'
-        return ax.bar(left=x, height=height, width=width,
-               bottom=bottom, label=label, color=color,  
-               **kwargs)
+    def barplot(self, canvas, x, height, width=0.8, bottom=0, color=Tango.colorsHex['mediumBlue'], label=None, **kwargs):
+        figure, _, _ = canvas
+        if 'barmode' in kwargs:
+            figure.layout['barmode'] = kwargs.pop('barmode')
+        return Bar(x=x, y=height, marker=Marker(color=color), name=label)
         
-    def xerrorbar(self, ax, X, Y, error, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, **kwargs):
-        if not('linestyle' in kwargs or 'ls' in kwargs):
-            kwargs['ls'] = 'none'
+    def xerrorbar(self, ax, X, Y, error, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, error_kwargs=None, **kwargs):
+        error_kwargs = error_kwargs or {}
+        if (error.shape[0] == 2) and (error.ndim == 2):
+            error_kwargs.update(dict(array=error[1], arrayminus=error[0], symmetric=False))
+        else:
+            error_kwargs.update(dict(array=error, symmetric=True))
         if Z is not None:
-            return ax.errorbar(X, Y, Z, xerr=error, ecolor=color, label=label, **kwargs)
-        return ax.errorbar(X, Y, xerr=error, ecolor=color, label=label, **kwargs)
-    
-    def yerrorbar(self, ax, X, Y, error, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, **kwargs):
-        if not('linestyle' in kwargs or 'ls' in kwargs):
-            kwargs['ls'] = 'none'
+            return Scatter3d(x=X, y=Y, z=Z, mode='markers', error_x=ErrorX(color=color, **error_kwargs or {}), marker=Marker(size='0'), name=label, **kwargs)
+        return Scatter(x=X, y=Y, mode='markers', error_x=ErrorX(color=color, **error_kwargs or {}), marker=Marker(size='0'), name=label, **kwargs)
+
+    def yerrorbar(self, ax, X, Y, error, Z=None, color=Tango.colorsHex['mediumBlue'], label=None, error_kwargs=None, **kwargs):
+        error_kwargs = error_kwargs or {}
+        if (error.shape[0] == 2) and (error.ndim == 2):
+            error_kwargs.update(dict(array=error[1], arrayminus=error[0], symmetric=False))
+        else:
+            error_kwargs.update(dict(array=error, symmetric=True))
         if Z is not None:
-            return ax.errorbar(X, Y, Z, yerr=error, ecolor=color, label=label, **kwargs)
-        return ax.errorbar(X, Y, yerr=error, ecolor=color, label=label, **kwargs)
+            return Scatter3d(x=X, y=Y, z=Z, mode='markers', error_y=ErrorY(color=color, **error_kwargs or {}), marker=Marker(size='0'), name=label, **kwargs)
+        return Scatter(x=X, y=Y, mode='markers', error_y=ErrorY(color=color, **error_kwargs or {}), marker=Marker(size='0'), name=label, **kwargs)
     
     def imshow(self, ax, X, extent=None, label=None, vmin=None, vmax=None, **imshow_kwargs):
         if 'origin' not in imshow_kwargs:
@@ -167,9 +161,8 @@ class PlotlyPlots(AbstractPlottingLibrary):
         return ax.imshow(X, label=label, extent=extent, vmin=vmin, vmax=vmax, **imshow_kwargs)
 
     def imshow_interact(self, ax, plot_function, extent=None, label=None, resolution=None, vmin=None, vmax=None, **imshow_kwargs):
-        if 'origin' not in imshow_kwargs:
-            imshow_kwargs['origin'] = 'lower'
-        return ImshowController(ax, plot_function, extent, resolution=resolution, vmin=vmin, vmax=vmax, **imshow_kwargs)
+        # TODO stream interaction?
+        super(PlotlyPlots, self).imshow_interact(ax, plot_function)
 
     def annotation_heatmap(self, ax, X, annotation, extent=None, label=None, imshow_kwargs=None, **annotation_kwargs):
         imshow_kwargs = imshow_kwargs or {}
diff --git a/GPy/testing/plotting_tests.py b/GPy/testing/plotting_tests.py
index 88aa32c6..4df227ec 100644
--- a/GPy/testing/plotting_tests.py
+++ b/GPy/testing/plotting_tests.py
@@ -127,7 +127,9 @@ def test_threed():
     m.optimize()
     m.plot_data(projection='3d')
     m.plot_mean(projection='3d')
-    for do_test in _image_comparison(baseline_images=['gp_3d_{}'.format(sub) for sub in ["data", "mean"]], extensions=extensions):
+    m.plot_samples(projection='3d', samples=1)
+    m.plot_samples(projection='3d', plot_raw=False, samples=1)
+    for do_test in _image_comparison(baseline_images=['gp_3d_{}'.format(sub) for sub in ["data", "mean", "samples", "samples_lik"]], extensions=extensions):
         yield (do_test, )
 
 def test_sparse():
@@ -164,9 +166,9 @@ def test_sparse_classification():
     Y = f+np.random.normal(0, .1, f.shape)
     m = GPy.models.SparseGPClassification(X, Y>Y.mean())
     m.optimize()
-    m.plot(plot_raw=False, apply_link=False)
-    m.plot(plot_raw=True, apply_link=False)
-    m.plot(plot_raw=True, apply_link=True)
+    m.plot(plot_raw=False, apply_link=False, samples_likelihood=3)
+    m.plot(plot_raw=True, apply_link=False, samples=3)
+    m.plot(plot_raw=True, apply_link=True, samples=3)
     for do_test in _image_comparison(baseline_images=['sparse_gp_class_{}'.format(sub) for sub in ["likelihood", "raw", 'raw_link']], extensions=extensions):
         yield (do_test, )
 
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 c188bb4e..2b4eb552 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.png b/GPy/testing/plotting_tests/baseline/sparse_gp_class_raw_link.png
index 453941fc..7a969dc7 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/travis_tests.py b/travis_tests.py
index 06d54de6..2d71de97 100644
--- a/travis_tests.py
+++ b/travis_tests.py
@@ -38,5 +38,5 @@ matplotlib.rcParams[u'figure.figsize'] = (4,3)
 matplotlib.rcParams[u'text.usetex'] = False
 
 import nose 
-nose.main('GPy', defaultTest='GPy/testing')  
+nose.main('GPy', defaultTest='GPy/testing/plotting_tests.py')