from alexlib development branch - #59
Merged
Merged
Conversation
- analysis.companion_indices: argsort/searchsorted instead of dense broadcast - interpolation.InverseDistanceWeighter.__call__: drop per-point copy and loop - io.iter_trajectories_ptvis: searchsorted-based row lookup - scene.iter_segments: intersect1d/in1d; drop unreachable return in trajectory_tags - trajectory: np.bool->bool; take_snapshot direct array indexing
- InverseDistanceWeighter.__call__: broadcast weights@data, drop per-point matched_data copy and Python loop (1.1). - neighb_dists: replace per-point loop with dists[use_parts].reshape (1.2). - rbf_interp: batched np.linalg.solve over the trajectory dimension for the boolean-mask path; index-array (scene) path keeps the original loop (1.4). - select_neighbs kept dense: test_compare_idw_rbf pins its exact np.argsort tie-breaking, so the planned KD-tree path is not applied (1.3 intentionally skipped). - fix typo 'No tracers im frame' -> 'in frame'.
- io.trajectories_table: one read filtered by time range, grouped by trajid.
- Scene.iter_trajectories: one read grouped by trajid, yielded in
self._trids order (same as before).
- Scene._iter_frame_arrays: one read grouped by time; yields one array per
frame in range (empty for missing frames, preserving iter_segments);
fix latent TypeError in cond handling ('&'.join(query_string, cond)).
- AnalysedScene._iter_frame_arrays: same batching; collect uses np.isin/
np.intersect1d instead of per-row set membership; fix Py3 it.izip.
- pairs.particle_pairs: KD-tree nearest neighbour (was dense argmin); dict lookups for trajectory membership; hoisted secondary start/end times into trajectories_in_frame (3.1, 3.3). - smoothing.savitzky_golay: replace per-component np.convolve loop with a single windowed matmul (sliding_window_view); drop deprecated np.mat (3.4). - tests/test_smoothing.py: polynomial exactness + old-vs-new allclose regression guard.
- Add pytest-benchmark dev dependency. - pyproject.toml: testpaths=['tests'] so benchmarks/ is excluded from the default run but runnable via 'pytest benchmarks/'. - benchmarks/bench_interp.py, bench_io.py, bench_pairs.py seeded with np.random.default_rng(42). - benchmarks/BASELINE.md records the post-optimization numbers.
- select_neighbs now dispatches: KD-tree path for num_neighbs mode when n > num_neighbs and no boundary tie exists; dense fallback otherwise. - Boundary-tie detection: recompute exact distances for the k=num_neighbs+1 candidates (O(m·k)) and check for tie at the selection boundary. If any row has a tie, fall back to dense O(m·n) to preserve exact np.argsort tie-breaking (required by test_compare_idw_rbf). - _select_neighbs_kdtree helper scatters exact distances into the dense (m,n) return format using a penalty-enabled argsort that respects the inf->0 round-trip for forbidden entries (self/companion). - _select_neighbs_dense: extracted existing dense body into a helper. - RBF's tracer_dists (tracer×tracer distance matrix) always uses the dense path because rbf_interp extracts kernel submatrices at arbitrary tracer-index pairs that are not guaranteed to be mutual neighbours in the tracer_tracer distance matrix. Benchmark (n=5000, m=2000, num_neighbs=4): test_idw_call 659 ms -> 140 ms (4.7x) test_neighb_dists 596 ms -> 37 ms (16x)
…t out-of-bounds error
… vtk export into flowtracks (v1.2.0)
…erage Demonstrates that shift_phase/phase_average already generalize to any periodic-flow dataset (LV, aorta, ...) whose acquisitions start at different points in the cycle: per-set shift_phase before phase_average recovers the true peak that naive averaging would smear.
Walks through the shift_phase + phase_average workflow from test_phase_align_synthetic_sets.py: peak selection (argmax), cyclic rolling of each set's cycle onto a reference phase (with a polar view of the roll), and naive vs. shift-aligned phase averaging.
No logic change: marimo reformatted multi-line return tuples/decorator args and dropped the unused bump binding from a cell's return tuple when the notebook was opened in an edit session.
For every test/script/notebook that reads data/particles.h5 or data/tracers.h5, add a companion path that converts the same data to Zarr and compares results end to end: - benchmarks/bench_io.py: zarr save/read benchmarks + a correctness check that Zarr round-trips to the same trajectories as HDF5. - examples/marimo_hdf5_vs_zarr_scene_analysis.py: companion to hdf5_scene_analysis.ipynb, comparing Scene(h5) vs zarr-round-tripped trajectories (positions/velocities/accel) — exact match. - examples/marimo_linking_trajectories_zarr.py: companion to linking_trajectories.ipynb, running the identical trajectory-welding algorithm on HDF5- and Zarr-backed input — identical welded output. - tests/test_zarr_pipeline_end_to_end.py: real (non-mocked) zarr trajectories feeding eulerian_grid directly through to a zarr output.
flowtracks/nhist.py imported matplotlib unconditionally at module level even though neither nhist() nor nhist_scipy() actually call into it (the plot= param is dead code outside the __main__ demo block) — that import was the only thing breaking `pytest` collection on GitHub Actions (ModuleNotFoundError: matplotlib), since CI's minimal install (pytest numpy scipy tables + pip install -e .) never pulled it in. - flowtracks/nhist.py: move the matplotlib import into the __main__ demo block, where it's actually used. - pyproject.toml: declare matplotlib as a real dependency anyway, since flowtracks/graphics.py does import it unconditionally at module level; also tighten requires-python to >=3.11 to match zarr>=3.0.0's actual floor (non-yanked zarr releases require Python>=3.11) and update classifiers accordingly. - setup.py: drop the stale, incomplete install_requires/python_requires/ classifiers/version that duplicated (and had drifted from) pyproject.toml's [project] table, which setuptools already treats as authoritative — one source of truth instead of a silent override. - .github/workflows/python-package.yml: drop Python 3.10 from the test matrix and bump the release-build Python to 3.11, matching the new requires-python floor.
…s target-id The correspondences-group fallback grouped 3D points across frames by row[3], the per-camera 2D-target array index for that frame. That index resets every frame and carries no persistent particle identity, so trajectories built from it stitched together unrelated particles frame to frame (multi-cm jumps well outside the tracker's own velocity bound). Add a linkage-based path that walks the tracker's own prev/next chain (openptv2's linkage/<name>/frame_NNNNN groups, same shape as ptv_is.#), mirroring iter_trajectories_ptvis' identity propagation. Verified against a real dataset: max inter-frame speed dropped from 87 mm/frame to 6.2 mm/frame, consistent with the tracker's configured velocity bound.
…lisions openptv2's linkage/<name>/frame_NNNNN prev arrays are not guaranteed collision-free: two particles in the same frame can claim the same predecessor (observed with the current default tracker, priority_segment_3d). The previous linkage-walk fix inherited trajid via plain fancy indexing, so a collision silently merged two distinct particles under one id - and once merged, that id kept re-colliding, snowballing into a single "trajectory" with several hundred points across a 100-frame run. Now any prev claim that is not unique within its frame breaks the chain (both claimants start new trajectories) instead of merging. Verified: no trajectory in a 100-frame run now exceeds 100 points, and per-frame speeds stay physically consistent (~0.1-0.3 mm/frame median).
…nment matrix stitch_trajectories compared every trajectory pair with per-pair Python- level np.linalg.norm calls, and solved a dense NxN Hungarian assignment each merge pass. Both are O(N^2) (the assignment solve worse), which made the function effectively unusable on a real tracking run's trajectory count (tens of thousands) - it never finished within any reasonable time. Bucket ends/starts by frame and evaluate one (end-frame, start-frame) bucket pair at a time as a single vectorised numpy op instead of one norm() call per candidate pair; solve the assignment only over the rows/ columns that actually have a candidate instead of a dense NxN matrix. Same candidate set, same costs, same matches - only how they get computed changes. Verified on a 21k-trajectory run: seconds instead of hanging indefinitely.
Adds flowtracks.zarr_scene.ZarrScene, implementing the same public surface
as scene.Scene (keys, shapes, trajectory_ids, trajectory_by_id,
iter_trajectories, iter_frames, frame_by_time, iter_segments, collect,
bounding_box, set_frame_range, frame_range, trajectory_tags) but reading a
Zarr store's trajectories/{pos,vel,time,trajid,accel} arrays instead of a
PyTables HDF5 /particles table. Duck-type compatible, not a subclass, since
the backing stores are fundamentally different -- existing callers that
already accept a Scene-like object (eulerian.py's `scene.collect(...)` fast
path) work unmodified against either.
Reads both flowtracks' own save_zarr_trajectories() output and openptv2's
RunStore.seal() output (which writes the same five arrays plus an optional
traj/{trajid,first,last,length} index -- the same triple as Scene's own
/bounds table) with no conversion step in either direction, since both sides
were designed to the same on-disk contract.
collect()'s `where` filtering is reimplemented directly in numpy instead of
composing a PyTables read_where() query string -- there's no query engine to
target, so this is strictly simpler than porting gen_query_string's string
DSL, not a re-implementation of it.
Adds open_scene(path, frame_range=None) to scene.py: dispatches to Scene or
ZarrScene via io.infer_format(), so callers don't need to know the format
up front.
Cross-validated: every ZarrScene method produces identical output to Scene
given the same trajectories written to both an HDF5 file and a Zarr store
(tests/test_zarr_scene.py, 13 tests). That cross-validation surfaced two
pre-existing bugs in Scene itself, fixed as part of this change:
- iter_segments() used np.in1d, removed in modern numpy; now np.isin.
- collect() with a `where` filter but no frame range set produced an
invalid PyTables expression ("& (...)", leading stray "&") because
self._frame_limit ("" when unset) was joined in unconditionally.
Full suite: 116 passed (103 existing + 13 new), 0 failures.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01MwY1Vo51jzAQDdv21mQaso
There was a problem hiding this comment.
Pull request overview
This PR brings in a large modernization/feature branch: adds Zarr-backed trajectory reading (including a ZarrScene compatible with the existing Scene API), expands the xarray-based post-analysis pipeline and ParaView/VTK export path, and introduces multiple vectorization/performance improvements with a broad new test/benchmark suite.
Changes:
- Add Zarr-native scene/trajectory support and end-to-end tests exercising Zarr → Eulerian grid → post-analysis → Zarr output.
- Introduce performance refactors (vectorized smoothing, KD-tree neighbor selection/pairing, batched RBF solves) with equivalence/benchmark coverage.
- Update packaging/CI (pyproject metadata, console entrypoint, version bump) and add new utilities/examples/docs.
Reviewed changes
Copilot reviewed 59 out of 59 changed files in this pull request and generated 4 comments.
Show a summary per file
| File | Description |
|---|---|
| tests/test_zarr_trajectories.py | New tests guarding linkage-vs-cache precedence and frame-gap behavior. |
| tests/test_zarr_scene.py | Cross-validates ZarrScene behavior vs Scene and openptv2 RunStore layout. |
| tests/test_zarr_pipeline_end_to_end.py | End-to-end Zarr pipeline test feeding eulerian_grid and reloading output. |
| tests/test_zarr_io.py | Zarr format inference, round-trip, and correspondence-reading tests. |
| tests/test_turbulent_statistics.py | Adds shift_fields idempotency test (and related helpers). |
| tests/test_streamlined_zarr_paraview.py | Comprehensive tests for CF metadata, Zarr/NetCDF equivalence, VTK export, streamlined pipeline. |
| tests/test_stitching.py | New tests for trajectory stitching logic and rejection conditions. |
| tests/test_speedup_equivalence.py | Equivalence tests for optimizations (RBF, Zarr I/O, tags, ptv_is parsing). |
| tests/test_smoothing.py | Regression tests validating vectorized Savitzky–Golay against legacy logic. |
| tests/test_sample_vtkcode.py | Tests for VTK writer utilities and numerical helpers; guards against divide warnings. |
| tests/test_post_analysis_xr.py | Tests for xarray post-analysis functions and recipe runs (NetCDF/Zarr). |
| tests/test_post_analysis_edge_cases.py | Edge-case tests for binning, phase shifting, stats weighting, export paths. |
| tests/test_phase_average_xr.py | Tests for phase averaging/fluctuations and recipe-driven NetCDF output. |
| tests/test_phase_align_synthetic_sets.py | Synthetic demonstration of shift-align then average preserving peaks. |
| tests/test_nhist.py | Tests comparing histogram implementations and empty-input behavior. |
| tests/test_eulerian_grid.py | Tests for Eulerian binning shapes/means and min_count masking. |
| tests/test_derived_fields.py | Validates derived fields vs legacy formulas; masks; VTK round-trip; NetCDF attrs. |
| tests/test_binning_equivalence.py | Validates vectorized binning reproduces legacy per-frame algorithm. |
| tests/helpers.py | Shared helpers for synthetic HDF5 inputs and fake scene objects. |
| tests/conftest.py | Adds shared pytest fixtures (repo/src paths, config/grid loaders). |
| setup.py | Moves most metadata to pyproject.toml; adds console script entrypoint; filters example files. |
| scripts/run_postptv_analysis.py | Adds a runnable analysis script wiring together conversion + post-analysis outputs. |
| requriments.txt | Updates dependency list (adds vtk). |
| pyproject.toml | Version bump + modern dependency set + pytest configuration. |
| peek_zarr.py | Adds a CLI-like utility to inspect Zarr trajectory stores. |
| future_plan.md | Adds roadmap/spec planning document for future development. |
| flowtracks/zarr_scene.py | Implements Zarr-backed Scene-compatible reader. |
| flowtracks/vtk_export.py | Adds/updates VTK export utilities and vectorized tensor metric computation. |
| flowtracks/trajectory.py | Fixes dtype usage and snapshot extraction logic. |
| flowtracks/stitching.py | Adds trajectory stitching implementation (gap bridging, assignment-based matching). |
| flowtracks/smoothing.py | Vectorizes Savitzky–Golay filtering and removes deprecated np.mat. |
| flowtracks/scene.py | Performance refactors for iteration/collection; adds open_scene dispatch. |
| flowtracks/pipeline.py | Adds a streamlined in-memory pipeline wrapper and keeps legacy stage wrappers. |
| flowtracks/phase_average.py | Introduces xarray-native phase averaging/fluctuations and a recipe runner. |
| flowtracks/pairs.py | Speeds up pairing via KD-tree and pre-indexing; hoists trajectory start/end times. |
| flowtracks/nhist.py | Adds scipy-backed histogram helper and refactors Matlab-like histogram implementation. |
| flowtracks/interpolation.py | Adds KD-tree neighbor selection path and batched RBF solve optimizations. |
| flowtracks/combine.py | Adds postptv-combine CLI entrypoint for combining realizations and exporting outputs. |
| flowtracks/analysis.py | Speeds up companion index mapping via sorting/searchsorted. |
| flowtracks/an_scene.py | Refactors analyzed-scene frame iteration and intersection logic for performance. |
| flowtracks/init.py | Updates package metadata and exports new public API surface. |
| examples/visualize_lv_results.py | Adds reader benchmarking + visualization example script. |
| examples/run_lv_pipeline.py | Adds a full example pipeline runner with logging and multiple outputs. |
| examples/marimo_zarr_dashboard.py | Adds an interactive marimo dashboard for Zarr exploration. |
| examples/marimo_phase_averaging_demo.py | Adds an interactive marimo demo for phase alignment + averaging. |
| examples/marimo_linking_trajectories_zarr.py | Adds an interactive marimo demo comparing linking on HDF5 vs Zarr. |
| examples/marimo_hdf5_vs_zarr_scene_analysis.py | Adds a marimo demo comparing scene outputs across backends. |
| examples/batch_Lagrangian_to_Eulerian.py | Adds/updates a thin wrapper preserving legacy HDF5 contract via new xarray binning. |
| benchmarks/bench_pairs.py | Adds pairing performance benchmark. |
| benchmarks/bench_io.py | Adds HDF5 vs Zarr I/O performance benchmarks. |
| benchmarks/bench_interp.py | Adds interpolation performance benchmarks. |
| benchmarks/BASELINE.md | Adds benchmark baseline documentation. |
| .github/workflows/python-publish.yml | Allows manual workflow dispatch for publishing. |
| .github/workflows/python-package.yml | Updates CI python versions and install steps; enables manual dispatch. |
Suppressed comments (1)
flowtracks/scene.py:253
_iter_frame_arrays()has the same empty-input bug asiter_trajectories(): whenread_where(read_cond)returns 0 rows,groupscontains an empty array andg['time'][0]raises IndexError. Filter out empty groups when buildingby_timeso frames with no particles still yield empty arrays as intended.
arr = self._table.read_where(read_cond)
order = np.argsort(arr['time'], kind='stable')
arr = arr[order]
bounds = np.flatnonzero(np.diff(arr['time'])) + 1
groups = np.split(arr, bounds)
by_time = {int(g['time'][0]): g for g in groups}
empty = arr[0:0]
💡 Add a code-review agent skill or configure MCP servers for context-aware, tailored reviews. Learn more in the docs.
Comment on lines
+205
to
+210
| order = np.argsort(arr['trajid'], kind='stable') | ||
| arr = arr[order] | ||
| bounds = np.flatnonzero(np.diff(arr['trajid'])) + 1 | ||
| groups = np.split(arr, bounds) | ||
| by_trid = {int(g['trajid'][0]): g for g in groups} | ||
| empty = arr[0:0] |
Comment on lines
+70
to
+76
| arr = self._table.read_where(read_cond) | ||
| order = np.argsort(arr['time'], kind='stable') | ||
| arr = arr[order] | ||
| bounds = np.flatnonzero(np.diff(arr['time'])) + 1 | ||
| groups = np.split(arr, bounds) | ||
| by_time = {int(g['time'][0]): g for g in groups} | ||
| empty = arr[0:0] |
Comment on lines
+54
to
+57
| for set_name, vel in (("s1", (1.0, 0.0, 0.0)), ("s2", (2.0, 0.0, 0.0))): | ||
| zarr_path = tmp_path / f"{set_name}_traj.zarr" | ||
| save_zarr_trajectories(_make_trajectories(vel, seed=hash(set_name) % 1000), zarr_path) | ||
| ds_sets[set_name] = eulerian_grid( |
Comment on lines
+20
to
+23
| @pytest.fixture(scope='session') | ||
| def config_dict(): | ||
| with open(SRC_DIR / 'config.yaml') as f: | ||
| return yaml.safe_load(f) |
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
No description provided.