diff --git a/changelog.d/vtk-algorithm-cache-observers.md b/changelog.d/vtk-algorithm-cache-observers.md new file mode 100644 index 00000000..37e86f3a --- /dev/null +++ b/changelog.d/vtk-algorithm-cache-observers.md @@ -0,0 +1,11 @@ +### 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 + `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/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 98e13571..80f77a14 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", @@ -151,17 +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. +# 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. -_datasets = {} # address -> ((modification time, time step), {device: DataSet}) +#: 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() @@ -170,11 +188,13 @@ def _key(vtk_data): def _remember(vtk_data, data): """Keep ``data`` as the Tack dataset of ``vtk_data`` as it is now.""" - 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)) - _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): @@ -197,10 +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 - 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: @@ -703,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]) @@ -1710,4 +1730,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 2e490f5c..5154a1ee 100644 --- a/packages/tack-vis/tests/test_vtk_algorithms.py +++ b/packages/tack-vis/tests/test_vtk_algorithms.py @@ -242,15 +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() - assert address not in ta._datasets + 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(): @@ -1042,3 +1062,119 @@ 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_takes_its_entry(): + """A dataset an algorithm read, freed by Python's collector in a cycle made after + 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 + + 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 + entries.append(weakref.ref(ta._as_dataset(data))) + holder = Holder() + holder.cycle = holder + holder.data = data + del data, holder + gc.collect() + assert all(entry() is None for entry in entries) + + +# ── 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)