Repository: holoviz/holoviews
Pull request #6336 (2024-07-22 04:59:18)

fix a problem where linked selections were resulting in repeated columns
While working with `link_selections` I ran into a problem where trying to link two plots that came from the same dataset would result in an error because at some point during the creation of the selection expression the columns would get duplicated. Here's an MRE that illustrates this problem:
```python
import holoviews as hv
from holoviews.selection import link_selections
import numpy as np
import pandas as pd
import panel as pn

hv.extension('bokeh')
pn.extension()

# some fake data
data = np.random.default_rng(seed=42).normal(size=(100, 3))
idx = np.arange(100)
cols = ['x', 'y', 'z']
df = pd.DataFrame(data, columns=cols)
df['id'] = idx

# want to link two plots across the `id` column
ls = link_selections.instance(index_cols=['id'])

img1 = hv.Points(df, kdims=['x', 'y'], vdims=['id']).opts(tools=['box_select'])
img2 = hv.Points(df, kdims=['x', 'z'], vdims=['id']).opts(tools=['box_select'])

pn.Row(ls(img1), ls(img2)).servable()
```

This results in the following slightly truncated error:

```
File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/element/selection.py", line 334, in _get_selection_expr_for_stream_value
    expr, _, _ = self._get_index_selection(kwargs['index'], index_cols)
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/element/selection.py", line 39, in _get_index_selection
    vals = dim(index_dim).apply(ds.iloc[index], expanded=False)
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/data/interface.py", line 33, in __getitem__
    res = self._perform_getitem(self.dataset, index)
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/data/interface.py", line 98, in _perform_getitem
    return dataset.clone(data, kdims=kdims, vdims=vdims, datatype=datatype)
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/data/__init__.py", line 1203, in clone
    return super().clone(data, shared_data, new_type, *args, **overrides)
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/dimension.py", line 561, in clone
    return clone_type(data, *args, **{k:v for k,v in settings.items()
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/data/__init__.py", line 329, in __init__
    initialized = Interface.initialize(type(self), data, kdims, vdims,
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/data/interface.py", line 253, in initialize
    (data, dims, extra_kws) = interface.init(eltype, data, kdims, vdims)
  File "/Users/jerry/Development/project/venv/lib/python3.10/site-packages/holoviews/core/data/pandas.py", line 83, in init
    raise DataError('Dimensions may not reference duplicated DataFrame '
holoviews.core.data.interface.DataError: Dimensions may not reference duplicated DataFrame columns (found duplicate 'id' columns). If you want to plot a column against itself simply declare two dimensions with the same name.
```

I traced this problem down into the `_get_index_selection` function, where the dataset is cloned on line 35:
```python
ds = self.clone(kdims=index_cols, new_type=Dataset)
```
Because cloning will take all the vdims if that parameter is left unspecified, if there is overlap between kdims and vdims, the selection on line 38 will throw an error because the columns are duplicated in the resulting dataframe:
```python
vals = dim(index_dim).apply(ds.iloc[index], expanded=False)
```

This PR makes the slight change that vdims are explicitly set to be cloned only if they do not overlap with the kdims. Then the selection is made on both rows and columns to ensure that each column appears only once.

This fixes the problem I observed in the MRE and doesn't seem to break anything else. It's my first attempt to contribute to HoloViews so I would appreciate any pointers on how to properly test this or anything else it might affect.

jerry-kobold: I'm wondering if anyone has any thoughts about this proposal? Just a sketch of whether this is on the right track or not.
grapesmoker: Just curious if anyone has any thoughts on this PR, even a hint about whether this is a sensible thing to attempt or not.
ahuang11: Thanks for fixing this! I'm not terribly familiar with this part of the code base so I can't comment on whether it's the right way or not, but I think other maintainers would appreciate it if you could add a UI test (something like https://github.com/holoviz/holoviews/blob/main/holoviews/tests/ui/bokeh/test_callback.py#L351-L378 basically copying your minimal example, and adding a few assert statements) and perhaps a short inline comment about what you're doing.

Test failures seem unrelated.
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/holoviz/holoviews/pull/6336?dropdown=coverage&src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Project coverage is 88.48%. Comparing base [(`a28189d`)](https://app.codecov.io/gh/holoviz/holoviews/commit/a28189de812f17074934bfe0990f95a2fb037dd8?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz) to head [(`1ce8ac6`)](https://app.codecov.io/gh/holoviz/holoviews/commit/1ce8ac63c06ffb1b15299e124b52dbd55100a95b?dropdown=coverage&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz).

<details><summary>Additional details and impacted files</summary>


```diff
@@            Coverage Diff             @@
##             main    #6336      +/-   ##
==========================================
- Coverage   88.48%   88.48%   -0.01%     
==========================================
  Files         323      323              
  Lines       68162    68164       +2     
==========================================
  Hits        60313    60313              
- Misses       7849     7851       +2     
```

</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/holoviz/holoviews/pull/6336?dropdown=coverage&src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=holoviz).


Patch:
diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml
index f63292851..261f73fe2 100644
--- a/.github/workflows/build.yaml
+++ b/.github/workflows/build.yaml
@@ -12,7 +12,7 @@ on:
 
 defaults:
   run:
-    shell: bash -el {0}
+    shell: bash -e {0}
 
 env:
   PYTHON_VERSION: "3.11"
@@ -33,14 +33,14 @@ jobs:
     name: Pixi lock
     runs-on: ubuntu-latest
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_lock@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_lock@v0
 
   conda_build:
     name: Build Conda
     needs: [pixi_lock]
     runs-on: "ubuntu-latest"
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_install@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_install@v0
         with:
           environments: "build"
           download-data: false
@@ -59,6 +59,9 @@ jobs:
     runs-on: ubuntu-latest
     needs: [conda_build, waiting_room]
     if: github.event_name == 'push' && startsWith(github.ref, 'refs/tags/')
+    defaults:
+      run:
+        shell: bash -el {0}
     steps:
       - uses: actions/download-artifact@v4
         with:
@@ -88,7 +91,7 @@ jobs:
     needs: [pixi_lock]
     runs-on: "ubuntu-latest"
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_install@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_install@v0
         with:
           environments: "build"
           download-data: false
@@ -118,10 +121,6 @@ jobs:
         run: python -m pip install dist/*.whl
       - name: Import package
         run: python -c "import $PACKAGE; print($PACKAGE.__version__)"
-      - name: Install test dependencies
-        run: python -m pip install "$(ls dist/*whl)[tests]"
-      - name: Test package
-        run: python -m pytest --pyargs $PACKAGE --color=yes
 
   pip_publish:
     name: Publish PyPI
diff --git a/.github/workflows/docs.yaml b/.github/workflows/docs.yaml
index 660cbf350..858db8fd3 100644
--- a/.github/workflows/docs.yaml
+++ b/.github/workflows/docs.yaml
@@ -22,14 +22,14 @@ on:
 
 defaults:
   run:
-    shell: bash -el {0}
+    shell: bash -e {0}
 
 jobs:
   pixi_lock:
     name: Pixi lock
     runs-on: ubuntu-latest
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_lock@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_lock@v0
 
   docs_build:
     name: Build Documentation
@@ -41,7 +41,7 @@ jobs:
     env:
       GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_install@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_install@v0
         with:
           environments: docs
       - name: Build documentation
diff --git a/.github/workflows/nightly_lock.yaml b/.github/workflows/nightly_lock.yaml
index 9b93b85a2..00e670c88 100644
--- a/.github/workflows/nightly_lock.yaml
+++ b/.github/workflows/nightly_lock.yaml
@@ -13,7 +13,7 @@ jobs:
     runs-on: ubuntu-latest
     timeout-minutes: 5
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_lock@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_lock@v0
       - name: Upload lock-file to S3
         if: "!github.event.pull_request.head.repo.fork"
         env:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index edd89360a..37fa992a5 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -1,11 +1,8 @@
-# This is the configuration for pre-commit, a local framework for managing pre-commit hooks
-#   Check out the docs at: https://pre-commit.com/
-
 exclude: (\.min\.js$|\.svg$|\.html$)
-default_stages: [commit]
+default_stages: [pre-commit]
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.6.0
+    rev: v5.0.0
     hooks:
       - id: end-of-file-fixer
       - id: trailing-whitespace
@@ -20,7 +17,7 @@ repos:
       - id: check-json
       - id: detect-private-key
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.5.5
+    rev: v0.6.9
     hooks:
       - id: ruff
         files: holoviews/|scripts/
@@ -55,7 +52,7 @@ repos:
           - yaml
           - css
   - repo: https://github.com/hoxbro/taplo-pre-commit
-    rev: v0.7.0
+    rev: v0.9.3
     hooks:
       - id: taplo
         args:
diff --git a/doc/releases.md b/doc/releases.md
index 8c59ab64f..7354f92d3 100644
--- a/doc/releases.md
+++ b/doc/releases.md
@@ -256,7 +256,7 @@ Enhancements:
 - Drop requirements for `OrderedDict` ([#5867](https://github.com/holoviz/holoviews/pull/5867), [#5890](https://github.com/holoviz/holoviews/pull/5890), [#5925](https://github.com/holoviz/holoviews/pull/5925))
 - Allow to link to an `Overlay` ([#5881](https://github.com/holoviz/holoviews/pull/5881))
 - Use contourpy for contour calculations ([#5910](https://github.com/holoviz/holoviews/pull/5910))
-- Use browser information to set `pixel_density` in `ResampleOperation2D` ([#5947](https://github.com/holoviz/holoviews/pull/5947))
+- Use browser information to set `pixel_ratio` in `ResampleOperation2D` ([#5947](https://github.com/holoviz/holoviews/pull/5947))
 - Avoid bounce back of events for `Range{X,Y,XY}` streams ([#5946](https://github.com/holoviz/holoviews/pull/5946))
 
 Bug fixes:
diff --git a/holoviews/__init__.py b/holoviews/__init__.py
index 1ca017cdf..a8b3f5868 100644
--- a/holoviews/__init__.py
+++ b/holoviews/__init__.py
@@ -73,9 +73,9 @@ In a notebook or ipython environment the usual
 To ask the community go to https://discourse.holoviz.org/.
 To report issues go to https://github.com/holoviz/holoviews.
 """
+import builtins
 import os
 import sys
-import warnings
 
 import param
 
@@ -114,12 +114,7 @@ from .util.warnings import (  # noqa: F401
     HoloviewsUserWarning,
 )
 
-# Suppress warnings generated by NumPy in matplotlib
-# Expected to be fixed in next matplotlib release
-warnings.filterwarnings("ignore",
-                        message="elementwise comparison failed; returning scalar instead")
-
-if "IPython" in sys.modules:
+if hasattr(builtins, "__IPYTHON__"):
     from .ipython import notebook_extension
     extension = notebook_extension
 else:
@@ -181,7 +176,7 @@ def help(obj, visualization=True, ansi=True, backend=None,
         pydoc.help(obj)
 
 
-del os, sys, rcfile, warnings
+del builtins, os, rcfile, sys
 
 def __getattr__(name):
     if name == "annotate":
diff --git a/holoviews/core/data/cudf.py b/holoviews/core/data/cudf.py
index 0cd8f6d0e..f316e5e6e 100644
--- a/holoviews/core/data/cudf.py
+++ b/holoviews/core/data/cudf.py
@@ -205,7 +205,8 @@ class cuDFInterface(PandasInterface):
                 with warnings.catch_warnings():
                     warnings.filterwarnings('ignore', r'invalid value encountered')
                     if sel.start is not None:
-                        new_masks.append(sel.start <= arr)
+                        # Comparison has to be in this order due to issues with datetime comparison (see #6407)
+                        new_masks.append(arr >= sel.start)
                     if sel.stop is not None:
                         new_masks.append(arr < sel.stop)
                 if not new_masks:
diff --git a/holoviews/core/data/dask.py b/holoviews/core/data/dask.py
index 9c0476b06..4d0034381 100644
--- a/holoviews/core/data/dask.py
+++ b/holoviews/core/data/dask.py
@@ -218,7 +218,7 @@ class DaskInterface(PandasInterface):
         for coord in indices:
             if any(isinstance(c, float) and np.isnan(c) for c in coord):
                 continue
-            if len(coord) == 1:
+            if len(coord) == 1 and not util.PANDAS_GE_220:
                 coord = coord[0]
             group = group_type(groupby.get_group(coord), **group_kwargs)
             data.append((coord, group))
diff --git a/holoviews/core/data/ibis.py b/holoviews/core/data/ibis.py
index fd3c11ac0..1023067a7 100644
--- a/holoviews/core/data/ibis.py
+++ b/holoviews/core/data/ibis.py
@@ -27,6 +27,10 @@ def ibis4():
 def ibis5():
     return ibis_version() >= Version("5.0")
 
+@lru_cache
+def ibis9_5():
+    return ibis_version() >= Version("9.5")
+
 
 class IbisInterface(Interface):
 
@@ -457,7 +461,8 @@ class IbisInterface(Interface):
         data = dataset.data
         if all(util.isscalar(s) or len(s) == 1 for s in samples):
             items = [s[0] if isinstance(s, tuple) else s for s in samples]
-            return data[data[dims[0].name].isin(items)]
+            subset = data[dims[0].name].isin(items)
+            return data.filter(subset) if ibis9_5() else data[subset]
 
         predicates = None
         for sample in samples:
diff --git a/holoviews/core/data/xarray.py b/holoviews/core/data/xarray.py
index 4daefd778..9da2170c8 100644
--- a/holoviews/core/data/xarray.py
+++ b/holoviews/core/data/xarray.py
@@ -93,7 +93,7 @@ class XArrayInterface(GridInterface):
         if isinstance(data, xr.DataArray):
             kdim_len = len(kdim_param.default) if kdims is None else len(kdims)
             vdim_len = len(vdim_param.default) if vdims is None else len(vdims)
-            if vdim_len > 1 and kdim_len == len(data.dims)-1 and data.shape[-1] == vdim_len:
+            if kdim_len == len(data.dims)-1 and data.shape[-1] == vdim_len:
                 packed = True
             elif vdims:
                 vdim = vdims[0]
@@ -125,7 +125,33 @@ class XArrayInterface(GridInterface):
                     vdims = [vdim]
                     data = data.to_dataset(name=vdim.name)
 
-        if not isinstance(data, (xr.Dataset, xr.DataArray)):
+        if isinstance(data, (xr.Dataset, xr.DataArray)):
+            # Started to warn in xarray 2023.12.0:
+            # The return type of `Dataset.dims` will be changed to return a
+            # set of dimension names in future, in order to be more consistent
+            # with `DataArray.dims`. To access a mapping from dimension names to
+            # lengths, please use `Dataset.sizes`.
+            data_info =  data.sizes if hasattr(data, "sizes") else data.dims
+            if not data.coords:
+                data = data.assign_coords(**{k: range(v) for k, v in data_info.items()})
+            if vdims is None:
+                vdims = list(data.data_vars)
+            if kdims is None:
+                xrdims = list(data_info)
+                xrcoords = list(data.coords)
+                kdims = [name for name in data.indexes.keys()
+                         if isinstance(data[name].data, np.ndarray)]
+                kdims = sorted(kdims, key=lambda x: (xrcoords.index(x) if x in xrcoords else float('inf'), x))
+                if packed:
+                    kdims = kdims[:-1]
+                elif set(xrdims) != set(kdims):
+                    virtual_dims = [xd for xd in xrdims if xd not in kdims]
+                    for c in data.coords:
+                        if c not in kdims and set(data[c].dims) == set(virtual_dims):
+                            kdims.append(c)
+            kdims = [retrieve_unit_and_label(kd) for kd in kdims]
+            vdims = [retrieve_unit_and_label(vd) for vd in vdims]
+        else:
             if kdims is None:
                 kdims = kdim_param.default
             if vdims is None:
@@ -176,32 +202,6 @@ class XArrayInterface(GridInterface):
                         arr = xr.DataArray(arr, coords=coords, **xr_kwargs)
                     arrays[vdim.name] = arr
                 data = xr.Dataset(arrays)
-        else:
-            # Started to warn in xarray 2023.12.0:
-            # The return type of `Dataset.dims` will be changed to return a
-            # set of dimension names in future, in order to be more consistent
-            # with `DataArray.dims`. To access a mapping from dimension names to
-            # lengths, please use `Dataset.sizes`.
-            data_info =  data.sizes if hasattr(data, "sizes") else data.dims
-            if not data.coords:
-                data = data.assign_coords(**{k: range(v) for k, v in data_info.items()})
-            if vdims is None:
-                vdims = list(data.data_vars)
-            if kdims is None:
-                xrdims = list(data_info)
-                xrcoords = list(data.coords)
-                kdims = [name for name in data.indexes.keys()
-                         if isinstance(data[name].data, np.ndarray)]
-                kdims = sorted(kdims, key=lambda x: (xrcoords.index(x) if x in xrcoords else float('inf'), x))
-                if packed:
-                    kdims = kdims[:-1]
-                elif set(xrdims) != set(kdims):
-                    virtual_dims = [xd for xd in xrdims if xd not in kdims]
-                    for c in data.coords:
-                        if c not in kdims and set(data[c].dims) == set(virtual_dims):
-                            kdims.append(c)
-            kdims = [retrieve_unit_and_label(kd) for kd in kdims]
-            vdims = [retrieve_unit_and_label(vd) for vd in vdims]
 
         not_found = []
         for d in kdims:
@@ -446,7 +446,11 @@ class XArrayInterface(GridInterface):
         Given a dataset object and data in the appropriate format for
         the interface, return a simple scalar.
         """
-        if not cls.packed(dataset) and len(data.data_vars) == 1:
+        if cls.packed(dataset):
+            array = data.squeeze()
+            if len(array.shape) == 0:
+                return array.item()
+        elif len(data.data_vars) == 1:
             array = data[dataset.vdims[0].name].squeeze()
             if len(array.shape) == 0:
                 return array.item()
diff --git a/holoviews/core/util.py b/holoviews/core/util.py
index 590b7ab6e..cd25fdd6a 100644
--- a/holoviews/core/util.py
+++ b/holoviews/core/util.py
@@ -44,6 +44,7 @@ pandas_version = Version(pd.__version__)
 
 NUMPY_GE_200 = numpy_version >= Version("2")
 PANDAS_GE_210 = pandas_version >= Version("2.1")
+PANDAS_GE_220 = pandas_version >= Version("2.2")
 
 # Types
 generator_types = (zip, range, types.GeneratorType)
diff --git a/holoviews/element/graphs.py b/holoviews/element/graphs.py
index 30578b86c..467a8632f 100644
--- a/holoviews/element/graphs.py
+++ b/holoviews/element/graphs.py
@@ -269,7 +269,7 @@ argument to specify a selection specification""")
 
         index_dim = self.nodes.kdims[2].name
         dimensions = self.kdims+self.vdims
-        node_selection = {index_dim: v for k, v in selection.items()  # noqa: RUF011
+        node_selection = {index_dim: v for k, v in selection.items()  # noqa: B035
                           if k in self.kdims}
         if selection_expr:
             mask = selection_expr.apply(self.nodes, compute=False, keep_index=True)
diff --git a/holoviews/element/selection.py b/holoviews/element/selection.py
index a6c8ff43a..349c3ee66 100644
--- a/holoviews/element/selection.py
+++ b/holoviews/element/selection.py
@@ -32,10 +32,12 @@ class SelectionIndexExpr:
         self._index_skip = True
         if not index:
             return None, None, None
-        ds = self.clone(kdims=index_cols, new_type=Dataset)
+        clone_vdims = [vdim.name for vdim in self.vdims if vdim.name not in index_cols]
+        cols = clone_vdims + index_cols
+        ds = self.clone(kdims=index_cols, vdims=clone_vdims, new_type=Dataset)
         if len(index_cols) == 1:
             index_dim = index_cols[0]
-            vals = dim(index_dim).apply(ds.iloc[index], expanded=False)
+            vals = dim(index_dim).apply(ds.iloc[index, cols], expanded=False)
             if vals.dtype.kind == 'O' and all(isinstance(v, np.ndarray) for v in vals):
                 vals = [v for arr in vals for v in util.unique_iterator(arr)]
             expr = dim(index_dim).isin(list(util.unique_iterator(vals)))
@@ -43,7 +45,7 @@ class SelectionIndexExpr:
             get_shape = dim(self.dataset.get_dimension(index_cols[0]), np.shape)
             index_cols = [dim(self.dataset.get_dimension(c), np.ravel) for c in index_cols]
             vals = dim(index_cols[0], util.unique_zip, *index_cols[1:]).apply(
-                ds.iloc[index], expanded=True, flat=True
+                ds.iloc[index, cols], expanded=True, flat=True
             )
             contains = dim(index_cols[0], util.lzip, *index_cols[1:]).isin(vals, object=True)
             expr = dim(contains, np.reshape, get_shape)
diff --git a/holoviews/ipython/__init__.py b/holoviews/ipython/__init__.py
index 7727ee55c..c71b81c63 100644
--- a/holoviews/ipython/__init__.py
+++ b/holoviews/ipython/__init__.py
@@ -2,6 +2,7 @@ import os
 from unittest import SkipTest
 
 import param
+from bokeh.settings import settings as bk_settings
 from IPython.core.completer import IPCompleter
 from IPython.display import HTML, publish_display_data
 from param import ipython as param_ext
@@ -164,6 +165,7 @@ class notebook_extension(extension):
             Store.set_display_hook('html+js', LabelledData, pprint_display)
             Store.set_display_hook('png', LabelledData, png_display)
             Store.set_display_hook('svg', LabelledData, svg_display)
+            bk_settings.simple_ids.set_value(False)
             notebook_extension._loaded = True
 
         css = ''
diff --git a/holoviews/operation/datashader.py b/holoviews/operation/datashader.py
index b24d28c7f..9dc6c010a 100644
--- a/holoviews/operation/datashader.py
+++ b/holoviews/operation/datashader.py
@@ -69,9 +69,10 @@ from ..element.util import connect_tri_edges_pd
 from ..streams import PointerXY
 from .resample import LinkableOperation, ResampleOperation2D
 
-ds_version = Version(ds.__version__)
-ds15 = ds_version >= Version('0.15.1')
-ds16 = ds_version >= Version('0.16.0')
+DATASHADER_VERSION = Version(ds.__version__).release
+DATASHADER_GE_0_14_0 = DATASHADER_VERSION >= (0, 14, 0)
+DATASHADER_GE_0_15_1 = DATASHADER_VERSION >= (0, 15, 1)
+DATASHADER_GE_0_16_0 = DATASHADER_VERSION >= (0, 16, 0)
 
 
 class AggregationOperation(ResampleOperation2D):
@@ -114,7 +115,7 @@ class AggregationOperation(ResampleOperation2D):
 
     @classmethod
     def _get_aggregator(cls, element, agg, add_field=True):
-        if ds15:
+        if DATASHADER_GE_0_15_1:
             agg_types = (rd.count, rd.any, rd.where)
         else:
             agg_types = (rd.count, rd.any)
@@ -379,7 +380,7 @@ class aggregate(LineAggregationOperation):
                         x_range=x_range, y_range=y_range)
 
         agg_kwargs = {}
-        if self.p.line_width and glyph == 'line' and ds_version >= Version('0.14.0'):
+        if self.p.line_width and glyph == 'line' and DATASHADER_GE_0_14_0:
             agg_kwargs['line_width'] = self.p.line_width
 
         dfdata = PandasInterface.as_dframe(data)
@@ -388,9 +389,9 @@ class aggregate(LineAggregationOperation):
         if sel_fn:
             if isinstance(params["vdims"], (Dimension, str)):
                 params["vdims"] = [params["vdims"]]
-            sum_agg = ds.summary(**{str(params["vdims"][0]): agg_fn, "index": ds.where(sel_fn)})
+            sum_agg = ds.summary(**{str(params["vdims"][0]): agg_fn, "__index__": ds.where(sel_fn)})
             agg = self._apply_datashader(dfdata, cvs_fn, sum_agg, agg_kwargs, x, y)
-            _ignore = [*params["vdims"], "index"]
+            _ignore = [*params["vdims"], "__index__"]
             sel_vdims = [s for s in agg if s not in _ignore]
             params["vdims"] = [*params["vdims"], *sel_vdims]
         else:
@@ -405,8 +406,7 @@ class aggregate(LineAggregationOperation):
 
         if isinstance(agg, xr.Dataset) or agg.ndim == 2:
             # Replacing x and y coordinates to avoid numerical precision issues
-            eldata = agg if ds_version > Version('0.5.0') else (xs, ys, agg.data)
-            return self.p.element_type(eldata, **params)
+            return self.p.element_type(agg, **params)
         else:
             params['vdims'] = list(map(str, agg.coords[agg_fn.column].data))
             return ImageStack(agg, **params)
@@ -421,14 +421,14 @@ class aggregate(LineAggregationOperation):
             )
             agg = cvs_fn(dfdata, x.name, y.name, agg_fn, **agg_kwargs)
 
-        is_where_index = ds15 and isinstance(agg_fn, ds.where) and isinstance(agg_fn.column, rd.SpecialColumn)
-        is_summary_index = isinstance(agg_fn, ds.summary) and "index" in agg
+        is_where_index = DATASHADER_GE_0_15_1 and isinstance(agg_fn, ds.where) and isinstance(agg_fn.column, rd.SpecialColumn)
+        is_summary_index = isinstance(agg_fn, ds.summary) and "__index__" in agg
         if is_where_index or is_summary_index:
             if is_where_index:
                 data = agg.data
-                agg = agg.to_dataset(name="index")
+                agg = agg.to_dataset(name="__index__")
             else:  # summary index
-                data = agg.index.data
+                data = agg["__index__"].data
             neg1 = data == -1
             for col in dfdata.columns:
                 if col in agg.coords:
@@ -708,7 +708,7 @@ class spikes_aggregate(LineAggregationOperation):
                         x_range=x_range, y_range=y_range)
 
         agg_kwargs = {}
-        if ds_version >= Version('0.14.0'):
+        if DATASHADER_GE_0_14_0:
             agg_kwargs['line_width'] = self.p.line_width
 
         rename_dict = {k: v for k, v in rename_dict.items() if k != v}
@@ -770,14 +770,12 @@ class geom_aggregate(AggregationOperation):
 
         if agg.ndim == 2:
             # Replacing x and y coordinates to avoid numerical precision issues
-            eldata = agg if ds_version > Version('0.5.0') else (xs, ys, agg.data)
-            return self.p.element_type(eldata, **params)
+            return self.p.element_type(agg, **params)
         else:
             layers = {}
             for c in agg.coords[agg_fn.column].data:
                 cagg = agg.sel(**{agg_fn.column: c})
-                eldata = cagg if ds_version > Version('0.5.0') else (xs, ys, cagg.data)
-                layers[c] = self.p.element_type(eldata, **params)
+                layers[c] = self.p.element_type(cagg, **params)
             return NdOverlay(layers, kdims=[element.get_dimension(agg_fn.column)])
 
 
@@ -788,7 +786,7 @@ class segments_aggregate(geom_aggregate, LineAggregationOperation):
 
     def _aggregate(self, cvs, df, x0, y0, x1, y1, agg_fn):
         agg_kwargs = {}
-        if ds_version >= Version('0.14.0'):
+        if DATASHADER_GE_0_14_0:
             agg_kwargs['line_width'] = self.p.line_width
 
         return cvs.line(df, [x0, x1], [y0, y1], agg_fn, axis=1, **agg_kwargs)
@@ -878,9 +876,6 @@ class regrid(AggregationOperation):
 
 
     def _process(self, element, key=None):
-        if ds_version <= Version('0.5.0'):
-            raise RuntimeError('regrid operation requires datashader>=0.6.0')
-
         # Compute coords, anges and size
         x, y = element.kdims
         coords = tuple(element.dimension_values(d, expanded=False) for d in [x, y])
@@ -1040,7 +1035,7 @@ class trimesh_rasterize(aggregate):
         precompute = self.p.precompute
         if interp == 'linear': interp = 'bilinear'
         wireframe = False
-        if (not (element.vdims or (isinstance(element, TriMesh) and element.nodes.vdims))) and ds_version <= Version('0.6.9'):
+        if (not (element.vdims or (isinstance(element, TriMesh) and element.nodes.vdims))) and DATASHADER_VERSION <= (0, 6, 9):
             self.p.aggregator = ds.any() if isinstance(agg, ds.any) or agg == 'any' else ds.count()
             return aggregate._process(self, element, key)
         elif ((not interp and (isinstance(agg, (ds.any, ds.count)) or
@@ -1101,11 +1096,11 @@ class quadmesh_rasterize(trimesh_rasterize):
     """
 
     def _precompute(self, element, agg):
-        if ds_version <= Version('0.7.0'):
+        if DATASHADER_VERSION <= (0, 7, 0):
             return super()._precompute(element.trimesh(), agg)
 
     def _process(self, element, key=None):
-        if ds_version <= Version('0.7.0'):
+        if DATASHADER_VERSION <= (0, 7, 0):
             return super()._process(element, key)
 
         if element.interface.datatype != 'xarray':
@@ -1286,21 +1281,39 @@ class shade(LinkableOperation):
             ydensity = element.ydensity
             bounds = element.bounds
 
+        # Convert to xarray if not already
+        if element.interface.datatype != 'xarray':
+            element = element.clone(datatype=['xarray'])
+
         kdims = element.kdims
         if isinstance(element, ImageStack):
             vdim = element.vdims
             array = element.data
-            if hasattr(array, "to_array"):
-                array = array.to_array("z")
-            array = array.transpose(*[kdim.name for kdim in kdims], ...)
+            # If data is a xarray Dataset it has to be converted to a
+            # DataArray, either by selecting the singular value
+            # dimension or by adding a z-dimension
+            kdims = [kdim.name for kdim in kdims]
+            if not element.interface.packed(element):
+                if len(vdim) == 1:
+                    array = array[vdim[0].name]
+                else:
+                    array = array.to_array("z")
+                    # If data is 3D then we have one extra constant dimension
+                    if array.ndim > 3:
+                        drop = [d for d in array.dims if d not in kdims+["z"]]
+                        array = array.squeeze(dim=drop)
+            array = array.transpose(*kdims, ...)
         else:
             vdim = element.vdims[0].name
             array = element.data[vdim]
 
+        # Dask is not supported by shade so materialize it
+        array = array.compute()
+
         shade_opts = dict(
             how=self.p.cnorm, min_alpha=self.p.min_alpha, alpha=self.p.alpha
         )
-        if ds_version >= Version('0.14.0'):
+        if DATASHADER_GE_0_14_0:
             shade_opts['rescale_discrete_levels'] = self.p.rescale_discrete_levels
 
         # Compute shading options depending on whether
@@ -1334,8 +1347,8 @@ class shade(LinkableOperation):
 
         if self.p.clims:
             shade_opts['span'] = self.p.clims
-        elif ds_version > Version('0.5.0') and self.p.cnorm != 'eq_hist':
-            shade_opts['span'] = element.range(vdim)
+        elif self.p.cnorm != 'eq_hist':
+            shade_opts['span'] = (array.min().item(), array.max().item())
 
         params = dict(get_param_values(element), kdims=kdims,
                       bounds=bounds, vdims=RGB.vdims[:],
@@ -1391,7 +1404,7 @@ class geometry_rasterize(LineAggregationOperation):
             data, col = self._precomputed[element._plot_id]
         else:
             if (('spatialpandas' not in element.interface.datatype) and
-                (not (ds16 and 'geodataframe' in element.interface.datatype))):
+                (not (DATASHADER_GE_0_16_0 and 'geodataframe' in element.interface.datatype))):
                 element = element.clone(datatype=['spatialpandas'])
             data = element.data
             col = element.interface.geo_column(data)
@@ -1406,7 +1419,7 @@ class geometry_rasterize(LineAggregationOperation):
         if isinstance(element, Polygons):
             agg = cvs.polygons(data, **agg_kwargs)
         elif isinstance(element, Path):
-            if self.p.line_width and ds_version >= Version('0.14.0'):
+            if self.p.line_width and DATASHADER_GE_0_14_0:
                 agg_kwargs['line_width'] = self.p.line_width
             agg = cvs.line(data, **agg_kwargs)
         elif isinstance(element, Points):
@@ -1609,7 +1622,7 @@ class SpreadingOperation(LinkableOperation):
     to make sparse plots more visible.
     """
 
-    how = param.ObjectSelector(default='source' if ds_version <= Version('0.11.1') else None,
+    how = param.ObjectSelector(default='source' if DATASHADER_VERSION <= (0, 11, 1) else None,
             objects=[None, 'source', 'over', 'saturate', 'add', 'max', 'min'], doc="""
         The name of the compositing operator to use when combining
         pixels. Default of None uses 'over' operator for RGB elements
@@ -1784,9 +1797,10 @@ class inspect_mask(Operation):
     operation.
     """
 
-    pixels = param.Integer(default=3, doc="""
+    pixels = param.ClassSelector(default=3, class_=(int, tuple), doc="""
        Size of the mask that should match the pixels parameter used in
-       the associated inspection operation.""")
+       the associated inspection operation. Pixels can be provided as
+       integer or x/y-tuple to perform asymmetric masking.""")
 
     streams = param.ClassSelector(default=[PointerXY], class_=(dict, list))
     x = param.Number(default=0)
@@ -1795,9 +1809,13 @@ class inspect_mask(Operation):
     @classmethod
     def _distance_args(cls, element, x_range, y_range,  pixels):
         ycount, xcount =  element.interface.shape(element, gridded=True)
+        if isinstance(pixels, tuple):
+            xpixels, ypixels = pixels
+        else:
+            xpixels = ypixels = pixels
         x_delta = abs(x_range[1] - x_range[0]) / xcount
         y_delta = abs(y_range[1] - y_range[0]) / ycount
-        return (x_delta*pixels, y_delta*pixels)
+        return (x_delta*xpixels, y_delta*ypixels)
 
     def _process(self, raster, key=None):
         if isinstance(raster, RGB):
@@ -1820,10 +1838,11 @@ class inspect(Operation):
     type.
     """
 
-    pixels = param.Integer(default=3, doc="""
+    pixels = param.ClassSelector(default=3, class_=(int, tuple), doc="""
        Number of pixels in data space around the cursor point to search
        for hits in. The hit within this box mask that is closest to the
-       cursor's position is displayed.""")
+       cursor's position is displayed. Pixels can be provided as
+       integer or x/y-tuple to perform asymmetric masking.""")
 
     null_value = param.Number(default=0, doc="""
        Value of raster which indicates no hits. For instance zero for
@@ -1856,9 +1875,9 @@ class inspect(Operation):
                                                y=PointerXY.param.y),
                                   class_=(dict, list))
 
-    x = param.Number(default=0, doc="x-position to inspect.")
+    x = param.Number(default=None, doc="x-position to inspect.")
 
-    y = param.Number(default=0, doc="y-position to inspect.")
+    y = param.Number(default=None, doc="y-position to inspect.")
 
     _dispatch = {}
 
@@ -1909,33 +1928,40 @@ class inspect_base(inspect):
 
     def _process(self, raster, key=None):
         self._validate(raster)
-        if isinstance(raster, RGB):
-            raster = raster[..., raster.vdims[-1]]
-        x_range, y_range = raster.range(0), raster.range(1)
-        xdelta, ydelta = self._distance_args(raster, x_range, y_range, self.p.pixels)
         x, y = self.p.x, self.p.y
-        val = raster[x-xdelta:x+xdelta, y-ydelta:y+ydelta].reduce(function=np.nansum)
-        if np.isnan(val):
-            val = self.p.null_value
-
-        if ((self.p.value_bounds and
-             not (self.p.value_bounds[0] < val < self.p.value_bounds[1]))
-             or val == self.p.null_value):
-            result = self._empty_df(raster.dataset)
-        else:
-            masked = self._mask_dataframe(raster, x, y, xdelta, ydelta)
-            result = self._sort_by_distance(raster, masked, x, y)
+        if x is not None and y is not None:
+            if isinstance(raster, RGB):
+                raster = raster[..., raster.vdims[-1]]
+            x_range, y_range = raster.range(0), raster.range(1)
+            xdelta, ydelta = self._distance_args(raster, x_range, y_range, self.p.pixels)
+            val = raster[x-xdelta:x+xdelta, y-ydelta:y+ydelta].reduce(function=np.nansum)
+            if np.isnan(val):
+                val = self.p.null_value
+
+            if ((self.p.value_bounds and
+                 not (self.p.value_bounds[0] < val < self.p.value_bounds[1]))
+                or val == self.p.null_value):
+                result = self._empty_df(raster.dataset)
+            else:
+                masked = self._mask_dataframe(raster, x, y, xdelta, ydelta)
+                result = self._sort_by_distance(raster, masked, x, y)
 
-        self.hits = result
+            self.hits = result
+        else:
+            self.hits = result = self._empty_df(raster.dataset)
         df = self.p.transform(result)
         return self._element(raster, df.iloc[:self.p.max_indicators])
 
     @classmethod
-    def _distance_args(cls, element, x_range, y_range,  pixels):
-        ycount, xcount =  element.interface.shape(element, gridded=True)
+    def _distance_args(cls, element, x_range, y_range, pixels):
+        ycount, xcount = element.interface.shape(element, gridded=True)
+        if isinstance(pixels, tuple):
+            xpixels, ypixels = pixels
+        else:
+            xpixels = ypixels = pixels
         x_delta = abs(x_range[1] - x_range[0]) / xcount
         y_delta = abs(y_range[1] - y_range[0]) / ycount
-        return (x_delta*pixels, y_delta*pixels)
+        return (x_delta*xpixels, y_delta*ypixels)
 
     @classmethod
     def _empty_df(cls, dataset):
@@ -1990,9 +2016,21 @@ class inspect_points(inspect_base):
         ds = raster.dataset.clone(df)
         xs, ys = (ds.dimension_values(kd) for kd in raster.kdims)
         dx, dy = xs - x, ys - y
-        distances = pd.Series(dx*dx + dy*dy)
-        return df.iloc[distances.argsort().values]
-
+        xtype, ytype = dx.dtype.kind, dy.dtype.kind
+        if xtype in 'Mm':
+            dx = dx.astype('int64')
+        if ytype in 'Mm':
+            dy = dx.astype('int64')
+        # If coordinate types don't match normalize
+        # coordinate space to ensure that distance
+        # in both direction is handled the same.
+        if xtype != ytype and len(dx) and len(dy):
+            with warnings.catch_warnings():
+                warnings.filterwarnings('ignore', r'invalid value encountered in (divide)')
+                dx = (dx - dx.min()) / (dx.max() - dx.min())
+                dy = (dy - dy.min()) / (dy.max() - dy.min())
+        distances = pd.Series(dx**2 + dy**2)
+        return df.iloc[distances.fillna(0).argsort().values]
 
 
 class inspect_polygons(inspect_base):
@@ -2026,7 +2064,7 @@ class inspect_polygons(inspect_base):
                 xs.append((np.min(gxs)+np.max(gxs))/2)
                 ys.append((np.min(gys)+np.max(gys))/2)
         dx, dy = np.array(xs) - x, np.array(ys) - y
-        distances = pd.Series(dx*dx + dy*dy)
+        distances = pd.Series(dx**2 + dy**2)
         return df.iloc[distances.argsort().values]
 
 
diff --git a/holoviews/operation/element.py b/holoviews/operation/element.py
index 131a14a08..7d302f6dd 100644
--- a/holoviews/operation/element.py
+++ b/holoviews/operation/element.py
@@ -830,7 +830,7 @@ class histogram(Operation):
 
         # Mask data
         if is_ibis_expr(data):
-            from ..core.data.ibis import ibis5
+            from ..core.data.ibis import ibis5, ibis9_5
 
             mask = data.notnull()
             if self.p.nonzero:
@@ -840,7 +840,7 @@ class histogram(Operation):
             else:
                 # to_projection removed in ibis 5.0.0
                 data = data.to_projection()
-            data = data[mask]
+            data = data.filter(mask) if ibis9_5() else data[mask]
             no_data = not len(data.head(1).execute())
             data = data[dim.name]
         else:
diff --git a/holoviews/operation/resample.py b/holoviews/operation/resample.py
index c97cefd15..bb482a79c 100644
--- a/holoviews/operation/resample.py
+++ b/holoviews/operation/resample.py
@@ -183,6 +183,11 @@ class ResampleOperation2D(ResampleOperation1D):
             if y and element.get_dimension_type(y[0]) in datetime_types:
                 ytype = 'datetime'
 
+        # Adjust width and height depending on pixel ratio
+        pixel_ratio = self._get_pixel_ratio()
+        width = int(width * pixel_ratio)
+        height = int(height * pixel_ratio)
+
         # Compute highest allowed sampling density
         xspan = xend - xstart
         yspan = yend - ystart
@@ -198,12 +203,11 @@ class ResampleOperation2D(ResampleOperation1D):
             yunit, height = 0, 0
         else:
             yunit = float(yspan)/height
-        xs, ys = (np.linspace(xstart+xunit/2., xend-xunit/2., width),
-                  np.linspace(ystart+yunit/2., yend-yunit/2., height))
 
-        pixel_ratio = self._get_pixel_ratio()
-        width = int(width * pixel_ratio)
-        height = int(height * pixel_ratio)
+        xs, ys = (
+            np.linspace(xstart+xunit/2., xend-xunit/2., width),
+            np.linspace(ystart+yunit/2., yend-yunit/2., height)
+        )
         return ((xstart, xend), (ystart, yend)), (xs, ys), (width, height), (xtype, ytype)
 
     def _get_pixel_ratio(self):
diff --git a/holoviews/plotting/bokeh/annotation.py b/holoviews/plotting/bokeh/annotation.py
index c97c44c95..64f9fbdc4 100644
--- a/holoviews/plotting/bokeh/annotation.py
+++ b/holoviews/plotting/bokeh/annotation.py
@@ -45,6 +45,11 @@ class _SyntheticAnnotationPlot(ColorbarPlot):
             raise ImportError(msg)
         super().__init__(element, **kwargs)
 
+    def _get_axis_dims(self, element):
+        if isinstance(element, (HLines, HSpans)):
+            return None, element.kdims[0], None
+        return element.kdims[0], None, None
+
     def _init_glyph(self, plot, mapping, properties):
         self._plot_methods = {"single": self._methods[self.invert_axes]}
         return super()._init_glyph(plot, mapping, properties)
diff --git a/holoviews/plotting/bokeh/callbacks.py b/holoviews/plotting/bokeh/callbacks.py
index 5e6011aaf..0386a20c1 100644
--- a/holoviews/plotting/bokeh/callbacks.py
+++ b/holoviews/plotting/bokeh/callbacks.py
@@ -50,6 +50,7 @@ from ...streams import (
     Lasso,
     MouseEnter,
     MouseLeave,
+    MultiAxisTap,
     PanEnd,
     PlotReset,
     PlotSize,
@@ -781,6 +782,54 @@ class TapCallback(PopupMixin, PointerXYCallback):
         return value
 
 
+class MultiAxisTapCallback(TapCallback):
+    """
+    Returns the mouse x/y-positions on tap event.
+    """
+
+    attributes = {'x': 'cb_obj.x', 'y': 'cb_obj.y'}
+
+    def _process_msg(self, msg):
+        x_range = self.plot.handles.get('x_range')
+        y_range = self.plot.handles.get('y_range')
+        extra_x = list(self.plot.handles.get('extra_x_ranges', {}).values())
+        extra_y = list(self.plot.handles.get('extra_y_ranges', {}).values())
+        xaxis = self.plot.handles.get('xaxis')
+        yaxis = self.plot.handles.get('yaxis')
+
+        # Compute x/y position relative to first axis
+        x, y = msg['x'], msg['y']
+        x0, x1 = x_range.start, x_range.end
+        y0, y1 = y_range.start, y_range.end
+        if isinstance(xaxis, DatetimeAxis):
+            x = convert_timestamp(x)
+            if isinstance(x0, (float, int)):
+                x0 = convert_timestamp(x0)
+            if isinstance(x1, (float, int)):
+                x1 = convert_timestamp(x1)
+        if isinstance(yaxis, DatetimeAxis):
+            y = convert_timestamp(y)
+            if isinstance(y0, (float, int)):
+                y0 = convert_timestamp(y0)
+            if isinstance(y1, (float, int)):
+                y1 = convert_timestamp(y1)
+        xs, ys = {x_range.name: x}, {y_range.name: y}
+        xspan, yspan = x1 - x0, y1 - y0
+        xfactor, yfactor = (x-x0) / xspan, (y-y0) / yspan
+
+        # Use computed factors to compute x/y position on other axes
+        for values, factor, ranges, axis in (
+            (xs, xfactor, extra_x, xaxis),
+            (ys, yfactor, extra_y, yaxis)
+        ):
+            for rng in ranges:
+                value = rng.start + (rng.end-rng.start) * factor
+                if isinstance(axis, DatetimeAxis) and isinstance(value, (float, int)):
+                    value = convert_timestamp(value)
+                values[rng.name] = value
+
+        return {'xs': xs, 'ys': ys}
+
 
 class SingleTapCallback(TapCallback):
     """
@@ -1142,13 +1191,25 @@ class Selection1DCallback(PopupMixin, Callback):
               }
               if (!xs || !ys) { return }
               for (const i of indices) {
-                const tx = xs[i]
+                let ix = xs[i]
+                let iy = ys[i]
+                let tx, ty
+                if (typeof ix === 'number') {
+                  tx = ix
+                  ty = iy
+                } else {
+                  while (ix.length && (typeof ix[0] !== 'number')) {
+                    ix = ix[0]
+                    iy = iy[0]
+                  }
+                  tx = Math.max(...ix)
+                  ty = Math.max(...iy)
+                }
                 if (!x || (tx > x)) {
-                  x = xs[i]
+                  x = tx
                 }
-                const ty = ys[i]
                 if (!y || (ty > y)) {
-                  y = ys[i]
+                  y = ty
                 }
               }
               if (x && y) {
@@ -1585,5 +1646,6 @@ Stream._callbacks['bokeh'].update({
     FreehandDraw: FreehandDrawCallback,
     PolyDraw    : PolyDrawCallback,
     PolyEdit    : PolyEditCallback,
-    SelectMode  : SelectModeCallback
+    SelectMode  : SelectModeCallback,
+    MultiAxisTap: MultiAxisTapCallback
 })
diff --git a/holoviews/plotting/bokeh/chart.py b/holoviews/plotting/bokeh/chart.py
index 17ef69983..dd76fcc68 100644
--- a/holoviews/plotting/bokeh/chart.py
+++ b/holoviews/plotting/bokeh/chart.py
@@ -880,7 +880,7 @@ class BarPlot(BarsMixin, ColorbarPlot, LegendPlot):
 
     def get_data(self, element, ranges, style):
         # Get x, y, group, stack and color dimensions
-        group_dim, stack_dim = None, None
+        group_dim, stack_dim, stack_order = None, None, None
         if element.ndims == 1:
             grouping = None
         elif self.stacked:
@@ -893,6 +893,8 @@ class BarPlot(BarsMixin, ColorbarPlot, LegendPlot):
             else:
                 stack_order = element.dimension_values(1, False)
             stack_order = list(stack_order)
+            stack_data = element.dimension_values(stack_dim)
+            stack_idx = stack_data == stack_data[0]
         else:
             grouping = 'grouped'
             group_dim = element.get_dimension(1)
@@ -921,7 +923,8 @@ class BarPlot(BarsMixin, ColorbarPlot, LegendPlot):
             grouped = {0: element}
             is_dt = isdatetime(xvals)
             if is_dt or xvals.dtype.kind not in 'OU':
-                xdiff = np.abs(np.diff(xvals))
+                xslice = stack_idx if stack_order else slice(None)
+                xdiff = np.abs(np.diff(xvals[xslice]))
                 diff_size = len(np.unique(xdiff))
                 if diff_size == 0 or (diff_size == 1 and xdiff[0] == 0):
                     xdiff = 1
diff --git a/holoviews/plotting/bokeh/element.py b/holoviews/plotting/bokeh/element.py
index 47d5950c7..b17662554 100644
--- a/holoviews/plotting/bokeh/element.py
+++ b/holoviews/plotting/bokeh/element.py
@@ -21,6 +21,7 @@ from bokeh.models import (
     GlyphRenderer,
     Legend,
     Renderer,
+    Span,
     Title,
     tools,
 )
@@ -53,6 +54,7 @@ from ...core.overlay import CompositeOverlay, NdOverlay
 from ...element import Annotation, Contours, Graph, Path, Tiles, VectorField
 from ...streams import Buffer, PlotSize, RangeXY
 from ...util.transform import dim
+from ...util.warnings import warn
 from ..plot import GenericElementPlot, GenericOverlayPlot
 from ..util import color_intervals, dim_axis_label, dim_range_key, process_cmap
 from .plot import BokehPlot
@@ -72,6 +74,7 @@ from .util import (
     bokeh32,
     bokeh34,
     bokeh35,
+    bokeh36,
     cds_column_replace,
     compute_layout_properties,
     date_to_integer,
@@ -208,7 +211,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         The scalebar_unit is only used if scalebar is True.""")
 
     scalebar_location = param.Selector(
-        default="bottom_right",
+        default=None,
         objects=[
             "top_left", "top_center", "top_right",
             "center_left", "center_center", "center_right",
@@ -218,6 +221,8 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         doc="""
             Location anchor for positioning scale bar.
 
+            Default to 'bottom_right', except if subcoordinate_y is True then it will default to 'right'.
+
             The scalebar_location is only used if scalebar is True.""")
 
     scalebar_label = param.String(
@@ -368,6 +373,8 @@ class ElementPlot(BokehPlot, GenericElementPlot):
 
         # Flag to check whether plot has been updated
         self._updated = False
+        # Counter to keep track of last stream update
+        self._stream_count = None
 
     def _hover_opts(self, element):
         if self.batched:
@@ -575,10 +582,11 @@ class ElementPlot(BokehPlot, GenericElementPlot):
             tool_list.append(tool)
 
         copied_tools = []
+        skip_models = (Span,)
         for tool in tool_list:
             if isinstance(tool, tools.Tool):
                 properties = {
-                    p: v.clone() if isinstance(v, Model) else v
+                    p: v.clone() if isinstance(v, Model) and not isinstance(v, skip_models) else v
                     for p, v in tool.properties_with_values(include_defaults=False).items()
                 }
                 tool = type(tool)(**properties)
@@ -802,20 +810,27 @@ class ElementPlot(BokehPlot, GenericElementPlot):
             if dim_range:
                 self._shared[shared_name] = True
 
+        # If we have a single dimension grab it so it can be set as the Range name
+        name = None
+        if dim:
+            name = dim.name
+        elif dims and len(dims) == 1:
+            name = dims[0].name
+
         if self._shared.get(shared_name) and not dim:
             pass
         elif categorical:
             axis_type = 'auto'
-            dim_range = FactorRange()
+            dim_range = FactorRange(name=name)
         elif None in [v0, v1] or any(
             True if isinstance(el, (str, bytes)+util.cftime_types)
             else not util.isfinite(el) for el in [v0, v1]
         ):
-            dim_range = range_type()
+            dim_range = range_type(name=name)
         elif issubclass(range_type, FactorRange):
-            dim_range = range_type(name=dim.name if dim else None)
+            dim_range = range_type(name=name)
         else:
-            dim_range = range_type(start=v0, end=v1, name=dim.name if dim else None)
+            dim_range = range_type(start=v0, end=v1, name=name)
 
         if not dim_range.tags and specs is not None:
             dim_range.tags.append(specs)
@@ -982,6 +997,13 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         fig.xaxis[0].update(**axis_props['x'])
         fig.yaxis[0].update(**axis_props['y'])
 
+        # Set up handlers to configure following behavior on streaming plots
+        if self.streaming:
+            fig.on_event('rangesupdate', self._disable_follow)
+            fig.on_event('reset', self._reset_follow)
+            code = "export default (_, cb_obj) => { cb_obj.origin.hold_render = false }"
+            fig.js_on_event('reset', CustomJS(code=code))
+
         # Do not add the extra axes to the layout if subcoordinates are used
         if self._subcoord_overlaid:
             return fig
@@ -1001,6 +1023,20 @@ class ElementPlot(BokehPlot, GenericElementPlot):
             fig.add_layout(ax, axis_position)
         return fig
 
+    def _disable_follow(self, event):
+        hold = self.state.hold_render
+        for stream in self.streaming:
+            stream.following = False
+        if not hold:
+            stream.trigger(self.streaming)
+            self.state.hold_render = True
+
+    def _reset_follow(self, event):
+        self.state.hold_render = False
+        for stream in self.streaming:
+            stream.following = True
+        stream.trigger(self.streaming)
+
     def _plot_properties(self, key, element):
         """
         Returns a dictionary of plot properties.
@@ -1312,7 +1348,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
                 self._shared.get(extra_y_range.name, False), log, streaming
             )
 
-    def _update_main_ranges(self, element, x_range, y_range, ranges):
+    def _update_main_ranges(self, element, x_range, y_range, ranges, subcoord=False):
         plot = self.handles['plot']
 
         l, b, r, t = None, None, None, None
@@ -1336,11 +1372,11 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         framewise = self.framewise
         streaming = (self.streaming and any(stream._triggering and stream.following
                                             for stream in self.streaming))
-        xupdate = ((not (self.model_changed(x_range) or self.model_changed(plot))
+        xupdate = not subcoord and ((not (self.model_changed(x_range) or self.model_changed(plot))
                     and (framewise or streaming))
                    or xfactors is not None)
-        yupdate = ((not (self.model_changed(x_range) or self.model_changed(plot))
-                    and (framewise or streaming) or yfactors is not None) and not self.subcoordinate_y)
+        yupdate = (not (self.model_changed(x_range) or self.model_changed(plot))
+                    and (framewise or streaming) or yfactors is not None)
 
         options = self._traverse_options(element, 'plot', ['width', 'height'], defaults=False)
         fixed_width = (self.frame_width or options.get('width'))
@@ -1452,7 +1488,27 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         if not self.drawn or xupdate:
             self._update_range(x_range, l, r, xfactors, self.invert_xaxis,
                                self._shared['x-main-range'], self.logx, streaming)
-        if not (self.drawn or self.subcoordinate_y) or yupdate:
+
+        # If subcoordinate_y is enabled we iterate over each of the
+        # subcoordinate ranges and let the subplot handle the update
+        if self.subcoordinate_y and yupdate and not subcoord:
+            updated = set()
+            for sp in (self.subplots or {}).values():
+                if isinstance(sp, GenericOverlayPlot):
+                    subcoord = False
+                    el_ranges = ranges
+                else:
+                    sp_range = sp.handles.get('y_range')
+                    if not sp_range or sp_range.name in updated:
+                        continue
+                    el_ranges = util.match_spec(sp.current_frame, ranges)
+                    updated.add(sp_range.name)
+                    subcoord = True
+                sp._update_main_ranges(
+                    sp.current_frame, sp.handles['x_range'], sp.handles['y_range'],
+                    el_ranges, subcoord=subcoord
+                )
+        elif (not self.drawn or yupdate) and (not self.subcoordinate_y or subcoord):
             self._update_range(
                 y_range, b, t, yfactors, self._get_tag(y_range, 'invert_yaxis'),
                 self._shared['y-main-range'], self.logy, streaming
@@ -2059,7 +2115,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         self._postprocess_hover(renderer, source)
 
         if self.scalebar:
-            self._draw_scalebar(plot)
+            self._draw_scalebar(plot=plot, renderer=renderer)
 
         zooms_subcoordy = self.handles.get('zooms_subcoordy')
         if zooms_subcoordy is not None:
@@ -2349,7 +2405,7 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         return any(self.lookup_options(frame, 'norm').options.get('framewise')
                    for frame in current_frames)
 
-    def _draw_scalebar(self, plot):
+    def _draw_scalebar(self, *, plot, renderer):
         """Draw scalebar on the plot
 
         This will draw a scalebar on the plot. See the documentation for
@@ -2357,10 +2413,14 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         `scalebar_opts`, and `scalebar_unit` for more information.
 
         Requires Bokeh 3.4
+
+        For scalebar on a subcoordinate_y plot Bokeh 3.6 is needed.
         """
 
         if not bokeh34:
             raise RuntimeError("Scalebar requires Bokeh >= 3.4.0")
+        elif not bokeh36 and self._subcoord_overlaid:
+            warn("Scalebar with subcoordinate_y requires Bokeh >= 3.6.0", RuntimeWarning)
 
         from bokeh.models import Metric, ScaleBar
 
@@ -2371,17 +2431,34 @@ class ElementPlot(BokehPlot, GenericElementPlot):
         else:
             base_unit = unit
 
-        _default_scalebar_opts = {"background_fill_alpha": 0.8}
-        opts = dict(_default_scalebar_opts, **self.scalebar_opts)
+        if self._subcoord_overlaid:
+            srange = renderer.coordinates.y_source
+            orientation = "vertical"
+            # Integer is used for the location as `major_label_overrides` overrides the
+            # label with {0: labelA, 1: labelB}
+            location = (self.scalebar_location or "right", len(plot.renderers) - 1)
+            default_scalebar_opts = {
+                'bar_line_width': 3,
+                'label_location': 'left',
+                'background_fill_color': 'white',
+                'background_fill_alpha': 0.6,
+                'length_sizing': 'adaptive',
+                # Adding location so people can overwrite the default
+                'location': location,
+            }
+        else:
+            srange = plot.x_range if self.scalebar_range == "x" else plot.y_range
+            orientation = "horizontal" if self.scalebar_range == "x" else "vertical"
+            location = self.scalebar_location or "bottom_right"
+            default_scalebar_opts = {"background_fill_alpha": 0.8, "location": location}
 
         scale_bar = ScaleBar(
-            range=plot.x_range if self.scalebar_range == "x" else plot.y_range,
-            orientation="horizontal" if self.scalebar_range == "x" else "vertical",
+            range=srange,
+            orientation=orientation,
             unit=unit,
             dimensional=Metric(base_unit=base_unit),
-            location=self.scalebar_location,
             label=self.scalebar_label,
-            **opts,
+            **dict(default_scalebar_opts, **self.scalebar_opts),
         )
         self.handles['scalebar'] = scale_bar
         plot.add_layout(scale_bar)
diff --git a/holoviews/plotting/bokeh/links.py b/holoviews/plotting/bokeh/links.py
index 5cf61d738..36ac59b92 100644
--- a/holoviews/plotting/bokeh/links.py
+++ b/holoviews/plotting/bokeh/links.py
@@ -13,7 +13,7 @@ from ..links import (
     VertexTableLink,
 )
 from ..plot import GenericElementPlot, GenericOverlayPlot
-from .util import bokeh34
+from .util import bokeh34, bokeh35
 
 
 class LinkCallback:
@@ -182,6 +182,23 @@ class RangeToolLinkCallback(LinkCallback):
                     ax.reset_end = end
 
         tool = RangeTool(**axes)
+
+        if bokeh35:
+            use_handles = getattr(link, 'use_handles', True)
+            start_gesture = getattr(link, 'start_gesture', 'tap')
+            inverted = getattr(link, 'inverted', True)
+
+            tool.overlay.use_handles = use_handles
+            tool.start_gesture = start_gesture
+            tool.overlay.inverted = inverted
+
+            if use_handles:
+                tool.overlay.handles.all.hover_fill_color = "grey"
+                tool.overlay.handles.all.hover_fill_alpha = 0.25
+                tool.overlay.handles.all.hover_line_alpha = 0
+                tool.overlay.handles.all.fill_alpha = 0.1
+                tool.overlay.handles.all.line_alpha = 0.25
+
         source_plot.state.add_tools(tool)
         if toolbars:
             toolbars[0].tools.append(tool)
diff --git a/holoviews/plotting/bokeh/plot.py b/holoviews/plotting/bokeh/plot.py
index 77658f43d..109b7067e 100644
--- a/holoviews/plotting/bokeh/plot.py
+++ b/holoviews/plotting/bokeh/plot.py
@@ -177,10 +177,15 @@ class BokehPlot(DimensionedPlot, CallbackPlot):
         if (self.streaming and self.streaming[0].data is self.current_frame.data
             and self._stream_data and not empty):
             stream = self.streaming[0]
-            if stream._triggering:
-                data = {k: v[-stream._chunk_length:] for k, v in data.items()}
-                source.stream(data, stream.length)
-            return
+            if stream._count == ((self._stream_count or 0) + 1):
+                if stream._triggering and stream.following:
+                    data = {k: v[-stream._chunk_length:] for k, v in data.items()}
+                    source.stream(data, stream.length)
+                    self._stream_count = stream._count
+                return
+            elif not stream.following:
+                return
+            self._stream_count = stream._count
 
         if cds_column_replace(source, data):
             source.data = data
@@ -244,10 +249,12 @@ class BokehPlot(DimensionedPlot, CallbackPlot):
             if plot.subplots:
                 plot.subplots.clear()
 
-            if isinstance(plot, GenericElementPlot):
-                for callback in plot.callbacks:
-                    streams += callback.streams
-                    callback.cleanup()
+            if not isinstance(plot, (GenericElementPlot, GenericOverlayPlot)):
+                continue
+
+            for callback in plot.callbacks:
+                streams += callback.streams
+                callback.cleanup()
 
             for stream in set(streams):
                 stream._subscribers = [
@@ -256,7 +263,6 @@ class BokehPlot(DimensionedPlot, CallbackPlot):
                     get_method_owner(subscriber) not in plots
                 ]
 
-
     def _fontsize(self, key, label='fontsize', common=True):
         """
         Converts integer fontsizes to a string specifying
diff --git a/holoviews/plotting/bokeh/util.py b/holoviews/plotting/bokeh/util.py
index bf69b888a..02d8cba1a 100644
--- a/holoviews/plotting/bokeh/util.py
+++ b/holoviews/plotting/bokeh/util.py
@@ -67,6 +67,7 @@ bokeh32 = bokeh_version >= Version("3.2")
 bokeh33 = bokeh_version >= Version("3.3")
 bokeh34 = bokeh_version >= Version("3.4")
 bokeh35 = bokeh_version >= Version("3.5")
+bokeh36 = bokeh_version >= Version("3.6")
 
 TOOL_TYPES = {
     'pan': tools.PanTool,
diff --git a/holoviews/plotting/links.py b/holoviews/plotting/links.py
index 0d528c4a2..25289f020 100644
--- a/holoviews/plotting/links.py
+++ b/holoviews/plotting/links.py
@@ -115,8 +115,17 @@ class RangeToolLink(Link):
     intervalsy = param.Tuple(default=None, length=2, doc="""
         (min, max) intervals for the y-axis""")
 
-    _requires_target = True
+    use_handles = param.Boolean(default=True, doc="""
+        Whether to display handles. Only available from Bokeh 3.5 onwards.""")
+
+    start_gesture = param.Selector(default='tap', objects=['pan', 'tap', 'none'],
+         doc="Gesture to start a range selection. Only available from Bokeh 3.5 onwards.")
 
+    inverted = param.Boolean(default=True, doc="""
+         Whether to invert the highlighting of the range selection.
+         Only available from Bokeh 3.5 onwards.""")
+
+    _requires_target = True
 
 class DataLink(Link):
     """
diff --git a/holoviews/plotting/mixins.py b/holoviews/plotting/mixins.py
index 524308a1c..8da7c7bf7 100644
--- a/holoviews/plotting/mixins.py
+++ b/holoviews/plotting/mixins.py
@@ -174,7 +174,7 @@ class BarsMixin:
         else:
             y0, y1 = ranges[vdim]['combined']
 
-        x0, x1 = (l, r) if util.isnumeric(l) and len(element.kdims) == 1 else ('', '')
+        x0, x1 = (l, r) if (util.isnumeric(l) or isinstance(l, util.datetime_types)) else ('', '')
         if range_type == 'data':
             return (x0, y0, x1, y1)
 
diff --git a/holoviews/plotting/mpl/chart.py b/holoviews/plotting/mpl/chart.py
index 5da620b85..917b4a631 100644
--- a/holoviews/plotting/mpl/chart.py
+++ b/holoviews/plotting/mpl/chart.py
@@ -936,9 +936,12 @@ class BarPlot(BarsMixin, ColorbarPlot, LegendPlot):
 
         cats = None
         style_dim = None
+        xslice = slice(None)
         if sdim:
             cats = values['stack']
             style_dim = sdim
+            stack_data = element.dimension_values(sdim)
+            xslice = stack_data == stack_data[0]
         elif cdim:
             cats = values['category']
             style_dim = cdim
@@ -954,14 +957,15 @@ class BarPlot(BarsMixin, ColorbarPlot, LegendPlot):
         is_dt = isdatetime(xvals)
         continuous = True
         if is_dt or xvals.dtype.kind not in 'OU' and not (cdim or len(element.kdims) > 1):
-            xdiff_vals = date2num(xvals) if is_dt else xvals
-            xdiff = np.abs(np.diff(xdiff_vals))
-            if len(np.unique(xdiff)) == 1:
-                # if all are same
+            xvals = xvals[xslice]
+            xdiff = np.abs(np.diff(xvals))
+            diff_size = len(np.unique(xdiff))
+            if diff_size == 0 or (diff_size == 1 and xdiff[0] == 0):
                 xdiff = 1
             else:
                 xdiff = np.min(xdiff)
-            width = (1 - self.bar_padding) * xdiff
+            width = (1 - self.bar_padding) * (date2num(xdiff) / 1000 if is_dt else xdiff)
+            data_width = (1 - self.bar_padding) * xdiff
         else:
             xdiff = len(values.get('category', [None]))
             width = (1 - self.bar_padding) / xdiff
@@ -1050,16 +1054,23 @@ class BarPlot(BarsMixin, ColorbarPlot, LegendPlot):
             axis.legend(title=title, **legend_opts)
 
         x_range = ranges[gdim.label]["data"]
-        if continuous and not is_dt:
-            if style.get('align', 'center') == 'center':
-                left_multiplier = 0.5
-                right_multiplier = 0.5
-            else:
-                left_multiplier = 0
-                right_multiplier = 1
+        if style.get('align', 'center') == 'center':
+            left_multiplier = 0.5
+            right_multiplier = 0.5
+        else:
+            left_multiplier = 0
+            right_multiplier = 1
+        if continuous:
+            ranges[gdim.label]["data"] = (
+                x_range[0] - data_width * left_multiplier,
+                x_range[1] + data_width * right_multiplier,
+            )
+        else:
+            locs = [item[0] for item in xticks]
+            xmin, xmax = min(locs), max(locs)
             ranges[gdim.label]["data"] = (
-                x_range[0] - width * left_multiplier,
-                x_range[1] + width * right_multiplier
+                - width * left_multiplier + xmin,
+                width * right_multiplier + xmax
             )
         return bars, xticks if not continuous else None, ax_dims
 
diff --git a/holoviews/plotting/mpl/plot.py b/holoviews/plotting/mpl/plot.py
index f166544d8..d8208709c 100644
--- a/holoviews/plotting/mpl/plot.py
+++ b/holoviews/plotting/mpl/plot.py
@@ -120,12 +120,18 @@ class MPLPlot(DimensionedPlot):
         Allows labeling the subaxes in each plot with various formatters
         including {Alpha}, {alpha}, {numeric} and {roman}.""")
 
+    sublabel_offset = param.Number(default=0, bounds=(0, None), doc="""
+        Allows offsetting the sublabel index.""")
+
     sublabel_position = param.NumericTuple(default=(-0.35, 0.85), doc="""
          Position relative to the plot for placing the optional subfigure label.""")
 
     sublabel_size = param.Number(default=18, doc="""
          Size of optional subfigure label.""")
 
+    sublabel_skip = param.List(default=None, item_type=int, doc="""
+        List of elements to skip when labeling subplots. Numbering starts at 1.""")
+
     projection = param.Parameter(default=None, doc="""
         The projection of the plot axis, default of None is equivalent to
         2D plot, '3d' and 'polar' are also supported by matplotlib by default.
@@ -207,20 +213,29 @@ class MPLPlot(DimensionedPlot):
 
 
     def _subplot_label(self, axis):
+        if self.sublabel_skip and self.layout_num in self.sublabel_skip:
+            return
         layout_num = self.layout_num if self.subplot else 1
-        if self.sublabel_format and not self.adjoined and layout_num > 0:
+        if self.sublabel_skip:
+            if any(n < 1 for n in self.sublabel_skip):
+                raise ValueError('sublabel_skip values must be greater than 0')
+            sublabel_num = len(set(range(layout_num)) - set(self.sublabel_skip))
+        else:
+            sublabel_num = layout_num
+        sublabel_num += self.sublabel_offset
+        if self.sublabel_format and not self.adjoined and sublabel_num > 0:
             from matplotlib.offsetbox import AnchoredText
             labels = {}
             if '{Alpha}' in self.sublabel_format:
-                labels['Alpha'] = int_to_alpha(layout_num-1)
+                labels['Alpha'] = int_to_alpha(sublabel_num-1)
             elif '{alpha}' in self.sublabel_format:
-                labels['alpha'] = int_to_alpha(layout_num-1, upper=False)
+                labels['alpha'] = int_to_alpha(sublabel_num-1, upper=False)
             elif '{numeric}' in self.sublabel_format:
-                labels['numeric'] = self.layout_num
+                labels['numeric'] = sublabel_num
             elif '{Roman}' in self.sublabel_format:
-                labels['Roman'] = int_to_roman(layout_num)
+                labels['Roman'] = int_to_roman(sublabel_num)
             elif '{roman}' in self.sublabel_format:
-                labels['roman'] = int_to_roman(layout_num).lower()
+                labels['roman'] = int_to_roman(sublabel_num).lower()
             at = AnchoredText(self.sublabel_format.format(**labels), loc=3,
                               bbox_to_anchor=self.sublabel_position, frameon=False,
                               prop=dict(size=self.sublabel_size, weight='bold'),
diff --git a/holoviews/plotting/mpl/tabular.py b/holoviews/plotting/mpl/tabular.py
index 122d57156..e779c5c6e 100644
--- a/holoviews/plotting/mpl/tabular.py
+++ b/holoviews/plotting/mpl/tabular.py
@@ -77,8 +77,7 @@ class TablePlot(ElementPlot):
                     cell_text = element.pprint_cell(adjusted_row, col)
                     if len(cell_text) > self.max_value_len:
                         cell_text = cell_text[:(self.max_value_len-3)]+'...'
-                if len(cell_text) + 2 > cell_widths[col]:
-                    cell_widths[col] = len(cell_text) + 2
+                cell_widths[col] = max(len(cell_text) + 2, cell_widths[col])
 
     def _cell_value(self, element, row, col):
         summarize = element.rows > self.max_rows
diff --git a/holoviews/plotting/plot.py b/holoviews/plotting/plot.py
index eb963afea..444dd6812 100644
--- a/holoviews/plotting/plot.py
+++ b/holoviews/plotting/plot.py
@@ -195,7 +195,7 @@ class Plot(param.Parameterized):
         """
         plots = self.traverse(lambda x: x, [Plot])
         for plot in plots:
-            if not isinstance(plot, (GenericCompositePlot, GenericElementPlot, GenericOverlayPlot)):
+            if not isinstance(plot, (GenericElementPlot, GenericOverlayPlot)):
                 continue
             for stream in set(plot.streams):
                 stream._subscribers = [
diff --git a/holoviews/streams.py b/holoviews/streams.py
index c8c37656d..e129c8b43 100644
--- a/holoviews/streams.py
+++ b/holoviews/streams.py
@@ -565,8 +565,6 @@ class Buffer(Pipe):
             data.stream.sink(self.send)
             self.sdf = data
 
-        if index and isinstance(example, pd.DataFrame):
-            example = example.reset_index()
         params['data'] = example
         super().__init__(**params)
         self.length = length
@@ -656,9 +654,6 @@ class Buffer(Pipe):
         """
         data = kwargs.get('data')
         if data is not None:
-            if (isinstance(data, pd.DataFrame) and
-                list(data.columns) != list(self.data.columns) and self._index):
-                data = data.reset_index()
             self.verify(data)
             kwargs['data'] = self._concat(data)
             self._count += 1
@@ -667,7 +662,7 @@ class Buffer(Pipe):
 
     @property
     def hashkey(self):
-        return {'hash': self._count}
+        return {'hash': (self._count, self._memoize_counter)}
 
 
 
@@ -1358,6 +1353,18 @@ class Tap(PointerXY):
            Pointer position along the y-axis in data coordinates""")
 
 
+class MultiAxisTap(LinkedStream):
+    """
+    The x/y-positions of a tap or click in data coordinates.
+    """
+
+    xs = param.Dict(default=None, constant=True, doc="""
+           Pointer positions along the x-axes in data coordinates""")
+
+    ys = param.Dict(default=None, constant=True, doc="""
+           Pointer positions along the y-axes in data coordinates""")
+
+
 class DoubleTap(PointerXY):
     """
     The x/y-position of a double-tap or -click in data coordinates.
diff --git a/pixi.toml b/pixi.toml
index 8c9d9ad02..c10edd0c5 100644
--- a/pixi.toml
+++ b/pixi.toml
@@ -1,6 +1,6 @@
 [project]
 name = "holoviews"
-channels = ["pyviz/label/dev", "bokeh", "conda-forge"]
+channels = ["pyviz/label/dev", "bokeh/label/rc", "conda-forge"]
 platforms = ["linux-64", "osx-arm64", "osx-64", "win-64"]
 
 [tasks]
@@ -18,28 +18,28 @@ test-310 = ["py310", "test-core", "test", "example", "test-example", "test-unit-
 test-311 = ["py311", "test-core", "test", "example", "test-example", "test-unit-task"]
 test-312 = ["py312", "test-core", "test", "example", "test-example", "test-unit-task"]
 test-ui = ["py312", "test-core", "test", "test-ui"]
-test-core = ["py312", "test-core", "test-unit-task"]
-test-gpu = ["py311", "test-core", "test", "test-gpu"]
+test-core = ["py313", "test-core", "test-unit-task"]
+test-gpu = ["py312", "test-core", "test", "test-gpu"]
 docs = ["py311", "example", "doc"]
 build = ["py311", "build"]
 lint = ["py311", "lint"]
 
 [dependencies]
+nomkl = "*"
+pip = "*"
+# Required
 bokeh = ">=3.1"
 colorcet = "*"
-numpy = ">=1.0"
+numpy = ">=1.21"
 packaging = "*"
-pandas = ">=0.20.0"
+pandas = ">=1.3"
 panel = ">=1.0"
 param = ">=2.0,<3.0"
-pip = "*"
 pyviz_comms = ">=2.1"
-# Recommended
-matplotlib-base = ">=3"
-plotly = ">=4.0"
 
 [feature.py39.dependencies]
 python = "3.9.*"
+panel = "1.4.*"
 
 [feature.py310.dependencies]
 python = "3.10.*"
@@ -53,16 +53,22 @@ bokeh_sampledata = "*"
 python = "3.12.*"
 bokeh_sampledata = "*"
 
+[feature.py313.dependencies]
+python = "3.13.*"
+bokeh_sampledata = "*"
+
 [feature.example.dependencies]
 cftime = "*"
 dask-core = "*"
 datashader = ">=0.11.1"
 ffmpeg = "*"
 ipython = ">=5.4.0"
+matplotlib-base = ">=3"
 netcdf4 = "*"
 networkx = "*"
 notebook = "*"
 pillow = "*"
+plotly = ">=4.0"
 pooch = "*"
 pyarrow = "*"
 scikit-image = "*"
@@ -94,7 +100,9 @@ datashader = ">=0.11.1"
 ffmpeg = "*"
 ibis-sqlite = "*"
 nbconvert = "*"
+matplotlib-base = ">=3"
 networkx = "*"
+plotly = ">=4.0"
 pillow = "*"
 scipy = ">=1.10" # Python 3.9 + Windows downloads 1.9
 selenium = "*"
@@ -113,7 +121,7 @@ test-example = { cmd = 'pytest -n logical --dist loadscope --nbval-lax examples'
 nbval = "*"
 
 [feature.test-ui]
-channels = ["pyviz/label/dev", "microsoft", "conda-forge"]
+channels = ["microsoft"]
 
 [feature.test-ui.dependencies]
 playwright = { version = "*", channel = "microsoft" }
@@ -127,7 +135,7 @@ cmd = 'pytest holoviews/tests/ui --ui --browser chromium'
 depends_on = ["_install-ui"]
 
 [feature.test-gpu]
-channels = ["pyviz/label/dev", "rapidsai", "conda-forge"]
+channels = ["rapidsai"]
 platforms = ["linux-64"]
 
 [feature.test-gpu.activation.env]
@@ -135,7 +143,7 @@ NUMBA_CUDA_LOW_OCCUPANCY_WARNINGS = "0"
 
 [feature.test-gpu.dependencies]
 cuda-version = "12.2.*"
-cudf = "24.06.*"
+cudf = "24.10.*"
 cupy = "*"
 cuspatial = "*"
 librmm = { version = "*", channel = "rapidsai" }
diff --git a/pyproject.toml b/pyproject.toml
index a95d6e4d5..914351eaa 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -21,6 +21,7 @@ classifiers = [
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Operating System :: OS Independent",
     "Intended Audience :: Science/Research",
     "Intended Audience :: Developers",
@@ -109,7 +110,6 @@ filterwarnings = [
     "ignore:datetime.datetime.utcfromtimestamp():DeprecationWarning:bokeh", # https://github.com/bokeh/bokeh/issues/13125
     "ignore:datetime.datetime.utcnow():DeprecationWarning:bokeh", # https://github.com/bokeh/bokeh/issues/13125
     # 2024-01: Pandas 2.2
-    "ignore:When grouping with a length-1 list::dask.dataframe.groupby", # https://github.com/dask/dask/issues/10572
     "ignore:\\s*Pyarrow will become a required dependency of pandas:DeprecationWarning", # Will go away by itself in Pandas 3.0
     "ignore:Passing a (SingleBlockManager|BlockManager) to (Series|GeoSeries|DataFrame|GeoDataFrame) is deprecated:DeprecationWarning", # https://github.com/holoviz/spatialpandas/issues/137
     # 2024-04
@@ -121,6 +121,7 @@ filterwarnings = [
     "ignore:unclosed file <_io.TextIOWrapper name='/dev/null' mode='w':ResourceWarning", # OK
     # 2024-07
     "ignore:The (non_)?interactive_bk attribute was deprecated in Matplotlib 3.9", # OK - Only happening in debug mode
+    "ignore: Jitify is performing a one-time only warm::cupy", # OK
 ]
 
 [tool.coverage]


Test patch:
diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml
index 2f5e4f208..9faaa2ac4 100644
--- a/.github/workflows/test.yaml
+++ b/.github/workflows/test.yaml
@@ -30,7 +30,7 @@ concurrency:
 
 defaults:
   run:
-    shell: bash -el {0}
+    shell: bash -e {0}
 
 env:
   DISPLAY: ":99.0"
@@ -107,7 +107,7 @@ jobs:
     name: Pixi lock
     runs-on: ubuntu-latest
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_lock@pixi
+      - uses: holoviz-dev/holoviz_tasks/pixi_lock@v0
         with:
           cache: ${{ github.event.inputs.cache == 'true' || github.event.inputs.cache == '' }}
 
@@ -115,29 +115,25 @@ jobs:
     name: unit:${{ matrix.environment }}:${{ matrix.os }}
     needs: [pre_commit, setup, pixi_lock]
     runs-on: ${{ matrix.os }}
+    if: needs.setup.outputs.code_change == 'true'
     strategy:
       fail-fast: false
       matrix: ${{ fromJson(needs.setup.outputs.matrix) }}
     timeout-minutes: 120
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_install@pixi
-        if: needs.setup.outputs.code_change == 'true'
+      - uses: holoviz-dev/holoviz_tasks/pixi_install@v0
         with:
           environments: ${{ matrix.environment }}
       - name: Check packages latest version
-        if: needs.setup.outputs.code_change == 'true'
         run: |
           pixi run -e ${{ matrix.environment }} check-latest-packages bokeh panel param datashader
       - name: Test Unit
-        if: needs.setup.outputs.code_change == 'true'
         run: |
           pixi run -e ${{ matrix.environment }} test-unit $COV
       - name: Test Examples
-        if: needs.setup.outputs.code_change == 'true'
         run: |
           pixi run -e ${{ matrix.environment }} test-example
       - uses: codecov/codecov-action@v4
-        if: needs.setup.outputs.code_change == 'true'
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
 
@@ -145,6 +141,7 @@ jobs:
     name: ui:${{ matrix.environment }}:${{ matrix.os }}
     needs: [pre_commit, setup, pixi_lock]
     runs-on: ${{ matrix.os }}
+    if: needs.setup.outputs.code_change == 'true'
     strategy:
       fail-fast: false
       matrix:
@@ -154,12 +151,10 @@ jobs:
     env:
       PANEL_LOG_LEVEL: info
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_install@pixi
-        if: needs.setup.outputs.code_change == 'true'
+      - uses: holoviz-dev/holoviz_tasks/pixi_install@v0
         with:
           environments: ${{ matrix.environment }}
       - name: Test UI
-        if: needs.setup.outputs.code_change == 'true'
         run: |
           # Create a .uicoveragerc file to set the concurrency library to greenlet
           # https://github.com/microsoft/playwright-python/issues/313
@@ -173,7 +168,6 @@ jobs:
           path: ./ui_screenshots
           if-no-files-found: ignore
       - uses: codecov/codecov-action@v4
-        if: needs.setup.outputs.code_change == 'true'
         with:
           token: ${{ secrets.CODECOV_TOKEN }}
 
@@ -181,6 +175,7 @@ jobs:
     name: core:${{ matrix.environment }}:${{ matrix.os }}
     needs: [pre_commit, setup, pixi_lock]
     runs-on: ${{ matrix.os }}
+    if: needs.setup.outputs.code_change == 'true'
     strategy:
       fail-fast: false
       matrix:
@@ -188,15 +183,22 @@ jobs:
         environment: ["test-core"]
     timeout-minutes: 120
     steps:
-      - uses: holoviz-dev/holoviz_tasks/pixi_install@pixi
-        if: needs.setup.outputs.code_change == 'true'
+      - uses: holoviz-dev/holoviz_tasks/pixi_install@v0
         with:
           environments: ${{ matrix.environment }}
       - name: Check packages latest version
-        if: needs.setup.outputs.code_change == 'true'
         run: |
           pixi run -e ${{ matrix.environment }} check-latest-packages numpy pandas bokeh panel param
       - name: Test Unit
-        if: needs.setup.outputs.code_change == 'true'
         run: |
           pixi run -e ${{ matrix.environment }} test-unit
+
+  result_test_suite:
+    name: result:test
+    needs: [unit_test_suite, ui_test_suite, core_test_suite]
+    if: always()
+    runs-on: ubuntu-latest
+    steps:
+      - name: check for failures
+        if: contains(needs.*.result, 'failure') || contains(needs.*.result, 'cancelled')
+        run: echo job failed && exit 1
diff --git a/examples/conftest.py b/examples/conftest.py
index a6ba9f708..b8af528b5 100644
--- a/examples/conftest.py
+++ b/examples/conftest.py
@@ -19,6 +19,8 @@ PD2 = Version(pd.__version__) >= Version("2.0")
 collect_ignore_glob = [
     # Needs selenium, phantomjs, firefox, and geckodriver to save a png picture
     "user_guide/Plotting_with_Bokeh.ipynb",
+    # Streaming data use streamz which is no longer maintained
+    "user_guide/16-Streaming_Data.ipynb",
     # Possible timeout error
     "user_guide/17-Dashboards.ipynb",
     # Give file not found
@@ -36,12 +38,6 @@ if PD2 and system == "Windows":
     ]
 
 
-# 2023-10-25, flaky on CI with timeout
-if system == "Darwin":
-    collect_ignore_glob += [
-        "user_guide/16-Streaming_Data.ipynb",
-    ]
-
 # 2024-01-15: See https://github.com/holoviz/holoviews/issues/6069
 if system == "Windows":
     collect_ignore_glob += [
diff --git a/holoviews/tests/element/test_selection.py b/holoviews/tests/element/test_selection.py
index 6c2e89101..05f8ff341 100644
--- a/holoviews/tests/element/test_selection.py
+++ b/holoviews/tests/element/test_selection.py
@@ -30,6 +30,7 @@ from holoviews.element import (
 )
 from holoviews.element.comparison import ComparisonTestCase
 from holoviews.element.selection import spatial_select_columnar
+from holoviews.util.transform import dim
 
 from ..utils import dask_switcher
 
@@ -61,6 +62,29 @@ ds_available = skipIf(ds is None, 'datashader not available')
 dd_available = pytest.mark.skipif(dd is None, reason='dask.dataframe not available')
 
 
+class TestIndexExpr(ComparisonTestCase):
+
+    def setUp(self):
+        import holoviews.plotting.bokeh # noqa
+        super().setUp()
+        self._backend = Store.current_backend
+        Store.set_current_backend('bokeh')
+
+    def tearDown(self):
+        Store.current_backend = self._backend
+
+    def test_index_selection_on_id_column(self):
+        # tests issue in https://github.com/holoviz/holoviews/pull/6336
+        x, y = np.random.randn(2, 100)
+        idx = np.arange(100)
+
+        points = Points(
+            {'x': x, 'y': y, 'id': idx}, kdims=['x', 'y'], vdims=['id'], datatype=['dataframe']
+        )
+        sel, _, _ = points._get_index_selection([3, 7], ['id'])
+        assert sel == dim('id').isin([3, 7])
+
+
 class TestSelection1DExpr(ComparisonTestCase):
 
     def setUp(self):
diff --git a/holoviews/tests/ipython/test_notebooks.py b/holoviews/tests/ipython/test_notebooks.py
index d23e5d0b6..97c27d9dd 100644
--- a/holoviews/tests/ipython/test_notebooks.py
+++ b/holoviews/tests/ipython/test_notebooks.py
@@ -6,6 +6,7 @@ import os
 import pytest
 
 pytest.importorskip("nbconvert")
+pytest.importorskip("IPython")
 
 import nbconvert
 import nbformat
diff --git a/holoviews/tests/operation/test_datashader.py b/holoviews/tests/operation/test_datashader.py
index 5cf1dceae..823b75f34 100644
--- a/holoviews/tests/operation/test_datashader.py
+++ b/holoviews/tests/operation/test_datashader.py
@@ -6,7 +6,6 @@ import numpy as np
 import pandas as pd
 import pytest
 from numpy import nan
-from packaging.version import Version
 
 from holoviews import (
     RGB,
@@ -38,29 +37,30 @@ from holoviews.streams import Tap
 from holoviews.util import render
 
 try:
-    import dask.dataframe as dd
     import datashader as ds
-    import xarray as xr
-
-    from holoviews.operation.datashader import (
-        AggregationOperation,
-        aggregate,
-        datashade,
-        directly_connect_edges,
-        ds_version,
-        dynspread,
-        inspect,
-        inspect_points,
-        inspect_polygons,
-        rasterize,
-        regrid,
-        shade,
-        spread,
-        stack,
-    )
 except ImportError:
     raise SkipTest('Datashader not available')
 
+import dask.dataframe as dd
+import xarray as xr
+
+from holoviews.operation.datashader import (
+    DATASHADER_VERSION,
+    AggregationOperation,
+    aggregate,
+    datashade,
+    directly_connect_edges,
+    dynspread,
+    inspect,
+    inspect_points,
+    inspect_polygons,
+    rasterize,
+    regrid,
+    shade,
+    spread,
+    stack,
+)
+
 try:
     import spatialpandas
 except ImportError:
@@ -121,6 +121,25 @@ class DatashaderAggregateTests(ComparisonTestCase):
                          vdims=[Dimension('Count', nodata=0)])
         self.assertEqual(img, expected)
 
+    def test_aggregate_points_empty(self):
+        points = Points([])
+        img = aggregate(points, dynamic=False,  x_range=(0, 1), y_range=(0, 1),
+                        width=2, height=2, pixel_ratio=1)
+        expected = Image(([0.25, 0.75], [0.25, 0.75], [[0, 0], [0, 0]]),
+                         vdims=[Dimension('Count', nodata=0)])
+        self.assertEqual(img, expected)
+
+    def test_aggregate_points_empty_with_pixel_ratio(self):
+        points = Points([])
+        img = aggregate(points, dynamic=False,  x_range=(0, 1), y_range=(0, 1),
+                        width=2, height=2, pixel_ratio=2)
+        expected = Image((
+            [0.125, 0.375, 0.625, 0.875],
+            [0.125, 0.375, 0.625, 0.875],
+            [[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]
+        ), vdims=[Dimension('Count', nodata=0)])
+        self.assertEqual(img, expected)
+
     def test_aggregate_points_count_column(self):
         points = Points([(0.2, 0.3, np.nan), (0.4, 0.7, 22), (0, 0.99,np.nan)], vdims='z')
         img = aggregate(points, dynamic=False,  x_range=(0, 1), y_range=(0, 1),
@@ -788,7 +807,7 @@ class DatashaderAggregateTests(ComparisonTestCase):
 class DatashaderCatAggregateTests(ComparisonTestCase):
 
     def setUp(self):
-        if ds_version < Version('0.11.0'):
+        if DATASHADER_VERSION < (0, 11, 0):
             raise SkipTest('Regridding operations require datashader>=0.11.0')
 
     def test_aggregate_points_categorical(self):
@@ -808,6 +827,22 @@ class DatashaderCatAggregateTests(ComparisonTestCase):
         actual = img.data
         assert (expected.data.to_array("z").values == actual.T.values).all()
 
+    def test_aggregate_points_categorical_one_category(self):
+        points = Points([(0.2, 0.3, 'A'), (0.4, 0.7, 'A'), (0, 0.99, 'A')], vdims='z')
+        img = aggregate(points, dynamic=False,  x_range=(0, 1), y_range=(0, 1),
+                        width=2, height=2, aggregator=ds.by('z', ds.count()))
+        x = np.array([0.25, 0.75])
+        y = np.array([0.25, 0.75])
+        a = np.array([[1, 2], [0, 0]])
+        xrds = xr.DataArray(
+            a,
+            dims=('x', 'y'),
+            coords={"x": x, "y": y}
+        )
+        expected = ImageStack(xrds, kdims=["x", "y"], vdims=["a"])
+        actual = img.data
+        assert (expected.data.to_array("z").values == actual.T.values).all()
+
     def test_aggregate_points_categorical_mean(self):
         points = Points([(0.2, 0.3, 'A', 0.1), (0.4, 0.7, 'B', 0.2), (0, 0.99, 'C', 0.3)], vdims=['cat', 'z'])
         img = aggregate(points, dynamic=False,  x_range=(0, 1), y_range=(0, 1),
@@ -883,10 +918,6 @@ class DatashaderRegridTests(ComparisonTestCase):
     Tests for datashader aggregation
     """
 
-    def setUp(self):
-        if ds_version <= Version('0.5.0'):
-            raise SkipTest('Regridding operations require datashader>=0.6.0')
-
     def test_regrid_mean(self):
         img = Image((range(10), range(5), np.arange(10) * np.arange(5)[np.newaxis].T))
         regridded = regrid(img, width=2, height=2, dynamic=False)
@@ -962,7 +993,7 @@ class DatashaderRasterizeTests(ComparisonTestCase):
     """
 
     def setUp(self):
-        if ds_version <= Version('0.6.4'):
+        if DATASHADER_VERSION <= (0, 6, 4):
             raise SkipTest('Regridding operations require datashader>=0.7.0')
 
         self.simplexes = [(0, 1, 2), (3, 2, 1)]
@@ -1279,12 +1310,12 @@ def test_rasterize_where_agg_no_column(point_plot, agg_input_fn, index_col):
     rast_input = dict(dynamic=False,  x_range=(-1, 1), y_range=(-1, 1), width=2, height=2)
     img = rasterize(point_plot, aggregator=agg_fn, **rast_input)
 
-    assert list(img.data) == ["index", "s", "val", "cat"]
+    assert list(img.data) == ["__index__", "s", "val", "cat"]
     assert list(img.vdims) == ["val", "s", "cat"]  # val first and no index
 
     # N=100 in point_data is chosen to have a big enough sample size
     # so that the index are not the same for the different agg_input_fn
-    np.testing.assert_array_equal(img.data["index"].data.flatten(), index_col)
+    np.testing.assert_array_equal(img.data["__index__"].data.flatten(), index_col)
 
     img_simple = rasterize(point_plot, aggregator=agg_input_fn("val"), **rast_input)
     np.testing.assert_array_equal(img_simple["val"], img["val"])
@@ -1322,7 +1353,7 @@ def test_rasterize_selector(point_plot, sel_fn):
     img = rasterize(point_plot, selector=sel_fn("val"), **rast_input)
 
     # Count is from the aggregator
-    assert list(img.data) == ["Count", "index", "s", "val", "cat"]
+    assert list(img.data) == ["Count", "__index__", "s", "val", "cat"]
     assert list(img.vdims) == ["Count", "s", "val", "cat"]  # no index
 
     # The output for the selector should be equal to the output for the aggregator using
@@ -1366,7 +1397,7 @@ class DatashaderSpreadTests(ComparisonTestCase):
         self.assertEqual(spreaded, RGB(arr))
 
     def test_spread_img_1px(self):
-        if ds_version < Version('0.12.0'):
+        if DATASHADER_VERSION < (0, 12, 0):
             raise SkipTest('Datashader does not support DataArray yet')
         arr = np.array([[0, 0, 0], [0, 0, 0], [1, 1, 1]]).T
         spreaded = spread(Image(arr))
@@ -1413,7 +1444,7 @@ class DatashaderStackTests(ComparisonTestCase):
 class GraphBundlingTests(ComparisonTestCase):
 
     def setUp(self):
-        if ds_version <= Version('0.7.0'):
+        if DATASHADER_VERSION <= (0, 7, 0):
             raise SkipTest('Regridding operations require datashader>=0.7.0')
         self.source = np.arange(8)
         self.target = np.zeros(8)
@@ -1423,14 +1454,25 @@ class GraphBundlingTests(ComparisonTestCase):
         direct = directly_connect_edges(self.graph)._split_edgepaths
         self.assertEqual(direct, self.graph.edgepaths)
 
+
 class InspectorTests(ComparisonTestCase):
     """
     Tests for inspector operations
     """
     def setUp(self):
         points = Points([(0.2, 0.3), (0.4, 0.7), (0, 0.99)])
-        self.pntsimg = rasterize(points, dynamic=False,
-                             x_range=(0, 1), y_range=(0, 1), width=4, height=4)
+        self.pntsimg = rasterize(
+            points, dynamic=False, height=4, width=4,
+            x_range=(0, 1), y_range=(0, 1)
+        )
+        date_pts = Points([
+            (np.datetime64('2024-09-25 11:00'), 0.3),
+            (np.datetime64('2024-09-25 11:01'), 0.7),
+            (np.datetime64('2024-09-25 11:04'), 0.99)])
+        self.datesimg = rasterize(
+            date_pts, dynamic=False, height=4, width=4,
+            x_range=(np.datetime64('2024-09-25 11:00'), np.datetime64('2024-09-25 11:04')), y_range=(0, 1)
+        )
         if spatialpandas is None:
             return
 
@@ -1445,7 +1487,6 @@ class InspectorTests(ComparisonTestCase):
     def tearDown(self):
         Tap.x, Tap.y = None, None
 
-
     def test_inspect_points_or_polygons(self):
         if spatialpandas is None:
             raise SkipTest('Polygon inspect tests require spatialpandas')
@@ -1502,6 +1543,18 @@ class InspectorTests(ComparisonTestCase):
         self.assertEqual(points.streams[0].x, 0.2)
         self.assertEqual(points.streams[0].y, 0.3)
 
+    def test_points_with_dates_inspection_1px_mask(self):
+        points = inspect_points(self.datesimg, max_indicators=3, dynamic=False, pixels=1,
+                                x=np.datetime64('2024-09-25 11:01'), y=-0.1)
+        self.assertEqual(points.dimension_values('x'), np.array([]))
+        self.assertEqual(points.dimension_values('y'), np.array([]))
+
+    def test_points_with_dates_inspection_2px_mask(self):
+        points = inspect_points(self.datesimg, max_indicators=3, dynamic=False, pixels=2,
+                                x=np.datetime64('2024-09-25 11:01'), y=-0.1)
+        self.assertEqual(points.dimension_values('x'), np.array([np.datetime64('2024-09-25 11:00')]))
+        self.assertEqual(points.dimension_values('y'), np.array([0.3]))
+
     def test_polys_inspection_1px_mask_hit(self):
         if spatialpandas is None:
             raise SkipTest('Polygon inspect tests require spatialpandas')
@@ -1521,7 +1574,6 @@ class InspectorTests(ComparisonTestCase):
         self.assertEqual(inspector.hits.iloc[0].geometry,
                          spatialpandas.geometry.polygon.Polygon(data))
 
-
     def test_polys_inspection_1px_mask_miss(self):
         if spatialpandas is None:
             raise SkipTest('Polygon inspect tests require spatialpandas')
diff --git a/holoviews/tests/plotting/bokeh/test_annotationplot.py b/holoviews/tests/plotting/bokeh/test_annotationplot.py
index 4520435b2..7ce298693 100644
--- a/holoviews/tests/plotting/bokeh/test_annotationplot.py
+++ b/holoviews/tests/plotting/bokeh/test_annotationplot.py
@@ -267,6 +267,21 @@ class TestHVLinesPlot(TestBokehPlot):
         assert list(source.data) == ["y"]
         assert (source.data["y"] == [0, 1, 2, 5.5]).all()
 
+    def test_hlines_plot_multi_y(self):
+        hlines = (
+            HLines({"y1": [1, 2, 3]}, 'y1') * HLines({'y2': [3, 4, 5]}, 'y2')
+        ).opts(multi_y=True)
+        plot = bokeh_renderer.get_plot(hlines)
+        sp1, sp2 = plot.subplots.values()
+        y1_range = sp1.handles['y_range']
+        assert y1_range.name == 'y1'
+        assert y1_range.start == 1
+        assert y1_range.end == 3
+        y2_range = sp2.handles['y_range']
+        assert y2_range.name == 'y2'
+        assert y2_range.start == 3
+        assert y2_range.end == 5
+
     def test_hlines_xlabel_ylabel(self):
         hlines = HLines(
             {"y": [0, 1, 2, 5.5], "extra": [-1, -2, -3, -44]}, vdims=["extra"]
diff --git a/holoviews/tests/plotting/bokeh/test_barplot.py b/holoviews/tests/plotting/bokeh/test_barplot.py
index 272b32103..777d81826 100644
--- a/holoviews/tests/plotting/bokeh/test_barplot.py
+++ b/holoviews/tests/plotting/bokeh/test_barplot.py
@@ -1,6 +1,11 @@
 import numpy as np
 import pandas as pd
-from bokeh.models import CategoricalColorMapper, LinearAxis, LinearColorMapper
+from bokeh.models import (
+    CategoricalColorMapper,
+    DatetimeAxis,
+    LinearAxis,
+    LinearColorMapper,
+)
 
 from holoviews.core.overlay import NdOverlay, Overlay
 from holoviews.element import Bars
@@ -313,6 +318,29 @@ class TestBarPlot(TestBokehPlot):
         plot = bokeh_renderer.get_plot(bars)
         np.testing.assert_almost_equal(plot.handles["glyph"].width, 69120000.0)
 
+    def test_bars_continuous_datetime_timezone_in_overlay(self):
+        # See: https://github.com/holoviz/holoviews/issues/6364
+        bars = Bars((pd.date_range("1/1/2000", periods=10, tz="UTC"), np.random.rand(10)))
+        overlay = Overlay([bars])
+        plot = bokeh_renderer.get_plot(overlay)
+        assert isinstance(plot.handles["xaxis"], DatetimeAxis)
+
+    def test_bars_continuous_datetime_stacked(self):
+        # See: https://github.com/holoviz/holoviews/issues/6288
+        data = pd.DataFrame({
+            "x": pd.to_datetime([
+                    "2017-01-01T00:00:00",
+                    "2017-01-01T00:00:00",
+                    "2017-01-01T01:00:00",
+                    "2017-01-01T01:00:00",
+                ]),
+            "cat": ["A", "B", "A", "B"],
+            "y": [1, 2, 3, 4],
+        })
+        bars = Bars(data, ["x", "cat"], ["y"]).opts(stacked=True)
+        plot = bokeh_renderer.get_plot(bars)
+        assert isinstance(plot.handles["xaxis"], DatetimeAxis)
+
     def test_bars_not_continuous_data_list(self):
         bars = Bars([("A", 1), ("B", 2), ("C", 3)])
         plot = bokeh_renderer.get_plot(bars)
@@ -356,3 +384,10 @@ class TestBarPlot(TestBokehPlot):
         )
         plot = bokeh_renderer.get_plot(bars)
         assert plot.handles["glyph"].width == 0.8
+
+    def test_bar_stacked_stack_variable_sorted(self):
+        # Check that if the stack dim is ordered
+        df = pd.DataFrame({"a": [*range(50), *range(50)], "b": sorted("ab" * 50), "c": range(100)})
+        bars = Bars(df, kdims=["a", "b"], vdims=["c"]).opts(stacked=True)
+        plot = bokeh_renderer.get_plot(bars)
+        assert plot.handles["glyph"].width == 0.8
diff --git a/holoviews/tests/plotting/bokeh/test_elementplot.py b/holoviews/tests/plotting/bokeh/test_elementplot.py
index c465ac46b..50c34c507 100644
--- a/holoviews/tests/plotting/bokeh/test_elementplot.py
+++ b/holoviews/tests/plotting/bokeh/test_elementplot.py
@@ -21,7 +21,7 @@ from holoviews import opts
 from holoviews.core import Dimension, DynamicMap, HoloMap, NdOverlay, Overlay
 from holoviews.core.util import dt_to_int
 from holoviews.element import Curve, HeatMap, Image, Labels, Scatter
-from holoviews.plotting.bokeh.util import bokeh34
+from holoviews.plotting.bokeh.util import bokeh34, bokeh36
 from holoviews.plotting.util import process_cmap
 from holoviews.streams import PointDraw, Stream
 from holoviews.util import render
@@ -885,6 +885,31 @@ class TestScalebarPlot:
         scalebar_icon = [tool for tool in toolbar.tools if tool.description == "Toggle ScaleBar"]
         assert len(scalebar_icon) == 1
 
+    @pytest.mark.skipif(not bokeh36, reason="requires Bokeh >= 3.6")
+    @pytest.mark.parametrize("enabled1", [True, False])
+    @pytest.mark.parametrize("enabled2", [True, False])
+    @pytest.mark.parametrize("enabled3", [True, False])
+    def test_scalebar_with_subcoordinate_y(self, enabled1, enabled2, enabled3):
+        from bokeh.models import ScaleBar
+
+        enabled = [enabled1, enabled2, enabled3]
+        curve1 = Curve([1, 2, 3], label='c1').opts(scalebar=enabled1, subcoordinate_y=True)
+        curve2 = Curve([1, 2, 3], label='c2').opts(scalebar=enabled2, subcoordinate_y=True)
+        curve3 = Curve([1, 2, 3], label='c3').opts(scalebar=enabled3, subcoordinate_y=True)
+        curves = curve1 * curve2 * curve3
+
+        plot = bokeh_renderer.get_plot(curves).handles["plot"]
+        coordinates = [r.coordinates for r in plot.renderers][::-1]
+        sb = (c for c in plot.center if isinstance(c, ScaleBar))
+        scalebars = [next(sb) if e else None for e in enabled]
+        assert sum(map(bool, scalebars)) == sum(enabled)
+
+        for coordinate, scalebar, idx in zip(coordinates, scalebars, "123"):
+            assert coordinate.y_source.name == f"c{idx}"
+            if scalebar is None:
+                continue
+            assert coordinate.y_source is scalebar.range
+
 
 class TestColorbarPlot(LoggingComparisonTestCase, TestBokehPlot):
 
diff --git a/holoviews/tests/plotting/bokeh/test_overlayplot.py b/holoviews/tests/plotting/bokeh/test_overlayplot.py
index 8f1e989eb..d2e65d1c4 100644
--- a/holoviews/tests/plotting/bokeh/test_overlayplot.py
+++ b/holoviews/tests/plotting/bokeh/test_overlayplot.py
@@ -15,7 +15,7 @@ from holoviews.element import (
     VLine,
 )
 from holoviews.plotting.bokeh.util import property_to_dict
-from holoviews.streams import Stream, Tap
+from holoviews.streams import Pipe, Stream, Tap
 from holoviews.util import Dynamic
 
 from ...utils import LoggingComparisonTestCase
@@ -307,6 +307,43 @@ class TestOverlayPlot(LoggingComparisonTestCase, TestBokehPlot):
             assert sp.handles['y_range'].start == 0
             assert sp.handles['y_range'].end == 27
 
+    def test_ndoverlay_subcoordinate_y_ranges_update(self):
+        def lines(data):
+            return NdOverlay({
+                'A': Curve(data, 'x', ('A', 'y1')).opts(subcoordinate_y=True, framewise=True),
+                'B': Curve(data, 'x', ('B', 'y2')).opts(subcoordinate_y=True, framewise=True),
+                'C': Curve(data, 'x', ('C', 'y3')).opts(subcoordinate_y=True, framewise=True),
+            })
+        data = {
+            'x': np.arange(10),
+            'A': np.arange(10),
+            'B': np.arange(10)*2,
+            'C': np.arange(10)*3
+        }
+        stream = Pipe(data=data)
+        dmap = DynamicMap(lines, streams=[stream])
+        plot = bokeh_renderer.get_plot(dmap)
+
+        assert plot.state.y_range.start == -0.5
+        assert plot.state.y_range.end == 2.5
+        for sp in plot.subplots.values():
+            y_range = sp.handles['y_range']
+            assert y_range.start == 0
+            assert y_range.end == data[y_range.name].max()
+
+        new_data = {
+            'x': np.arange(10),
+            'A': data['A'] + 1,
+            'B': data['B'] + 2,
+            'C': data['C'] + 3
+        }
+        stream.event(data=new_data)
+        for sp in plot.subplots.values():
+            y_range = sp.handles['y_range']
+            ydata = new_data[y_range.name]
+            assert y_range.start == ydata.min()
+            assert y_range.end == ydata.max()
+
     def test_overlay_subcoordinate_y_ranges(self):
         data = {
             'x': np.arange(10),
@@ -327,6 +364,43 @@ class TestOverlayPlot(LoggingComparisonTestCase, TestBokehPlot):
             assert sp.handles['y_range'].start == 0
             assert sp.handles['y_range'].end == 27
 
+    def test_overlay_subcoordinate_y_ranges_update(self):
+        def lines(data):
+            return Overlay([
+                Curve(data, 'x', ('A', 'y1'), label='A').opts(subcoordinate_y=True, framewise=True),
+                Curve(data, 'x', ('B', 'y2'), label='B').opts(subcoordinate_y=True, framewise=True),
+                Curve(data, 'x', ('C', 'y3'), label='C').opts(subcoordinate_y=True, framewise=True),
+            ])
+        data = {
+            'x': np.arange(10),
+            'A': np.arange(10),
+            'B': np.arange(10)*2,
+            'C': np.arange(10)*3
+        }
+        stream = Pipe(data=data)
+        dmap = DynamicMap(lines, streams=[stream])
+        plot = bokeh_renderer.get_plot(dmap)
+
+        assert plot.state.y_range.start == -0.5
+        assert plot.state.y_range.end == 2.5
+        for sp in plot.subplots.values():
+            y_range = sp.handles['y_range']
+            assert y_range.start == 0
+            assert y_range.end == data[y_range.name].max()
+
+        new_data = {
+            'x': np.arange(10),
+            'A': data['A'] + 1,
+            'B': data['B'] + 2,
+            'C': data['C'] + 3
+        }
+        stream.event(data=new_data)
+        for sp in plot.subplots.values():
+            y_range = sp.handles['y_range']
+            ydata = new_data[y_range.name]
+            assert y_range.start == ydata.min()
+            assert y_range.end == ydata.max()
+
 
 class TestLegends(TestBokehPlot):
 
diff --git a/holoviews/tests/plotting/bokeh/test_plot.py b/holoviews/tests/plotting/bokeh/test_plot.py
index b5c1a5978..03c9761f2 100644
--- a/holoviews/tests/plotting/bokeh/test_plot.py
+++ b/holoviews/tests/plotting/bokeh/test_plot.py
@@ -2,19 +2,23 @@ import numpy as np
 import pyviz_comms as comms
 from bokeh.models import (
     ColumnDataSource,
+    CrosshairTool,
     CustomJS,
     HoverTool,
     LinearColorMapper,
     LogColorMapper,
+    Span,
 )
 from param import concrete_descendents
 
 from holoviews import Curve
 from holoviews.core.element import Element
 from holoviews.core.options import Store
+from holoviews.core.spaces import DynamicMap
 from holoviews.element.comparison import ComparisonTestCase
 from holoviews.plotting.bokeh.callbacks import Callback
 from holoviews.plotting.bokeh.element import ElementPlot
+from holoviews.streams import Pipe
 
 bokeh_renderer = Store.renderers['bokeh']
 
@@ -80,11 +84,60 @@ class TestBokehPlot(ComparisonTestCase):
                     self.assertIn(lookup[2:-1], cds.data)
 
         # Ensure all the glyph renderers have a hover tool
-        print(renderers, hover)
         for renderer in renderers:
             self.assertTrue(any(renderer in h.renderers for h in hover))
 
 
+def test_element_plot_stream_cleanup():
+    stream = Pipe()
+
+    dmap = DynamicMap(Curve, streams=[stream])
+
+    plot = bokeh_renderer.get_plot(dmap)
+
+    assert len(stream._subscribers) == 1
+
+    plot.cleanup()
+
+    assert not stream._subscribers
+
+
+def test_overlay_plot_stream_cleanup():
+    stream1 = Pipe()
+    stream2 = Pipe()
+
+    dmap1 = DynamicMap(Curve, streams=[stream1])
+    dmap2 = DynamicMap(Curve, streams=[stream2])
+
+    plot = bokeh_renderer.get_plot(dmap1 * dmap2)
+
+    assert len(stream1._subscribers) == 4
+    assert len(stream2._subscribers) == 4
+
+    plot.cleanup()
+
+    assert not stream1._subscribers
+    assert not stream2._subscribers
+
+
+def test_layout_plot_stream_cleanup():
+    stream1 = Pipe()
+    stream2 = Pipe()
+
+    dmap1 = DynamicMap(Curve, streams=[stream1])
+    dmap2 = DynamicMap(Curve, streams=[stream2])
+
+    plot = bokeh_renderer.get_plot(dmap1 + dmap2)
+
+    assert len(stream1._subscribers) == 2
+    assert len(stream2._subscribers) == 2
+
+    plot.cleanup()
+
+    assert not stream1._subscribers
+    assert not stream2._subscribers
+
+
 def test_sync_two_plots():
     curve = lambda i: Curve(np.arange(10) * i, label="ABC"[i])
     plot1 = curve(0) * curve(1)
@@ -124,3 +177,17 @@ def test_sync_three_plots():
                 assert v[0].code == "dst.muted = src.muted"
                 assert isinstance(v[1], CustomJS)
                 assert v[1].code == "dst.muted = src.muted"
+
+
+def test_span_not_cloned_crosshair():
+    # See https://github.com/holoviz/holoviews/issues/6386
+    height = Span(dimension="height")
+    cht = CrosshairTool(overlay=height)
+
+    layout = Curve([]).opts(tools=[cht]) + Curve([]).opts(tools=[cht])
+
+    (fig1, *_), (fig2, *_) = bokeh_renderer.get_plot(layout).handles["plot"].children
+    tool1 = next(t for t in fig1.tools if isinstance(t, CrosshairTool))
+    tool2 = next(t for t in fig2.tools if isinstance(t, CrosshairTool))
+
+    assert tool1.overlay is tool2.overlay
diff --git a/holoviews/tests/plotting/matplotlib/test_barplot.py b/holoviews/tests/plotting/matplotlib/test_barplot.py
index 5e4f85cd1..7b9d66f50 100644
--- a/holoviews/tests/plotting/matplotlib/test_barplot.py
+++ b/holoviews/tests/plotting/matplotlib/test_barplot.py
@@ -45,7 +45,7 @@ class TestBarPlot(LoggingComparisonTestCase, TestMPLPlot):
         bars = Bars([("A", 1), ("B", 2), ("C", 3)])
         plot = mpl_renderer.get_plot(bars)
         ax = plot.handles["axis"]
-        np.testing.assert_almost_equal(ax.get_xlim(), (-0.54, 2.54))
+        np.testing.assert_almost_equal(ax.get_xlim(), (-0.4, 2.4))
         assert ax.patches[0].get_width() == 0.8
         np.testing.assert_equal(ax.get_xticks(), [0, 1, 2])
         np.testing.assert_equal(
@@ -69,7 +69,7 @@ class TestBarPlot(LoggingComparisonTestCase, TestMPLPlot):
         plot = mpl_renderer.get_plot(bars)
         ax = plot.handles["axis"]
 
-        np.testing.assert_almost_equal(ax.get_xlim(), (-0.3233333, 3.8566667))
+        np.testing.assert_almost_equal(ax.get_xlim(), (-0.1333333,  3.6666667))
         assert ax.patches[0].get_width() == 0.26666666666666666
         ticklabels = ax.get_xticklabels()
         expected = [
@@ -113,7 +113,7 @@ class TestBarPlot(LoggingComparisonTestCase, TestMPLPlot):
         plot = mpl_renderer.get_plot(bars)
         ax = plot.handles["axis"]
 
-        np.testing.assert_almost_equal(ax.get_xlim(), (-0.59, 3.59))
+        np.testing.assert_almost_equal(ax.get_xlim(), (-0.4, 3.4))
         assert ax.patches[0].get_width() == 0.8
         ticklabels = ax.get_xticklabels()
         expected = [
@@ -136,7 +136,7 @@ class TestBarPlot(LoggingComparisonTestCase, TestMPLPlot):
         plot = mpl_renderer.get_plot(bars)
         ax = plot.handles["axis"]
 
-        np.testing.assert_almost_equal(ax.get_xlim(), (-0.34,  2.74))
+        np.testing.assert_almost_equal(ax.get_xlim(), (-0.2,  2.6))
         assert ax.patches[0].get_width() == 0.4
         assert len(ax.get_xticks()) > 3
 
diff --git a/holoviews/tests/plotting/matplotlib/test_layoutplot.py b/holoviews/tests/plotting/matplotlib/test_layoutplot.py
index 97ed2ab41..5c32d7395 100644
--- a/holoviews/tests/plotting/matplotlib/test_layoutplot.py
+++ b/holoviews/tests/plotting/matplotlib/test_layoutplot.py
@@ -58,3 +58,25 @@ class TestLayoutPlot(LoggingComparisonTestCase, TestMPLPlot):
         cp1, cp2 = plot.traverse(lambda x: x, [CurvePlot])
         self.assertTrue(cp1.handles['axis'].get_ylim(), (1, 3))
         self.assertTrue(cp2.handles['axis'].get_ylim(), (10, 30))
+
+    def test_layout_sublabel_offset(self):
+        from holoviews.plotting.mpl import CurvePlot
+        layout = Curve([]) + Curve([]) + Curve([]) + Curve([])
+        layout.opts(sublabel_offset=1)
+        plot = mpl_renderer.get_plot(layout)
+        cps = plot.traverse(lambda x: x, [CurvePlot])
+        assert cps[0].handles["sublabel"].get_text() == "B"
+        assert cps[1].handles["sublabel"].get_text() == "C"
+        assert cps[2].handles["sublabel"].get_text() == "D"
+        assert cps[3].handles["sublabel"].get_text() == "E"
+
+    def test_layout_sublabel_skip(self):
+        from holoviews.plotting.mpl import CurvePlot
+        layout = Curve([]) + Curve([]) + Curve([]) + Curve([])
+        layout.opts(sublabel_skip=[1, 3])
+        plot = mpl_renderer.get_plot(layout)
+        cps = plot.traverse(lambda x: x, [CurvePlot])
+        assert "sublabel" not in cps[0].handles
+        assert cps[1].handles["sublabel"].get_text() == "A"
+        assert "sublabel" not in cps[2].handles
+        assert cps[3].handles["sublabel"].get_text() == "B"
diff --git a/holoviews/tests/plotting/plotly/test_plot.py b/holoviews/tests/plotting/plotly/test_plot.py
index 3a033f884..0f71a9952 100644
--- a/holoviews/tests/plotting/plotly/test_plot.py
+++ b/holoviews/tests/plotting/plotly/test_plot.py
@@ -2,16 +2,68 @@ import plotly.graph_objs as go
 import pyviz_comms as comms
 from param import concrete_descendents
 
-from holoviews.core import Store
+from holoviews.core import DynamicMap, Store
+from holoviews.element import Curve
 from holoviews.element.comparison import ComparisonTestCase
 from holoviews.plotting.plotly.element import ElementPlot
 from holoviews.plotting.plotly.util import figure_grid
+from holoviews.streams import Pipe
 
 from .. import option_intersections
 
 plotly_renderer = Store.renderers['plotly']
 
 
+def test_element_plot_stream_cleanup():
+    stream = Pipe()
+
+    dmap = DynamicMap(Curve, streams=[stream])
+
+    plot = plotly_renderer.get_plot(dmap)
+
+    assert len(stream._subscribers) == 1
+
+    plot.cleanup()
+
+    assert not stream._subscribers
+
+
+def test_overlay_plot_stream_cleanup():
+    stream1 = Pipe()
+    stream2 = Pipe()
+
+    dmap1 = DynamicMap(Curve, streams=[stream1])
+    dmap2 = DynamicMap(Curve, streams=[stream2])
+
+    plot = plotly_renderer.get_plot(dmap1 * dmap2)
+
+    assert len(stream1._subscribers) == 4
+    assert len(stream2._subscribers) == 4
+
+    plot.cleanup()
+
+    assert not stream1._subscribers
+    assert not stream2._subscribers
+
+
+def test_layout_plot_stream_cleanup():
+    stream1 = Pipe()
+    stream2 = Pipe()
+
+    dmap1 = DynamicMap(Curve, streams=[stream1])
+    dmap2 = DynamicMap(Curve, streams=[stream2])
+
+    plot = plotly_renderer.get_plot(dmap1 + dmap2)
+
+    assert len(stream1._subscribers) == 2
+    assert len(stream2._subscribers) == 2
+
+    plot.cleanup()
+
+    assert not stream1._subscribers
+    assert not stream2._subscribers
+
+
 class TestPlotlyPlot(ComparisonTestCase):
 
     def setUp(self):
diff --git a/holoviews/tests/test_streams.py b/holoviews/tests/test_streams.py
index 2bc957f67..3f7124a85 100644
--- a/holoviews/tests/test_streams.py
+++ b/holoviews/tests/test_streams.py
@@ -938,7 +938,7 @@ class TestBufferDataFrameStream(ComparisonTestCase):
     def test_init_buffer_dframe_with_index(self):
         data = pd.DataFrame({'x': np.array([1]), 'y': np.array([2])})
         buff = Buffer(data)
-        self.assertEqual(buff.data, data.reset_index())
+        self.assertEqual(buff.data, data)
 
     def test_buffer_dframe_send(self):
         data = pd.DataFrame({'x': np.array([0]), 'y': np.array([1])})
@@ -952,7 +952,7 @@ class TestBufferDataFrameStream(ComparisonTestCase):
         buff = Buffer(data)
         buff.send(pd.DataFrame({'x': np.array([1]), 'y': np.array([2])}))
         dframe = pd.DataFrame({'x': np.array([0, 1]), 'y': np.array([1, 2])}, index=[0, 0])
-        self.assertEqual(buff.data.values, dframe.reset_index().values)
+        pd.testing.assert_frame_equal(buff.data, dframe)
 
     def test_buffer_dframe_larger_than_length(self):
         data = pd.DataFrame({'x': np.array([0]), 'y': np.array([1])})
@@ -980,7 +980,7 @@ class TestBufferDataFrameStream(ComparisonTestCase):
         data = pd.DataFrame({'a': [1, 2, 3]})
         buff = Buffer(data)
         buff.clear()
-        self.assertEqual(buff.data, data.iloc[:0, :].reset_index())
+        pd.testing.assert_frame_equal(buff.data, data.iloc[:0, :])
 
 
 class Sum(Derived):
diff --git a/holoviews/tests/ui/bokeh/test_callback.py b/holoviews/tests/ui/bokeh/test_callback.py
index 725df2869..1fe6ba082 100644
--- a/holoviews/tests/ui/bokeh/test_callback.py
+++ b/holoviews/tests/ui/bokeh/test_callback.py
@@ -1,11 +1,12 @@
 import numpy as np
+import pandas as pd
 import panel as pn
 import pytest
 
 import holoviews as hv
 from holoviews import Curve, DynamicMap, Scatter
 from holoviews.plotting.bokeh.util import bokeh34
-from holoviews.streams import BoundsX, BoundsXY, BoundsY, Lasso, RangeXY
+from holoviews.streams import BoundsX, BoundsXY, BoundsY, Lasso, MultiAxisTap, RangeXY
 
 from .. import expect, wait_until
 
@@ -164,6 +165,54 @@ def test_multi_axis_rangexy(serve_hv):
     ), page)
 
 
+@pytest.mark.usefixtures("bokeh_backend")
+def test_multi_axis_tap(serve_hv):
+    c1 = Curve(np.arange(10).cumsum(), vdims='y1')
+    c2 = Curve(np.arange(20).cumsum(), vdims='y2')
+
+    overlay = (c1 * c2).opts(multi_y=True)
+
+    s = MultiAxisTap(source=overlay)
+
+    page = serve_hv(overlay)
+
+    hv_plot = page.locator('.bk-events')
+
+    expect(hv_plot).to_have_count(1)
+
+    hv_plot.click()
+
+    def test():
+        assert s.xs == {'x': 11.560240963855422}
+        assert s.ys == {'y1': 18.642857142857146, 'y2': 78.71428571428572}
+
+    wait_until(test, page)
+
+
+@pytest.mark.usefixtures("bokeh_backend")
+def test_multi_axis_tap_datetime(serve_hv):
+    c1 = Curve((pd.date_range('2024-01-01', '2024-01-10'), np.arange(10).cumsum()), vdims='y1')
+    c2 = Curve((pd.date_range('2024-01-01', '2024-01-20'), np.arange(20).cumsum()), vdims='y2')
+
+    overlay = (c1 * c2).opts(multi_y=True)
+
+    s = MultiAxisTap(source=overlay)
+
+    page = serve_hv(overlay)
+
+    hv_plot = page.locator('.bk-events')
+
+    expect(hv_plot).to_have_count(1)
+
+    hv_plot.click()
+
+    def test():
+        assert s.xs == {'x': np.datetime64('2024-01-12T13:26:44.819277')}
+        assert s.ys == {'y1': 18.13070539419087, 'y2': 76.551867219917}
+
+    wait_until(test, page)
+
+
 @pytest.mark.usefixtures("bokeh_backend")
 def test_bind_trigger(serve_hv):
     # Regression test for https://github.com/holoviz/holoviews/issues/6013
@@ -213,7 +262,43 @@ def test_stream_popup(serve_hv):
     hv_plot.click()
     expect(hv_plot).to_have_count(1)
 
-    locator = page.locator("#tap")
+    locator = page.locator(".markdown")
+    expect(locator).to_have_count(1)
+
+
+@skip_popup
+@pytest.mark.usefixtures("bokeh_backend")
+def test_stream_popup_polygons_tap(serve_hv):
+    def popup_form(name):
+        return f"# {name}"
+
+    points = hv.Polygons([(0, 0), (0, 1), (1, 1), (1, 0)]).opts(tools=["tap"])
+    hv.streams.Tap(source=points, popup=popup_form("Tap"))
+
+    page = serve_hv(points)
+    hv_plot = page.locator('.bk-events')
+    hv_plot.click()
+    expect(hv_plot).to_have_count(1)
+
+    locator = page.locator(".markdown")
+    expect(locator).to_have_count(1)
+
+
+@skip_popup
+@pytest.mark.usefixtures("bokeh_backend")
+def test_stream_popup_polygons_selection1d(serve_hv):
+    def popup_form(name):
+        return f"# {name}"
+
+    points = hv.Polygons([(0, 0), (0, 1), (1, 1), (1, 0)]).opts(tools=["tap"])
+    hv.streams.Selection1D(source=points, popup=popup_form("Tap"))
+
+    page = serve_hv(points)
+    hv_plot = page.locator('.bk-events')
+    hv_plot.click()
+    expect(hv_plot).to_have_count(1)
+
+    locator = page.locator(".markdown")
     expect(locator).to_have_count(1)