Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Use contourpy for contour calculations #5910

Merged
merged 15 commits into from
Oct 16, 2023
Merged

Use contourpy for contour calculations #5910

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
15 commits
Select commit Hold shift + click to select a range
1bef1ae
Use contourpy for contour calculations
ianthomas23 Sep 28, 2023
a398dad
Add contourpy to tests_core dependencies
ianthomas23 Oct 3, 2023
4fe9e2f
Add minimal MaxNLocator implementation
ianthomas23 Oct 3, 2023
bdac077
Add MaxNLocator unit tests
ianthomas23 Oct 4, 2023
b33418e
Correctly test contour calculations with datetime dimensions
ianthomas23 Oct 9, 2023
0d653a5
Raise error if attempting filled contours using datetime spatial coor…
ianthomas23 Oct 9, 2023
fb4eb3f
Only require matplotlib for contours if using datatime data
ianthomas23 Oct 13, 2023
662d4fb
Improved testing
ianthomas23 Oct 13, 2023
4d5efe7
Update holoviews/operation/element.py
ianthomas23 Oct 13, 2023
682c3d1
Request correct number contour levels in MaxNLocator
ianthomas23 Oct 13, 2023
4f9af4f
Correctly propagate empty contours
ianthomas23 Oct 13, 2023
279ab20
Review comments
ianthomas23 Oct 13, 2023
7dcf8fe
Merge branch 'main' into contours_via_contourpy
hoxbro Oct 13, 2023
939456e
Sort imports
hoxbro Oct 13, 2023
ba342a5
Don't create empty contour geometries
ianthomas23 Oct 16, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
136 changes: 89 additions & 47 deletions holoviews/operation/element.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
from ..element.raster import RGB, Image
from ..element.util import categorical_aggregate2d # noqa (API import)
from ..streams import RangeXY
from ..util.locator import MaxNLocator

column_interfaces = [ArrayInterface, DictInterface, PandasInterface]

Expand Down Expand Up @@ -560,13 +561,9 @@ class contours(Operation):

def _process(self, element, key=None):
try:
from matplotlib.axes import Axes
from matplotlib.contour import QuadContourSet
from matplotlib.dates import date2num, num2date
from matplotlib.figure import Figure
from contourpy import FillType, LineType, contour_generator
except ImportError:
raise ImportError("contours operation requires matplotlib.")
extent = element.range(0) + element.range(1)[::-1]
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is no longer used as I don't think it is required. There is explicit code to set the x and y arrays and pass them through, so the extent kwarg would have been ignored.

raise ImportError("contours operation requires contourpy.")

xs = element.dimension_values(0, True, flat=False)
ys = element.dimension_values(1, True, flat=False)
Expand All @@ -586,6 +583,15 @@ def _process(self, element, key=None):
# if any data is a datetime, transform to matplotlib's numerical format
data_is_datetime = tuple(isdatetime(arr) for k, arr in enumerate(data))
if any(data_is_datetime):
if any(data_is_datetime[:2]) and self.p.filled:
raise RuntimeError("Datetime spatial coordinates are not supported "
"for filled contour calculations.")

try:
from matplotlib.dates import date2num, num2date
except ImportError:
raise ImportError("contours operation using datetimes requires matplotlib.") from None

data = tuple(
date2num(d) if is_datetime else d
for d, is_datetime in zip(data, data_is_datetime)
Expand All @@ -598,61 +604,97 @@ def _process(self, element, key=None):
contour_type = Contours
vdims = element.vdims[:1]

kwargs = {}
levels = self.p.levels
zmin, zmax = element.range(2)
if isinstance(self.p.levels, int):
if isinstance(levels, int):
if zmin == zmax:
contours = contour_type([], [xdim, ydim], vdims)
return (element * contours) if self.p.overlaid else contours
data += (levels,)
else:
# The +1 is consistent with Matplotlib's use of MaxNLocator for contours.
locator = MaxNLocator(levels + 1)
hoxbro marked this conversation as resolved.
Show resolved Hide resolved
levels = locator.tick_values(zmin, zmax)
else:
kwargs = {'levels': levels}
levels = np.array(levels)

if data_is_datetime[2]:
levels = date2num(levels)

fig = Figure()
ax = Axes(fig, [0, 0, 1, 1])
contour_set = QuadContourSet(ax, *data, filled=self.p.filled,
extent=extent, **kwargs)
levels = np.array(contour_set.get_array())
crange = levels.min(), levels.max()
if self.p.filled:
levels = levels[:-1] + np.diff(levels)/2.
vdims = [vdims[0].clone(range=crange)]

cont_gen = contour_generator(
*data,
line_type=LineType.ChunkCombinedOffset,
fill_type=FillType.ChunkCombinedOffsetOffset,
)

def points_to_datetime(points):
# transform x/y coordinates back to datetimes
xs, ys = np.split(points, 2, axis=1)
if data_is_datetime[0]:
xs = np.array(num2date(xs))
if data_is_datetime[1]:
ys = np.array(num2date(ys))
return np.concatenate((xs, ys), axis=1)

paths = []
empty = np.array([[np.nan, np.nan]])
for level, cset in zip(levels, contour_set.collections):
exteriors = []
interiors = []
for geom in cset.get_paths():
interior = []
polys = geom.to_polygons(closed_only=False)
for ncp, cp in enumerate(polys):
if any(data_is_datetime[0:2]):
# transform x/y coordinates back to datetimes
xs, ys = np.split(cp, 2, axis=1)
if data_is_datetime[0]:
xs = np.array(num2date(xs))
if data_is_datetime[1]:
ys = np.array(num2date(ys))
cp = np.concatenate((xs, ys), axis=1)
if ncp == 0:
exteriors.append(cp)
if self.p.filled:
empty = np.array([[np.nan, np.nan]])
for lower_level, upper_level in zip(levels[:-1], levels[1:]):
filled = cont_gen.filled(lower_level, upper_level)
# Only have to consider last index 0 as we are using contourpy without chunking
if (points := filled[0][0]) is None:
continue

exteriors = []
interiors = []
if any(data_is_datetime[0:2]):
points = points_to_datetime(points)

offsets = filled[1][0]
outer_offsets = filled[2][0]

# Loop through exterior polygon boundaries.
for jstart, jend in zip(outer_offsets[:-1], outer_offsets[1:]):
if exteriors:
exteriors.append(empty)
else:
interior.append(cp)
if len(polys):
exteriors.append(points[offsets[jstart]:offsets[jstart + 1]])

# Loop over the (jend-jstart-1) interior boundaries.
interior = [points[offsets[j]:offsets[j + 1]] for j in range(jstart+1, jend)]
interiors.append(interior)
if not exteriors:
continue
geom = {
element.vdims[0].name:
num2date(level) if data_is_datetime[2] else level,
(xdim, ydim): np.concatenate(exteriors[:-1])
}
if self.p.filled and interiors:
geom['holes'] = interiors
paths.append(geom)
level = (lower_level + upper_level) / 2
geom = {
element.vdims[0].name:
num2date(level) if data_is_datetime[2] else level,
(xdim, ydim): np.concatenate(exteriors) if exteriors else [],
}
if interiors:
geom['holes'] = interiors
paths.append(geom)
else:
for level in levels:
lines = cont_gen.lines(level)
# Only have to consider last index 0 as we are using contourpy without chunking
if (points := lines[0][0]) is None:
continue

if any(data_is_datetime[0:2]):
points = points_to_datetime(points)

offsets = lines[1][0]
if offsets is not None and len(offsets) > 2:
# Casting offsets to int64 to avoid possible numpy UFuncOutputCastingError
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a concern. It isn't needed for test_operation.py but is for one of the KDE tests in test_statsoperations.py.

Copy link
Member

@hoxbro hoxbro Oct 2, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What was the type of the KDE test?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

With a clean environment this is now needed to avoid the "UFunc add" error in many tests (test_operation.py as well as test_statsoperations.py). The complaining line is in numpy function_base.py when adding extra values to these uint32 offsets. The values added could in theory be negative (although they are always positive in practice) and hence it doesn't like the unsigned-ness of uint32. int64 is the next largest signed integer, so it seems a sensible choice.

offsets = offsets[1:-1].astype(np.int64)
points = np.insert(points, offsets, np.nan, axis=0)
geom = {
element.vdims[0].name:
num2date(level) if data_is_datetime[2] else level,
(xdim, ydim): points if points is not None else [],
}
paths.append(geom)
contours = contour_type(paths, label=element.label, kdims=element.kdims, vdims=vdims)
if self.p.overlaid:
contours = element * contours
Expand Down
179 changes: 171 additions & 8 deletions holoviews/tests/operation/test_operation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -94,20 +94,126 @@ def test_image_gradient(self):
self.assertEqual(op_img, img.clone(np.array([[3.162278, 3.162278], [3.162278, 3.162278]]), group='Gradient'))

def test_image_contours(self):
img = Image(np.array([[0, 1, 0], [3, 4, 5.], [6, 7, 8]]))
img = Image(np.array([[0, 1, 0], [0, 1, 0]]))
op_contours = contours(img, levels=[0.5])
contour = Contours([[(-0.166667, 0.333333, 0.5), (-0.333333, 0.277778, 0.5),
(np.nan, np.nan, 0.5), (0.333333, 0.3, 0.5),
(0.166667, 0.333333, 0.5)]],
# Note multiple lines which are nan-separated.
contour = Contours([[(-0.166667, 0.25, 0.5), (-0.1666667, -0.25, 0.5),
(np.nan, np.nan, 0.5), (0.1666667, -0.25, 0.5),
(0.1666667, 0.25, 0.5)]],
vdims=img.vdims)
self.assertEqual(op_contours, contour)

def test_image_contours_empty(self):
img = Image(np.array([[0, 1, 0], [0, 1, 0]]))
# Contour level outside of data limits
op_contours = contours(img, levels=[23.0])
contour = Contours([], vdims=img.vdims)
self.assertEqual(op_contours, contour)

def test_image_contours_auto_levels(self):
z = np.array([[0, 1, 0], [3, 4, 5.], [6, 7, 8]])
img = Image(z)
for nlevels in range(3, 20):
op_contours = contours(img, levels=nlevels)
levels = [item['z'] for item in op_contours.data]
assert len(levels) <= nlevels + 2
assert np.min(levels) <= z.min()
assert np.max(levels) < z.max()

def test_image_contours_no_range(self):
img = Image(np.zeros((2, 2)))
op_contours = contours(img, levels=2)
contour = Contours([], vdims=img.vdims)
self.assertEqual(op_contours, contour)

def test_image_contours_x_datetime(self):
x = np.array(['2023-09-01', '2023-09-03', '2023-09-05'], dtype='datetime64')
y = [14, 15]
z = np.array([[0, 1, 0], [0, 1, 0]])
img = Image((x, y, z))
op_contours = contours(img, levels=[0.5])
# Note multiple lines which are nan-separated.
tz = dt.timezone.utc
expected_x = np.array(
[dt.datetime(2023, 9, 2, tzinfo=tz), dt.datetime(2023, 9, 2, tzinfo=tz), np.nan,
dt.datetime(2023, 9, 4, tzinfo=tz), dt.datetime(2023, 9, 4, tzinfo=tz)],
dtype=object)

# Separately compare nans and datetimes
x = op_contours.dimension_values('x')
mask = np.array([True, True, False, True, True]) # Mask ignoring nans
np.testing.assert_array_equal(x[mask], expected_x[mask])
np.testing.assert_array_equal(x[~mask].astype(float), expected_x[~mask].astype(float))

np.testing.assert_array_almost_equal(op_contours.dimension_values('y').astype(float),
[15, 14, np.nan, 14, 15])
np.testing.assert_array_almost_equal(op_contours.dimension_values('z'), [0.5]*5)

def test_image_contours_y_datetime(self):
x = [14, 15, 16]
y = np.array(['2023-09-01', '2023-09-03'], dtype='datetime64')
z = np.array([[0, 1, 0], [0, 1, 0]])
img = Image((x, y, z))
op_contours = contours(img, levels=[0.5])
# Note multiple lines which are nan-separated.
np.testing.assert_array_almost_equal(op_contours.dimension_values('x').astype(float),
[14.5, 14.5, np.nan, 15.5, 15.5])

tz = dt.timezone.utc
expected_y = np.array(
[dt.datetime(2023, 9, 3, tzinfo=tz), dt.datetime(2023, 9, 1, tzinfo=tz), np.nan,
dt.datetime(2023, 9, 1, tzinfo=tz), dt.datetime(2023, 9, 3, tzinfo=tz)],
dtype=object)

# Separately compare nans and datetimes
y = op_contours.dimension_values('y')
mask = np.array([True, True, False, True, True]) # Mask ignoring nans
np.testing.assert_array_equal(y[mask], expected_y[mask])
np.testing.assert_array_equal(y[~mask].astype(float), expected_y[~mask].astype(float))

np.testing.assert_array_almost_equal(op_contours.dimension_values('z'), [0.5]*5)

def test_image_contours_xy_datetime(self):
x = np.array(['2023-09-01', '2023-09-03', '2023-09-05'], dtype='datetime64')
y = np.array(['2023-10-07', '2023-10-08'], dtype='datetime64')
z = np.array([[0, 1, 0], [0, 1, 0]])
img = Image((x, y, z))
op_contours = contours(img, levels=[0.5])
# Note multiple lines which are nan-separated.

tz = dt.timezone.utc
expected_x = np.array(
[dt.datetime(2023, 9, 2, tzinfo=tz), dt.datetime(2023, 9, 2, tzinfo=tz), np.nan,
dt.datetime(2023, 9, 4, tzinfo=tz), dt.datetime(2023, 9, 4, tzinfo=tz)],
dtype=object)
expected_y = np.array(
[dt.datetime(2023, 10, 8, tzinfo=tz), dt.datetime(2023, 10, 7, tzinfo=tz), np.nan,
dt.datetime(2023, 10, 7, tzinfo=tz), dt.datetime(2023, 10, 8, tzinfo=tz)],
dtype=object)

# Separately compare nans and datetimes
x = op_contours.dimension_values('x')
mask = np.array([True, True, False, True, True]) # Mask ignoring nans
np.testing.assert_array_equal(x[mask], expected_x[mask])
np.testing.assert_array_equal(x[~mask].astype(float), expected_x[~mask].astype(float))

y = op_contours.dimension_values('y')
np.testing.assert_array_equal(y[mask], expected_y[mask])
np.testing.assert_array_equal(y[~mask].astype(float), expected_y[~mask].astype(float))

np.testing.assert_array_almost_equal(op_contours.dimension_values('z'), [0.5]*5)

def test_image_contours_z_datetime(self):
z = np.array([['2023-09-10', '2023-09-10'], ['2023-09-10', '2023-09-12']], dtype='datetime64')
img = Image(z)
op_contours = contours(img, levels=[np.datetime64('2023-09-11')])
np.testing.assert_array_almost_equal(op_contours.dimension_values('x'), [0.25, 0.0])
np.testing.assert_array_almost_equal(op_contours.dimension_values('y'), [0.0, -0.25])
expected_z = np.array([
dt.datetime(2023, 9, 11, 0, 0, tzinfo=dt.timezone.utc),
dt.datetime(2023, 9, 11, 0, 0, tzinfo=dt.timezone.utc)], dtype=object)
np.testing.assert_array_equal(op_contours.dimension_values('z'), expected_z)

def test_qmesh_contours(self):
qmesh = QuadMesh(([0, 1, 2], [1, 2, 3], np.array([[0, 1, 0], [3, 4, 5.], [6, 7, 8]])))
op_contours = contours(qmesh, levels=[0.5])
Expand Down Expand Up @@ -146,13 +252,70 @@ def test_qmesh_curvilinear_edges_contours(self):
self.assertEqual(op_contours, contour)

def test_image_contours_filled(self):
img = Image(np.array([[0, 2, 0], [0, 2, 0]]))
# Two polygons (nan-separated) without holes
op_contours = contours(img, filled=True, levels=[0.5, 1.5])
data = [[(-0.25, -0.25, 1), (-0.08333333, -0.25, 1), (-0.08333333, 0.25, 1),
(-0.25, 0.25, 1), (-0.25, -0.25, 1), (np.nan, np.nan, 1), (0.08333333, -0.25, 1),
(0.25, -0.25, 1), (0.25, 0.25, 1), (0.08333333, 0.25, 1), (0.08333333, -0.25, 1)]]
polys = Polygons(data, vdims=img.vdims[0].clone(range=(0.5, 1.5)))
self.assertEqual(op_contours, polys)

def test_image_contours_filled_with_hole(self):
img = Image(np.array([[0, 0, 0], [0, 1, 0.], [0, 0, 0]]))
# Single polygon with hole
op_contours = contours(img, filled=True, levels=[0.25, 0.75])
data = [[(-0.25, 0.0, 0.5), (0.0, -0.25, 0.5), (0.25, 0.0, 0.5), (0.0, 0.25, 0.5),
(-0.25, 0.0, 0.5)]]
polys = Polygons(data, vdims=img.vdims[0].clone(range=(0.25, 0.75)))
self.assertEqual(op_contours, polys)
expected_holes = [[[np.array([[0.0, -0.08333333], [-0.08333333, 0.0], [0.0, 0.08333333],
[0.08333333, 0.0], [0.0, -0.08333333]])]]]
np.testing.assert_array_almost_equal(op_contours.holes(), expected_holes)

def test_image_contours_filled_multi_holes(self):
img = Image(np.array([[0, 0, 0, 0, 0], [0, 1, 0, 1, 0], [0, 0, 0, 0, 0]]))
# Single polygon with two holes
op_contours = contours(img, filled=True, levels=[-0.5, 0.5])
data = [[(-0.4, -0.3333333, 0), (-0.2, -0.3333333, 0), (0, -0.3333333, 0),
(0.2, -0.3333333, 0), (0.4, -0.3333333, 0), (0.4, 0, 0), (0.4, 0.3333333, 0),
(0.2, 0.3333333, 0), (0, 0.3333333, 0), (-0.2, 0.3333333, 0), (-0.4, 0.3333333, 0),
(-0.4, 0, 0), (-0.4, -0.3333333, 0)]]
polys = Polygons(data, vdims=img.vdims[0].clone(range=(-0.5, 0.5)))
self.assertEqual(op_contours, polys)
expected_holes = [[[np.array([[-0.2, -0.16666667], [-0.3, 0], [-0.2, 0.16666667], [-0.1, 0],
[-0.2, -0.16666667]]),
np.array([[0.2, -0.16666667], [0.1, 0], [0.2, 0.16666667], [0.3, 0],
[0.2, -0.16666667]])]]]
np.testing.assert_array_almost_equal(op_contours.holes(), expected_holes)

def test_image_contours_filled_empty(self):
img = Image(np.array([[0, 1, 0], [3, 4, 5.], [6, 7, 8]]))
op_contours = contours(img, filled=True, levels=[2, 2.5])
data = [[(0., 0.166667, 2.25), (0.333333, 0.166667, 2.25), (0.333333, 0.2, 2.25), (0., 0.222222, 2.25),
(-0.333333, 0.111111, 2.25), (-0.333333, 0.055556, 2.25), (0., 0.166667, 2.25)]]
polys = Polygons(data, vdims=img.vdims[0].clone(range=(2, 2.5)))
# Contour level outside of data limits
op_contours = contours(img, filled=True, levels=[20.0, 23.0])
polys = Polygons([], vdims=img.vdims[0].clone(range=(20.0, 23.0)))
self.assertEqual(op_contours, polys)

def test_image_contours_filled_auto_levels(self):
z = np.array([[0, 1, 0], [3, 4, 5], [6, 7, 8]])
img = Image(z)
for nlevels in range(3, 20):
op_contours = contours(img, filled=True, levels=nlevels)
levels = [item['z'] for item in op_contours.data]
assert len(levels) <= nlevels + 1
hoxbro marked this conversation as resolved.
Show resolved Hide resolved
delta = 0.5*(levels[1] - levels[0])
assert np.min(levels) <= z.min() + delta
assert np.max(levels) >= z.max() - delta

def test_image_contours_filled_x_datetime(self):
x = np.array(['2023-09-01', '2023-09-05', '2023-09-09'], dtype='datetime64')
y = np.array([6, 7])
z = np.array([[0, 2, 0], [0, 2, 0]])
img = Image((x, y, z))
msg = r'Datetime spatial coordinates are not supported for filled contour calculations.'
with pytest.raises(RuntimeError, match=msg):
_ = contours(img, filled=True, levels=[0.5, 1.5])

def test_points_histogram(self):
points = Points([float(i) for i in range(10)])
op_hist = histogram(points, num_bins=3, normed=True)
Expand Down
Loading