@@ -259,7 +259,7 @@ class Learner2D(BaseLearner):

                 self._vdim = len(value)

             except TypeError:

                 self._vdim = 1

-        return self._vdim if self._vdim is not None else 1

+        return self._vdim or 1

     @property

     def bounds_are_done(self):

@@ -385,10 +385,34 @@ class Learner2D(BaseLearner):

                 self._stack[p] = np.inf

     def plot(self, n=None, tri_alpha=0):

+        """Plot the Learner2D's current state.

+

+        This plot function interpolates the data on a regular grid.

+        The gridspacing is evaluated by checking the size of the smallest

+        triangle.

+

+        Parameters

+        ----------

+        n : int

+            Number of points in x and y. If None (default) this number is

+            evaluated by looking at the size of the smallest triangle.

+        tri_alpha : float

+            The opacity (0 <= tri_alpha <= 1) of the triangles overlayed on

+            top of the image. By default the triangulation is not visible.

+

+        Returns

+        -------

+        plot : holoviews.Overlay or holoviews.HoloMap

+            A `holoviews.Overlay` of `holoviews.Image * holoviews.EdgePaths`.

+            If the `learner.function` returns a vector output, a

+            `holoviews.HoloMap` of the `holoviews.Overlay`s wil be returned.

+

+        Notes

+        -----

+        The plot object that is returned if `learner.function` returns a

+        vector *cannot* be used with the live_plotting functionality.

+        """

         hv = ensure_holoviews()

-        if self.vdim > 1:

-            raise NotImplementedError('holoviews currently does not support',

-                                      '3D surface plots in bokeh.')

         x, y = self.bounds

         lbrt = x[0], y[0], x[1], y[1]

@@ -404,7 +428,12 @@ class Learner2D(BaseLearner):

             x = y = np.linspace(-0.5, 0.5, n)

             z = ip(x[:, None], y[None, :] * self.aspect_ratio).squeeze()

-            im = hv.Image(np.rot90(z), bounds=lbrt)

+            if self.vdim > 1:

+                ims = {i: hv.Image(np.rot90(z[:, :, i]), bounds=lbrt)

+                       for i in range(z.shape[-1])}

+                im = hv.HoloMap(ims)

+            else:

+                im = hv.Image(np.rot90(z), bounds=lbrt)

             if tri_alpha:

                 points = self._unscale(ip.tri.points[ip.tri.vertices])