From 9e358325991c0552c6836761d6dea7d119112c76 Mon Sep 17 00:00:00 2001 From: Berk Geveci Date: Wed, 4 Nov 2026 15:32:03 -0500 Subject: [PATCH 1/3] VTK algorithms' dataset cache: observers kept reachable, entries dropped at the next use The cache of Tack datasets per VTK dataset drops an entry from a DeleteEvent observer, a function only the VTK object referred to. When Python's collector freed such a dataset in a reference cycle made after it was read (a failing test's traceback holding its frames), it took the observer for garbage too and cleared its globals and closure first; the deletion then called it, and Python segfaulted (_PyEval_EvalFrameDefault from vtkPythonCommand::Execute, from vtkObject::ObjectFinalize, from the collector). VTK 9.7.1 and master both; a VTK-only reproduction is VTK bugs 19 in the vault. The observers are now kept in a module-level dict until they fire, and do no more than note the address; the cache drops the entry, and the observer, at its next use, outside whatever freed the dataset. The regression test makes the cycle after the dataset was read, as the crash needs; it segfaulted before. --- changelog.d/vtk-algorithm-cache-observers.md | 8 +++++ .../src/tack/interop/vtk/algorithms.py | 35 +++++++++++++++++-- .../tack-vis/tests/test_vtk_algorithms.py | 29 ++++++++++++++- 3 files changed, 69 insertions(+), 3 deletions(-) create mode 100644 changelog.d/vtk-algorithm-cache-observers.md diff --git a/changelog.d/vtk-algorithm-cache-observers.md b/changelog.d/vtk-algorithm-cache-observers.md new file mode 100644 index 00000000..9edde338 --- /dev/null +++ b/changelog.d/vtk-algorithm-cache-observers.md @@ -0,0 +1,8 @@ +### Fixed + +- Python crashed in its garbage collector when it freed, in a reference cycle (a + traceback's frames, say), a VTK dataset one of Tack's VTK algorithms had read: the + observer that drops the dataset's cache entry, referenced only by the dataset, was + cleared as garbage before the deletion called it (a VTK bug, in 9.7.1 and master). + The observers are kept reachable until they fire, and only note the deletion; the + cache drops the entry at its next use. diff --git a/packages/tack-vis/src/tack/interop/vtk/algorithms.py b/packages/tack-vis/src/tack/interop/vtk/algorithms.py index 98e13571..ed095af8 100644 --- a/packages/tack-vis/src/tack/interop/vtk/algorithms.py +++ b/packages/tack-vis/src/tack/interop/vtk/algorithms.py @@ -160,6 +160,32 @@ _datasets = {} # address -> ((modification time, time step), {device: DataSet}) +#: Addresses of VTK datasets deleted since the cache was last used. A DeleteEvent only +#: notes one, and its entry goes at the next use (``_drop_deleted``), outside whatever +#: freed the dataset. +_deleted = [] + +#: Each DeleteEvent observer, by its dataset's address, kept here until it has fired. +#: Referenced only by its VTK object, an observer is garbage when Python's collector +#: frees that object in a cycle (a traceback's frame holding it): the collector clears +#: the function's globals and closure, then the deletion calls it, and Python crashed. +_observers = {} + + +def _note_deleted(address): + """A DeleteEvent observer noting ``address`` deleted.""" + def note(caller, event): + _deleted.append(address) + return note + + +def _drop_deleted(): + """Drop the entries, and observers, of the VTK datasets deleted since the last use.""" + while _deleted: + address = _deleted.pop() + _datasets.pop(address, None) + _observers.pop(address, None) + def _key(vtk_data): from tack.interop.vtk import _time_key @@ -170,10 +196,14 @@ def _key(vtk_data): def _remember(vtk_data, data): """Keep ``data`` as the Tack dataset of ``vtk_data`` as it is now.""" + _drop_deleted() address = vtk_data.GetAddressAsString("vtkObject") if address not in _datasets: - # The address, not the caller: a DeleteEvent's caller is None. - vtk_data.AddObserver("DeleteEvent", lambda caller, event: _datasets.pop(address, None)) + # The address, not the caller: a DeleteEvent's caller is None. Only noted + # (see ``_deleted``); a later object at this address gets its own observer, + # since the noted address is dropped before it is remembered. + _observers[address] = observer = _note_deleted(address) + vtk_data.AddObserver("DeleteEvent", observer) _datasets[address] = (_key(vtk_data), {None: data}) @@ -197,6 +227,7 @@ def _as_dataset(vtk_data, device=None): data = vtk_to_dataset(vtk_data, dtype=None, copy=False) return data if device is None else data.to_device(device) return data + _drop_deleted() entry = _datasets.get(vtk_data.GetAddressAsString("vtkObject")) if entry is None or entry[0] != _key(vtk_data): _remember(vtk_data, vtk_to_dataset(vtk_data, dtype=None, copy=False)) diff --git a/packages/tack-vis/tests/test_vtk_algorithms.py b/packages/tack-vis/tests/test_vtk_algorithms.py index 2e490f5c..ca8fa417 100644 --- a/packages/tack-vis/tests/test_vtk_algorithms.py +++ b/packages/tack-vis/tests/test_vtk_algorithms.py @@ -250,7 +250,10 @@ def test_an_entry_goes_with_its_vtk_dataset(): assert address in ta._datasets del grid gc.collect() - assert address not in ta._datasets + # The deletion is noted, and the entry goes at the cache's next use. + assert address in ta._deleted + ta._drop_deleted() + assert address not in ta._datasets and address not in ta._observers def test_a_field_filters_output_keeps_its_inputs_topology(): @@ -1042,3 +1045,27 @@ def ends(data): furthest_halves = ends(on_halves)[:, 0].max() assert furthest_whole > 0.0 # the paths do cross x = 0 np.testing.assert_allclose(furthest_halves, furthest_whole, rtol=1e-6) + + +def test_a_dataset_the_collector_frees_leaves_its_observer_callable(): + """A dataset an algorithm read, freed by Python's collector in a cycle made after + it was read (a traceback's frames hold one so): its DeleteEvent observer, referenced + only by it, was cleared as garbage first, and calling it crashed Python (VTK bugs + 19). The observers are kept reachable, so this collects without crashing.""" + import gc + + class Holder: + pass + + addresses = [] + for _ in range(3): + data = _wavelet((-4, 4, -4, 4, -4, 4))() + assert ta.Contour(field="RTData", value=150.0)(data).number_of_cells > 0 + addresses.append(data.GetAddressAsString("vtkObject")) + holder = Holder() + holder.cycle = holder + holder.data = data + del data, holder + gc.collect() + ta._drop_deleted() + assert not any(address in ta._observers for address in addresses) From 9f80f4965ce9ec8331dcfda32aae24ad8979df86 Mon Sep 17 00:00:00 2001 From: Berk Geveci Date: Wed, 4 Nov 2026 15:44:30 -0500 Subject: [PATCH 2/3] VTK algorithms: Statistics, Entropy and Integrate, one answer for all the input The summaries reduce everything they are given -- every dataset of a partitioned dataset or a tree, on every MPI rank (the run's world) -- to one answer, as VTK's vtkDescriptiveStatistics and vtkIntegrateAttributes reduce composite data; every rank gets it. - Statistics: a vtkTable, a row per field and a vector field's components one by one (ExtractComponent), columns field, where, component, td.statistics' values under Tack's names, and the unit. - Entropy: the same rows, with the entropy in bits over `bins` bins. - Integrate: vtkIntegrateAttributes' output, one vertex at the measure-weighted centroid (the integral of position over the measure), each point field's integral as point data and each cell field's as cell data, the measure as Length, Area or Volume by the cells' dimension. - Fields: `fields`, SetInputArrayToProcess, or every field. Held to td.statistics and td.entropy called directly, Integrate to vtkIntegrateAttributes (centroid, volume, a point and a vector cell field), pieces and trees to the whole (each shared point counted once where GhostCells or GenerateGlobalIds mark its owner), a surface's measure to its area, and 3 MPI ranks to the whole. --- changelog.d/vtk-algorithms-summaries.md | 9 + docs/visualization/vtk.md | 21 +- .../src/tack/interop/vtk/algorithms.py | 258 ++++++++++++++++++ .../tack-vis/tests/_vtk_algorithms_mpi.py | 8 + .../tack-vis/tests/test_vtk_algorithms.py | 92 +++++++ 5 files changed, 385 insertions(+), 3 deletions(-) create mode 100644 changelog.d/vtk-algorithms-summaries.md diff --git a/changelog.d/vtk-algorithms-summaries.md b/changelog.d/vtk-algorithms-summaries.md new file mode 100644 index 00000000..7e33929d --- /dev/null +++ b/changelog.d/vtk-algorithms-summaries.md @@ -0,0 +1,9 @@ +### Added + +- `Statistics`, `Entropy` and `Integrate`, VTK algorithms of `td.statistics`, + `td.entropy` and `td.integrate` that reduce all their input (every dataset of a + partitioned dataset or a tree, on every MPI rank) to one answer, as VTK's + `vtkDescriptiveStatistics` and `vtkIntegrateAttributes` do. `Statistics` and + `Entropy` give a `vtkTable`, a row per field and a vector field's components one + by one; `Integrate` gives `vtkIntegrateAttributes`' one vertex at the centroid, + each field's integral on it and the measure as `Length`, `Area` or `Volume`. diff --git a/docs/visualization/vtk.md b/docs/visualization/vtk.md index ed978763..9086b97a 100644 --- a/docs/visualization/vtk.md +++ b/docs/visualization/vtk.md @@ -178,9 +178,9 @@ and what the algorithm supplies itself (the device, a locator, the communicator): `Contour(field="p", value=0.5, merge_points=False)`. Every algorithm of the [reference](../vis-guide/reference/algorithms.md) has -one, named as the Tack function in CamelCase, but for the summaries (to come) -and the finite-volume operators that take values on faces (`net_flux`, -`divergence`, `perot`), which VTK keeps nowhere. "Fields" means the input's +one, named as the Tack function in CamelCase, but for the finite-volume +operators that take values on faces (`net_flux`, `divergence`, `perot`), which +VTK keeps nowhere. "Fields" means the input's type, shallow-copied, with the new arrays added; "faces" means every face as a polygon (a 2D mesh's edges as lines), or every edge as a line, the new values their cell data. @@ -203,6 +203,8 @@ their cell data. | Fields from fields | `Gradients`, `FlowQuantities`, `Project`, `ValuesAtCenters`, `Interpolate`, `Magnitude`, `Dot`, `Cross`, `Composite`, `LogValues` | fields | | | `Convert` | the input's structure | | Finite volumes | `Jump`, `UpwindFlux` | faces | +| Summaries | `Statistics`, `Entropy` | `vtkTable`, a row per field (and component) | +| | `Integrate` | `vtkUnstructuredGrid` of one vertex, as `vtkIntegrateAttributes` | | Probing, sampling and particles | `Probe` (the source on port 1), `SphInterpolate` (the particles on port 1) | the input's structure, with the source's fields | | | `Voxelize` | `vtkImageData` (the `grid` parameter's) | | | `Streamlines`, `Advect` (seeds on port 1, or `seeds=`) | `vtkUnstructuredGrid` | @@ -286,6 +288,19 @@ empty, and the algorithm's `error` attribute keeps the exception for Python callers. `error` is `None` after a run that succeeded, so a script checks it after `Update()` or a pipeline call. +The summaries reduce all their input to one answer, as VTK's +`vtkDescriptiveStatistics` and `vtkIntegrateAttributes` do: every dataset of a +partitioned dataset or a tree, on every MPI rank, and every rank gets the +answer. A point that pieces share counts once where its copies' owner is +marked: run `GhostCells` or `GenerateGlobalIds` first (pieces cut without +either, such as `vtkExtractVOI`'s, count it in each). `Statistics` and +`Entropy` take `fields` (by default every field) and give a row per field, a +vector field's components one by one; the statistics keep Tack's names +(`kurtosis` is not the excess). `Integrate` puts its one point at the input's +centroid, each field's integral on it, and the measure as `Length`, `Area` or +`Volume`; integrals are by quadrature through each cell's map, so curved and +higher-order cells exactly. + ## Write your own VTK algorithm Subclass `TackAlgorithm` and give it `compute(data)`, a Tack dataset from a diff --git a/packages/tack-vis/src/tack/interop/vtk/algorithms.py b/packages/tack-vis/src/tack/interop/vtk/algorithms.py index ed095af8..c5c05e84 100644 --- a/packages/tack-vis/src/tack/interop/vtk/algorithms.py +++ b/packages/tack-vis/src/tack/interop/vtk/algorithms.py @@ -74,6 +74,7 @@ "Dot", "EdgeLengths", "Elevation", + "Entropy", "ExtractCells", "ExtractEdges", "ExtractExtent", @@ -91,6 +92,7 @@ "Gradients", "HaloCatalog", "ImplicitValues", + "Integrate", "Interpolate", "Jump", "LogValues", @@ -112,6 +114,7 @@ "SphInterpolate", "Spherical", "Split", + "Statistics", "Streamlines", "Subdivide", "SurfaceNormals", @@ -1741,4 +1744,259 @@ class GhostPoints(_OverPieces): FUNCTION = "ghost_points" +# ── Reductions: one answer for all the input ──────────────────────── + +class _Reduction(TackAlgorithm): + """A filter whose answer is one result for all its input: every dataset in it (a + partitioned dataset's pieces, a tree's datasets) on every rank, reduced together, + as VTK's ``vtkDescriptiveStatistics`` and ``vtkIntegrateAttributes`` reduce + composite data. Every rank gets the result. ``reduce(pieces)`` gives the VTK + output from a ``PartitionedDataSet`` of every dataset in the input.""" + + OUTPUT = "vtkTable" + + fields = Parameter(None, "the fields: names, or (name, 'points' / 'cells'); by " + "default every field. A vector field is reduced a " + "component at a time") + + def RequestDataObject(self, request, inInfo, outInfo): + import vtkmodules.vtkCommonDataModel as dm + + kind = getattr(dm, self.OUTPUT) + output = outInfo.GetInformationObject(0) + current = output.Get(vtkDataObject.DATA_OBJECT()) + if current is None or not current.IsA(self.OUTPUT): + output.Set(vtkDataObject.DATA_OBJECT(), kind()) + return 1 + + def RequestData(self, request, inInfo, outInfo): + from tack.dataset import PartitionedDataSet + + self._run_world = self._world() + input = self.GetInputData(inInfo, 0, 0) + if isinstance(input, (*_TREES, vtkPartitionedDataSet)): + leaves = [] + iterator = input.NewIterator() + iterator.InitTraversal() + while not iterator.IsDoneWithTraversal(): + leaves.append(iterator.GetCurrentDataObject()) + iterator.GoToNextItem() + else: + leaves = [input] + pieces = PartitionedDataSet([_as_dataset(leaf, self.device) for leaf in leaves], + world=self._run_world) + self.GetOutputData(outInfo, 0).ShallowCopy(self.reduce(pieces)) + return 1 + + def _selected(self, pieces): + """``(label, where, key, component)`` for each value to reduce: the chosen fields + (``fields``, or ``SetInputArrayToProcess``'s, or every field), a vector field's + components one by one (``component`` None for a scalar).""" + import tack.dataset as td + from tack.arrays import width_of + + first = next(iter(pieces), None) + if first is None: + return [] + chosen = self.fields + if chosen is None: + chosen = [a for a in map(self._input_array, range(self._arrays_set())) + if a is not None] or None + if chosen is None: + chosen = [(name, "points" if _on_points(field.space) else "cells") + for name, field in first.fields.entries()] + elif isinstance(chosen, (str, tuple)): + chosen = [chosen] + selected = [] + for value in chosen: + key = _field_argument(value) + field = first.fields[key] + name = key if isinstance(key, str) else key[0] + where = "points" if _on_points(field.space) else "cells" + width = width_of(field.values) + components = [None] if not width else list(range(width)) + selected += [(name, where, (name, td.H1 if where == "points" else td.Constant), c) + for c in components] + return selected + + def _arrays_set(self): + arrays = self.GetInformation().Get(vtkAlgorithm.INPUT_ARRAYS_TO_PROCESS()) + return 0 if arrays is None else arrays.GetNumberOfInformationObjects() + + +def _on_points(space): + from tack.dataset import H1, Values + + return isinstance(space, H1) or (isinstance(space, Values) and space.on == "points") + + +def _component_pieces(pieces, key, component): + """``pieces`` with the field ``key``'s ``component`` (a scalar field's self for None) + as the field ``__value__``.""" + from tack.arrays import ExtractComponent + from tack.dataset import DataSet, Field, PartitionedDataSet + + out = [] + for piece in pieces: + field = piece.fields[key] + values = field.values if component is None else ExtractComponent(field.values, + component) + out.append(DataSet(piece.topology, piece.geometry, + fields={"__value__": Field(field.space, values, field.unit)}, + ghosts=piece.ghosts, global_ids=piece.global_ids, + **piece._frame())) + return PartitionedDataSet(out, world=pieces.world) + + +def _table(columns): + """A ``vtkTable`` of ``columns``, ``{name: list}``: strings as a string column, + integers as ``i64``, the rest as doubles (``None`` as NaN).""" + from vtkmodules.vtkCommonCore import vtkDoubleArray, vtkStringArray, vtkTypeInt64Array + from vtkmodules.vtkCommonDataModel import vtkTable + + table = vtkTable() + for name, values in columns.items(): + if all(isinstance(v, str) for v in values) and values: + array = vtkStringArray() + for v in values: + array.InsertNextValue(v) + elif all(isinstance(v, int) and not isinstance(v, bool) for v in values) and values: + array = vtkTypeInt64Array() + for v in values: + array.InsertNextValue(v) + else: + array = vtkDoubleArray() + for v in values: + array.InsertNextValue(float("nan") if v is None else float(v)) + array.SetName(name) + table.AddColumn(array) + return table + + +class Statistics(_Reduction): + """``td.statistics`` of each field (each component of a vector field) over all the + input: a ``vtkTable``, a row per field and component, with columns ``field``, + ``where``, ``component`` (-1 for a scalar), the statistics as Tack names them -- + ``n``, ``min``, ``max``, ``sum``, ``mean``, ``m2``, ``m3``, ``m4``, + ``sample_variance``, ``population_variance``, ``sample_stddev``, + ``population_stddev``, ``skewness``, ``kurtosis`` (not the excess) and + ``weight`` -- and ``unit``, the field's (``""`` unknown). VTK's + ``vtkDescriptiveStatistics``.""" + + NAMES = ("n", "min", "max", "sum", "mean", "m2", "m3", "m4", "sample_variance", + "population_variance", "sample_stddev", "population_stddev", "skewness", + "kurtosis", "weight") + + weights = Parameter(None, "weights: 'volume' (each cell's measure, a point's share " + "of its cells'), or a field; by default none") + + def reduce(self, pieces): + import tack.dataset as td + + columns = {"field": [], "where": [], "component": [], + **{name: [] for name in self.NAMES}, "unit": []} + for name, where, key, component in self._selected(pieces): + stats = td.statistics(_component_pieces(pieces, key, component), "__value__", + weights=self.weights) + columns["field"].append(name) + columns["where"].append(where) + columns["component"].append(-1 if component is None else component) + for statistic in self.NAMES: + columns[statistic].append(getattr(stats, statistic)) + columns["unit"].append(stats.unit or "") + return _table(columns) + + +class Entropy(_Reduction): + """``td.entropy`` of each field (each component of a vector field) over all the input, + in bits, its values in ``bins`` equal bins over their range: a ``vtkTable`` with + columns ``field``, ``where``, ``component`` (-1 for a scalar) and ``entropy``. + Viskores' ``Entropy``.""" + + bins = Parameter(10, "the number of equal bins over the values' range") + weights = Parameter(None, "weights: 'volume', or a field; by default none") + + def reduce(self, pieces): + import tack.dataset as td + + columns = {"field": [], "where": [], "component": [], "entropy": []} + for name, where, key, component in self._selected(pieces): + columns["field"].append(name) + columns["where"].append(where) + columns["component"].append(-1 if component is None else component) + columns["entropy"].append(td.entropy(_component_pieces(pieces, key, component), + "__value__", bins=self.bins, + weights=self.weights)) + return _table(columns) + + +class Integrate(_Reduction): + """``td.integrate`` over all the input, as VTK's ``vtkIntegrateAttributes`` gives it: + a ``vtkUnstructuredGrid`` of one vertex at the input's centroid (the integral of + position over the measure), each point field's integral as point data and each + cell field's as cell data, and the measure as the cell array ``Length``, ``Area`` + or ``Volume``, by the cells' dimension. Integrals are by quadrature through each + cell's own map, so curved and higher-order cells exactly.""" + + OUTPUT = "vtkUnstructuredGrid" + + fields = Parameter(None, "the fields to integrate: names, or (name, 'points' / " + "'cells'); by default every field") + + def reduce(self, pieces): + import numpy as np + from vtkmodules.vtkCommonCore import vtkDoubleArray, vtkPoints + from vtkmodules.vtkCommonDataModel import vtkUnstructuredGrid + + import tack.dataset as td + from tack.interop.vtk import _label + + selected = {} + for name, where, key, _ in self._selected(pieces): + selected.setdefault((name, where), key) + first = next(iter(pieces), None) + origin = np.asarray(first.origin if first is not None else (0.0, 0.0, 0.0), float) + out = [] + for piece in pieces: + fields = {f"{where}:{name}": piece.fields[key] + for (name, where), key in selected.items()} + position = piece.geometry + if np.any(np.asarray(piece.origin, float) != origin): + # Pieces on other origins: their positions relative to the first's. + values = piece.world_positions() - origin + position = td.Field(piece.geometry.space, + tack.array(values.astype(np.float32 if values.dtype + == np.float32 else + np.float64))) + fields["__position__"] = position + out.append(td.DataSet(piece.topology, piece.geometry, fields=fields, + ghosts=piece.ghosts, global_ids=piece.global_ids, + **piece._frame())) + integrals = td.integrate(td.PartitionedDataSet(out, world=pieces.world), + fields=list(fields) if out else None) + measure = integrals.measure + centre = origin + (np.asarray(integrals.values["__position__"], float) / measure + if measure else 0.0) + grid = vtkUnstructuredGrid() + points = vtkPoints() + points.InsertNextPoint(*centre) + grid.SetPoints(points) + grid.InsertNextCell(1, 1, [0]) # VTK_VERTEX + for (name, where), _ in selected.items(): + value = np.atleast_1d(np.asarray(integrals.values[f"{where}:{name}"], float)) + array = vtkDoubleArray() + array.SetNumberOfComponents(value.size) + array.InsertNextTuple(value.tolist()) + array.SetName(name) + (grid.GetPointData() if where == "points" else grid.GetCellData()).AddArray(array) + _label(array, integrals.units.get(f"{where}:{name}")) + dimension = max((piece.topology.dimension for piece in pieces), default=3) + total = vtkDoubleArray() + total.InsertNextValue(measure) + total.SetName({1: "Length", 2: "Area"}.get(dimension, "Volume")) + _label(total, integrals.measure_unit) + grid.GetCellData().AddArray(total) + return grid + + del tack diff --git a/packages/tack-vis/tests/_vtk_algorithms_mpi.py b/packages/tack-vis/tests/_vtk_algorithms_mpi.py index 6029609d..fd73f6fb 100644 --- a/packages/tack-vis/tests/_vtk_algorithms_mpi.py +++ b/packages/tack-vis/tests/_vtk_algorithms_mpi.py @@ -61,6 +61,14 @@ def main(): assert owned == cells, (owned, cells) assert size == 1 or np.count_nonzero(marks), "no ghost cells from the neighbours" + # Reductions: every rank gets the whole's, each shared point counted once. + stats = (ta.GenerateGlobalIds(**options) >> ta.Statistics(fields=["RTData"], + **options))(pieces) + count = stats.GetColumnByName("n").GetValue(0) + assert count == points, (count, points) + volume = ta.Integrate(**options)(pieces).GetCellData().GetArray("Volume").GetValue(0) + assert volume == 16.0 ** 3, volume + merged = ta.MergePieces(**options)(pieces) if rank == 0: assert merged.number_of_points == points, (merged.number_of_points, points) diff --git a/packages/tack-vis/tests/test_vtk_algorithms.py b/packages/tack-vis/tests/test_vtk_algorithms.py index ca8fa417..edaecddf 100644 --- a/packages/tack-vis/tests/test_vtk_algorithms.py +++ b/packages/tack-vis/tests/test_vtk_algorithms.py @@ -1069,3 +1069,95 @@ class Holder: gc.collect() ta._drop_deleted() assert not any(address in ta._observers for address in addresses) + + +# ── Reductions: Statistics, Entropy, Integrate ────────────────────── + +def _row(table, row=0): + return {table.GetColumnName(i): table.GetValue(row, i).ToString() + for i in range(table.GetNumberOfColumns())} + + +def test_statistics_are_tacks_in_a_table(): + image = _wavelet()() + table = ta.Statistics(fields=["RTData"])(image) + assert isinstance(table, vtk.vtkTable) and table.GetNumberOfRows() == 1 + expected = td.statistics(vtk_to_dataset(image, dtype=None), "RTData") + row = _row(table) + assert (row["field"], row["where"], row["component"]) == ("RTData", "points", "-1") + for name in ta.Statistics.NAMES: + column = table.GetColumnByName(name) + assert column.GetValue(0) == pytest.approx(getattr(expected, name), rel=1e-12), name + + +def test_a_vector_field_is_reduced_a_component_at_a_time(): + image = (_wavelet() >> ta.Gradients(field="RTData", at="points"))() + table = ta.Statistics(fields=[("gradient", "points")])(image) + assert table.GetNumberOfRows() == 3 + assert [table.GetColumnByName("component").GetValue(i) for i in range(3)] == [0, 1, 2] + values = np.asarray(image.point_data["gradient"]) + np.testing.assert_allclose([table.GetColumnByName("mean").GetValue(i) for i in range(3)], + values.astype(np.float64).mean(axis=0), rtol=1e-6) + + +def test_every_field_by_default_and_by_set_input_array_to_process(): + image = (_wavelet() >> ta.Elevation())() + table = ta.Statistics()(image) + assert sorted(table.GetColumnByName("field").GetValue(i) + for i in range(table.GetNumberOfRows())) == ["RTData", "elevation"] + chosen = ta.Statistics() + chosen.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, + "elevation") + assert _row(chosen(image))["field"] == "elevation" + + +def test_entropy_is_tacks(): + image = _wavelet()() + table = ta.Entropy(fields=["RTData"], bins=16)(image) + expected = td.entropy(vtk_to_dataset(image, dtype=None), "RTData", bins=16) + assert table.GetColumnByName("entropy").GetValue(0) == pytest.approx(expected, rel=1e-12) + + +@pytest.mark.parametrize("marks", [ta.GenerateGlobalIds, ta.GhostCells]) +def test_a_reduction_of_pieces_and_trees_is_the_wholes(marks): + """The halves share a plane of points: with its copies' owner marked (global ids or + ghost cells), each is counted once, and the halves' statistics are the whole's.""" + whole, halves = _halves() + expected = ta.Statistics(fields=["RTData"])(whole) + tree = vtk.vtkPartitionedDataSetCollection() + tree.SetNumberOfPartitionedDataSets(1) + tree.SetPartitionedDataSet(0, marks()(halves)) + for table in [(marks() >> ta.Statistics(fields=["RTData"]))(halves), + ta.Statistics(fields=["RTData"])(tree)]: + assert table.GetColumnByName("n").GetValue(0) == expected.GetColumnByName( + "n").GetValue(0) + # Over pieces, sums of f32 values and moments combine in another order. + for name in ta.Statistics.NAMES: + assert table.GetColumnByName(name).GetValue(0) == pytest.approx( + expected.GetColumnByName(name).GetValue(0), rel=1e-6), name + + +def test_integrate_is_vtks_integrate_attributes(): + image = (_wavelet() >> ta.Gradients(field="RTData", name="g"))() + ours = ta.Integrate()(image) + theirs = vtk.vtkIntegrateAttributes()(image) + assert isinstance(ours, vtk.vtkUnstructuredGrid) and ours.number_of_points == 1 + np.testing.assert_allclose(ours.points.data.GetTuple3(0), theirs.points.data.GetTuple3(0), + atol=1e-9) + assert ours.cell_data["Volume"].GetValue(0) == pytest.approx( + theirs.cell_data["Volume"].GetValue(0), rel=1e-12) + assert ours.point_data["RTData"].GetValue(0) == pytest.approx( + theirs.point_data["RTData"].GetValue(0), rel=1e-6) + np.testing.assert_allclose(ours.cell_data["g"].GetTuple3(0), + theirs.cell_data["g"].GetTuple3(0), rtol=1e-6) + _, halves = _halves() + pieces = ta.Integrate()(halves) + assert pieces.cell_data["Volume"].GetValue(0) == pytest.approx( + theirs.cell_data["Volume"].GetValue(0), rel=1e-12) + + +def test_a_surfaces_measure_is_its_area(): + surface = (_wavelet() >> ta.ExtractSurface())() + out = ta.Integrate()(surface) + assert "Area" in out.cell_data.keys() + assert out.cell_data["Area"].GetValue(0) == pytest.approx(6 * 16 * 16, rel=1e-12) From c9f494e88c7fef9f01297432a930d25f1e9422bd Mon Sep 17 00:00:00 2001 From: Berk Geveci Date: Wed, 4 Nov 2026 16:08:24 -0500 Subject: [PATCH 3/3] VTK algorithms' dataset cache: an attribute of the VTK dataset's wrapper, no observer The cache of Tack datasets per VTK dataset kept entries by address and dropped them from a DeleteEvent observer, which needed a registry to keep the observers reachable (VTK bugs 19: the collector cleared one as garbage, then the deletion called it) and deferred drops. The entry is now an attribute of the dataset's Python wrapper: VTK keeps a wrapper's attributes while its C++ object lives, on every later wrapper of it (GetInputData's in the next filter), and frees them with the object, so an entry lives exactly as long as its dataset with no observer. Measured (yavin, the wheel and a VTK master build): on master a wrapper costs 0.6 us however many objects keep attributes (10,000), and Contour over 1,000 and 5,000 boxes runs as fast with the attribute as with the observer (951 against 943 ms, 4,803 against 4,817); VTK 9.7.1 sweeps kept attributes at every wrapper (50 us with 10,000), which it does for the attributes the pythonic API gives every wrapper anyway, and the two caches cost the same there too. An entry is stale when the dataset's modification time or time step differs from when it was cached, as before. The comment says what that sees: Modified() (an array's, the points'), arrays added or removed; not values written in place without it, nor InsertNextCell, which VTK's own pipeline would not run again for either. Tests: the entry outlives the wrapper and goes with the dataset, goes in the collector's cycle (the crash's case), and is stale after Modified() or a new time step but not after a write without it. --- changelog.d/vtk-algorithm-cache-observers.md | 11 ++- .../src/tack/interop/vtk/algorithms.py | 80 ++++++++----------- .../tack-vis/tests/test_vtk_algorithms.py | 51 ++++++++---- 3 files changed, 74 insertions(+), 68 deletions(-) diff --git a/changelog.d/vtk-algorithm-cache-observers.md b/changelog.d/vtk-algorithm-cache-observers.md index 9edde338..37e86f3a 100644 --- a/changelog.d/vtk-algorithm-cache-observers.md +++ b/changelog.d/vtk-algorithm-cache-observers.md @@ -2,7 +2,10 @@ - Python crashed in its garbage collector when it freed, in a reference cycle (a traceback's frames, say), a VTK dataset one of Tack's VTK algorithms had read: the - observer that drops the dataset's cache entry, referenced only by the dataset, was - cleared as garbage before the deletion called it (a VTK bug, in 9.7.1 and master). - The observers are kept reachable until they fire, and only note the deletion; the - cache drops the entry at its next use. + `DeleteEvent` observer that dropped the dataset's cache entry, referenced only by + the dataset, was cleared as garbage before the deletion called it (a VTK bug, in + 9.7.1 and master). The cache entry is now an attribute of the dataset's Python + wrapper, which VTK keeps while the dataset lives and frees with it: no observer. + An entry is converted again when the dataset's modification time or time step + changes; after changing a dataset in place, call `Modified()`, as VTK's own + filters need. diff --git a/packages/tack-vis/src/tack/interop/vtk/algorithms.py b/packages/tack-vis/src/tack/interop/vtk/algorithms.py index c5c05e84..80f77a14 100644 --- a/packages/tack-vis/src/tack/interop/vtk/algorithms.py +++ b/packages/tack-vis/src/tack/interop/vtk/algorithms.py @@ -154,43 +154,32 @@ # A Tack dataset keeps what it derives from its topology -- locators, faces, point # links, groups -- on the topology, so filters that read one VTK dataset (several # on one input, or one run again after a parameter changed) build that once when -# they share one Tack dataset. Kept by the VTK object's address, modification time -# and time step: the modification time is a global counter, so the same object -# modified, or a later one at the same address, does not match. An entry goes when -# its VTK object is deleted, and holds nothing that keeps that object alive. Beside -# the dataset sharing VTK's memory, an entry keeps its copy on each device a filter -# asked for (``device=``), so the data moves there once. - -_datasets = {} # address -> ((modification time, time step), {device: DataSet}) - -#: Addresses of VTK datasets deleted since the cache was last used. A DeleteEvent only -#: notes one, and its entry goes at the next use (``_drop_deleted``), outside whatever -#: freed the dataset. -_deleted = [] - -#: Each DeleteEvent observer, by its dataset's address, kept here until it has fired. -#: Referenced only by its VTK object, an observer is garbage when Python's collector -#: frees that object in a cycle (a traceback's frame holding it): the collector clears -#: the function's globals and closure, then the deletion calls it, and Python crashed. -_observers = {} - - -def _note_deleted(address): - """A DeleteEvent observer noting ``address`` deleted.""" - def note(caller, event): - _deleted.append(address) - return note - +# they share one Tack dataset. The entry is an attribute of the VTK dataset's +# Python wrapper (``_CACHE``): VTK keeps a wrapper's attributes while its C++ object +# lives, on every later wrapper of it (the one GetInputData gives in the next +# filter), and frees them with the object, so an entry lives exactly as long as its +# dataset, with no observer to call (a DeleteEvent observer crashed Python's +# collector, VTK bugs 19 in the vault). On VTK master a wrapper costs the same however +# many objects keep attributes (0.6 us); VTK 9.7.1 sweeps them at every wrapper (50 us +# with 10,000), as it does the attributes the pythonic API gives every wrapper. +# +# An entry is stale when the dataset's modification time or time step differs from +# what it was cached at (``_key``). The modification time follows VTK's contract: it +# sees a change made with ``Modified()`` (an array's or the points'), arrays added or +# removed, not values written in place without it, nor ``InsertNextCell``. Tack shares +# VTK's memory, so such values reach its fields, but what it derived from the old +# state (a locator, the faces) is kept: call ``Modified()`` after changing a dataset +# in place, as VTK's own filters need. Beside the dataset sharing VTK's memory, an +# entry keeps its copy on each device a filter asked for (``device=``), so the data +# moves there once. -def _drop_deleted(): - """Drop the entries, and observers, of the VTK datasets deleted since the last use.""" - while _deleted: - address = _deleted.pop() - _datasets.pop(address, None) - _observers.pop(address, None) +#: The attribute of a VTK dataset's wrapper that holds its entry: ``((modification +#: time, time step), {device: DataSet})``. +_CACHE = "_tack_cache" def _key(vtk_data): + """What an entry is valid for: the dataset's modification time and time step.""" from tack.interop.vtk import _time_key info = vtk_data.GetInformation() @@ -199,15 +188,13 @@ def _key(vtk_data): def _remember(vtk_data, data): """Keep ``data`` as the Tack dataset of ``vtk_data`` as it is now.""" - _drop_deleted() - address = vtk_data.GetAddressAsString("vtkObject") - if address not in _datasets: - # The address, not the caller: a DeleteEvent's caller is None. Only noted - # (see ``_deleted``); a later object at this address gets its own observer, - # since the noted address is dropped before it is remembered. - _observers[address] = observer = _note_deleted(address) - vtk_data.AddObserver("DeleteEvent", observer) - _datasets[address] = (_key(vtk_data), {None: data}) + setattr(vtk_data, _CACHE, (_key(vtk_data), {None: data})) + + +def _entry(vtk_data): + """``vtk_data``'s entry, or ``None`` (none, or a stale one).""" + entry = getattr(vtk_data, _CACHE, None) + return entry if entry is not None and entry[0] == _key(vtk_data) else None def _as_dataset(vtk_data, device=None): @@ -230,11 +217,10 @@ def _as_dataset(vtk_data, device=None): data = vtk_to_dataset(vtk_data, dtype=None, copy=False) return data if device is None else data.to_device(device) return data - _drop_deleted() - entry = _datasets.get(vtk_data.GetAddressAsString("vtkObject")) - if entry is None or entry[0] != _key(vtk_data): + entry = _entry(vtk_data) + if entry is None: _remember(vtk_data, vtk_to_dataset(vtk_data, dtype=None, copy=False)) - entry = _datasets[vtk_data.GetAddressAsString("vtkObject")] + entry = getattr(vtk_data, _CACHE) copies = entry[1] data = copies[None] if device is None or data.device == device: @@ -737,7 +723,7 @@ def _set_unit_output(self, output, result): output.CompositeShallowCopy(result) return output.ShallowCopy(result) - kept = _datasets.get(result.GetAddressAsString("vtkObject")) + kept = _entry(result) if kept is not None: # The output is a new object over the same arrays: the same Tack dataset. _remember(output, kept[1][None]) diff --git a/packages/tack-vis/tests/test_vtk_algorithms.py b/packages/tack-vis/tests/test_vtk_algorithms.py index edaecddf..5154a1ee 100644 --- a/packages/tack-vis/tests/test_vtk_algorithms.py +++ b/packages/tack-vis/tests/test_vtk_algorithms.py @@ -242,18 +242,35 @@ def test_a_changed_dataset_is_converted_again(): def test_an_entry_goes_with_its_vtk_dataset(): + """The entry is an attribute of the dataset's wrapper: it outlives the wrapper while + the dataset lives (a fresh wrapper has it) and goes with the dataset.""" import gc + import weakref - grid = _rotation() - ta._as_dataset(grid) - address = grid.GetAddressAsString("vtkObject") - assert address in ta._datasets - del grid + holder = vtk.vtkPartitionedDataSet() + holder.append(_rotation()) + data = weakref.ref(ta._as_dataset(holder.GetPartition(0))) + gc.collect() + assert data() is not None and ta._as_dataset(holder.GetPartition(0)) is data() + holder.SetPartition(0, None) gc.collect() - # The deletion is noted, and the entry goes at the cache's next use. - assert address in ta._deleted - ta._drop_deleted() - assert address not in ta._datasets and address not in ta._observers + assert data() is None + + +def test_an_entry_is_stale_after_modified_or_a_new_time(): + """``Modified()`` (here an array's) or another time step converts again; values + written in place without ``Modified()`` reach the fields, which share VTK's memory, + but the entry stays, as VTK's own pipeline would not run again either.""" + grid = _rotation() + data = ta._as_dataset(grid) + np.asarray(grid.GetPointData().GetArray("v"))[0] = 7.0 + assert ta._as_dataset(grid) is data + assert data.fields["v"].values.numpy()[0, 0] == 7.0 + grid.GetPointData().GetArray("v").Modified() + moved = ta._as_dataset(grid) + assert moved is not data + grid.GetInformation().Set(vtk.vtkDataObject.DATA_TIME_STEP(), 3.0) + assert ta._as_dataset(grid) is not moved def test_a_field_filters_output_keeps_its_inputs_topology(): @@ -1047,28 +1064,28 @@ def ends(data): np.testing.assert_allclose(furthest_halves, furthest_whole, rtol=1e-6) -def test_a_dataset_the_collector_frees_leaves_its_observer_callable(): +def test_a_dataset_the_collector_frees_takes_its_entry(): """A dataset an algorithm read, freed by Python's collector in a cycle made after - it was read (a traceback's frames hold one so): its DeleteEvent observer, referenced - only by it, was cleared as garbage first, and calling it crashed Python (VTK bugs - 19). The observers are kept reachable, so this collects without crashing.""" + it was read (as a traceback's frames hold one): a DeleteEvent observer, referenced + only by it, was cleared as garbage first and calling it crashed Python (VTK bugs + 19). The entry is the wrapper's attribute, no observer: it goes with the dataset.""" import gc + import weakref class Holder: pass - addresses = [] + entries = [] for _ in range(3): data = _wavelet((-4, 4, -4, 4, -4, 4))() assert ta.Contour(field="RTData", value=150.0)(data).number_of_cells > 0 - addresses.append(data.GetAddressAsString("vtkObject")) + entries.append(weakref.ref(ta._as_dataset(data))) holder = Holder() holder.cycle = holder holder.data = data del data, holder gc.collect() - ta._drop_deleted() - assert not any(address in ta._observers for address in addresses) + assert all(entry() is None for entry in entries) # ── Reductions: Statistics, Entropy, Integrate ──────────────────────