diff --git a/design/README.md b/design/README.md index 90565b927..e52897a7a 100644 --- a/design/README.md +++ b/design/README.md @@ -57,3 +57,19 @@ and section 11 for integration and manual acceptance. The implementation and user documentation are complete on the unreleased phy 2.2 branch. Final large-dataset smoke testing and save/reopen validation remain. + +## Long-recording curation (proposed) + +The [long-recording curation Phy-side design contract](large-recording-curation.md) +describes bounded local datasets for ~1000-hour sorting datasets, the three spike-index +spaces, optional upstream waveform evidence, full-duration metrics, replay safety, +and the implementation stages. The accompanying synthetic 1000-hour +fixture and reference tests are in `phy/apps/tests/large_recording_fixtures.py` +and `phy/apps/tests/test_large_recording_contract.py`. + +This is a proposal and a reference test contract, **not yet** a `prepare-large` +implementation or an upstream export-format agreement. + +The [integration audit and plan](large-recording-integration-plan.md) refines +the contract with correctness findings, module boundaries, ordered commits in one PR and +acceptance gates. Follow its gates before enabling editable bounded sessions. diff --git a/design/large-recording-curation.md b/design/large-recording-curation.md new file mode 100644 index 000000000..26f6d482e --- /dev/null +++ b/design/large-recording-curation.md @@ -0,0 +1,348 @@ +# Long-recording curation in Phy 2.x: Phy-side design contract + +**Status:** implementation present; acceptance validation is ongoing. This +specification is for Phy-side behavior; it does not authorize changes to the +AIND/Code Ocean export pipeline. The current branch implements the bounded +workflow behind an explicit manifest while retaining standard behavior for +ordinary datasets. See the [integration plan](large-recording-integration-plan.md) +and [validation record](large-recording-validation.md) for implementation scope, +checks, and remaining acceptance work. + +The [integration audit and plan](large-recording-integration-plan.md) refines +this proposal with assignment revisions, transactional persistence, capability +gates and ordered commit acceptance criteria for one implementation PR. Its enablement gates take precedence +over the original implementation sequence in Section 9. + +## 1. Objective and boundaries + +Curate chronic Neuropixels sorting datasets of roughly 1000 hours without opening hundreds +of millions or billions of spikes in the existing in-memory Phy clustering model. +The workstation downloads the **complete sorter output once**; downloading the +complete raw voltage recording is **not required**. All preparation, interactive +curation, full-data replay, and regeneration run locally: + +```text +complete, immutable sorting (local disk) + | prepare-large: streamed/batch + v +bounded, temporally distributed Phy dataset + full-duration summaries + | existing Template GUI + Unit Stability + v +semantic curation actions (not only sampled row assignments) + | apply-curation: streamed/batch against complete sorting + v +new complete curated sorting (original untouched) + | prepare-large again, potentially different evidence windows + +-----------------------> new bounded GUI dataset +``` + +The **complete spike-level sorting stays outside the interactive GUI**. +The bounded dataset must reproduce useful current Phy workflows rather than +pretend to be a statistically representative replacement for every spike. +Full-duration plots and rates come from separate complete-sorting summaries. + +Out of scope: cloud services, rewriting phylib for out-of-core clustering, +longitudinal unit-identity graphs (phy3), automatic change-point detection, and +implementing upstream waveform extraction in this repository. + +## 2. The three index spaces (essential contract) + +All index values below are **zero-based row positions**, not timestamps, sample +numbers, SpikeInterface unit IDs, or stable IDs assigned by a sorter. + +| Symbol | Meaning | Where it lives | +| --- | --- | --- | +| `g` | Full-sorting spike row index, in immutable source order | Original `spike_times.npy` / `spike_clusters.npy` | +| `b` | Bounded spike row index used by the GUI | Bounded `spike_times.npy` / `spike_clusters.npy` | +| `w` | Row in the stored waveform/channel tensors | `_phy_spikes_subset.waveforms.npy` / `.channels.npy` | + +**Canonical mapping:** the bounded package writes +`_phy_spike_global_ids.npy`, a one-dimensional strictly increasing `int64` +array of length `N_bounded`. It is `G[b] = g`. The underscore prefix is +intentional: phylib treats ordinary `spike_*.npy` files as per-spike +attributes, so `spike_global_ids.npy` would be needlessly interpreted as one. + +**Dataset-local waveform IDs:** `_phy_spikes_subset.spikes.npy[w] = b` in the +bounded package. In the *upstream full export*, the same filename instead +contains `g` because that file accompanies the full sorting. The file **never** +means an intrinsic immutable/global ID; it indexes whichever Phy dataset is in +its directory. + +Example: + +```text +full indices selected by prepare-large, G: + [14, 19, 27, 117, 203, 839] +bounded rows b: + [ 0, 1, 2, 3, 4, 5] +upstream waveform full indices: + [14, 117, 839, 1902] +intersection with G: + [14, 117, 839] +bounded _phy_spikes_subset.spikes.npy: + [ 0, 3, 5] +``` + +Recover full source row numbers for waveforms by composition: + +```python +G = np.load('_phy_spike_global_ids.npy', mmap_mode='r') +B = np.load('_phy_spikes_subset.spikes.npy', mmap_mode='r') +full_waveform_rows = G[B] # [14, 117, 839] +``` + +`prepare-large` must preserve source row order, subset all per-spike arrays +consistently, retain absolute spike sample times, and remap matching waveform +rows **along with the corresponding channel/waveform data**. A sorted +intersection / `np.searchsorted` implementation needs no billion-element +inverse lookup array. Sorting by sample time is not a substitute for preserving +row identity; use the source export's exact row order (including ties). + +## 3. Source and optional upstream waveform inputs + +The source is a locally accessible standard Phy/SpikeInterface export (or an +adapter yielding equivalent arrays). The original spike index `g` refers to +**the exact exported row ordering**, not an assumed row in an unrelated +SpikeInterface spike vector. Source arrays and their provenance must be +validated before processing; never infer a row mapping from spike timestamps. + +The existing, unmodified Phy waveform files are: + +| Filename | Meaning | +| --- | --- | +| `_phy_spikes_subset.spikes.npy` | `int64[N_waveforms]`: source full spike row indices, sorted/unique | +| `_phy_spikes_subset.channels.npy` | `int32[N_waveforms, C_sparse]`: model channel indices, with `-1` padding; `channel_map.npy` maps these to raw-file columns | +| `_phy_spikes_subset.waveforms.npy` | numeric `[N_waveforms, N_samples, C_sparse]`: aligned waveform data, zero on padded channels | + +The upstream AIND pipeline **may later** supply the following proposed files; +these names/schema are *not yet finalized with its developers*: + +| Optional file | Initial proposal | +| --- | --- | +| `_phy_spikes_subset.window_bounds.npy` | `int64[W, 2]` recording sample intervals, `[start, end)` | +| `_phy_spikes_subset.window_types.npy` | `uint8[W]`: `0` survey, `1` rescue | +| `_phy_spikes_subset.window_strata.npy` | `int32[W]`: survey stratum index, `-1` for rescue | +| `_phy_spikes_subset.manifest.json` | Small versioned scalar metadata: source identity, time base, extractor geometry, parameters, seed, array semantics | + +There are deliberately **no structured NumPy dtypes**, mixed-purpose `N x 5` +matrices, or massive JSON window arrays. The upstream manifest is useful for +planning downstream windows but is not required for fixture-based development. +The presence of a window denotes *potential* waveform coverage, not one stored +waveform for every spike within that interval: per-unit waveform caps still +apply. The three established waveform files remain authoritative for which +individual spikes actually have stored waveforms. + +In the absence of the proposed window sidecars, `prepare-large` may use the +existing three files and their spike-index/time distribution as available +evidence, with reduced guarantees about raw-I/O locality and survey windows. +A dataset with neither raw traces nor pre-extracted waveforms must not claim +that Waveform/Feature View evidence is available. + +The full export's identifier/fingerprint should combine an immutable sorting +asset ID when available, full spike count, and a digest of a canonical +row-ordered immutable spike array (`spike_times.npy` is one candidate). The +upstream reservoir must be verified against that *exact* source. A mutable +`spike_clusters.npy` hash must not be the sole identity used across curation +iterations. Byte-level hashes in examples are illustrative; the definitive +AIND provenance method is a pending joint decision. + +Coordinate identity and assignment revision are separate: the latter binds the +initial full cluster assignments and metadata for replay. Each curated output +records its parent revision and journal digest. Validate relevant immutable +feature, acquisition and time-base provenance as well as spike rows; equal +timestamp arrays alone do not establish source identity. Production hashing +must stream file contents rather than read a complete array into memory. + +## 4. Bounded-package Phy-side output + +`prepare-large` creates a **new directory**, not an in-place rewrite of the +original. At minimum it contains ordinary, internally consistent Phy files +(`params.py`, `spike_times.npy`, `spike_templates.npy`, `spike_clusters.npy`, +channel/template arrays, and present per-spike attributes), plus: + +| File | Role | +| --- | --- | +| `_phy_spike_global_ids.npy` | Required `int64[N_bounded]` sorted mapping `b -> g`, no duplicates | +| `_phy_curation_manifest.json` | Versioned local contract; immutable source ID/fingerprint, complete/bounded counts, durations, preparation algorithm/parameters/seed, evidence limitations and file names | +| `_phy_spikes_subset.{spikes,channels,waveforms}.npy` | Existing Phy-format stored waveforms, *remapped to bounded row indices* (if available) | +| Full-duration metrics | Separate time-binned evidence, described in Section 6; do not pretend proxy metrics are full metrics | + +Unselected spike rows are absent from the GUI arrays, **not deleted from the +complete sorting**. Absolute `spike_times.npy` values remain in acquisition +sample coordinates and require a dtype capable of representing 1000 h at the +source sample rate (typically `int64`, not `int32`). For 30 kHz, 1000 h is +108 billion samples. Keep recording start/end, session/block boundaries, valid +acquisition exposure, and gaps distinct from periods of zero firing. + +Do not rely on phylib's fallback of estimating `model.duration` from the last +spike when raw traces are absent. The future GUI metrics/time-axis integration +must take the authoritative duration/coverage metadata instead. Standard +`params.py` still identifies the bounded Phy dataset; the manifest augments it. + +## 5. Preparation algorithm: scientific and memory constraints + +A first version should: + +1. Scan the **full** local sorting in a bounded-memory batch pass; calculate + full-duration per-unit time bins and coverage/exposure. +2. Partition *valid acquisition time* into coarse strata. Select short shared + contiguous survey windows at deterministic seeded offsets across sessions/ + blocks (illustrative defaults: one 30 s window per 1 h stratum). +3. Include **all units and all spikes within each chosen window** as initial + bounded candidates; preserve source indices and absolute times. +4. Enforce configurable **total spike / RAM budgets**. Add shared rescue windows + for underrepresented units, prioritizing benefit to several units and + temporally separated evidence; retain bounded deterministic selection. +5. Where a waveform reservoir exists, favor windows with usable waveform + coverage and report per-unit counts and uncovered regions. Selecting a + window does not guarantee all its spikes have stored waveforms. +6. Write an atomic, versioned bounded package. Verify array-length alignment, + sorted global IDs, waveform/feature alignment, provenance, coverage and + source immutability before marking it usable. + +These choices are **heuristics**, not a claim of unbiased sampling, complete +physiological representativeness, or equivalent CCG/ISI statistics. Never label +metrics from the bounded subset as full-data measurements. Windows identified +manually in a future Unit Stability View may motivate a regenerated package; +automated anomaly detection is out of scope. + +The upstream `20_000` waveform cap per unit is **not** a mandatory local budget: +for 1000 units x 20k waveforms x 90 samples x 32 channels x 4 bytes it would +approach **230 GB** uncompressed. Extraction, transfer, storage and bounded +GUI consumption require independently configurable *total* byte/spike budgets. + +## 6. Unit Stability evidence from the full recording + +A Unit Stability View should work on ordinary 1-10 h recordings too. For a +selected cluster it plots time-binned **full-sorting** spike count / firing rate, +amplitude median/quantiles, depth/position and presence/coverage when those +quantities exist in the source. Use an explicit time axis, valid acquisition +exposure per bin, data-availability flags, and recording gaps. Counts divided +by wall-clock bin width are wrong for acquisition gaps. No missing metric should +be silently synthesized. + +Initial storage may use regular homogeneous `.npy` arrays plus a small manifest, +with explicit row-to-cluster-ID mapping and bin edges; exact metric schema and +plot integration belong to the Unit Stability PR. At each full-data curation +iteration, merged/split clusters must receive **correctly recomputed or +incrementally updated** summaries, not merely the stale parent-unit curves. + +Current `BaseController.get_mean_firing_rate()` divides the loaded cluster's +`n_spikes` by `model.duration`; this would severely undercount a bounded unit. +When a bounded package is detected, its full-count/exposure-aware value must +replace that calculation wherever presented. The same provenance distinction +applies to the Cluster View, correlograms, ISI and other proxy-based statistics. + +Window statistics need explicit observation intervals. Do not connect ISIs or +CCG pairs across independent segments or unsupported gaps. CCG normalization +must use a declared lag-exposure estimator for those intervals; initially hide +unsupported normalized diagnostics. Whole-window sampling, rescue selection +and incomplete waveform coverage do not establish unbiased full-data quality +metrics or chronic biological unit identity. + +## 7. Semantic curation and full-data replay + +The current `Clustering` operates on arrays of loaded spike rows, creates new +cluster IDs during `assign/split`, and saves those local assignments. They +cannot by themselves define the fate of billions of unobserved spikes. + +Before mutations in a bounded GUI session, record **versioned semantic events** +with preconditions, affected source clusters, deterministic operation parameters, +resolved cluster IDs and provenance. Track edits through existing undo/redo and +save. Apply them in order against the *complete* local sorting to write a new +curated output (never overwrite the original). + +| Operation | Initial replay policy | +| --- | --- | +| Merge complete units | Reapply to every source spike in the units; replayable | +| Cluster quality/metadata labels | Replay per resulting cluster identity; replayable | +| Threshold split | Replay only if the predicate/feature definition is serialized and an equivalent full-data per-spike attribute exists; otherwise disallow | +| Lasso or arbitrary selected spike-ID split | Not replayable for unseen spikes by default; **block in bounded mode** unless a valid full-data rule is supplied | +| Undo/redo | Journal must represent the committed final event sequence and allow exact save/reopen, not duplicate reverted edits | + +Unbounded/ordinary Phy behavior must remain unchanged. It is not acceptable to +silently treat a selection of bounded spike IDs as a valid full-unit split. +Existing `curation.json` merge proposals and `curation_review.json` are separate +producer/review artifacts; the new semantic event journal must not overwrite +either file. Re-running `prepare-large` after `apply-curation` may reuse the same +source spike/waveform bytes because curation changes assignments, not spike +identity. + +For the first editable release, support merges and labels only. Keep all splits +blocked until their specific replay rules pass end-to-end tests. Allocate IDs +from the complete unit catalog, including units absent from the subset. Record +resolved inherited metadata, rather than recomputing it from full parent counts. +Commit journal, assignments and metadata as one validated save generation; +reopen must recover a coherent generation and allocator state. Validate operation +intent before mutation, including plugin and console entry points. + +## 8. Test/reference fixture and future acceptance + +`phy/apps/tests/large_recording_fixtures.py` deterministically writes a minimal +32-spike, **1000 h** full Phy-format sorting and a bounded eight-spike package, +including upstream waveforms for eight different source rows and remapped +bounded waveforms for five. The fixture intentionally omits the raw binary. +`phy/apps/tests/test_large_recording_contract.py` asserts: + +- exact `int64` sample indexing well past the 32-bit limit, source-order and + time preservation; +- full -> bounded -> waveform mapping, including copied waveform/channel rows; +- existing Phy waveform filenames, row/dtype/shape and `-1` padding invariants; +- simple window array dtypes, sample-coordinate bounds and shared row alignment; +- source fingerprint agreement, separate source-vs-bounded identity; +- empty/disjoint intersections and existing phylib load compatibility. + +**These are reference contract tests only:** they do not claim that +`prepare-large`, Unit Stability or semantic replay already exists. Subsequent +PRs should reuse the fixture but replace the reference remapping with tests of +the production command. Add mismatched-source, reordered-row, excessive-memory, +partially covered units, gap, sparse feature and replay/undo corruption tests +before enabling bounded curation on real datasets. + +## 9. Ownership and implementation order + +1. **This PR:** Phy-side design, source-coordinate/format reference fixture and + tests only. No GUI changes, no new CLI, no upstream exporter changes. +2. **Unit Stability:** full-data summaries and selected-unit GUI view, initially + usable on normal recordings. Keep user-visible full/partial provenance clear. +3. **Semantic curation journal:** integration with `Supervisor`, event semantics, + undo/save/reopen and blocking unsafe splits in bounded mode. +4. **`prepare-large`:** separate streaming batch CLI, input validation, window + selection, output/remapped waveform files and bounded-manifest generation. +5. **`apply-curation`:** validate/replay semantic actions against full local + sorting, write a new curated full sorting, verify, then regenerate a bounded + package. Test deterministic repeatability and failures atomically. + +The integration plan now organizes these stages as commits within a single PR, +rather than separate PRs. **Do not ship an editable bounded GUI without replay safety**. + +Agent starting points: `AGENTS.md`, `design/README.md`, +`phy/apps/__init__.py` (Click CLI), `phy/apps/template/gui.py` (controller), +`phy/apps/base.py` (views/rates/save), `phy/cluster/{clustering,supervisor}.py` +(actions/history), `phy/apps/template/tests/test_gui.py`, and +`phy/apps/tests/test_large_recording_contract.py`. The existing waveform +reader is in the separate `phylib/io/model.py` dependency. + +## 10. Decisions not frozen by this Phy-side PR + +The upstream team should review the exact optional window-sidecar names, AIND +sorting asset identity/fingerprinting, segment/time/recording-block semantics, +waveform extraction alignment, total storage budget, and whether windows or a +more general sampling-plan format fits its pipeline best. Phy should accept the +existing three-file waveform format regardless. Values such as `1 h`, `30 s`, +`2,000` and `20,000` are preliminary, not schema requirements. + +Code Ocean handles **creation of a reusable evidence reservoir**; Phy handles +**bounded selection, GUI coordinate remapping, full-duration summaries, +semantic curation and replay**. Neither side should infer an absent waveform +or silently redefine source spike IDs. + +## Implementation references + +- [Phy `Clustering.assign`](../phy/cluster/clustering.py) and + [Supervisor](../phy/cluster/supervisor.py) for current row-based edits. +- [Template GUI](../phy/apps/template/gui.py) and + [Waveform/Feature mixins](../phy/apps/base.py) for standard stored-waveform + support and full/proxy metric safeguards. +- [Existing stored-waveform GUI test](../phy/apps/template/tests/test_gui.py). +- [SpikeInterface `export_to_phy` source (AIND inspected commit)](https://github.com/SpikeInterface/spikeinterface/blob/896ea840da72b4a5f06dcf8870ad820d20fe7031/src/spikeinterface/exporters/to_phy.py). diff --git a/design/large-recording-integration-plan.md b/design/large-recording-integration-plan.md new file mode 100644 index 000000000..de9e7043d --- /dev/null +++ b/design/large-recording-integration-plan.md @@ -0,0 +1,283 @@ +# Long recording audit and integration plan + +## Assessment + +The downloaded contract bundle is a useful foundation. Keep the bounded GUI, +immutable full sorting, explicit source-row mapping, and full-data replay model. +Do not enable editing until replay and persistence agree with the GUI. A 1000-hour +timestamp fixture validates coordinate range; it does not validate capacity or +scientific representativeness. + +This review targets phy `7e439491` and the downloaded +`phy_large_recording_contract_bundle.zip`. The accompanying contract remains a +proposal. This plan refines its implementation order and supplies additional +requirements; the phases below govern enablement. + +## Findings that affect correctness + +1. **Identity needs two layers.** A spike-coordinate fingerprint identifies + immutable rows and their time/channel/feature provenance. An assignment + revision identifies the starting full cluster assignments and metadata. + Replay must validate both. A times-only hash cannot distinguish recordings + with identical timestamps or changes to feature values. Include the acquisition + identity, time base, relevant immutable array digests and schema in coordinate + identity; hash assignments and metadata separately for each revision. Record + parent revision, journal digest, output revision and preparation identity. + Compute hashes incrementally, without `read_bytes()` on production arrays. + +2. **Cluster IDs must come from the full namespace.** + [Clustering](../phy/cluster/clustering.py) initializes its allocator from the + maximum loaded cluster ID. A subset can omit the highest full ID. Seed the + allocator from a full unit catalog and persist its high-water mark. Record + explicit output IDs in events and reject collisions during replay. Units with + no bounded evidence remain in the catalog and output, with missing evidence + visible to the user. + +3. **Save is a transaction.** + [BaseController](../phy/apps/base.py) currently saves assignments and metadata + separately; [Supervisor.save](../phy/cluster/supervisor.py) emits further save + events and clears dirty state. Adding another journal file to that chain does + not make it crash consistent. Bounded mode needs immutable save generations + containing assignments, metadata, journal and identity, published by a single + commit pointer. Reopen only complete validated generations. Preserve the last + valid generation after interruption; keep dirty state after a failed commit. + Define whether undo survives reopen. Recommended v1: reopen establishes a new + undo baseline while retaining the committed semantic history for replay. + +4. **Record resolved metadata.** `Clustering._do_merge` selects the largest + parent using loaded counts, and Supervisor inherits its metadata. Full and + sampled counts can choose different parents. Serialize the resolved group and + custom labels in merge events, including an explicit policy for conflicting + labels. Replay must not infer this choice from full counts. + +5. **Capture intent before mutation.** + [AmplitudeView.on_request_split](../phy/cluster/views/amplitude.py) returns + selected row IDs and clears the threshold. A listener after clustering cannot + recover the predicate. Route bounded edits through a command service before + any mutation. Gate direct assign/split, plugins, console operations, shortcuts, + merge proposals and threshold providers at the mutation boundary. V1 supports + whole-unit merges and labels; block every split until its full-data semantics + are implemented and tested. + +6. **Statistics require observation intervals.** Preserve shared contiguous + windows and segment boundaries. ISIs must not connect windows or acquisition + gaps. CCGs must exclude pairs crossing unsupported boundaries, and their + normalization must match observed lag exposure; the existing full-duration + baseline in `BaseController._get_correlograms_rate` is not valid for a window + subset. Initially suppress unsupported normalization and diagnostic claims. + Label window statistics as sampled evidence. Rescue selection and waveform + availability introduce bias; they do not justify extrapolating quality, + contamination, refractory violations or chronic unit identity to the full run. + +7. **Define time and channel coordinates.** V1 accepts one validated monotonic + sample axis with a fixed sampling rate and explicit valid acquisition + intervals. Reject unsupported segment layouts rather than concatenate silently. + Preserve tied timestamps and source row order. Sparse waveform channel IDs + index model channels; `channel_map.npy` maps these to raw-file columns. They + are not interchangeable with hardware or SpikeInterface channel IDs. Test a + nonidentity map, sparse padding, whitening/scaling and waveform alignment. + Keep sample coordinates as integers and CPU times as float64; subtract the + viewport origin before converting fine-scale plotting positions to float32. + +8. **All preparation products need budgets.** A dense units-by-time-bin summary, + exact quantile calculation, waveform gather, or full index inversion can + exceed RAM despite streamed input. Specify batch, output, summary, cache and + GUI memory limits, plus temporary disk estimates. Distinguish peak working + memory from file-backed pages. Whole shared windows exceeding the spike budget + must be shortened, rejected or replaced with documented deterministic rules; + do not silently thin spikes within a window and retain its statistical claims. + +9. **Persist the final window plan.** Store retained intervals, segment/block + identity, survey/rescue classification and valid exposure in the bounded + package. A seed and parameters do not identify trimmed or replaced windows. + Count acquisition exposure once when composing merged-unit summaries. + +## Module boundaries + +Use a small `phy/apps/large_recording/` package, with pure services callable +without Qt and thin Click/controller adapters. Refine filenames while coding; +keep these responsibilities separate: + +| Responsibility | Contract | +| --- | --- | +| Source adapter and validation | Stream canonical arrays; classify per-spike, template, channel and cluster data explicitly; reject ambiguous unsupported arrays | +| Manifest and identity | Validate versions, paths, coordinate identity, assignment revisions, exposure and capabilities | +| Preparation | Deterministic window plan, bounded output, streamed waveform remapping, budget report | +| Summaries | Full counts/exposure, availability, declared exact or approximate estimators and multiresolution queries | +| Semantic commands and journal | Operation validation, explicit IDs/metadata, history cursor and transactional commit | +| Replay | Pure event interpretation, chunked output assignments, invariants and lineage | +| GUI adapter | Capability gating, evidence labels, bounded model/controller integration | + +Do not load the complete source through TemplateModel or Clustering in batch +commands. Do not hide full-data scans in a selected-unit GUI callback. Limit +changes to phylib to demonstrated loader/coordinate fixes, with separate tests +and dependency requirements. Preserve existing ordinary-dataset behavior. + +## One implementation PR with ordered commits and acceptance gates + +Implement the complete workflow on `feature/long-recording-curation` in the +main phy checkout. Use the stages below as reviewable commit groups within one +PR. Intermediate capability gates support development; merge the PR only after +the complete prepare, curate, save, replay and regenerate workflow passes its +acceptance checks. Predicate splits require their own full-data equivalence +checks within the same PR if included in its supported functionality. + +### 1. Contract and reference fixtures + +Start with the revised contract, this plan and format fixtures. +Fix the bundle's missing `dtype` in `params.py`, allow equal sample times, and +reject malformed reference row IDs. Run the real released-phylib loader test +without accepting a skip. Its success proves this small file layout loads, +not that all GUI workflows or long-recording statistics are supported. + +### 2. Source validation and summary engine + +Implement identities, assignment revisions, array inventory, exposure and +capability validation independently of Qt. Start with full counts and rates; +add amplitude/depth only with defined units and available source attributes. +Do not add per-bin quantiles until their exact/approximate algorithm, merge +behavior and memory costs are specified. Medians cannot be merged by averaging +parent medians; invalidate or recompute unavailable derived summaries. + +Gate: exact tiny-oracle counts; missing exposure; gaps; source/feature mismatch; +reordered rows; column-vector input arrays; unknown versions; unsupported +segments; bounded peak memory as source size grows; actionable disk estimates. + +### 3. Read-only preparation and inspection + +Implement `prepare-large` with deterministic shared windows, rescue policies, +per-unit evidence reports and explicit budgets. Write into a temporary sibling +directory and publish only after complete validation. Rebuild paths in params, +subset every registered per-spike array, retain template-indexed structures, and +avoid copying stale `.phy` caches and derived cluster metrics. Unknown attributes +must be explicitly included, excluded with a report, or rejected. + +Add a bounded-manifest GUI adapter with authoritative duration, full/proxy +counts and unsupported-view safeguards. Until phase 6, require read-only bounded +sessions, including save and plugin mutation entry points. + +Gate: production command vs independently specified expected arrays; shuffled +channel mapping; actual `get_waveforms()` values; sparse PC feature alignment; +missing/partial waveform coverage; tied samples; zero evidence units; empty +selection policy; deterministic budgets; interrupted publication; no source +writes; GUI open/close and time-axis precision tests. + +### 4. Semantic commands and durable saves + +Implement merge/label command validation, full-catalog ID allocation, resolved +metadata, undo/redo active history and transactional generations. Keep the +feature disabled for editable bounded sessions while replay is incomplete. +Existing `curation.json` and `curation_review.json` retain their separate roles. + +Gate: merge chains; conflicting labels; undo then branch; redo; reopen baseline; +omitted high-ID units; proposal-driven merges; rejected direct split/assign; +failed-save dirty state; fault injection at every commit boundary; assignments +equal the journal's projection onto bounded rows. + +### 5. Full-data replay + +Implement `apply-curation` against an exact source assignment revision. For v1 +merges, resolve the final source-unit mapping once, then stream assignments into +a new full dataset. Preserve row count/order, times and immutable arrays; write +metadata and lineage coherently. Reject repeated application to the wrong parent +revision. Publish only a verified complete output. Reprepare from that revision +with the same coordinate identity and fresh revision-dependent summaries. + +Gate: compare every assignment with an independent tiny oracle; preserve all +unseen spikes; replay chains and labels exactly; reject mismatched parent or +corrupted journal; interruption never exposes a complete-looking partial output; +same source/events produce the same semantic result; source bytes unchanged. + +### 6. Enable bounded merges and labels with Unit Stability + +Expose editable bounded mode only when phases 3–5 pass together. Integrate a +Unit Stability view that distinguishes acquisition exposure, evidence windows, +available measurements and unresolved identity. Add bounded status and operation +capabilities to the existing workflow. Count/exposure summaries can update for +merges; quantiles and other nonadditive metrics require recomputation or explicit +unavailability. Keep all splits blocked. + +Gate: prepare → inspect → merge/label → undo/redo → save/reopen → replay → +reprepare → verify. Include ordinary 1–10 h datasets to prevent regressions. +Run `make test-full`, `make doc-check`, and `uv build` when CLI/package changes +are included. Complete manual smoke testing with an actual chronic dataset. + +### 7. Optional predicate splits + +Add one operation at a time after the initial workflow is accepted. Store feature +identity/units, predicate operator, threshold, missing-value policy, parent +revision and explicit child IDs. The same predicate must evaluate identically +on GUI evidence and all source spikes. Stored waveform-derived or local PCA +amplitudes do not automatically define a full-data predicate. Lasso splits +remain unavailable without a separately validated full-data rule. + +## Standard recording compatibility gate + +Recording length must never select a new behavior implicitly. Activate bounded +curation only through an explicit validated bounded manifest. Without that +manifest, retain the standard loader, selector, metrics, mutation capabilities, +save format and history behavior. A malformed bounded manifest must produce an +actionable error rather than silently enter standard editable mode. Optional +Unit Stability on ordinary recordings must not replace their existing views, +change their default layout or introduce full-source scans at GUI startup. + +Keep the semantic command, replay and transactional save implementation behind +the bounded adapter wherever possible. Any necessary shared-code change requires +focused standard-mode regression coverage in addition to bounded-mode tests. +Do not impose bounded split restrictions or journal files on ordinary datasets. + +Before merging the single PR, run `make test-full` and verify ordinary datasets +with raw traces, stored waveforms only, and sparse features. Cover standard +merge and label inheritance, lasso/threshold/template splits, undo/redo and +branching, save/reopen assignments and metadata, merge proposals, plugin-driven +edits, channel mapping, sampled views and existing rate/ISI/CCG behavior. Assert +that ordinary saves do not create bounded journals or manifests. Compare startup +memory and selection latency on the same ordinary datasets before and after the +change; any meaningful regression must be explained and resolved before merge. +Existing tests and these checks establish compatibility evidence; the plan +itself does not guarantee compatibility before implementation. + +## Capacity and scientific acceptance + +Use three kinds of evidence: tiny exact-oracle datasets for correctness; +generated high-density on-disk datasets for resource behavior; real chronic +recordings for usability and scientific acceptance. Duration, spike count, +number of units, array widths and available waveform coverage are separate axes. +Measure preparation/replay throughput, peak process memory, temporary disk, +GUI startup and selected-unit update latency on documented hardware. Establish +acceptance budgets from those measurements before advertising capacity. + +Have a domain curator review known drift, rare/bursting units, transient +artifacts, gaps and intentionally incorrect merges. Compare selected evidence +with withheld contiguous windows and full-data diagnostics where available. +Record what the curator can and cannot infer; finite selected evidence cannot +guarantee biological identity throughout 1000 hours. Software integrity and +scientific evidence sufficiency are separate release criteria. + +## Producer decisions + +Agree with the producer team on source asset identity, full unit IDs, authoritative +duration/exposure, segment mapping, waveform time/channel/amplitude semantics, +reservoir coverage and total storage costs. Existing three-file waveforms remain +supported; proposed window sidecars must remain optional until jointly agreed. +The inspected [SpikeInterface exporter revision](https://github.com/SpikeInterface/spikeinterface/blob/896ea840da72b4a5f06dcf8870ad820d20fe7031/src/spikeinterface/exporters/to_phy.py) +supports one segment and distinguishes Phy channel indices from exported original +channel IDs. Do not assume that arbitrary multi-session inputs already satisfy +this contract or that a deployed exporter matches this inspected revision. + +## Implementation and validation status + +The branch now includes batched source validation and preparation, full-source +summaries, bounded GUI curation and checkpoints, semantic replay, Unit Stability, +and dedicated long-recording tests. The feature is implemented, but its +acceptance remains incomplete: this test and synthetic benchmark evidence does +not replace real chronic-recording review by domain curators. + +On 2026-10-13, `make test-full` completed with 745 passed and 1 skipped. Three +additional bounded GUI lifecycle regressions pass in +`phy/apps/large_recording/tests/test_gui_lifecycle.py`. +[`design/large-recording-validation.md`](large-recording-validation.md) records +the synthetic capacity runs and their limits. The separate 32-spike, 1000-hour +fixture verifies sample-coordinate handling and file mapping; it is not a +billion-spike capacity test or evidence of chronic biological identity. Curator +acceptance on real chronic recordings remains outstanding. diff --git a/design/large-recording-validation.md b/design/large-recording-validation.md new file mode 100644 index 000000000..2dcd3e761 --- /dev/null +++ b/design/large-recording-validation.md @@ -0,0 +1,77 @@ +# Long-recording implementation validation + +This record summarizes the current implementation checks and their limits. The +results validate the measured code paths and fixtures described here. They do +not establish capacity for every source layout or scientific suitability for +chronic curation. + +## Automated checks + +On 2026-10-13, `make test-full` completed with **756 passed and 1 +skipped** (496 non-application tests and 260 application tests). Final channel +geometry and nonidentity waveform-mapping additions passed a further focused +source/preparation run of **33 tests**. Lint, format checking, package build, +and the strict documentation build also passed. The dedicated bounded GUI lifecycle file passed all three +regressions: reentrant observer rollback with the GUI attached, merge-proposal +review through undo/redo/save/reopen, and pending full metrics for threshold +children in the default bounded views. + +The reference long-recording fixture is deliberately sparse: 32 spikes spread +over a 1000-hour sample clock. It checks long `int64` coordinates, tied times, +source-row mapping, waveform remapping, and loader compatibility. It is not a +billion-spike capacity test, a realistic chronic recording, or evidence of +biological identity across a recording. + +## Synthetic preparation and replay runs + +The reproducible benchmark is `uv run python +tools/benchmark_long_recording.py`. It writes synthetic source arrays in +batches, with no raw voltage data, and runs preparation and replay in separate +processes. Both runs below used Python 3.13.3, NumPy 2.4.3, Linux +`x86_64`/glibc 2.39, 30 kHz sampling, 1000 units, a 1000-hour sample axis, +100,000-row batches, 30-second windows, 3600-second strata, a 100,000-spike +budget, a 512 MiB output budget, seed 0, and 3600-second summary bins. The +oracle independently streamed every source and output assignment row and +matched the expected merge in both cases. + +| Synthetic source spikes | Source files | Bounded package | Bounded spikes | Curated output | Prepare time / peak RSS | Replay time / peak RSS | Oracle rows | +| ---: | ---: | ---: | ---: | ---: | ---: | ---: | ---: | +| 250,000 | 4,032,892 B | 8,625,460 B | 2,389 | 5,038,086 B | 1.902 s / 215,773,184 B | 0.114 s / 214,847,488 B | 250,000 | +| 2,000,000 | 32,032,892 B | 9,157,120 B | 16,676 | 40,038,086 B | 1.227 s / 214,433,792 B | 0.292 s / 215,285,760 B | 2,000,000 | + +Each measurement is a single run, not a throughput curve. Preparation and +replay durations exclude child-process startup; peak RSS is the process high +water mark. The source is a regular synthetic distribution with tiny templates +and no stored waveforms, raw data, irregular acquisition gaps, or real feature +loads. These figures therefore do not predict every preparation cost or GUI +interaction on a chronic dataset. The full JSON reports are retained at +`/tmp/phy-long-recording-benchmark-250k.json` and +`/tmp/phy-long-recording-benchmark-2m.json` on the validation workstation. + +## Ordinary-dataset smoke comparison + +A separate three-run comparison used fresh subprocesses, configs, and caches in +the same Python environment. It compared baseline `a30d5c65` with current +`f212f069e5b5b19469cfa73781dac08de093d9a1` on a 12,000-spike, two-cluster +dataset with 1,024 stored waveform rows of shape `[40, 4]` and dense template +features. Counts and output dimensions matched exactly. + +| Measure (median of 3) | Baseline | Current | +| --- | ---: | ---: | +| Controller startup | 0.03879 s | 0.02822 s | +| Waveform provider | 0.00353 s | 0.00293 s | +| Feature provider | 0.00372 s | 0.00302 s | + +This small case showed no regression; the apparent gains are within the range +where run-to-run noise matters. The comparison did not measure Qt rendering, +raw-voltage I/O, or large-dataset selection churn. Its report is retained at +`/tmp/phy-standard-recording-comparison.json` on the validation workstation. + +## Remaining acceptance work + +Real chronic-recording review by domain curators remains outstanding. It must +cover drift, rare and bursting units, artifacts, gaps, waveform coverage, +intentionally incorrect merges, withheld contiguous windows, and full-data +diagnostics where available. Software integrity checks and synthetic throughput +do not establish that the evidence is scientifically sufficient or that a +cluster represents one biological unit throughout 1000 hours. diff --git a/docs/api.md b/docs/api.md index a5f0f2506..52733d154 100644 --- a/docs/api.md +++ b/docs/api.md @@ -10117,6 +10117,13 @@ Component that brings manual clustering facilities to a GUI: Whether automatic navigation and similar-cluster selection skip noise and MUA clusters. The default is True. +* `clustering_factory : callable, optional` + Construct the clustering model. Defaults to Clustering; bounded curation + supplies a guarded implementation with the same constructor interface. + +* `cluster_meta_factory : callable, optional` + Construct metadata from groups and labels. Defaults to create_cluster_meta. + **Events** When this component is attached to a GUI, the following events are emitted: @@ -13197,6 +13204,15 @@ Create a trace view. --- +#### TemplateController.create_unit_stability_view + + +**`TemplateController.create_unit_stability_view(self)`** + +Optionally inspect selected-unit activity over the loaded recording. + +--- + #### TemplateController.create_waveform_view diff --git a/docs/changelog.md b/docs/changelog.md index 2842a69f3..f18f7b05f 100644 --- a/docs/changelog.md +++ b/docs/changelog.md @@ -11,10 +11,27 @@ been included in a stable release. The current entries cover all user-visible changes committed since 23 July 2026; test-only commits are represented by the behavior they verify rather than listed separately. +### Long recordings + +- Added `prepare-large` and `apply-curation` for bounded shared-window evidence, + full-duration exposure-aware summaries, semantic merges/labels/source-amplitude + splits, and replay into a new complete sorting. +- Added Unit Stability, precise window-relative amplitude/feature inspection, + sampled interval/correlogram views, and explicit full-versus-evidence counts. + Standard datasets retain their existing curation workflow and default layout. +- Bounded saves use validated atomic generations and reject conflicting saves + from stale sessions. Unsupported selected-row splits are blocked. +- Sparse cluster IDs no longer require dense allocations during clustering. +- Bounded amplitude plots label seconds from window start and display the + window's recording-time range above the plot. + ### Documentation - Documented uv-first installation of exact previous phy releases in separate environments, including how to intentionally replace a tool installation. +- Added a long-recording workflow guide covering bounded evidence preparation, + supported curation, full-sorting replay, and limits on interpreting sampled + measurements. ### Added diff --git a/docs/cli.md b/docs/cli.md index 0bb61a764..16b554a1f 100644 --- a/docs/cli.md +++ b/docs/cli.md @@ -25,6 +25,35 @@ GUI commands accept: --clear-cache remove and rebuild the dataset .phy cache ``` +## Long recordings + +Prepare bounded evidence, curate it with `template-gui`, then replay supported +actions into a new full sorting. See the [long-recording guide](long-recordings.md) +for input requirements, interpretation limits, and the complete workflow. + +```bash +phy prepare-large /data/recording/params.py /data/recording-bounded \ + --duration 3600000 \ + --valid-intervals /data/recording/valid-intervals.npy \ + --max-spikes 2000000 \ + --max-bytes 536870912 \ + --window-seconds 30 \ + --stratum-seconds 3600 \ + --seed 0 \ + --bin-seconds 3600 \ + --batch-size 100000 + +phy template-gui /data/recording-bounded/params.py + +phy apply-curation /data/recording/params.py /data/recording-bounded \ + /data/recording-curated --batch-size 100000 +``` + +`--duration` is acquisition duration in seconds. `--valid-intervals` reads an +integer `.npy` array of `[start, end)` sample pairs. `--max-bytes` caps the +prepared package size, not process memory. The output directories must be new; +the original sorting remains the replay parent. + ## Raw trace viewer ```bash diff --git a/docs/images/long-recording-evidence.png b/docs/images/long-recording-evidence.png new file mode 100644 index 000000000..294d03924 Binary files /dev/null and b/docs/images/long-recording-evidence.png differ diff --git a/docs/images/long-recording-overview.png b/docs/images/long-recording-overview.png new file mode 100644 index 000000000..cfd729bf7 Binary files /dev/null and b/docs/images/long-recording-overview.png differ diff --git a/docs/long-recordings.md b/docs/long-recordings.md new file mode 100644 index 000000000..b5bf8d424 --- /dev/null +++ b/docs/long-recordings.md @@ -0,0 +1,183 @@ +# Curating long recordings + +Phy's long-recording workflow prepares a bounded evidence dataset for interactive +review, then replays supported curation actions against the complete local +sorting. Preparation reads source arrays in batches; it does not load the full +sorting into the GUI. The source dataset is retained unchanged, and replay +writes a separate output dataset. + +This workflow is intended for a Phy sorting with one fixed-rate, monotonic sample +clock. Multi-segment recordings are not supported. Supply acquisition duration +and valid acquisition intervals when they are not already recorded in the +dataset metadata. Duration describes exposure; Phy cannot infer acquisition +gaps from periods with no spikes. Intervals use half-open sample coordinates +`[start, end)`. + +The screenshots below show the actual bounded GUI views using a deterministic +synthetic eight-hour sorting. The overview plots full-duration rates and +available source means; gray marks the synthetic acquisition gap and cyan marks +the retained evidence windows. The evidence view places windowed amplitudes at +top left, sampled ISI at top right, sampled correlograms at middle right, and +features below. They illustrate the interface and its sampling context; they +are not biological data or scientific examples. + +![Full-duration Unit Stability with a synthetic acquisition gap](images/long-recording-overview.png) + +![Windowed source amplitudes, features, sampled ISI and correlogram](images/long-recording-evidence.png) + +Recreate these images from source with Qt's X11 platform and software OpenGL +(the script writes a temporary synthetic source and bounded package): + +```bash +QT_QPA_PLATFORM=xcb LIBGL_ALWAYS_SOFTWARE=1 \ + xvfb-run -a -s '-screen 0 2400x1600x24 +extension GLX' \ + uv run python tools/capture_long_recording_views.py +``` + +## Prepare bounded evidence + +Run `prepare-large` with the source `params.py` and a new output directory: + +```bash +phy prepare-large /data/recording/params.py /data/recording-bounded \ + --duration 3600000 \ + --valid-intervals /data/recording/valid-intervals.npy \ + --max-spikes 2000000 \ + --max-bytes 536870912 \ + --window-seconds 30 \ + --stratum-seconds 3600 \ + --seed 0 \ + --bin-seconds 3600 \ + --batch-size 100000 +``` + +`--duration` is in seconds and must agree with declared acquisition metadata. +The `.npy` valid-interval file contains an integer array of shape `[W, 2]` with +ordered, nonoverlapping sample pairs. +Without it, the full declared duration is treated as valid exposure. Do not use +that default when the recording contains acquisition gaps. + +`--max-bytes` limits the bounded output package size, including its summaries and +window plan. It is not a process-memory or peak-RSS limit. The source unit IDs +must fit nonnegative `int32`; sample coordinates and full-source spike row IDs +use `int64`. + +Preparation chooses shared time windows under the spike and byte budgets. The +window length, strata, and seed make the plan reproducible; inspect the saved +window plan in the output package when documenting or repeating an analysis. +Windows that exceed the budget are omitted, so a high-density interval can +reduce retained evidence. Preparation does not shorten an over-budget window or +silently thin spikes inside it. A survey of shared windows is supplemented by +deterministic rescue evidence for units missed by the survey. This improves the +chance of seeing rare units, but makes the evidence selection unequal: rescue +windows must not be treated as a probability sample. + +Preparation records exact full-recording spike counts and rates using valid +acquisition exposure. If source scalar amplitudes or depths are available, it +also records per-bin finite-value means and availability counts. It does not +invent missing measurements or convert source values to physical units. Stored +waveforms are included only where available in the source; raw traces are not +required when stored waveform data are sufficient for the views being used. +Producer-supplied per-unit metrics and annotations remain inherited metadata; +Phy does not recalculate those values for merged or split units. Use the +full-count, rate, amplitude-mean, and depth-mean summaries produced by this +workflow for the measurements they explicitly cover. + +## Review and curate + +Open the bounded output with the normal GUI command: + +```bash +phy template-gui /data/recording-bounded/params.py +``` + +The bounded manifest explicitly selects this mode. Bounded mode has its own +configuration namespace and default view layout; it does not replace the normal +Phy layout or restore ordinary saved dock geometry. The GUI keeps full-recording +counts and rates distinct from counts in selected evidence rows. Unit Stability +shows time-binned counts/rates, available source scalar-amplitude means, +evidence windows, and acquisition gaps. Missing measurements and units without +evidence remain visible as unavailable. + +Evidence Amplitude shows raw source scalar amplitudes from the currently +selected evidence window. The Feature view also uses window-relative seconds. +The amplitude plot labels this axis as **Seconds from window start** and shows +the window's recording-time range above the plot. +Both views label the current window and reset their time origin when it changes. +To inspect another retained window, choose **View > Inspect evidence window** +and enter its zero-based index. Unit Stability keeps the full-recording time +context visible while you inspect individual windows. + +The Sampled ISI view excludes intervals crossing retained-window boundaries. +The Sampled Correlogram view shows raw pair counts from retained windows; it +does not apply a full-recording rate baseline. Both views are explicitly labeled +as sampled evidence. Its central bin counts unordered within-unit pairs and +combines both temporal orientations for cross-unit pairs whose absolute delay +is less than one bin, so tied sample times do not depend on source row order. +For exact pair counts, each unit may have at most 46,340 spikes in an evidence +window; the view reports an error for a larger per-unit window. Prepare shorter +windows if this limit is reached. + +Bounded curation supports whole-unit merges, cluster labels, and a source +amplitude threshold action. Use the GUI action **Split by source amplitude** +with one unit selected and a threshold in the source amplitude units. The action +uses the source `amplitudes.npy` scalar for every spike during replay; nonfinite +values stay with the remainder. Child units show full-data metrics as pending +until the replayed sorting is re-prepared, because their complete counts and +trends have not yet been calculated. Merge and label actions are journaled with +the bounded assignments and metadata in an atomic checkpoint generation. +Reopening loads the last complete generation. A reopened session starts a new +undo baseline. Each save checks that the checkpoint has not changed since this +session opened. If another session has saved first, the stale save fails; reopen +the bounded dataset before continuing to avoid overwriting that generation. + +Lasso, feature, and waveform splits remain unavailable in bounded mode because +the selected rows do not define how unseen source spikes should be assigned. +Standard, unbounded datasets retain their existing curation behavior. + +## Unit Stability on ordinary recordings + +Unit Stability can also be opened as an optional view for an ordinary dataset. +It evaluates only the selected unit and bins its loaded spikes over the recording +duration. Since ordinary datasets do not supply valid acquisition intervals to +this view, its rate assumes continuous acquisition; the plot labels that +assumption. The ordinary default view layout remains unchanged. Use the bounded +manifest workflow above when known acquisition gaps must be represented +explicitly. + +## Apply to the complete sorting + +After saving the bounded session, replay its journal against the exact source +assignment revision used during preparation: + +```bash +phy apply-curation /data/recording/params.py /data/recording-bounded \ + /data/recording-curated --batch-size 100000 +``` + +Replay checks the source identity and parent assignment revision, streams the +complete assignment array in batches, and writes a new full sorting. It rejects +a journal prepared from a different source revision. Keep the original source +and bounded package; they establish the lineage of the curated output. To review +the curated full sorting with bounded evidence again, run `prepare-large` on +the curated output to produce a fresh package and revision-specific summaries. +This replay and re-preparation step also resolves full-data counts and trends +for threshold-split children. + +## Interpreting evidence + +Full counts and exposure-based rates describe the complete sorting. Window plots, +waveforms, feature scatterplots, and other bounded measurements describe only +the retained evidence. Survey and rescue selection can overrepresent some +periods or units, and waveform availability can vary across units and time. +These views do not establish full-recording quality, contamination, refractory +violations, or continuous biological identity. Use them to guide review, not as +scientific claims about unseen spikes. Long-duration coordinates preserve the +source sample clock; silence alone is never evidence that the acquisition was +inactive. + +For reproducible synthetic capacity measurements, run +`uv run python tools/benchmark_long_recording.py`. The benchmark reports +preparation and replay time and peak memory in separate processes and checks +replay against a streamed oracle. It is capped at two million synthetic spikes; +its results are not a substitute for measurements on real chronic datasets. diff --git a/mkdocs.yml b/mkdocs.yml index bf8a4ad12..f59e59417 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -17,6 +17,7 @@ nav: - How-to and customization: - Configuration and state: configuration.md - Performance and spike sampling: performance.md + - Curating long recordings: long-recordings.md - Change keyboard shortcuts: keyboard_customization.md - Plugin introduction: customization.md - Plugin examples: plugins.md diff --git a/phy/apps/__init__.py b/phy/apps/__init__.py index db51602ea..c6767d88c 100644 --- a/phy/apps/__init__.py +++ b/phy/apps/__init__.py @@ -174,6 +174,57 @@ def cli_template_describe(ctx, params_path): template_describe(params_path) +# ------------------------------------------------------------------------------ +# Long-recording batch commands +# ------------------------------------------------------------------------------ + + +@phycli.command('prepare-large') +@click.argument('params-path', type=click.Path(exists=True, dir_okay=False, path_type=Path)) +@click.argument('output-dir', type=click.Path(path_type=Path)) +@click.option('--duration', 'duration_seconds', type=click.FloatRange(min=0, min_open=True)) +@click.option('--valid-intervals', type=click.Path(exists=True, dir_okay=False, path_type=Path)) +@click.option('--max-spikes', type=click.IntRange(min=1), default=2000000, show_default=True) +@click.option('--max-bytes', type=click.IntRange(min=1), default=536870912, show_default=True) +@click.option('--window-seconds', type=click.FloatRange(min=0, min_open=True), default=30.0) +@click.option('--stratum-seconds', type=click.FloatRange(min=0, min_open=True), default=3600.0) +@click.option('--seed', type=click.IntRange(min=0), default=0) +@click.option('--bin-seconds', type=click.FloatRange(min=0, min_open=True), default=3600.0) +@click.option('--batch-size', type=click.IntRange(min=1), default=100000) +def cli_prepare_large(params_path, output_dir, valid_intervals=None, **kwargs): + """Prepare bounded shared-window evidence and complete-sorting summaries. + + VALID_INTERVALS is an optional int64 .npy array of acquisition sample [start, end) pairs. + Duration/exposure must be authoritative; it is never inferred from spike silence. + """ + import numpy as np + + from .large_recording.prepare import prepare_large + + try: + intervals = np.load(valid_intervals, allow_pickle=False) if valid_intervals else None + output = prepare_large(params_path, output_dir, valid_intervals=intervals, **kwargs) + except (ValueError, OSError) as error: + raise click.ClickException(str(error)) from error + click.echo(f'Prepared bounded evidence: {output}') + + +@phycli.command('apply-curation') +@click.argument('source-params', type=click.Path(exists=True, dir_okay=False, path_type=Path)) +@click.argument('bounded-dir', type=click.Path(exists=True, file_okay=False, path_type=Path)) +@click.argument('output-dir', type=click.Path(path_type=Path)) +@click.option('--batch-size', type=click.IntRange(min=1), default=100000) +def cli_apply_curation(source_params, bounded_dir, output_dir, batch_size): + """Replay committed bounded curation into a new complete sorting.""" + from .large_recording.replay import apply_curation + + try: + output = apply_curation(source_params, bounded_dir, output_dir, batch_size=batch_size) + except (ValueError, OSError) as error: + raise click.ClickException(str(error)) from error + click.echo(f'Published complete curated sorting: {output}') + + # ------------------------------------------------------------------------------ # Kwik GUI # ------------------------------------------------------------------------------ @@ -224,6 +275,11 @@ def cli_alf_convert(ctx, subdirs, out_dir): from phylib.io.merge import Merger from phylib.io.model import load_model + if any((Path(path) / '_phy_curation_manifest.json').exists() for path in subdirs): + raise click.ClickException( + 'Apply bounded curation to the complete sorting before ALF conversion.' + ) + out_dir = Path(out_dir) if len(subdirs) >= 2: @@ -253,6 +309,11 @@ def template_extract_waveforms( """Extract spike waveforms.""" from phylib.io.model import load_model + if (Path(params_path).parent / '_phy_curation_manifest.json').exists(): + raise click.ClickException( + 'Extract waveforms on the original source and regenerate bounded evidence; bounded package assets are immutable.' + ) + model = load_model(params_path) model.save_spikes_subset_waveforms( max_n_spikes_per_template=n_spikes_per_cluster, max_n_channels=nc diff --git a/phy/apps/base.py b/phy/apps/base.py index 593ddfcff..63c8a7fc2 100644 --- a/phy/apps/base.py +++ b/phy/apps/base.py @@ -2449,7 +2449,11 @@ def create_gui(self, default_views=None, **kwargs): subtitle=str(self.dir_path), config_dir=self.config_dir, local_path=self.cache_dir / 'state.json', - default_state_path=Path(inspect.getfile(self.__class__)).parent / 'static/state.json', + default_state_path=getattr( + self, + 'default_state_path', + Path(inspect.getfile(self.__class__)).parent / 'static/state.json', + ), view_creator=self.view_creator, default_views=default_views, enable_threading=self._enable_threading, diff --git a/phy/apps/large_recording/__init__.py b/phy/apps/large_recording/__init__.py new file mode 100644 index 000000000..373145000 --- /dev/null +++ b/phy/apps/large_recording/__init__.py @@ -0,0 +1 @@ +"""Bounded evidence preparation and full-sorting curation replay.""" diff --git a/phy/apps/large_recording/adapter.py b/phy/apps/large_recording/adapter.py new file mode 100644 index 000000000..cc5a351ac --- /dev/null +++ b/phy/apps/large_recording/adapter.py @@ -0,0 +1,372 @@ +"""Manifest-gated template controller for bounded recording evidence.""" + +from copy import deepcopy +from pathlib import Path + +import numpy as np +from phylib.io.model import get_template_params +from phylib.utils import Bunch + +from phy.cluster._propositions import MergePropositionController, decode_review_mapping + +from ..base import _allocate_spike_counts, _sample_spikes_evenly +from ..template.gui import TemplateController +from .curation import BoundedSession +from .model import BoundedTemplateModel +from .statistics import window_correlograms, window_isi +from .supervisor import BoundedSupervisor +from .views import EvidenceAmplitudeView, SampledCorrelogramView, SampledISIView + + +class BoundedTemplateController(TemplateController): + gui_name = 'BoundedTemplateGUI' + # Bounded mode has its own state namespace and initial view list. Importing + # ordinary saved dock geometry would suppress the new evidence views. + default_state_path = None + default_views = ( + 'UnitStabilityView', + 'WaveformView', + 'EvidenceAmplitudeView', + 'FeatureView', + 'ProbeView', + 'SampledISIView', + 'SampledCorrelogramView', + ) + + def __init__(self, dir_path=None, model=None, **kwargs): + if model is not None: + raise ValueError( + 'Open bounded packages from their directory, without a supplied model.' + ) + self.bounded_session = BoundedSession(dir_path) + self.evidence_window_index = 0 + kwargs['enable_threading'] = False + kwargs.update(get_template_params(Path(dir_path) / 'params.py')) + kwargs.pop('dir_path', None) + super().__init__(dir_path=dir_path, **kwargs) + + def _create_model(self, dir_path=None, **kwargs): + model = BoundedTemplateModel(dir_path=dir_path, session=self.bounded_session, **kwargs) + if model.sample_rate != self.bounded_session.manifest['sample_rate']: + model.close() + raise ValueError('Bounded params sampling rate disagrees with the manifest.') + model.spike_clusters = self.bounded_session.assignments.copy() + model.metadata = deepcopy(self.bounded_session.metadata) + model.duration = self.bounded_session.manifest['num_samples'] / model.sample_rate + + def reject_save(*args, **kwargs): + raise ValueError('Save bounded curation through supervisor.save() as one transaction.') + + model.save_spike_clusters = reject_save + model.save_metadata = reject_save + model.save_spikes_subset_waveforms = reject_save + model.freeze_source() + return model + + def _cache_methods(self): + # Ordinary persistent feature caches omit the current evidence window. + # Bounded queries are already capped; evaluate them against live state. + return + + def close(self, close_model=True): + super().close(close_model=close_model) + if close_model: + self.bounded_session.close() + + def _set_supervisor(self): + session = self.bounded_session + propositions = self._load_merge_propositions() + self.supervisor = BoundedSupervisor( + spike_clusters=self.model.spike_clusters.copy(), + cluster_groups=self.model.metadata.get('group', {}), + cluster_metrics=self.cluster_metrics, + cluster_labels=self.model.metadata, + similarity=self.similarity_functions[self.similarity], + session=session, + merge_propositions=propositions, + ) + self.supervisor._validate_native() + + def _load_merge_propositions(self): + controller = super()._load_merge_propositions() + reviews = self.bounded_session.extra_documents.get('curation_review.json') + if controller is not None and reviews is not None: + source, mapping = decode_review_mapping(reviews) + if source.get('sha256') != self._merge_proposition_sha256: + raise ValueError('Committed reviews belong to different merge propositions.') + controller = MergePropositionController(controller.catalog._replace(reviews=mapping)) + return controller + + def _set_cluster_metrics(self): + super()._set_cluster_metrics() + self.cluster_metrics['full_count'] = self.get_full_count + self.cluster_metrics['evidence_count'] = lambda c: self.supervisor.n_spikes(c) + self.cluster_metrics['waveform_count'] = self.get_waveform_count + self.cluster_metrics['evidence'] = lambda c: ( + 'pending replay' if self.get_full_count(c) is None else 'sampled windows' + ) + + def _set_view_creator(self): + super()._set_view_creator() + for name in ( + 'ISIView', + 'CorrelogramView', + 'FiringRateView', + 'TraceView', + 'TraceImageView', + ): + self.view_creator.pop(name, None) + self.view_creator.pop('AmplitudeView', None) + self.view_creator['EvidenceAmplitudeView'] = self.create_evidence_amplitude_view + self.view_creator['UnitStabilityView'] = self.create_unit_stability_view + self.view_creator['SampledISIView'] = lambda: SampledISIView(cluster_stat=self._get_isi) + self.view_creator['SampledCorrelogramView'] = lambda: SampledCorrelogramView( + correlograms=self._window_correlograms, sample_rate=self.model.sample_rate + ) + + def _get_waveforms_dict(self): + functions = super()._get_waveforms_dict() + if not self.bounded_session.manifest['capabilities']['stored_waveforms']: + functions.pop('waveforms', None) + functions.pop('mean_waveforms', None) + return functions + + def create_feature_view(self): + if ( + self.model.features is None + and not self.bounded_session.manifest['capabilities']['stored_waveforms'] + ): + return + view = super().create_feature_view() + if view is not None: + view.ex_status = self._window_status() + return view + + def _source_units(self, cluster_id): + version = self.bounded_session.events + cached = getattr(self, '_lineage_cache', None) + if cached is not None and cached[0] is version: + return cached[1].get(int(cluster_id)) + units = {int(c): {int(c)} for c in self.bounded_session.unit_ids} + for event in self.bounded_session.events: + if event['operation'] == 'merge': + parents = [units.pop(c) for c in event['parents']] + units[event['target']] = ( + None if any(p is None for p in parents) else set.union(*parents) + ) + elif event['operation'] == 'threshold_split': + units.pop(event['parent']) + units.update(dict.fromkeys(event['children'])) + self._lineage_cache = version, units + return units.get(int(cluster_id)) + + def get_full_count(self, cluster_id): + units = self._source_units(cluster_id) + if units is None: + return None + if not hasattr(self, '_full_counts'): + self._full_counts = dict( + zip(self.bounded_session.unit_ids, self.bounded_session.manifest['counts']) + ) + counts = self._full_counts + return sum(int(counts[c]) for c in units) + + def get_waveform_count(self, cluster_id): + stored = self.model.spike_waveforms + if stored is None: + return 0 + if ( + not hasattr(self, '_waveform_membership_cache') + or self._waveform_membership_cache[0] is not self.bounded_session.events + ): + labels = self.supervisor.clustering.spike_clusters[stored.spike_ids] + ids, counts = np.unique(labels, return_counts=True) + self._waveform_membership_cache = self.bounded_session.events, dict(zip(ids, counts)) + return int(self._waveform_membership_cache[1].get(cluster_id, 0)) + + def get_mean_firing_rate(self, cluster_id): + count = self.get_full_count(cluster_id) + if count is None: + return None + intervals = np.asarray(self.bounded_session.manifest['valid_intervals']) + exposure = np.sum(intervals[:, 1] - intervals[:, 0]) / self.model.sample_rate + return count / exposure + + @property + def evidence_windows(self): + if not hasattr(self, '_evidence_windows'): + self._evidence_windows = np.load( + self.dir_path / self.bounded_session.manifest['windows_file'] + ) + return self._evidence_windows + + def _get_isi(self, cluster_id): + ids = self.supervisor.clustering.spikes_per_cluster.get(cluster_id, []) + samples = self.model.spike_samples[ids] + return Bunch(data=window_isi(samples, self.evidence_windows, self.model.sample_rate)) + + def _window_correlograms(self, cluster_ids, bin_size, window_size): + return window_correlograms( + self.model.spike_samples, + self.supervisor.clustering.spike_clusters, + self.evidence_windows, + cluster_ids, + self.model.sample_rate, + bin_size, + window_size, + ) + + def create_unit_stability_view(self): + from .stability import UnitStabilityView + + return UnitStabilityView(self) + + @property + def current_evidence_window(self): + return self.evidence_windows[self.evidence_window_index] + + def _window_ids(self, cluster_id=None): + start, end = self.current_evidence_window + if cluster_id is None: + first, last = np.searchsorted(self.model.spike_samples, [start, end]) + return np.arange(first, last, dtype=np.int64) + ids = self.supervisor.clustering.spikes_per_cluster.get( + cluster_id, np.array([], dtype=np.int64) + ) + times = self.model.spike_samples[ids] + first, last = np.searchsorted(times, [start, end]) + return ids[first:last] + + def _window_time(self, ids): + # Subtract int64 samples before floating conversion, even at hour 1000. + return ( + self.model.spike_samples[ids] - self.current_evidence_window[0] + ) / self.model.sample_rate + + def _get_feature_view_spike_ids(self, cluster_id=None, load_all=False): + ids = self._window_ids(cluster_id) + if self.model.features is None: + if self.model.spike_waveforms is None: + return np.array([], dtype=np.int64) + ids = np.intersect1d(ids, self.model.spike_waveforms.spike_ids) + if getattr(self.model, 'features_rows', None) is not None: + ids = np.intersect1d(ids, self.model.features_rows) + limit = self.n_spikes_features_background if cluster_id is None else self.n_spikes_features + return ids if load_all else _sample_spikes_evenly(ids, limit) + + def _get_feature_view_spike_times(self, cluster_id=None, load_all=False): + ids = self._get_feature_view_spike_ids(cluster_id, load_all=load_all) + if not len(ids): + return + start, end = self.current_evidence_window + return Bunch( + data=self._window_time(ids), + spike_ids=ids, + lim=(0.0, (end - start) / self.model.sample_rate), + ) + + def _evidence_amplitudes(self, cluster_ids, load_all=False): + active = [c for c in cluster_ids if c is not None] + eligible = [self._window_ids(c) for c in active] + counts = _allocate_spike_counts( + [len(ids) for ids in eligible], + per_cluster=None if load_all else self.n_spikes_amplitudes, + total=None if load_all else self.n_spikes_amplitudes_total, + ) + values = self.bounded_session.attributes['amplitudes.npy'].reshape(-1) + output = [] + if None in cluster_ids: + ids = self._window_ids() + ids = ids[~np.isin(self.supervisor.clustering.spike_clusters[ids], active)] + ids = _sample_spikes_evenly(ids, self.n_spikes_amplitudes_background) + ids = ids[np.isfinite(values[ids])] + output.append( + Bunch(amplitudes=values[ids], spike_ids=ids, spike_times=self._window_time(ids)) + ) + for ids, count in zip(eligible, counts): + ids = _sample_spikes_evenly(ids, int(count)) + ids = ids[np.isfinite(values[ids])] + output.append( + Bunch(amplitudes=values[ids], spike_ids=ids, spike_times=self._window_time(ids)) + ) + return output + + def create_evidence_amplitude_view(self): + array = self.bounded_session.attributes.get('amplitudes.npy') + if array is None or (array.ndim != 1 and array.shape[1:] != (1,)): + return + start, end = self.current_evidence_window + view = EvidenceAmplitudeView( + amplitudes={'source': self._evidence_amplitudes}, + amplitudes_type='source', + duration=(end - start) / self.model.sample_rate, + split_is_eligible=lambda: False, + ) + view.ex_status = self._window_status() + view.set_window_context( + self.evidence_window_index, + start / self.model.sample_rate, + end / self.model.sample_rate, + ) + return view + + def _window_status(self): + start, end = self.current_evidence_window + return ( + f'Evidence window {self.evidence_window_index}: {start / self.model.sample_rate:.3f}–' + f'{end / self.model.sample_rate:.3f} s; time axis is seconds from window start; nonfinite amplitudes omitted.' + ) + + def inspect_evidence_window(self, index): + index = int(index) + if not 0 <= index < len(self.evidence_windows): + raise ValueError('Evidence window index is outside the retained plan.') + self.evidence_window_index = index + gui = getattr(self.supervisor, 'gui', None) + if gui is not None: + start, end = self.current_evidence_window + for view in gui.list_views('EvidenceAmplitudeView', 'FeatureView'): + view.ex_status = self._window_status() + if isinstance(view, EvidenceAmplitudeView): + view.duration = (end - start) / self.model.sample_rate + view.set_window_context( + index, start / self.model.sample_rate, end / self.model.sample_rate + ) + view.canvas.panzoom.reset() + if view.cluster_ids: + view.plot() + view.update_status() + + def create_gui(self, *args, **kwargs): + gui = super().create_gui(*args, **kwargs) + gui.setWindowTitle( + f'phy — Bounded evidence ({self.model.n_spikes:,} / ' + f'{self.bounded_session.manifest["num_spikes"]:,} spikes)' + ) + + @self.supervisor.actions.add(prompt=True, name='Split by source amplitude') + def split_by_source_amplitude(threshold): + """Split one selected unit below a source amplitude threshold (source units).""" + ids = self.supervisor.selected + if len(ids) != 1: + raise ValueError('Select one unit for a source-amplitude threshold split.') + return self.supervisor.split_source_amplitude(ids[0], float(threshold)) + + @gui.view_actions.add(prompt=True, prompt_default=lambda: self.evidence_window_index) + def inspect_evidence_window(index): + """Inspect a retained evidence window by its zero-based index.""" + self.inspect_evidence_window(index) + + split_action = self.supervisor.actions.get('split') + if split_action is not None: + split_action.setEnabled(False) + split_action.setToolTip( + 'Selected-row splits cannot assign unseen source spikes; use a source-amplitude predicate.' + ) + predicate_action = self.supervisor.actions.get('Split by source amplitude') + if predicate_action is not None: + predicate_action.setEnabled( + 'amplitudes.npy' in self.bounded_session.capabilities['predicate_split_attributes'] + ) + + return gui diff --git a/phy/apps/large_recording/curation.py b/phy/apps/large_recording/curation.py new file mode 100644 index 000000000..6a45d7f70 --- /dev/null +++ b/phy/apps/large_recording/curation.py @@ -0,0 +1,793 @@ +"""Semantic bounded curation and crash-consistent, immutable save generations.""" + +import copy +import hashlib +import json +import os +import shutil +import uuid +from contextlib import contextmanager +from pathlib import Path + +import numpy as np + +from .source import hash_json + +MANIFEST = '_phy_curation_manifest.json' +BASE_ASSIGNMENTS = '_phy_spike_clusters_base.npy' +BOOLEAN_CAPABILITIES = ( + 'merge', + 'label', + 'split', + 'stored_waveforms', + 'raw_traces', + 'normalized_correlograms', +) + + +def _json_bytes(value): + return ( + json.dumps(value, sort_keys=True, separators=(',', ':'), allow_nan=False) + '\n' + ).encode() + + +def _digest(value): + return hashlib.sha256(_json_bytes(value)).hexdigest() + + +def _hash_file(path): + digest = hashlib.sha256() + with Path(path).open('rb') as stream: + while chunk := stream.read(1024 * 1024): + digest.update(chunk) + return digest.hexdigest() + + +def _read_json(path): + with Path(path).open(encoding='utf8') as stream: + return json.load(stream) + + +def _write_json(path, value): + with Path(path).open('wb') as stream: + stream.write(_json_bytes(value)) + stream.flush() + os.fsync(stream.fileno()) + + +def _fsync_directory(path): + if os.name != 'nt': + descriptor = os.open(path, os.O_RDONLY) + try: + os.fsync(descriptor) + finally: + os.close(descriptor) + + +def _safe_file(root, name): + if not isinstance(name, str) or Path(name).name != name or name in ('', '.', '..'): + raise ValueError('Manifest/checkpoint filename must be a plain local filename.') + path = Path(root) / name + if path.is_symlink(): + raise ValueError('Manifest/checkpoint files must not be symbolic links.') + return path + + +@contextmanager +def _commit_lock(root): + """Hold a process-safe OS lock; process exit releases it automatically.""" + with _safe_file(root, 'commit.lock').open('a+b') as stream: + if os.name == 'nt': + import msvcrt + + # Windows byte-range locking requires at least one byte in the file. + stream.seek(0, os.SEEK_END) + if stream.tell() == 0: + stream.write(b'\0') + stream.flush() + stream.seek(0) + msvcrt.locking(stream.fileno(), msvcrt.LK_LOCK, 1) + try: + yield + finally: + stream.seek(0) + msvcrt.locking(stream.fileno(), msvcrt.LK_UNLCK, 1) + else: + import fcntl + + fcntl.flock(stream.fileno(), fcntl.LOCK_EX) + try: + yield + finally: + fcntl.flock(stream.fileno(), fcntl.LOCK_UN) + + +def _current_pointer(root): + try: + return _read_json(root / 'current.json') + except FileNotFoundError: + return None + + +def _integer(value, name): + if isinstance(value, bool) or not isinstance(value, (int, np.integer)) or value < 0: + raise ValueError(f'{name} must be a nonnegative integer.') + return int(value) + + +def _validate_capabilities(manifest): + capabilities = manifest.get('capabilities', {}) + if not isinstance(capabilities, dict): + raise ValueError('Bounded capabilities must be a mapping.') + if ( + set(capabilities) + - set(BOOLEAN_CAPABILITIES) + - { + 'full_attributes', + 'predicate_split_attributes', + } + ): + raise ValueError('Unsupported bounded capability declaration.') + result = {} + for name in BOOLEAN_CAPABILITIES: + value = capabilities.get(name, False) + if type(value) is not bool: + raise ValueError(f'Bounded capability {name!r} must be a boolean.') + result[name] = value + if result['split']: + raise ValueError('General selected-row splits are unsupported in bounded mode.') + for name in ('full_attributes', 'predicate_split_attributes'): + values = capabilities.get(name, []) + if ( + not isinstance(values, list) + or any(not isinstance(value, str) for value in values) + or len(set(values)) != len(values) + ): + raise ValueError(f'Bounded capability {name!r} must be a unique string list.') + result[name] = list(values) + if set(result['predicate_split_attributes']) - {'amplitudes.npy'}: + raise ValueError( + 'Unsupported bounded predicate attribute; only raw amplitudes are supported.' + ) + return result + + +def _field_name(field): + if ( + not isinstance(field, str) + or not field + or field in ('.', '..') + or any(char in field for char in ('/', '\\', '\0', '\n', '\r', '\t')) + ): + raise ValueError('Metadata fields must be nonempty safe local field names.') + return field + + +def normalize_metadata(metadata): + """Return JSON-compatible metadata with integer cluster keys.""" + if not isinstance(metadata, dict): + raise ValueError('Cluster metadata must be a field-to-cluster mapping.') + result = {} + for field, values in metadata.items(): + _field_name(field) + if not isinstance(values, dict): + raise ValueError('Invalid cluster metadata field.') + result[field] = {} + for key, value in values.items(): + try: + cluster = int(key) + except (ValueError, TypeError) as error: + raise ValueError('Invalid metadata cluster ID.') from error + if str(cluster) != str(key) or cluster < 0: + raise ValueError('Invalid metadata cluster ID.') + _json_bytes(value) + result[field][cluster] = copy.deepcopy(value) + return result + + +def validate_predicate(predicate): + """Validate a threshold on a canonical full-source scalar attribute.""" + if not isinstance(predicate, dict) or set(predicate) != { + 'attribute', + 'operator', + 'threshold', + 'nonfinite', + }: + raise ValueError('Threshold requires attribute/operator/threshold/nonfinite.') + if predicate['attribute'] not in ('amplitudes.npy', 'spike_templates.npy'): + raise ValueError('Unsupported full-source threshold attribute.') + if predicate['operator'] not in ('<', '<=', '>', '>='): + raise ValueError('Unsupported threshold comparator.') + threshold = predicate['threshold'] + if isinstance(threshold, bool) or not isinstance(threshold, (int, float)): + raise ValueError('Threshold must be a finite number.') + if not np.isfinite(threshold) or predicate['nonfinite'] not in ('keep', 'reject'): + raise ValueError('Threshold needs a finite threshold and explicit nonfinite policy.') + return copy.deepcopy(predicate) + + +def evaluate_predicate(values, predicate): + values = np.asarray(values).reshape(-1) + finite = np.isfinite(values) + if predicate['nonfinite'] == 'reject' and not finite.all(): + raise ValueError('Nonfinite source values violate threshold policy.') + operators = {'<': np.less, '<=': np.less_equal, '>': np.greater, '>=': np.greater_equal} + return finite & operators[predicate['operator']](values, predicate['threshold']) + + +def interpret_events(unit_ids, metadata, events): + """Validate event intent and resolve the final catalog and metadata.""" + active = set(map(int, unit_ids)) + if not active or any(i < 0 for i in active): + raise ValueError('Full source unit catalog must contain nonnegative IDs.') + allocated = set(active) + metadata = normalize_metadata(metadata) + if not isinstance(events, list): + raise ValueError('Semantic journal must be an ordered event list.') + next_id = max(active) + 1 + for event in events: + if not isinstance(event, dict) or event.get('version') != 1: + raise ValueError('Unsupported semantic journal event version.') + operation = event.get('operation') + if operation == 'merge': + if set(event) != {'version', 'operation', 'parents', 'target', 'metadata'}: + raise ValueError('Malformed merge event.') + parents = [_integer(i, 'Parent ID') for i in event['parents']] + target = _integer(event['target'], 'Target ID') + if target > np.iinfo(np.int32).max: + raise ValueError('Target cluster ID exceeds the current Phy loader namespace.') + if len(parents) < 2 or len(set(parents)) != len(parents): + raise ValueError('Merge requires at least two distinct complete units.') + if not set(parents) <= active or target in allocated or target < next_id: + raise ValueError('Merge references inactive units or colliding target ID.') + resolved = event['metadata'] + if not isinstance(resolved, dict) or any(not isinstance(f, str) for f in resolved): + raise ValueError('Merge metadata must contain resolved field values.') + _json_bytes(resolved) + for field in resolved: + _field_name(field) + active.difference_update(parents) + active.add(target) + allocated.add(target) + next_id = target + 1 + for values in metadata.values(): + for parent in parents: + values.pop(parent, None) + for field, value in resolved.items(): + metadata.setdefault(field, {})[target] = copy.deepcopy(value) + elif operation == 'label': + if set(event) != {'version', 'operation', 'field', 'ids', 'value'}: + raise ValueError('Malformed label event.') + ids = [_integer(i, 'Label ID') for i in event['ids']] + field = event['field'] + _field_name(field) + if not ids or len(set(ids)) != len(ids) or not set(ids) <= active: + raise ValueError('Label references inactive units or invalid metadata field.') + _json_bytes(event['value']) + values = metadata.setdefault(field, {}) + for cluster in ids: + values[cluster] = copy.deepcopy(event['value']) + elif operation == 'threshold_split': + if set(event) != { + 'version', + 'operation', + 'parent', + 'children', + 'predicate', + 'metadata', + }: + raise ValueError('Malformed threshold split event.') + parent = _integer(event['parent'], 'Parent ID') + children = [_integer(i, 'Child ID') for i in event['children']] + if any(child > np.iinfo(np.int32).max for child in children): + raise ValueError('Child cluster ID exceeds the current Phy loader namespace.') + validate_predicate(event['predicate']) + if parent not in active or len(children) != 2 or len(set(children)) != 2: + raise ValueError('Threshold split requires an active parent and two children.') + if any(i in allocated or i < next_id for i in children): + raise ValueError('Threshold split child ID collision.') + resolved = normalize_metadata(event['metadata']) + if any(not set(v) <= set(children) for v in resolved.values()): + raise ValueError('Split metadata must refer only to explicit children.') + active.remove(parent) + active.update(children) + allocated.update(children) + next_id = max(children) + 1 + for values in metadata.values(): + values.pop(parent, None) + for field, values in resolved.items(): + metadata.setdefault(field, {}).update(values) + else: + raise ValueError(f'Unsupported semantic operation: {operation!r}.') + return active, metadata, next_id + + +def project_assignments(assignments, events, attributes=None): + """Interpret a validated journal against a bounded array or full-source batch.""" + output = np.asarray(assignments, dtype=np.int64).reshape(-1).copy() + attributes = attributes or {} + for event in events: + if event['operation'] == 'merge': + output[np.isin(output, event['parents'])] = event['target'] + elif event['operation'] == 'threshold_split': + predicate = event['predicate'] + attribute = predicate['attribute'] + if attribute not in attributes: + raise ValueError(f'Full-source threshold attribute is unavailable: {attribute}.') + mask = output == event['parent'] + values = np.asarray(attributes[attribute]).reshape(-1) + if len(values) != len(output): + raise ValueError('Threshold attribute is not a scalar per-spike array.') + selected = np.zeros(len(output), dtype=bool) + selected[mask] = evaluate_predicate(values[mask], predicate) + output[mask & ~selected] = event['children'][0] + output[mask & selected] = event['children'][1] + return output + + +class BoundedSession: + """Bounded semantic state; reopen establishes a new undo baseline. + + Snapshot/restore integrate with the controller's undo stack. Saving publishes + one pointer after every generation file has been flushed and verified. + """ + + def __init__(self, dataset_path): + self.global_ids = None + self.base_assignments = None + self.attributes = {} + try: + self._initialize(dataset_path) + except BaseException: + self.close() + raise + + def _initialize(self, dataset_path): + self.path = Path(dataset_path) + if self.path.is_file(): + self.path = self.path.parent + self.manifest = _read_json(self.path / MANIFEST) + manifest = self.manifest + required = { + 'source', + 'sample_rate', + 'num_samples', + 'valid_intervals', + 'num_spikes', + 'bounded_num_spikes', + 'unit_ids', + 'counts', + 'next_cluster_id', + 'global_spike_ids_file', + 'windows_file', + } + if manifest.get('format') != 'phy_bounded_curation' or manifest.get('format_version') != 1: + raise ValueError('Unsupported bounded curation manifest version.') + if not required <= set(manifest): + raise ValueError( + 'Incomplete bounded curation manifest; regenerate with prepare-large.' + ) + source = manifest['source'] + if not isinstance(source, dict) or not {'identity', 'revision', 'path'} <= set(source): + raise ValueError('Bounded manifest needs source identity and assignment revision.') + if not isinstance(source['identity'], dict) or not isinstance(source['revision'], str): + raise ValueError('Invalid source identity or revision.') + identity = source['identity'] + identity_fields = {'format_version', 'acquisition', 'time_base', 'arrays', 'sha256'} + if ( + not identity_fields <= set(identity) + or identity.get('format_version') != 1 + or identity['sha256'] + != hash_json({k: v for k, v in identity.items() if k != 'sha256'}) + or len(source['revision']) != 64 + or any(c not in '0123456789abcdef' for c in source['revision']) + ): + raise ValueError('Invalid source coordinate identity or assignment revision digest.') + n = _integer(manifest['bounded_num_spikes'], 'Bounded spike count') + full_n = _integer(manifest['num_spikes'], 'Full spike count') + samples = _integer(manifest['num_samples'], 'Recording sample count') + if not n or n > full_n or not samples or not np.isfinite(manifest['sample_rate']): + raise ValueError('Invalid bounded counts or recording duration.') + if manifest['sample_rate'] <= 0: + raise ValueError('Sampling rate must be positive.') + if identity['time_base'] != { + 'sample_rate': manifest['sample_rate'], + 'num_samples': samples, + 'valid_intervals': manifest['valid_intervals'], + }: + raise ValueError('Bounded time base disagrees with immutable source identity.') + self.capabilities = _validate_capabilities(manifest) + self.unit_ids = [_integer(i, 'Unit ID') for i in manifest['unit_ids']] + if any(unit > np.iinfo(np.int32).max for unit in self.unit_ids): + raise ValueError('Unit IDs exceed the current Phy loader namespace.') + if not self.unit_ids or len(set(self.unit_ids)) != len(self.unit_ids): + raise ValueError('Full unit catalog must contain distinct unit IDs.') + counts = np.asarray(manifest['counts']) + if ( + counts.shape != (len(self.unit_ids),) + or counts.dtype.kind not in 'iu' + or np.any(counts < 0) + or sum(map(int, counts)) != full_n + ): + raise ValueError('Full unit counts do not match the source spike count/catalog.') + self.base_next_cluster_id = _integer(manifest['next_cluster_id'], 'Next cluster ID') + if self.base_next_cluster_id <= max(self.unit_ids): + raise ValueError('Cluster allocator collides with the complete source catalog.') + intervals = np.asarray(manifest['valid_intervals']) + if ( + intervals.ndim != 2 + or intervals.shape[1] != 2 + or not len(intervals) + or intervals.dtype.kind not in 'iu' + ): + raise ValueError('Valid acquisition intervals must be integer [start, end) pairs.') + if ( + np.any(intervals[:, 0] < 0) + or np.any(intervals[:, 1] > samples) + or np.any(intervals[:, 1] <= intervals[:, 0]) + or np.any(intervals[1:, 0] < intervals[:-1, 1]) + ): + raise ValueError('Invalid or overlapping acquisition intervals.') + self.global_ids = np.load( + _safe_file(self.path, manifest['global_spike_ids_file']), mmap_mode='r' + ) + if self.global_ids.dtype != np.int64 or self.global_ids.shape != (n,): + raise ValueError('Global spike IDs must be a one-dimensional int64 mapping.') + if ( + np.any(self.global_ids < 0) + or np.any(self.global_ids >= full_n) + or np.any(self.global_ids[1:] <= self.global_ids[:-1]) + ): + raise ValueError('Global spike IDs must be sorted, unique source row indices.') + windows = np.load(_safe_file(self.path, manifest['windows_file']), mmap_mode='r') + if ( + windows.ndim != 2 + or windows.shape[1] != 2 + or not len(windows) + or windows.dtype != np.int64 + ): + raise ValueError('Evidence windows must be int64 sample intervals.') + if ( + np.any(windows[:, 0] < 0) + or np.any(windows[:, 1] > samples) + or np.any(windows[:, 1] <= windows[:, 0]) + or np.any(windows[1:, 0] < windows[:-1, 1]) + ): + raise ValueError('Invalid or overlapping evidence windows.') + for start, stop in windows: + if not any(left <= start and stop <= right for left, right in intervals): + raise ValueError('Evidence windows cross unsupported acquisition gaps.') + for name, digest in manifest.get('bounded_files', {}).items(): + if _hash_file(_safe_file(self.path, name)) != digest: + raise ValueError(f'Bounded package integrity check failed: {name}.') + hashed_files = manifest.get('bounded_files', {}) + if not { + BASE_ASSIGNMENTS, + manifest['global_spike_ids_file'], + manifest['windows_file'], + 'spike_times.npy', + } <= set(hashed_files): + raise ValueError('Bounded manifest lacks required immutable array digests.') + self.base_assignments = np.load(self.path / BASE_ASSIGNMENTS, mmap_mode='r') + if ( + self.base_assignments.shape != (n,) + or self.base_assignments.dtype.kind not in 'iu' + or not set(np.unique(self.base_assignments)) <= set(self.unit_ids) + ): + raise ValueError('Immutable bounded assignments do not match the full catalog.') + self.base_metadata = normalize_metadata( + manifest.get('base_metadata', manifest.get('metadata', {})) + ) + if any(not set(values) <= set(self.unit_ids) for values in self.base_metadata.values()): + raise ValueError('Base metadata refers to units outside the full source catalog.') + bounded_times = np.load(self.path / 'spike_times.npy', mmap_mode='r').reshape(-1) + if ( + bounded_times.shape != (n,) + or bounded_times.dtype.kind not in 'iu' + or np.any(bounded_times[1:] < bounded_times[:-1]) + ): + raise ValueError('Bounded spike sample coordinates have invalid alignment/order.') + window_ids = np.searchsorted(windows[:, 0], bounded_times, side='right') - 1 + if np.any(window_ids < 0) or np.any(bounded_times >= windows[window_ids, 1]): + raise ValueError('Bounded spikes lie outside the recorded evidence windows.') + self.attributes = {} + for name in ('amplitudes.npy', 'spike_templates.npy'): + if (self.path / name).exists(): + self.attributes[name] = np.load(self.path / name, mmap_mode='r') + for name in self.capabilities['predicate_split_attributes']: + array = self.attributes.get(name) + if ( + array is None + or name not in identity['arrays'] + or not (array.ndim == 1 or array.shape[1:] == (1,)) + or len(array) != n + or array.dtype.kind not in 'iuf' + ): + raise ValueError( + f'Bounded predicate attribute is not available as a source scalar: {name}.' + ) + self.assignments = np.asarray(self.base_assignments, dtype=np.int64).copy() + self.metadata = copy.deepcopy(self.base_metadata) + self.events = [] + self.extra_documents = {} + self.next_cluster_id = self.base_next_cluster_id + self.checkpoint_path = None + self._expected_pointer = None + pointer_path = self.path / '.phy-curation' / 'current.json' + if pointer_path.exists(): + self._load_checkpoint(pointer_path) + + def _validated_state(self, assignments, metadata, events, next_cluster_id): + self._validate_event_capabilities(events) + _, projected_metadata, minimum = interpret_events( + self.unit_ids, self.base_metadata, events + ) + projected = project_assignments(self.base_assignments, events, self.attributes) + values = np.asarray(assignments) + if values.dtype.kind not in 'iu' or values.shape != projected.shape: + raise ValueError('Checkpoint assignments have invalid shape or dtype.') + if not np.array_equal(values, projected): + raise ValueError('Checkpoint assignments disagree with semantic journal projection.') + if normalize_metadata(metadata) != projected_metadata: + raise ValueError('Checkpoint metadata disagrees with semantic journal projection.') + minimum = max(minimum, self.base_next_cluster_id) + next_cluster_id = _integer(next_cluster_id, 'Next cluster ID') + if next_cluster_id < minimum: + raise ValueError('Checkpoint allocator collides with journal IDs.') + return projected, projected_metadata, next_cluster_id + + def _load_checkpoint(self, pointer_path): + pointer = _read_json(pointer_path) + generation = _safe_file(pointer_path.parent / 'generations', pointer['generation']) + descriptor = _read_json(generation / 'checkpoint.json') + if _digest(descriptor) != pointer['sha256']: + raise ValueError('Checkpoint commit pointer digest is invalid.') + if descriptor.get('version') != 1 or descriptor.get('manifest_sha256') != _digest( + self.manifest + ): + raise ValueError('Checkpoint belongs to a different bounded manifest.') + for name, digest in descriptor['files'].items(): + if _hash_file(_safe_file(generation, name)) != digest: + raise ValueError(f'Checkpoint corruption: {name}.') + if not {'spike_clusters.npy', 'journal.json', 'metadata.json'} <= set(descriptor['files']): + raise ValueError('Incomplete checkpoint file inventory.') + events = _read_json(generation / 'journal.json') + for event in events: + if event.get('operation') == 'threshold_split': + event['metadata'] = normalize_metadata(event['metadata']) + metadata = _read_json(generation / 'metadata.json') + assignments = np.load(generation / 'spike_clusters.npy', mmap_mode='r') + values, metadata, next_id = self._validated_state( + assignments, metadata, events, descriptor['next_cluster_id'] + ) + self.assignments, self.metadata = values, metadata + self.events = events + self.next_cluster_id = next_id + self.checkpoint_path = generation + self._expected_pointer = copy.deepcopy(pointer) + self.extra_documents = { + name: _read_json(generation / name) for name in descriptor.get('extra_documents', []) + } + + def snapshot(self): + return { + 'assignments': self.assignments.copy(), + 'metadata': copy.deepcopy(self.metadata), + 'events': copy.deepcopy(self.events), + 'next_cluster_id': self.next_cluster_id, + } + + def journal_state(self): + """Return a compact history prefix without copying bounded spike arrays. + + Commands append event dictionaries and never modify earlier events. The + tuple shares those immutable event records between history entries. + """ + return {'events': tuple(self.events), 'next_cluster_id': self.next_cluster_id} + + def restore_journal(self, state): + """Restore an active history prefix while preserving allocator high-water.""" + events = list(state['events']) + self._validate_event_capabilities(events) + _, metadata, minimum = interpret_events(self.unit_ids, self.base_metadata, events) + next_id = _integer(state['next_cluster_id'], 'Next cluster ID') + if next_id < max(minimum, self.base_next_cluster_id): + raise ValueError('History allocator collides with semantic journal IDs.') + assignments = project_assignments(self.base_assignments, events, self.attributes) + self.assignments, self.metadata, self.events = assignments, metadata, events + self.next_cluster_id = max(self.next_cluster_id, next_id) + + def close(self): + """Release file-backed bounded readers independently of the GUI model.""" + for array in (self.global_ids, self.base_assignments, *self.attributes.values()): + mapping = getattr(array, '_mmap', None) + if mapping is not None and not mapping.closed: + mapping.close() + + def restore(self, snapshot): + assignments, metadata, next_id = self._validated_state(**snapshot) + self.assignments, self.metadata = assignments, metadata + self.events = copy.deepcopy(snapshot['events']) + self.next_cluster_id = next_id + + def _append(self, event): + self._validate_event_capabilities([event]) + events = self.events + [event] + _, metadata, next_id = interpret_events(self.unit_ids, self.base_metadata, events) + assignments = project_assignments(self.base_assignments, events, self.attributes) + self.assignments, self.metadata, self.events = assignments, metadata, events + self.next_cluster_id = max(self.next_cluster_id, next_id) + return copy.deepcopy(event) + + def _require_capability(self, operation, predicate=None): + if operation == 'threshold_split': + predicate = validate_predicate(predicate) + if predicate['attribute'] not in self.capabilities['predicate_split_attributes']: + raise ValueError( + 'Predicate split capability is unavailable in this bounded package.' + ) + elif not self.capabilities.get(operation, False): + raise ValueError( + f'{operation.capitalize()} capability is disabled in this bounded package.' + ) + + def _validate_event_capabilities(self, events): + for event in events: + if isinstance(event, dict): + operation = event.get('operation') + if operation in ('merge', 'label', 'threshold_split'): + self._require_capability(operation, event.get('predicate')) + + def preflight_merge(self, parents, target, metadata=None): + self._require_capability('merge') + event = { + 'version': 1, + 'operation': 'merge', + 'parents': [int(i) for i in parents], + 'target': _integer(target, 'Target ID'), + 'metadata': copy.deepcopy(metadata or {}), + } + if event['target'] < self.next_cluster_id: + raise ValueError('Merge target is below the complete catalog allocator.') + interpret_events(self.unit_ids, self.base_metadata, self.events + [event]) + return event + + def record_merge(self, parents, target, metadata=None): + return self._append(self.preflight_merge(parents, target, metadata)) + + def preflight_label(self, field, ids, value): + self._require_capability('label') + event = { + 'version': 1, + 'operation': 'label', + 'field': field, + 'ids': [int(i) for i in ids], + 'value': copy.deepcopy(value), + } + interpret_events(self.unit_ids, self.base_metadata, self.events + [event]) + return event + + def record_label(self, field, ids, value): + return self._append(self.preflight_label(field, ids, value)) + + def record_threshold_split(self, parent, children, predicate, metadata=None): + """Apply an explicitly enabled full-source scalar predicate.""" + self._require_capability('threshold_split', predicate) + children = [int(i) for i in children] + if any(i < self.next_cluster_id for i in children): + raise ValueError('Threshold children collide with the complete catalog allocator.') + return self._append( + { + 'version': 1, + 'operation': 'threshold_split', + 'parent': int(parent), + 'children': children, + 'predicate': validate_predicate(predicate), + 'metadata': copy.deepcopy(metadata or {}), + } + ) + + def preflight_split(self, parent, selected_id, remainder_id, predicate, metadata=None): + self._require_capability('threshold_split', predicate) + children = [int(remainder_id), int(selected_id)] + if any(i < self.next_cluster_id for i in children): + raise ValueError('Threshold children collide with the complete catalog allocator.') + event = { + 'version': 1, + 'operation': 'threshold_split', + 'parent': int(parent), + 'children': children, + 'predicate': validate_predicate(predicate), + 'metadata': copy.deepcopy(metadata or {}), + } + interpret_events(self.unit_ids, self.base_metadata, self.events + [event]) + project_assignments(self.base_assignments, self.events + [event], self.attributes) + return event + + def record_split(self, parent, selected_id, remainder_id, predicate, metadata=None): + return self._append( + self.preflight_split(parent, selected_id, remainder_id, predicate, metadata) + ) + + def save( + self, + assignments=None, + metadata=None, + events=None, + next_cluster_id=None, + extra_documents=None, + ): + """Validate, flush a generation, then atomically replace its commit pointer.""" + assignments = self.assignments if assignments is None else assignments + metadata = self.metadata if metadata is None else metadata + events = self.events if events is None else events + next_cluster_id = self.next_cluster_id if next_cluster_id is None else next_cluster_id + values, metadata, next_id = self._validated_state( + assignments, metadata, events, next_cluster_id + ) + documents = copy.deepcopy( + self.extra_documents if extra_documents is None else extra_documents + ) + if not isinstance(documents, dict): + raise ValueError('Extra checkpoint documents must be a filename-to-JSON mapping.') + for name, value in documents.items(): + _safe_file(self.path, name) + if name in ('checkpoint.json', 'spike_clusters.npy', 'journal.json', 'metadata.json'): + raise ValueError('Extra document collides with a checkpoint state file.') + _json_bytes(value) + root = self.path / '.phy-curation' + generations = root / 'generations' + generations.mkdir(parents=True, exist_ok=True) + _fsync_directory(root) + _fsync_directory(self.path) + with _commit_lock(root): + if _current_pointer(root) != self._expected_pointer: + raise ValueError( + 'Bounded curation checkpoint changed in another session; ' + 'reopen the dataset before saving to avoid losing its changes.' + ) + return self._save_generation(generations, values, metadata, events, next_id, documents) + + def _save_generation(self, generations, values, metadata, events, next_id, documents): + root = generations.parent + identifier = uuid.uuid4().hex + stage = generations / f'.{identifier}.tmp' + final = generations / identifier + stage.mkdir() + pointer_tmp = root / f'.current-{identifier}.tmp' + try: + np.save(stage / 'spike_clusters.npy', values) + # Windows _commit (os.fsync) requires a writable descriptor. + with (stage / 'spike_clusters.npy').open('r+b') as stream: + os.fsync(stream.fileno()) + _write_json(stage / 'journal.json', events) + _write_json(stage / 'metadata.json', metadata) + for name, value in documents.items(): + _write_json(stage / name, value) + descriptor = { + 'version': 1, + 'manifest_sha256': _digest(self.manifest), + 'next_cluster_id': next_id, + 'files': { + name: _hash_file(stage / name) + for name in ('spike_clusters.npy', 'journal.json', 'metadata.json', *documents) + }, + 'extra_documents': sorted(documents), + } + _write_json(stage / 'checkpoint.json', descriptor) + _fsync_directory(stage) + os.replace(stage, final) + _fsync_directory(generations) + pointer = {'generation': identifier, 'sha256': _digest(descriptor)} + _write_json(pointer_tmp, pointer) + os.replace(pointer_tmp, root / 'current.json') + _fsync_directory(root) + except BaseException: + shutil.rmtree(stage, ignore_errors=True) + pointer_tmp.unlink(missing_ok=True) + raise + self.assignments, self.metadata, self.events = values, metadata, copy.deepcopy(events) + self.next_cluster_id, self.checkpoint_path = next_id, final + self.extra_documents = documents + self._expected_pointer = pointer + return final diff --git a/phy/apps/large_recording/model.py b/phy/apps/large_recording/model.py new file mode 100644 index 000000000..9f0de0389 --- /dev/null +++ b/phy/apps/large_recording/model.py @@ -0,0 +1,65 @@ +"""Bounded template loading without dense caches indexed by cluster namespace.""" + +from types import MappingProxyType + +import numpy as np +from phylib.io.model import TemplateModel + + +class BoundedTemplateModel(TemplateModel): + """Load committed assignments and original templates as separate concepts. + + The controller derives live unit evidence from Supervisor assignments. + phylib's eager merged-cluster waveform cache allocates by the maximum + absolute ID and does not track subsequent GUI edits; bounded mode does not + expose that cache as current scientific evidence. + """ + + def __init__(self, *args, session, **kwargs): + self._bounded_session = session + super().__init__(*args, **kwargs) + self.sparse_clusters = None + + def freeze_source(self): + """Prevent public model edits from bypassing semantic curation.""" + names = set() + for name, value in vars(self).items(): + if not name.startswith('_') and isinstance(value, np.ndarray): + value.flags.writeable = False + names.add(name) + self.metadata = MappingProxyType( + {name: MappingProxyType(values) for name, values in self.metadata.items()} + ) + names.update( + ( + 'metadata', + 'save_spike_clusters', + 'save_metadata', + 'save_spikes_subset_waveforms', + 'sample_rate', + 'duration', + 'n_channels', + 'n_channels_dat', + 'n_samples_waveforms', + 'dtype', + 'offset', + ) + ) + self._protected_source_names = frozenset(names) + + def __setattr__(self, name, value): + if name in getattr(self, '_protected_source_names', ()): + raise ValueError('Bounded model source is read-only; use semantic curation commands.') + super().__setattr__(name, value) + + def _load_spike_clusters(self): + return self._bounded_session.assignments.copy() + + def get_merge_map(self): + return {}, np.array([], dtype=np.int64) + + def cluster_waveforms(self): + return None + + def get_cluster_mean_waveforms(self, *args, **kwargs): + raise ValueError('Use the controller waveform/template providers for live bounded units.') diff --git a/phy/apps/large_recording/prepare.py b/phy/apps/large_recording/prepare.py new file mode 100644 index 000000000..4f4b5e3dc --- /dev/null +++ b/phy/apps/large_recording/prepare.py @@ -0,0 +1,526 @@ +"""Publish a self-contained bounded evidence package from a full sorting.""" + +import csv +import json +import os +import shutil +import tempfile +from pathlib import Path + +import numpy as np + +from .source import ( + PER_SPIKE, + STATIC, + WAVEFORMS, + Source, + close_array, + hash_file, + hash_json, + publish_directory, +) +from .summaries import build_summaries, scalar_attributes, summary_size + + +def _row_bounds(times, windows): + return [ + ( + int(np.searchsorted(times, start, side='left')), + int(np.searchsorted(times, end, side='left')), + ) + for start, end in windows + ] + + +def _union(windows, block_starts=()): + result = [] + for start, end in sorted(windows): + if result and ( + start < result[-1][1] or (start == result[-1][1] and start not in block_starts) + ): + result[-1][1] = max(result[-1][1], int(end)) + else: + result.append([int(start), int(end)]) + return result + + +def _write_metadata(path, metadata): + for field, values in metadata.items(): + if ( + not isinstance(field, str) + or not field + or '/' in field + or '\\' in field + or field in ('.', '..') + ): + raise ValueError('Metadata field names must be safe filenames.') + with (path / f'cluster_{field}.tsv').open('w', encoding='utf-8', newline='') as stream: + writer = csv.writer(stream, delimiter='\t') + writer.writerow(['cluster_id', field]) + writer.writerows(sorted(values.items())) + + +def _selected_waveform_count(source, windows): + ids = source.arrays.get('_phy_spikes_subset.spikes.npy') + if ids is None: + return 0 + return sum( + int(np.searchsorted(ids, stop) - np.searchsorted(ids, start)) + for start, stop in _row_bounds(source.arrays['spike_times.npy'].reshape(-1), windows) + ) + + +def prepare_large( + params_path, + output_dir, + duration_seconds=None, + valid_intervals=None, + max_spikes=2000000, + max_bytes=536870912, + window_seconds=30, + stratum_seconds=3600, + seed=0, + batch_size=100000, + bin_seconds=3600, +): + """Choose complete shared windows, write bounded arrays, and publish atomically. + + Survey windows use a seeded offset within each acquisition stratum. A second + pass tries one first-spike rescue window for every unit without evidence. + Candidate windows that exceed either budget are explicitly omitted; spikes + within accepted windows are never independently thinned. + """ + with Source(params_path, duration_seconds, valid_intervals, batch_size) as source: + return _prepare_source( + source, + output_dir, + max_spikes, + max_bytes, + window_seconds, + stratum_seconds, + seed, + batch_size, + bin_seconds, + ) + + +def _prepare_source( + source, + output_dir, + max_spikes, + max_bytes, + window_seconds, + stratum_seconds, + seed, + batch_size, + bin_seconds, +): + output_dir = Path(output_dir).resolve() + missing_static = { + 'templates.npy', + 'channel_map.npy', + 'channel_positions.npy', + } - source.arrays.keys() + if missing_static: + raise ValueError(f'Bounded GUI preparation requires: {", ".join(sorted(missing_static))}.') + if output_dir.exists(): + raise ValueError(f'Output already exists: {output_dir}.') + if output_dir == source.path or source.path in output_dir.parents: + raise ValueError('Bounded output must be outside the source dataset directory.') + if int(max_spikes) != max_spikes or max_spikes <= 0 or max_bytes <= 0: + raise ValueError('max_spikes and max_bytes must be positive budgets.') + widths = [] + for label, seconds in ( + ('window_seconds', window_seconds), + ('stratum_seconds', stratum_seconds), + ): + width = float(seconds) * source.sample_rate + if ( + not np.isfinite(width) + or width < 1 + or not np.isclose(width, round(width), rtol=0, atol=1e-6) + ): + raise ValueError(f'{label} must specify a positive integer sample width.') + widths.append(int(round(width))) + width, stratum = widths + times = source.arrays['spike_times.npy'].reshape(-1) + assignments = source.arrays['spike_clusters.npy'].reshape(-1) + summary_bytes = summary_size(source, bin_seconds)[0] + producer_documents = [ + source.path / name + for name in ('curation.json', 'curation_review.json') + if (source.path / name).is_file() + ] + fixed_bytes = ( + sum(path.stat().st_size for path in producer_documents) + + summary_bytes + + sum((source.path / name).stat().st_size for name in source.arrays.keys() & STATIC) + ) + # Reserve headers, manifest, params, metadata, and window records explicitly. + metadata_bytes = sum( + len(str(value).encode('utf-8')) + 64 + for values in source.metadata.values() + for value in values.values() + ) + overhead = 32768 + metadata_bytes + 1024 * len(source.arrays) + fixed_bytes += overhead + row_bytes = 16 # int64 global row mapping and baseline assignment + for name in source.arrays.keys() & PER_SPIKE: + array = source.arrays[name] + row_bytes += array.dtype.itemsize * int(np.prod(array.shape[1:], dtype=np.int64)) + waveform_bytes = 0 + if source.arrays.keys() >= WAVEFORMS: + for name in WAVEFORMS: + array = source.arrays[name] + waveform_bytes += array.dtype.itemsize * int(np.prod(array.shape[1:], dtype=np.int64)) + if fixed_bytes > max_bytes: + raise ValueError( + f'Output budget exceeded before selecting evidence: estimated ' + f'{fixed_bytes} bytes; increase max_bytes or bin_seconds.' + ) + windows, plan, omitted = [], [], [] + rng = np.random.default_rng(seed) + + def accept(start, end, kind, segment, stratum_id): + candidate = _union([*windows, [start, end]], set(source.valid_intervals[:, 0])) + n = sum(stop - first for first, stop in _row_bounds(times, candidate)) + n_waveforms = _selected_waveform_count(source, candidate) + size = fixed_bytes + n * row_bytes + n_waveforms * waveform_bytes + len(candidate) * 256 + record = { + 'start': int(start), + 'end': int(end), + 'kind': kind, + 'segment': 0, + 'block': int(segment), + 'stratum': int(stratum_id), + } + if n > max_spikes or size > max_bytes: + omitted.append( + {**record, 'reason': 'spike_budget' if n > max_spikes else 'byte_budget'} + ) + return False + windows[:] = candidate + plan.append(record) + return True + + # Candidate enumeration is scalar and bounded by output-window metadata, + # rather than constructing recording-length arrays. + for segment, (first, last) in enumerate(source.valid_intervals): + for start in range(int(first), int(last), stratum): + end = min(start + stratum, int(last)) + length = min(width, end - start) + offset = int(rng.integers(0, end - start - length + 1)) + accept(start + offset, start + offset + length, 'survey', segment, start // stratum) + if (len(plan) + len(omitted)) * 256 + fixed_bytes > max_bytes: + raise ValueError( + 'Window plan metadata exceeds max_bytes; increase stratum_seconds.' + ) + evidence = set() + for first, stop in _row_bounds(times, windows): + for start in range(first, stop, batch_size): + evidence.update(np.unique(assignments[start : min(start + batch_size, stop)]).tolist()) + missing = set(source.unit_ids.tolist()) - evidence + for batch in source.iter_batches(): + if not missing: + break + units, positions = np.unique(batch['spike_clusters.npy'], return_index=True) + for unit, position in zip(units.tolist(), positions.tolist(), strict=True): + if unit not in missing: + continue + sample = int(batch['spike_times.npy'][position]) + segment = int(np.searchsorted(source.valid_intervals[:, 0], sample, side='right') - 1) + first, last = source.valid_intervals[segment] + start = max(int(first), min(sample - width // 2, int(last) - width)) + end = min(start + width, int(last)) + accept(start, end, 'rescue', segment, -1) + # One deterministic rescue attempt per unit, including rejected attempts. + missing.remove(unit) + row_bounds = _row_bounds(times, windows) + n_bounded = sum(stop - first for first, stop in row_bounds) + if not n_bounded: + raise ValueError( + 'No bounded evidence fits the budgets; increase max_spikes/max_bytes ' + 'or shorten window_seconds (tied spikes remain together).' + ) + output_dir.parent.mkdir(parents=True, exist_ok=True) + estimated_bytes = ( + fixed_bytes + + n_bounded * row_bytes + + _selected_waveform_count(source, windows) * waveform_bytes + + len(windows) * 256 + ) + free = shutil.disk_usage(output_dir.parent).free + if free < estimated_bytes: + raise OSError( + f'Insufficient disk space for bounded preparation: estimated {estimated_bytes:,} bytes, {free:,} available.' + ) + temporary = Path( + tempfile.mkdtemp(prefix=f'.{output_dir.name}.prepare-', dir=output_dir.parent) + ) + global_ids, base, destinations = None, None, {} + try: + global_ids = np.lib.format.open_memmap( + temporary / '_phy_spike_global_ids.npy', mode='w+', dtype=np.int64, shape=(n_bounded,) + ) + for name in source.arrays.keys() & PER_SPIKE: + array = source.arrays[name] + destinations[name] = np.lib.format.open_memmap( + temporary / name, mode='w+', dtype=array.dtype, shape=(n_bounded, *array.shape[1:]) + ) + base = np.lib.format.open_memmap( + temporary / '_phy_spike_clusters_base.npy', + mode='w+', + dtype=np.int64, + shape=(n_bounded,), + ) + cursor = 0 + bounded_counts = {int(unit): 0 for unit in source.unit_ids} + for first, stop in row_bounds: + for start in range(first, stop, batch_size): + end = min(start + batch_size, stop) + length = end - start + global_ids[cursor : cursor + length] = np.arange(start, end, dtype=np.int64) + base[cursor : cursor + length] = assignments[start:end] + for name, target in destinations.items(): + target[cursor : cursor + length] = source.arrays[name][start:end] + ids, counts = np.unique(assignments[start:end], return_counts=True) + for unit, count in zip(ids.tolist(), counts.tolist(), strict=True): + bounded_counts[unit] += count + cursor += length + for array in (global_ids, base, *destinations.values()): + array.flush() + for name in source.arrays.keys() & STATIC: + shutil.copyfile(source.path / name, temporary / name) + if source.arrays.keys() >= WAVEFORMS: + _copy_waveforms(source, temporary, global_ids, row_bounds) + waveform_counts = {int(unit): 0 for unit in source.unit_ids} + waveform_ids_path = temporary / '_phy_spikes_subset.spikes.npy' + if waveform_ids_path.exists(): + waveform_ids = np.load(waveform_ids_path, mmap_mode='r') + for start in range(0, len(waveform_ids), batch_size): + ids, counts = np.unique( + base[waveform_ids[start : start + batch_size]], return_counts=True + ) + for unit, count in zip(ids, counts): + waveform_counts[int(unit)] += int(count) + waveform_ids._mmap.close() + np.save(temporary / '_phy_curation_windows.npy', np.asarray(windows, dtype=np.int64)) + _write_metadata(temporary, source.metadata) + (temporary / '_phy_cluster_metadata.json').write_text( + json.dumps(source.metadata, sort_keys=True) + '\n', encoding='utf-8' + ) + # Standard loader stays self-contained; no original raw path is referenced. + params = { + key: source.params[key] + for key in ( + 'sample_rate', + 'n_channels_dat', + 'dtype', + 'offset', + 'hp_filtered', + 'acquisition_id', + 'depth_axis', + 'amplitude_units', + 'depth_units', + 'position_units', + 'waveform_units', + ) + if key in source.params + } + params.update( + dat_path=None, + sample_rate=source.sample_rate, + num_samples=source.num_samples, + valid_intervals=source.valid_intervals.tolist(), + ) + params.setdefault('dtype', 'int16') + params.setdefault('n_channels_dat', len(source.arrays.get('channel_map.npy', []))) + (temporary / 'params.py').write_text( + ''.join(f'{key} = {value!r}\n' for key, value in params.items()), encoding='utf-8' + ) + summaries = build_summaries(source, temporary, bin_seconds, max_bytes=max_bytes) + final_windows = [] + for start, end in windows: + contributors = [ + record for record in plan if record['start'] < end and record['end'] > start + ] + kinds = sorted({record['kind'] for record in contributors}) + final_windows.append( + { + 'start': start, + 'end': end, + 'segment': 0, + 'block': contributors[0]['block'], + 'kind': kinds[0] if len(kinds) == 1 else 'mixed', + 'kinds': kinds, + 'strata': sorted({record['stratum'] for record in contributors}), + 'exposure_seconds': (end - start) / source.sample_rate, + } + ) + manifest = { + 'format': 'phy_bounded_curation', + 'format_version': 1, + 'source': { + 'identity': source.identity, + 'revision': source.revision, + 'path': str(source.path), + }, + 'sample_rate': source.sample_rate, + 'num_samples': source.num_samples, + 'duration_seconds': source.num_samples / source.sample_rate, + 'valid_intervals': source.valid_intervals.tolist(), + 'num_spikes': source.n_spikes, + 'bounded_num_spikes': n_bounded, + 'unit_ids': source.unit_ids.tolist(), + 'counts': source.counts.tolist(), + 'bounded_counts': [bounded_counts[int(unit)] for unit in source.unit_ids], + 'waveform_counts': [waveform_counts[int(unit)] for unit in source.unit_ids], + 'next_cluster_id': source.next_cluster_id, + 'global_spike_ids_file': '_phy_spike_global_ids.npy', + 'base_assignments_file': '_phy_spike_clusters_base.npy', + 'windows_file': '_phy_curation_windows.npy', + 'windows': final_windows, + 'accepted_window_candidates': plan, + 'omitted_windows': omitted, + 'metadata': source.metadata, + 'base_metadata': source.metadata, + 'summaries': summaries, + **{key: value for key, value in summaries.items() if key.endswith('_file')}, + 'selection': { + 'algorithm': 'shared_windows_v1', + 'seed': int(seed), + 'window_seconds': window_seconds, + 'stratum_seconds': stratum_seconds, + }, + 'budgets': { + 'max_spikes': int(max_spikes), + 'max_bytes': int(max_bytes), + 'batch_size': int(batch_size), + 'summary_bytes': int(summary_bytes), + 'estimated_bytes': int( + fixed_bytes + + n_bounded * row_bytes + + _selected_waveform_count(source, windows) * waveform_bytes + + len(windows) * 256 + ), + }, + 'capabilities': { + 'merge': True, + 'label': True, + 'split': False, + 'full_attributes': sorted(scalar_attributes(source)), + 'predicate_split_attributes': ( + ['amplitudes.npy'] if 'amplitude' in scalar_attributes(source) else [] + ), + 'stored_waveforms': bool(sum(waveform_counts.values())), + 'raw_traces': False, + 'normalized_correlograms': False, + }, + 'source_array_inventory': { + name: { + 'role': 'per_spike' + if name in PER_SPIKE + else 'stored_waveform' + if name in WAVEFORMS + else 'static' + if name in STATIC + else 'excluded_window_sidecar', + 'shape': list(array.shape), + 'dtype': str(array.dtype), + 'bytes': int(array.nbytes), + } + for name, array in source.arrays.items() + }, + 'excluded_files': sorted( + p.name + for p in source.path.iterdir() + if p.name not in (PER_SPIKE | STATIC | WAVEFORMS) + and p.name + not in ( + 'params.py', + '_phy_cluster_metadata.json', + 'curation.json', + 'curation_review.json', + ) + and (not p.name.startswith('cluster_') or p.stem == 'cluster_info') + ), + } + for document in producer_documents: + shutil.copyfile(document, temporary / document.name) + bound_paths = [ + *temporary.glob('*.npy'), + temporary / 'params.py', + temporary / '_phy_cluster_metadata.json', + ] + manifest['bounded_files'] = {p.name: hash_file(p) for p in sorted(bound_paths)} + manifest['preparation_identity'] = hash_json(manifest) + (temporary / '_phy_curation_manifest.json').write_text( + json.dumps(manifest, indent=2, sort_keys=True) + '\n', encoding='utf-8' + ) + actual_bytes = sum(p.stat().st_size for p in temporary.iterdir()) + if actual_bytes > max_bytes: + raise ValueError(f'Final output exceeds max_bytes: {actual_bytes} bytes.') + for array in (global_ids, base, *destinations.values()): + close_array(array) + for path in temporary.iterdir(): + if path.is_file(): + # Windows _commit (os.fsync) requires a writable descriptor. + with path.open('r+b') as stream: + os.fsync(stream.fileno()) + _fsync_directory(temporary) + if output_dir.exists(): + raise ValueError(f'Output already exists: {output_dir}.') + from .curation import BoundedSession + + checked = BoundedSession(temporary) + checked.close() + publish_directory(temporary, output_dir) + _fsync_directory(output_dir.parent) + except BaseException: + for array in (global_ids, base, *destinations.values()): + close_array(array) + shutil.rmtree(temporary, ignore_errors=True) + raise + return output_dir + + +def _fsync_directory(path): + if os.name != 'nt': + descriptor = os.open(path, os.O_RDONLY) + try: + os.fsync(descriptor) + finally: + os.close(descriptor) + + +def _copy_waveforms(source, output, global_ids, row_bounds): + """Intersect sparse row IDs with selected spans in bounded working memory.""" + upstream_ids = source.arrays['_phy_spikes_subset.spikes.npy'] + spans = [ + (int(np.searchsorted(upstream_ids, first)), int(np.searchsorted(upstream_ids, stop))) + for first, stop in row_bounds + ] + n = sum(stop - first for first, stop in spans) + targets = {} + try: + for name in WAVEFORMS: + original = source.arrays[name] + targets[name] = np.lib.format.open_memmap( + output / name, mode='w+', dtype=original.dtype, shape=(n, *original.shape[1:]) + ) + cursor = 0 + for first, stop in spans: + for start in range(first, stop, source.batch_size): + end = min(start + source.batch_size, stop) + length = end - start + targets['_phy_spikes_subset.spikes.npy'][cursor : cursor + length] = ( + np.searchsorted(global_ids, upstream_ids[start:end]) + ) + for name in WAVEFORMS - {'_phy_spikes_subset.spikes.npy'}: + targets[name][cursor : cursor + length] = source.arrays[name][start:end] + cursor += length + for array in targets.values(): + array.flush() + finally: + for array in targets.values(): + close_array(array) diff --git a/phy/apps/large_recording/replay.py b/phy/apps/large_recording/replay.py new file mode 100644 index 000000000..96c39ba8c --- /dev/null +++ b/phy/apps/large_recording/replay.py @@ -0,0 +1,230 @@ +"""Stream semantic curation into a new complete sorting, never the source.""" + +import csv +import os +import shutil +import tempfile +from pathlib import Path + +import numpy as np + +from .curation import ( + BoundedSession, + _digest, + _fsync_directory, + _write_json, + interpret_events, + project_assignments, +) +from .source import PER_SPIKE, Source, close_array, publish_directory + + +def _merge_mapping(unit_ids, events): + """Resolve a whole-unit merge journal once, independently of spike count.""" + lineage = {int(unit): [int(unit)] for unit in unit_ids} + for event in events: + if event['operation'] == 'merge': + members = [] + for parent in event['parents']: + members.extend(lineage.pop(parent)) + lineage[event['target']] = members + mapping = {old: new for new, members in lineage.items() for old in members} + keys = np.array(sorted(mapping), dtype=np.int64) + targets = np.array([mapping[int(key)] for key in keys], dtype=np.int64) + return keys, targets + + +def _write_metadata(path, metadata): + for field, values in sorted(metadata.items()): + if not field or any(char in field for char in ('/', '\\', '\0')): + raise ValueError('Metadata fields must be safe local filenames.') + with (path / f'cluster_{field}.tsv').open('w', encoding='utf8', newline='') as stream: + writer = csv.writer(stream, delimiter='\t') + writer.writerow(['cluster_id', field]) + for cluster, value in sorted(values.items()): + if value is not None: + writer.writerow([cluster, value]) + stream.flush() + os.fsync(stream.fileno()) + _write_json(path / '_phy_cluster_metadata.json', metadata) + + +def apply_curation(source_params, bounded_dir, output_dir, batch_size=100000): + """Validate identity/revision and replay only the last committed generation. + + The output is published as one new directory after validation. Full source + arrays remain file-backed; only one assignment batch and the unit map reside + in memory. Arbitrary selected-row splits have no journal interpretation. + """ + session = BoundedSession(bounded_dir) + source = None + try: + if session.checkpoint_path is None: + raise ValueError('Save a bounded semantic checkpoint before applying curation.') + manifest = session.manifest + source = Source( + source_params, + duration_seconds=manifest['num_samples'] / manifest['sample_rate'], + valid_intervals=manifest['valid_intervals'], + batch_size=batch_size, + ) + return _apply_curation(source, session, output_dir, batch_size) + finally: + if source is not None: + source.close() + session.close() + + +def _apply_curation(source, session, output_dir, batch_size): + manifest = session.manifest + expected = manifest['source'] + if source.identity != expected['identity']: + raise ValueError('Source coordinate/feature identity does not match bounded preparation.') + if source.revision != expected['revision']: + raise ValueError('Source assignment revision does not match journal parent revision.') + if ( + source.unit_ids.tolist() != manifest['unit_ids'] + or source.n_spikes != manifest['num_spikes'] + ): + raise ValueError('Source catalog or spike count does not match bounded preparation.') + # Identity/revision bind the complete source. Validate the actual projection + # too, so a tampered bounded base cannot silently replay unrelated evidence. + source_assignments = source.arrays['spike_clusters.npy'].reshape(-1) + for start in range(0, len(session.global_ids), batch_size): + rows = session.global_ids[start : start + batch_size] + if not np.array_equal( + source_assignments[rows], + session.base_assignments[start : start + batch_size], + ): + raise ValueError('Immutable bounded assignments do not match source rows.') + for name in sorted(source.arrays.keys() & (PER_SPIKE - {'spike_clusters.npy'})): + bounded = np.load(session.path / name, mmap_mode='r', allow_pickle=False) + try: + if bounded.shape != (len(session.global_ids), *source.arrays[name].shape[1:]): + raise ValueError(f'Bounded attribute has invalid source alignment: {name}.') + for start in range(0, len(session.global_ids), batch_size): + rows = session.global_ids[start : start + batch_size] + if not np.array_equal( + source.arrays[name][rows], + bounded[start : start + batch_size], + equal_nan=True, + ): + raise ValueError(f'Bounded attribute does not match source rows: {name}.') + finally: + close_array(bounded) + active_units, metadata, _ = interpret_events(source.unit_ids, source.metadata, session.events) + output = Path(output_dir).resolve() + if output.exists(): + raise FileExistsError(f'Curated output already exists: {output}.') + if output == source.path or source.path in output.parents: + raise ValueError('Curated output must be outside the immutable source directory.') + output.parent.mkdir(parents=True, exist_ok=True) + estimated_bytes = sum( + (source.path / name).stat().st_size + for name in source.arrays + if name != 'spike_clusters.npy' + ) + estimated_bytes += source.n_spikes * 8 + 1024 * 1024 + free = shutil.disk_usage(output.parent).free + if free < estimated_bytes: + raise OSError( + f'Insufficient disk space for full replay: estimated at least {estimated_bytes:,} bytes, {free:,} available.' + ) + stage = Path(tempfile.mkdtemp(prefix=f'.{output.name}-', dir=output.parent)) + curated = None + assignments = None + try: + # Copy only validated canonical assets, not GUI caches or producer reviews. + for name in source.arrays: + if name != 'spike_clusters.npy': + shutil.copyfile(source.path / name, stage / name) + shutil.copyfile(source.params_path, stage / 'params.py') + upstream_manifest = source.path / '_phy_spikes_subset.manifest.json' + if upstream_manifest.exists(): + shutil.copyfile(upstream_manifest, stage / upstream_manifest.name) + # Preserve source params semantics and make external raw paths absolute. + # The full export can contain no raw traces; immutable array rows suffice. + dat_path = source.params.get('dat_path') + if dat_path: + raw_paths = dat_path if isinstance(dat_path, (list, tuple)) else [dat_path] + resolved = [str((source.path / raw).resolve()) for raw in raw_paths] + dat_path = resolved if isinstance(dat_path, (list, tuple)) else resolved[0] + with (stage / 'params.py').open('a', encoding='utf8') as stream: + stream.write(f'\nnum_samples = {source.num_samples!r}\n') + stream.write(f'valid_intervals = {source.valid_intervals.tolist()!r}\n') + stream.write(f'dat_path = {dat_path!r}\n') + assignments = np.lib.format.open_memmap( + stage / 'spike_clusters.npy', + mode='w+', + dtype=np.int64, + shape=source.arrays['spike_clusters.npy'].shape, + ) + flat = assignments.reshape(-1) + has_split = any(event['operation'] == 'threshold_split' for event in session.events) + keys, targets = ( + (None, None) if has_split else _merge_mapping(source.unit_ids, session.events) + ) + for batch in source.iter_batches(): + values = batch['spike_clusters.npy'] + if has_split: + projected = project_assignments(values, session.events, batch) + else: + positions = np.searchsorted(keys, values) + projected = targets[positions] + flat[batch['start'] : batch['stop']] = projected + assignments.flush() + assignments._mmap.close() + del flat, assignments + assignments = None + _write_metadata(stage, metadata) + _write_json(stage / '_phy_curation_journal.json', session.events) + # Re-reading verifies alignment, full catalog and immutable byte digests. + curated = Source(stage / 'params.py', batch_size=batch_size) + if curated.identity != source.identity or curated.n_spikes != source.n_spikes: + raise ValueError('Curated output changed immutable source coordinates.') + if curated.metadata != metadata: + raise ValueError('Curated output metadata failed exact round-trip validation.') + actual = curated.arrays['spike_clusters.npy'].reshape(-1) + for batch in source.iter_batches(): + if has_split: + expected_batch = project_assignments( + batch['spike_clusters.npy'], session.events, batch + ) + else: + expected_batch = targets[np.searchsorted(keys, batch['spike_clusters.npy'])] + if not np.array_equal(actual[batch['start'] : batch['stop']], expected_batch): + raise ValueError('Curated output assignments failed semantic replay verification.') + _write_json( + stage / '_phy_curation_lineage.json', + { + 'format': 'phy_curation_lineage', + 'format_version': 1, + 'coordinate_identity': source.identity, + 'parent_revision': source.revision, + 'journal_sha256': _digest(session.events), + 'revision': curated.revision, + 'preparation_sha256': _digest(manifest), + 'next_cluster_id': session.next_cluster_id, + 'unit_ids': sorted(active_units), + }, + ) + curated.close() + for path in stage.iterdir(): + if path.is_file(): + # Windows _commit (os.fsync) requires a writable descriptor. + with path.open('r+b') as stream: + os.fsync(stream.fileno()) + _fsync_directory(stage) + # Refuse replacement even when a racing caller creates the target. + if output.exists(): + raise FileExistsError(f'Curated output already exists: {output}.') + publish_directory(stage, output) + _fsync_directory(output.parent) + except BaseException: + if curated is not None: + curated.close() + if assignments is not None: + assignments._mmap.close() + shutil.rmtree(stage, ignore_errors=True) + raise + return output diff --git a/phy/apps/large_recording/source.py b/phy/apps/large_recording/source.py new file mode 100644 index 000000000..07bacf21f --- /dev/null +++ b/phy/apps/large_recording/source.py @@ -0,0 +1,570 @@ +"""Validated, file-backed source sorting on one monotonic sample clock.""" + +import csv +import ctypes +import errno +import hashlib +import json +import os +import runpy +import sys +from pathlib import Path + +import numpy as np + +PER_SPIKE = frozenset( + ( + 'spike_times.npy', + 'spike_clusters.npy', + 'spike_templates.npy', + 'amplitudes.npy', + 'pc_features.npy', + 'template_features.npy', + 'spike_positions.npy', + 'spike_depths.npy', + ) +) +STATIC = frozenset( + ( + 'templates.npy', + 'templates_ind.npy', + 'template_ind.npy', + 'channel_map.npy', + 'channel_positions.npy', + 'channel_shanks.npy', + 'channel_map_si.npy', + 'channel_groups.npy', + 'whitening_mat.npy', + 'whitening_mat_inv.npy', + 'similar_templates.npy', + 'pc_feature_ind.npy', + 'template_feature_ind.npy', + ) +) +WAVEFORMS = frozenset( + ( + '_phy_spikes_subset.spikes.npy', + '_phy_spikes_subset.channels.npy', + '_phy_spikes_subset.waveforms.npy', + ) +) +WINDOW_SIDECARS = frozenset( + ( + '_phy_spikes_subset.window_bounds.npy', + '_phy_spikes_subset.window_types.npy', + '_phy_spikes_subset.window_strata.npy', + ) +) + + +def hash_file(path, batch_bytes=1024 * 1024): + """Hash bytes incrementally, including the array header.""" + digest = hashlib.sha256() + with Path(path).open('rb') as stream: + for block in iter(lambda: stream.read(batch_bytes), b''): + digest.update(block) + return digest.hexdigest() + + +def hash_json(value): + """Hash canonical JSON, normalizing integer metadata dictionary keys.""" + + def normalize(item): + if isinstance(item, dict): + return {str(key): normalize(val) for key, val in item.items()} + if isinstance(item, (list, tuple)): + return [normalize(val) for val in item] + if isinstance(item, np.generic): + return item.item() + return item + + payload = json.dumps(normalize(value), sort_keys=True, separators=(',', ':'), allow_nan=False) + return hashlib.sha256(payload.encode('utf-8')).hexdigest() + + +def assignment_revision(path, metadata): + """Identify mutable assignments and resolved cluster metadata.""" + path = Path(path) + if path.is_dir(): + path = path / 'spike_clusters.npy' + return hash_json({'assignments': hash_file(path), 'metadata': metadata}) + + +def close_array(array): + """Release a file-backed array handle, including on Windows.""" + mapping = getattr(array, '_mmap', None) + if mapping is not None and not mapping.closed: + mapping.close() + + +def publish_directory(staged, destination): + """Atomically publish a complete directory without replacing any target. + + Ordinary POSIX rename can replace an empty destination directory. Native + exclusive rename prevents a racing publication from silently replacing it. + """ + if os.name == 'nt': + os.rename(staged, destination) + return + library = ctypes.CDLL(None, use_errno=True) + if sys.platform.startswith('linux'): + rename = getattr(library, 'renameat2', None) + if rename is None: + raise OSError(errno.ENOSYS, 'Atomic exclusive directory publication is unavailable.') + rename.argtypes = [ + ctypes.c_int, + ctypes.c_char_p, + ctypes.c_int, + ctypes.c_char_p, + ctypes.c_uint, + ] + rename.restype = ctypes.c_int + result = rename(-100, os.fsencode(staged), -100, os.fsencode(destination), 1) + elif sys.platform == 'darwin': + rename = getattr(library, 'renamex_np', None) + if rename is None: + raise OSError(errno.ENOSYS, 'Atomic exclusive directory publication is unavailable.') + rename.argtypes = [ctypes.c_char_p, ctypes.c_char_p, ctypes.c_uint] + rename.restype = ctypes.c_int + result = rename(os.fsencode(staged), os.fsencode(destination), 4) + else: + raise OSError(errno.ENOSYS, 'Atomic exclusive directory publication is unavailable.') + if result: + code = ctypes.get_errno() + raise OSError(code, os.strerror(code), str(destination)) + + +def _vector(array, name): + if array.ndim == 1 or (array.ndim == 2 and array.shape[1] == 1): + return array.reshape(-1) + raise ValueError(f'{name} must be a vector or column vector.') + + +def _metadata(path): + sidecar = path / '_phy_cluster_metadata.json' + if sidecar.exists(): + metadata = json.loads(sidecar.read_text(encoding='utf-8')) + if not isinstance(metadata, dict): + raise ValueError('Cluster metadata JSON must be a field mapping.') + result = {} + for field, values in metadata.items(): + if not isinstance(values, dict): + raise ValueError('Cluster metadata fields must map integer cluster IDs to values.') + result[field] = {int(key): value for key, value in values.items()} + if any(key < 0 for key in result[field]): + raise ValueError('Cluster metadata IDs must be nonnegative.') + return result + result = {} + for filename in sorted([*path.glob('cluster_*.tsv'), *path.glob('cluster_*.csv')]): + # cluster_info is generated GUI metrics, not a source label table. + if filename.stem == 'cluster_info': + continue + with filename.open(encoding='utf-8', newline='') as stream: + reader = csv.DictReader(stream, delimiter='\t' if filename.suffix == '.tsv' else ',') + if not reader.fieldnames or 'cluster_id' not in reader.fieldnames: + raise ValueError(f'{filename.name} requires a cluster_id column.') + for row in reader: + cluster_id = int(row['cluster_id']) + if cluster_id < 0: + raise ValueError('Cluster metadata IDs must be nonnegative.') + for field, value in row.items(): + if field != 'cluster_id': + result.setdefault(field, {})[cluster_id] = value + return result + + +class Source: + """Read source arrays in batches without constructing a GUI model. + + Duration must come from producer metadata or an explicit argument: the last + spike is never used as acquisition exposure. Intervals are half-open samples. + """ + + def __init__( + self, params_path, duration_seconds=None, valid_intervals=None, batch_size=100000 + ): + self.arrays = {} + try: + self._initialize(params_path, duration_seconds, valid_intervals, batch_size) + except BaseException: + self.close() + raise + + def _initialize( + self, params_path, duration_seconds=None, valid_intervals=None, batch_size=100000 + ): + self.params_path = Path(params_path).resolve() + self.path = self.params_path.parent + self.params = { + key: value + for key, value in runpy.run_path(str(self.params_path)).items() + if not key.startswith('__') + } + self.batch_size = int(batch_size) + if self.batch_size <= 0: + raise ValueError('batch_size must be positive.') + if 'dtype' in self.params: + try: + raw_dtype = np.dtype(self.params['dtype']) + except (TypeError, ValueError) as error: + raise ValueError('params.py dtype must be a valid raw numeric dtype.') from error + if not np.issubdtype(raw_dtype, np.number): + raise ValueError('params.py dtype must be a valid raw numeric dtype.') + self.sample_rate = float(self.params.get('sample_rate', 0)) + if not np.isfinite(self.sample_rate) or self.sample_rate <= 0: + raise ValueError('A positive finite sample_rate is required.') + upstream = self.path / '_phy_spikes_subset.manifest.json' + recording = {} + info = {} + if upstream.exists(): + info = json.loads(upstream.read_text(encoding='utf-8')) + if info.get('format') != 'phy_spikes_subset' or info.get('format_version') != 1: + raise ValueError('Unsupported waveform manifest format or version.') + recording = info.get('recording', {}) + upstream_rate = recording.get('sampling_frequency_hz', self.sample_rate) + if upstream_rate != self.sample_rate: + raise ValueError('Waveform manifest sampling rate disagrees with params.py.') + if self.params.get('num_segments', recording.get('num_segments', 1)) != 1: + raise ValueError('Only a single monotonic acquisition segment is supported.') + duration_samples = None + if duration_seconds is not None: + duration_samples = float(duration_seconds) * self.sample_rate + if ( + not np.isfinite(duration_samples) + or duration_samples <= 0 + or duration_samples > np.iinfo(np.int64).max + or not np.isclose(duration_samples, round(duration_samples), rtol=0, atol=1e-6) + ): + raise ValueError('duration_seconds must specify a positive integer sample count.') + duration_samples = int(round(duration_samples)) + if ( + 'num_samples' in self.params + and 'num_samples' in recording + and self.params['num_samples'] != recording['num_samples'] + ): + raise ValueError( + 'Producer sample counts disagree between params and waveform manifest.' + ) + declared_samples = self.params.get('num_samples', recording.get('num_samples')) + if declared_samples is None: + declared_samples = duration_samples + elif duration_samples is not None and int(declared_samples) != duration_samples: + raise ValueError('Explicit duration disagrees with acquisition metadata.') + if declared_samples is None: + raise ValueError( + 'Acquisition exposure is missing; provide duration_seconds or num_samples.' + ) + self.num_samples = int(declared_samples) + if ( + self.num_samples <= 0 + or self.num_samples > np.iinfo(np.int64).max + or self.num_samples != declared_samples + ): + raise ValueError('num_samples must be a positive int64 sample count.') + if ( + 'valid_intervals' in self.params + and 'valid_intervals' in recording + and not np.array_equal( + np.asarray(self.params['valid_intervals']), + np.asarray(recording['valid_intervals']), + ) + ): + raise ValueError( + 'Producer valid_intervals disagree between params and waveform manifest.' + ) + declared_intervals = self.params.get('valid_intervals', recording.get('valid_intervals')) + if valid_intervals is not None and declared_intervals is not None: + if not np.array_equal(np.asarray(valid_intervals), np.asarray(declared_intervals)): + raise ValueError('Explicit valid_intervals disagree with acquisition metadata.') + if valid_intervals is None: + valid_intervals = declared_intervals + intervals = [[0, self.num_samples]] if valid_intervals is None else valid_intervals + raw_intervals = np.asarray(intervals) + if raw_intervals.ndim != 2 or raw_intervals.shape[1] != 2 or not len(raw_intervals): + raise ValueError('valid_intervals must be nonempty sample [start, end) pairs.') + self.valid_intervals = raw_intervals.astype(np.int64) + if ( + not np.array_equal(raw_intervals, self.valid_intervals) + or np.any(self.valid_intervals[:, 0] < 0) + or np.any(self.valid_intervals[:, 1] > self.num_samples) + or np.any(self.valid_intervals[:, 1] <= self.valid_intervals[:, 0]) + or np.any(self.valid_intervals[1:, 0] < self.valid_intervals[:-1, 1]) + ): + raise ValueError( + 'valid_intervals must be ordered, disjoint integer acquisition samples.' + ) + names = {p.name for p in self.path.glob('*.npy')} + unknown = names - PER_SPIKE - STATIC - WAVEFORMS - WINDOW_SIDECARS + # Replay lineage and original optional sidecars are retained; bounded + # outputs deliberately cannot be used as full-source sorting inputs. + if unknown: + raise ValueError(f'Unregistered source arrays: {", ".join(sorted(unknown))}.') + self.arrays = { + name: np.load(self.path / name, mmap_mode='r', allow_pickle=False) + for name in sorted(names) + } + for required in ('spike_times.npy', 'spike_clusters.npy', 'spike_templates.npy'): + if required not in self.arrays: + raise ValueError(f'Missing required source array {required}.') + self.n_spikes = len(_vector(self.arrays['spike_times.npy'], 'spike_times.npy')) + self.metadata = _metadata(self.path) + self.validate() + counts = {} + for batch in self.iter_batches(): + ids, values = np.unique(batch['spike_clusters.npy'], return_counts=True) + for cluster_id, count in zip(ids.tolist(), values.tolist(), strict=True): + counts[cluster_id] = counts.get(cluster_id, 0) + count + for values in self.metadata.values(): + for cluster_id in values: + counts.setdefault(cluster_id, 0) + lineage_path = self.path / '_phy_curation_lineage.json' + lineage = {} + if lineage_path.exists(): + lineage = json.loads(lineage_path.read_text(encoding='utf-8')) + if ( + lineage.get('format') != 'phy_curation_lineage' + or lineage.get('format_version') != 1 + ): + raise ValueError('Unsupported curation lineage format or version.') + for unit in lineage.get('unit_ids', []): + if not isinstance(unit, int) or unit < 0: + raise ValueError('Lineage unit IDs must be nonnegative integers.') + counts.setdefault(unit, 0) + if any(unit < 0 or unit > np.iinfo(np.int32).max for unit in counts): + raise ValueError( + 'Unit IDs must fit nonnegative int32 for the current Phy loader; spike sample/row coordinates remain int64.' + ) + self.unit_ids = np.array(sorted(counts), dtype=np.int64) + minimum_next = int(self.unit_ids[-1]) + 1 if len(self.unit_ids) else 0 + self.next_cluster_id = lineage.get('next_cluster_id', minimum_next) + if not isinstance(self.next_cluster_id, int) or self.next_cluster_id < minimum_next: + raise ValueError('Lineage next_cluster_id collides with the source unit catalog.') + self.counts = np.array( + [counts[cluster_id] for cluster_id in self.unit_ids], dtype=np.int64 + ) + immutable = sorted(names - {'spike_clusters.npy'} - WINDOW_SIDECARS) + # Optional subset planning metadata describes a sampling policy, not a + # new immutable spike-coordinate identity. + acquisition = { + key: self.params[key] + for key in ( + 'acquisition_id', + 'n_channels_dat', + 'dtype', + 'offset', + 'depth_axis', + 'amplitude_units', + 'depth_units', + 'position_units', + 'waveform_units', + ) + if key in self.params + } + self.identity = { + 'format_version': 1, + 'acquisition': acquisition, + 'time_base': { + 'sample_rate': self.sample_rate, + 'num_samples': self.num_samples, + 'valid_intervals': self.valid_intervals.tolist(), + }, + 'arrays': {name: hash_file(self.path / name) for name in immutable}, + } + self.identity['sha256'] = hash_json(self.identity) + declared_source = info.get('source_sorting', {}) + if not isinstance(declared_source, dict): + raise ValueError('Waveform manifest source_sorting must be a mapping.') + if declared_source.get('num_spikes', self.n_spikes) != self.n_spikes: + raise ValueError('Waveform reservoir source spike count does not match this sorting.') + declared_times = declared_source.get('spike_times_sha256') + if ( + declared_times is not None + and declared_times != self.identity['arrays']['spike_times.npy'] + ): + raise ValueError( + 'Waveform reservoir source coordinate fingerprint does not match this sorting.' + ) + self.revision = assignment_revision(self.path, self.metadata) + if lineage: + if ( + lineage.get('coordinate_identity') != self.identity + or lineage.get('revision') != self.revision + ): + raise ValueError( + 'Curation lineage does not match immutable identity or assignments.' + ) + + def close(self): + """Release canonical array mappings after a batch command.""" + for array in self.arrays.values(): + close_array(array) + self.arrays.clear() + + def __enter__(self): + return self + + def __exit__(self, *args): + self.close() + + def iter_batches(self): + """Yield bounded row slices; vector fields use one-dimensional views.""" + for start in range(0, self.n_spikes, self.batch_size): + stop = min(start + self.batch_size, self.n_spikes) + batch = {'start': start, 'stop': stop} + for name in self.arrays.keys() & PER_SPIKE: + array = self.arrays[name] + if name in ('spike_times.npy', 'spike_clusters.npy', 'spike_templates.npy'): + array = _vector(array, name) + batch[name] = array[start:stop] + yield batch + + def validate(self): + """Check source row, clock, sparse channel, and feature alignment.""" + for name in self.arrays.keys() & PER_SPIKE: + if self.arrays[name].shape[0] != self.n_spikes: + raise ValueError(f'{name} is not aligned with source spikes.') + for name in ('spike_times.npy', 'spike_clusters.npy', 'spike_templates.npy'): + array = _vector(self.arrays[name], name) + if not np.issubdtype(array.dtype, np.integer): + raise ValueError(f'{name} must contain integer coordinates.') + previous = -1 + for batch in self.iter_batches(): + times = batch['spike_times.npy'] + if ( + np.any(times < previous) + or np.any(times[1:] < times[:-1]) + or np.any(times < 0) + or np.any(times >= self.num_samples) + ): + raise ValueError( + 'Spike sample times must be monotonic and within acquisition bounds.' + ) + if len(times): + previous = int(times[-1]) + interval_ids = np.searchsorted(self.valid_intervals[:, 0], times, side='right') - 1 + if np.any(interval_ids < 0) or np.any(times >= self.valid_intervals[interval_ids, 1]): + raise ValueError('Spike times fall outside valid acquisition intervals.') + for name in ('spike_clusters.npy', 'spike_templates.npy'): + if np.any(batch[name] < 0): + raise ValueError(f'{name} must contain nonnegative IDs.') + templates = self.arrays.get('templates.npy') + if templates is not None and np.any(batch['spike_templates.npy'] >= len(templates)): + raise ValueError('spike_templates.npy references an unavailable template.') + waveforms = self.arrays.keys() & WAVEFORMS + if waveforms and waveforms != WAVEFORMS: + raise ValueError('Stored waveforms require spikes, channels, and waveforms together.') + if waveforms: + ids = self.arrays['_phy_spikes_subset.spikes.npy'] + if ids.ndim != 1 or ids.dtype != np.dtype('int64'): + raise ValueError('Stored waveform row IDs must be a one-dimensional int64 array.') + channels = self.arrays['_phy_spikes_subset.channels.npy'] + values = self.arrays['_phy_spikes_subset.waveforms.npy'] + n_channels = len(self.arrays.get('channel_map.npy', [])) + if ( + values.ndim != 3 + or channels.shape != (len(ids), values.shape[2]) + or len(values) != len(ids) + or not np.issubdtype(channels.dtype, np.integer) + ): + raise ValueError('Stored waveform tensors and channels are not aligned.') + last = -1 + # Sparse-padding checks create temporary boolean arrays and gathers; + # cap their tensor bytes independently of the source row batch size. + row_bytes = max(1, values.dtype.itemsize * int(np.prod(values.shape[1:]))) + waveform_batch = min(self.batch_size, max(1, 8 * 1024 * 1024 // row_bytes)) + for start in range(0, len(ids), waveform_batch): + rows = ids[start : start + waveform_batch] + slots = channels[start : start + waveform_batch] + if ( + np.any(rows <= last) + or np.any(rows[1:] <= rows[:-1]) + or np.any(rows < 0) + or np.any(rows >= self.n_spikes) + ): + raise ValueError('Stored waveform row IDs must be sorted, unique source rows.') + if len(rows): + last = int(rows[-1]) + if np.any((slots < -1) | (slots >= n_channels)): + raise ValueError( + 'Sparse waveform channels must index model channels or use -1.' + ) + block = values[start : start + waveform_batch] + if np.any(np.moveaxis(block, 1, 2)[slots == -1] != 0): + raise ValueError('Sparse waveform padding channels must contain zeros.') + channel_map = self.arrays.get('channel_map.npy') + if channel_map is not None: + channels = _vector(channel_map, 'channel_map.npy') + if ( + not np.issubdtype(channels.dtype, np.integer) + or np.any(channels < 0) + or len(np.unique(channels)) != len(channels) + ): + raise ValueError('channel_map.npy must contain unique nonnegative raw column IDs.') + raw_channels = self.params.get('n_channels_dat') + if raw_channels is not None and ( + int(raw_channels) != raw_channels + or raw_channels <= 0 + or np.any(channels >= raw_channels) + ): + raise ValueError('channel_map.npy contains columns outside n_channels_dat.') + for name in ('channel_groups.npy', 'channel_shanks.npy', 'channel_map_si.npy'): + if name in self.arrays and len(_vector(self.arrays[name], name)) != len(channels): + raise ValueError(f'{name} must align with model channels.') + for name in ('whitening_mat.npy', 'whitening_mat_inv.npy'): + if name in self.arrays and self.arrays[name].shape != ( + len(channels), + len(channels), + ): + raise ValueError(f'{name} must contain a model-channel square matrix.') + templates = self.arrays.get('templates.npy') + if templates is not None: + # phylib uses singular template_ind; KS2's plural sidecar does + # not turn a dense tensor into a sparse tensor in the loader. + sparse = self.arrays.get('template_ind.npy') + if templates.ndim != 3: + raise ValueError('templates.npy must be a template/sample/channel tensor.') + if sparse is None and templates.shape[2] != len(channels): + raise ValueError('Dense templates must align with model channels.') + for index_name in ('template_ind.npy', 'templates_ind.npy'): + sparse = self.arrays.get(index_name) + if sparse is None: + continue + if sparse.shape != (len(templates), templates.shape[2]): + raise ValueError(f'{index_name} must align with sparse templates.') + if not np.issubdtype(sparse.dtype, np.integer): + raise ValueError(f'{index_name} must contain integer channel indices.') + if np.any((sparse < -1) | (sparse >= len(channels))): + raise ValueError(f'{index_name} contains invalid model channel IDs.') + positions = self.arrays.get('channel_positions.npy') + if positions is not None and ( + positions.shape != (len(channels), 2) or not np.all(np.isfinite(positions)) + ): + raise ValueError( + 'channel_positions.npy must contain finite x/y model channel coordinates.' + ) + for feature, indices in ( + ('pc_features.npy', 'pc_feature_ind.npy'), + ('template_features.npy', 'template_feature_ind.npy'), + ): + if feature in self.arrays and indices not in self.arrays: + raise ValueError(f'{feature} requires {indices}.') + if feature in self.arrays: + values = self.arrays[feature] + lookup = self.arrays[indices] + if lookup.ndim != 2 or values.ndim < 2 or values.shape[-1] != lookup.shape[1]: + raise ValueError(f'{feature} and {indices} sparse axes do not agree.') + templates = self.arrays.get('templates.npy') + if templates is not None and len(lookup) != len(templates): + raise ValueError(f'{indices} must align with templates.') + if not np.issubdtype(lookup.dtype, np.integer): + raise ValueError(f'{indices} must contain integer indices.') + limit = ( + len(self.arrays.get('channel_map.npy', [])) + if feature == 'pc_features.npy' + else len(lookup) + ) + for start in range(0, len(lookup), self.batch_size): + rows = lookup[start : start + self.batch_size] + if np.any((rows < -1) | (rows >= limit)): + raise ValueError(f'{indices} contains invalid sparse indices.') + return self diff --git a/phy/apps/large_recording/stability.py b/phy/apps/large_recording/stability.py new file mode 100644 index 000000000..131268fd1 --- /dev/null +++ b/phy/apps/large_recording/stability.py @@ -0,0 +1,176 @@ +"""Selected-unit full-duration evidence with explicit acquisition exposure.""" + +import numpy as np +from phylib.utils import Bunch, connect, unconnect + +from phy.cluster.views import ManualClusteringView +from phy.plot.plot import PlotCanvasMpl + + +class UnitStabilityView(ManualClusteringView): + """Full counts/rates, available source trends, gaps and retained evidence windows.""" + + plot_canvas_class = PlotCanvasMpl + max_n_clusters = 1 + _default_position = 'right' + ex_status = 'Full-duration evidence; sampling cannot establish chronic unit identity.' + + def __init__(self, controller): + super().__init__() + self.controller = controller + self.canvas.subplots(3, 1, sharex=True) + self._summary_arrays = {} + + def _array(self, name): + if name not in self._summary_arrays: + self._summary_arrays[name] = np.load(self.controller.dir_path / name, mmap_mode='r') + return self._summary_arrays[name] + + def get_data(self, cluster_id): + controller = self.controller + session = getattr(controller, 'bounded_session', None) + if session is None: + # Optional ordinary-data view: only the selected cluster is evaluated. + times = controller.get_spike_times(cluster_id) + edges = np.linspace(0, controller.model.duration, 101) + counts, _ = np.histogram(times, edges) + exposure = np.diff(edges) + return Bunch( + edges=edges, + counts=counts, + exposure=exposure, + windows=[], + intervals=[[0, controller.model.duration]], + trends={}, + provenance='Loaded recording; continuous acquisition assumed.', + ) + manifest = session.manifest + summary = manifest.get('summaries', manifest) + edges = self._array(summary['edges_file']) / controller.model.sample_rate + exposure = self._array(summary['exposure_file']) + units = controller._source_units(cluster_id) + if units is None: + return Bunch( + edges=edges, + counts=None, + exposure=exposure, + windows=[], + intervals=[], + trends={}, + provenance='Split evidence pending full-data replay.', + ) + lookup = {int(c): i for i, c in enumerate(session.unit_ids)} + rows = [lookup[c] for c in units] + # Bound temporaries by one unit's summary row, including merged units. + count_grid = self._array(summary['counts_file']) + counts = np.zeros(len(exposure), dtype=np.int64) + for row in rows: + counts += count_grid[row] + trends = {} + for attribute in ('amplitude', 'depth'): + sum_file, count_file = f'{attribute}_sum_file', f'{attribute}_count_file' + if sum_file not in summary: + continue + total, valid = np.zeros(len(exposure)), np.zeros(len(exposure), dtype=np.int64) + for row in rows: + total += self._array(summary[sum_file])[row] + valid += self._array(summary[count_file])[row] + mean = np.full(len(exposure), np.nan) + np.divide(total, valid, out=mean, where=valid > 0) + trends[attribute] = mean + return Bunch( + edges=edges, + counts=counts, + exposure=exposure, + trends=trends, + windows=controller.evidence_windows / controller.model.sample_rate, + intervals=np.asarray(manifest['valid_intervals']) / controller.model.sample_rate, + provenance='Complete source sorting; exact counts and finite-value means.', + ) + + def plot(self, **kwargs): + axes = list(self.canvas.iter_ax()) + for ax in axes: + ax.clear() + session = getattr(self.controller, 'bounded_session', None) + if not self.cluster_ids: + axes[0].text( + 0.02, + 0.6, + 'Select a unit to inspect full-duration evidence.', + transform=axes[0].transAxes, + ) + else: + cluster_id = self.cluster_ids[0] + data = self.get_data(cluster_id) + if data.counts is None: + axes[0].text(0.02, 0.6, data.provenance, transform=axes[0].transAxes) + else: + rates = np.full(len(data.exposure), np.nan) + np.divide(data.counts, data.exposure, out=rates, where=data.exposure > 0) + x = (data.edges[:-1] + data.edges[1:]) / 7200 + axes[0].plot(x, rates) + for ax, attribute in zip(axes[1:], ('amplitude', 'depth')): + if attribute in data.trends: + ax.plot(x, data.trends[attribute]) + else: + ax.text( + 0.02, + 0.6, + f'{attribute.title()} measurement unavailable.', + transform=ax.transAxes, + ) + for ax in axes: + for start, end in data.windows: + ax.axvspan(start / 3600, end / 3600, color='cyan', alpha=0.15) + previous = 0 + for start, end in data.intervals: + if start > previous: + ax.axvspan(previous / 3600, start / 3600, color='gray', alpha=0.4) + previous = end + if previous < data.edges[-1]: + ax.axvspan(previous / 3600, data.edges[-1] / 3600, color='gray', alpha=0.4) + axes[0].set_title(f'Unit {cluster_id} — {data.provenance}', fontsize=9) + axes[0].set_ylabel('Rate (Hz)') + axes[1].set_ylabel('Mean amplitude\n(source units)') + axes[2].set_ylabel('Mean depth\n(source units)') + axes[2].set_xlabel('Elapsed acquisition clock (hours)') + if session is not None: + missing = [ + c for c, n in zip(session.unit_ids, session.manifest['bounded_counts']) if not n + ] + if missing: + axes[2].text( + 0.02, + -0.42, + f'{len(missing)} source units without evidence: ' + f'{", ".join(map(str, missing[:12]))}' + (' …' if len(missing) > 12 else ''), + transform=axes[2].transAxes, + fontsize=8, + ) + self.canvas.figure.tight_layout() + self.canvas.update() + + def on_select(self, cluster_ids=None, **kwargs): + self.cluster_ids = tuple(cluster_ids or ())[:1] + self.plot(**kwargs) + + def attach(self, gui): + super().attach(gui) + + @self.actions.add(prompt=True) + def show_source_unit(cluster_id): + """Inspect a full source unit, including one without sampled evidence.""" + session = getattr(self.controller, 'bounded_session', None) + if session is not None and int(cluster_id) not in session.unit_ids: + raise ValueError('Unknown source unit ID.') + self.on_select([int(cluster_id)]) + + @connect(event='bounded_state_changed', sender=self.controller.supervisor) + def on_state(sender): + self.plot() + + @connect(event='close_view', sender=self) + def on_close(sender, *args): + unconnect(on_state, on_close) + self._summary_arrays.clear() diff --git a/phy/apps/large_recording/statistics.py b/phy/apps/large_recording/statistics.py new file mode 100644 index 000000000..04742de84 --- /dev/null +++ b/phy/apps/large_recording/statistics.py @@ -0,0 +1,67 @@ +"""Descriptive statistics restricted to retained contiguous observation windows.""" + +import numpy as np +from phylib.stats import correlograms + + +def window_isi(samples, windows, sample_rate): + intervals = [] + for start, end in windows: + first, last = np.searchsorted(samples, [start, end]) + intervals.append(np.diff(samples[first:last]) / sample_rate) + return np.concatenate(intervals) if intervals else np.array([], dtype=np.float64) + + +def window_correlograms( + samples, clusters, windows, cluster_ids, sample_rate, bin_size, window_size +): + """Exact window-local pair counts, without a full-duration rate baseline. + + Positive/negative lag bins follow phylib's floor-to-bin convention. The + central bin counts unordered within-unit pairs and cross-unit pairs with + absolute delays below one bin. Cross-unit central counts sum both temporal + orientations, making them independent of equal-time source row order. + """ + selected = np.isin(clusters, cluster_ids) + order = np.argsort(cluster_ids) + sorted_ids = np.asarray(cluster_ids)[order] + result = None + for start, end in windows: + first, last = np.searchsorted(samples, [start, end]) + mask = selected[first:last] + ticks = samples[first:last][mask] - start + if len(ticks) and ticks.max() >= 2**50: + raise ValueError( + 'Observation window is too long for exact correlogram sample conversion.' + ) + # phylib truncates float times * rate to integer samples. A common + # quarter-sample offset preserves every integer tick after conversion + # without changing any pair delay, including late-clock recordings. + times = (ticks.astype(np.float64) + 0.25) / sample_rate + # Avoid phylib's lookup allocation indexed by absolute cluster ID. + labels = order[np.searchsorted(sorted_ids, clusters[first:last][mask])] + if np.any(np.bincount(labels, minlength=len(cluster_ids)) > 46340): + raise ValueError( + 'Evidence window has too many spikes for exact correlogram counts; prepare shorter windows.' + ) + positive = correlograms( + times, + labels, + sample_rate=sample_rate, + cluster_ids=np.arange(len(cluster_ids)), + bin_size=bin_size, + window_size=window_size, + symmetrize=False, + ).astype(np.int64) + # Accumulate before symmetrization to preserve exact central counts and + # int64 totals across many windows. Self-spike pairs are never counted. + if result is None: + result = positive + else: + result += positive + if result is None: + raise ValueError('Correlograms require at least one observation window.') + center = result[..., 0] + result[..., 0].T + np.fill_diagonal(center, np.diag(result[..., 0])) + result[..., 0] = center + return np.dstack((result[..., 1:][..., ::-1].transpose(1, 0, 2), result)) diff --git a/phy/apps/large_recording/summaries.py b/phy/apps/large_recording/summaries.py new file mode 100644 index 000000000..2f498fab5 --- /dev/null +++ b/phy/apps/large_recording/summaries.py @@ -0,0 +1,122 @@ +"""Exact file-backed counts and additive source-measurement summaries.""" + +from contextlib import ExitStack +from pathlib import Path + +import numpy as np + +from .source import close_array + + +def scalar_attributes(source): + """Return available scalar measurements without inventing physical units.""" + attributes = {} + amplitude = source.arrays.get('amplitudes.npy') + if amplitude is not None and (amplitude.ndim == 1 or amplitude.shape[1:] == (1,)): + attributes['amplitude'] = amplitude.reshape(-1) + depth = source.arrays.get('spike_depths.npy') + if depth is not None and (depth.ndim == 1 or depth.shape[1:] == (1,)): + attributes['depth'] = depth.reshape(-1) + elif (positions := source.arrays.get('spike_positions.npy')) is not None: + axis = source.params.get('depth_axis') + if axis is not None: + if ( + not isinstance(axis, int) + or positions.ndim != 2 + or not 0 <= axis < positions.shape[1] + ): + raise ValueError('depth_axis must declare a valid spike_positions.npy coordinate.') + attributes['depth'] = positions[:, axis] + return attributes + + +def summary_size(source, bin_seconds=3600): + """Estimate exact array bytes before allocating the dense output grid.""" + step = float(bin_seconds) * source.sample_rate + if not np.isfinite(step) or step < 1 or not np.isclose(step, round(step), rtol=0, atol=1e-6): + raise ValueError('bin_seconds must specify a positive integer sample width.') + step = int(round(step)) + bins = (source.num_samples + step - 1) // step + size = len(source.unit_ids) * bins * 8 * (1 + 2 * len(scalar_attributes(source))) + size += (bins * 2 + 1) * 8 + return size, step, bins + + +def build_summaries(source, output_dir, bin_seconds=3600, max_bytes=536870912): + """Stream exact full counts, exposure, and finite-value sums/counts. + + Exposure is one shared vector. Summing parent-unit counts on merge must + divide by this vector once. Scalar means are sums/counts, never medians. + """ + size, step, bins = summary_size(source, bin_seconds) + if size + 1024 * (1 + 2 * len(scalar_attributes(source))) > max_bytes: + raise ValueError( + f'Summary budget exceeded: need at least {size} array bytes; ' + 'increase bin_seconds or max_bytes.' + ) + with ExitStack() as stack: + output_dir = Path(output_dir) + output_dir.mkdir(parents=True, exist_ok=True) + result = { + 'counts_file': '_phy_summary_counts.npy', + 'edges_file': '_phy_summary_edges.npy', + 'exposure_file': '_phy_summary_exposure.npy', + } + shape = (len(source.unit_ids), bins) + counts = np.lib.format.open_memmap( + output_dir / result['counts_file'], mode='w+', dtype=np.int64, shape=shape + ) + stack.callback(close_array, counts) + counts[:] = 0 + attributes = scalar_attributes(source) + grids = {} + for field in attributes: + for kind, dtype in (('sum', np.float64), ('count', np.int64)): + key = f'{field}_{kind}_file' + result[key] = f'_phy_summary_{field}_{kind}.npy' + array = np.lib.format.open_memmap( + output_dir / result[key], mode='w+', dtype=dtype, shape=shape + ) + stack.callback(close_array, array) + array[:] = 0 + grids[field, kind] = array + edges = np.lib.format.open_memmap( + output_dir / result['edges_file'], mode='w+', dtype=np.int64, shape=(bins + 1,) + ) + stack.callback(close_array, edges) + # Bounded chunks also avoid a recording-duration-sized temporary arange. + for start in range(0, bins + 1, source.batch_size): + stop = min(start + source.batch_size, bins + 1) + edges[start:stop] = np.arange(start, stop, dtype=np.int64) * step + edges[-1] = source.num_samples + exposure = np.lib.format.open_memmap( + output_dir / result['exposure_file'], mode='w+', dtype=np.float64, shape=(bins,) + ) + stack.callback(close_array, exposure) + exposure[:] = 0 + for start, end in source.valid_intervals: + first, last = int(start // step), int((end - 1) // step) + for offset in range(first, last + 1, source.batch_size): + stop = min(offset + source.batch_size, last + 1) + exposure[offset:stop] += ( + np.minimum(edges[offset + 1 : stop + 1], end) + - np.maximum(edges[offset:stop], start) + ) / source.sample_rate + for batch in source.iter_batches(): + rows = np.searchsorted(source.unit_ids, batch['spike_clusters.npy']) + columns = batch['spike_times.npy'] // step + np.add.at(counts, (rows, columns), 1) + for field, values in attributes.items(): + values = values[batch['start'] : batch['stop']] + finite = np.isfinite(values) + np.add.at(grids[field, 'sum'], (rows[finite], columns[finite]), values[finite]) + np.add.at(grids[field, 'count'], (rows[finite], columns[finite]), 1) + for array in (counts, edges, exposure, *grids.values()): + array.flush() + result['estimators'] = { + 'counts': 'exact', + 'exposure': 'acquisition_seconds', + **dict.fromkeys(attributes, 'finite_sum_and_count'), + } + result['measurement_units'] = {field: f'source {field} units' for field in attributes} + return result diff --git a/phy/apps/large_recording/supervisor.py b/phy/apps/large_recording/supervisor.py new file mode 100644 index 000000000..a5765232f --- /dev/null +++ b/phy/apps/large_recording/supervisor.py @@ -0,0 +1,459 @@ +"""Guarded native commands paired with semantic state, history and checkpoints.""" + +from contextlib import contextmanager +from copy import deepcopy +from functools import wraps +from types import MappingProxyType + +import numpy as np +from phylib.utils import emit + +from phy.cluster._history import History +from phy.cluster._utils import ClusterMeta +from phy.cluster.clustering import Clustering +from phy.cluster.supervisor import Supervisor + +from .curation import evaluate_predicate, interpret_events + + +class MutationPolicy: + """Only semantic commands may mutate native bounded state.""" + + def __init__(self): + self.depth = 0 + + def require(self): + if not self.depth: + raise ValueError( + 'Bounded curation requires a semantic command. Use supervisor.merge(), ' + 'supervisor.label(), or split_source_amplitude(); arbitrary splits are unavailable.' + ) + + @contextmanager + def command(self): + self.depth += 1 + try: + yield + finally: + self.depth -= 1 + + +class GuardedClustering(Clustering): + def __init__(self, *args, policy, **kwargs): + self.policy = policy + super().__init__(*args, **kwargs) + + @property + def spike_clusters(self): + result = self._spike_clusters.view() + result.flags.writeable = bool(self.policy.depth) + return result + + @property + def cluster_ids(self): + result = self._cluster_ids.view() + result.flags.writeable = False + return result + + @property + def spikes_per_cluster(self): + return MappingProxyType(self._spikes_per_cluster) + + def _update_cluster_ids(self, *args, **kwargs): + super()._update_cluster_ids(*args, **kwargs) + for ids in self._spikes_per_cluster.values(): + ids.flags.writeable = False + + def merge(self, *args, **kwargs): + self.policy.require() + return super().merge(*args, **kwargs) + + def assign(self, *args, **kwargs): + self.policy.require() + return super().assign(*args, **kwargs) + + def reset(self): + self.policy.require() + return super().reset() + + def undo(self): + self.policy.require() + return super().undo() + + def redo(self): + self.policy.require() + return super().redo() + + +class GuardedMeta(ClusterMeta): + def __init__(self, policy): + self.policy = policy + super().__init__() + + def set(self, *args, **kwargs): + self.policy.require() + return super().set(*args, **kwargs) + + def add_field(self, *args, **kwargs): + self.policy.require() + return super().add_field(*args, **kwargs) + + def from_dict(self, *args, **kwargs): + self.policy.require() + return super().from_dict(*args, **kwargs) + + def undo(self): + self.policy.require() + return super().undo() + + def redo(self): + self.policy.require() + return super().redo() + + +class SessionHistory(History): + def __init__(self, session): + self.session = session + super().__init__(session.journal_state()) + + def capture(self): + self.add(self.session.journal_state()) + + def undo(self): + self.back() + self.session.restore_journal(self.current_item) + + def redo(self): + self.forward() + self.session.restore_journal(self.current_item) + + +def semantic_command(method): + @wraps(method) + def execute(self, *args, **kwargs): + with self._transaction(): + return method(self, *args, **kwargs) + + return execute + + +class BoundedSupervisor(Supervisor): + """Native UI/history operations paired with a validated semantic journal.""" + + def __init__(self, *args, session, **kwargs): + self.session = session + self.policy = MutationPolicy() + self.select_actions = None + self.cluster_view = None + self.similarity_view = None + + def clustering_factory(*args, **kwargs): + return GuardedClustering(*args, policy=self.policy, **kwargs) + + def meta_factory(groups, labels): + meta = GuardedMeta(self.policy) + data = {} + for field, values in labels.items(): + for cluster, value in values.items(): + data.setdefault(cluster, {})[field] = value + with self.policy.command(): + meta.add_field('group') + meta.from_dict(data) + return meta + + # Generation state supersedes disposable clustering caches. + kwargs['context'] = None + kwargs['new_cluster_id'] = session.next_cluster_id + kwargs['clustering_factory'] = clustering_factory + kwargs['cluster_meta_factory'] = meta_factory + super().__init__(*args, **kwargs) + self.session_history = SessionHistory(session) + + @contextmanager + def _transaction(self): + """Roll native and semantic state back if an observer or command fails.""" + clustering = self.clustering + arrays = ( + clustering._spike_clusters.copy(), + dict(clustering._spikes_per_cluster), + clustering._cluster_ids.copy(), + ) + data, fields = deepcopy(self.cluster_meta._data), dict(self.cluster_meta._fields) + histories = [ + clustering._undo_stack, + self.cluster_meta._undo_stack, + self._global_history, + self.session_history, + ] + if self.merge_propositions is not None: + histories.append(self.merge_propositions._history) + history_states = [(list(h._history), h._index) for h in histories] + journal, selection, dirty = ( + self.session.journal_state(), + self.selection.snapshot(), + self._is_dirty, + ) + propositions = ( + self.merge_propositions.snapshot() if self.merge_propositions is not None else None + ) + workflow, columns = self._workflow_context(), list(self.columns) + try: + yield + except Exception: + clustering._spike_clusters, clustering._spikes_per_cluster, clustering._cluster_ids = ( + arrays + ) + self.cluster_meta._data, self.cluster_meta._fields = data, fields + for history, (items, index) in zip(histories, history_states): + history._history, history._index = items, index + self.session.restore_journal(journal) + self.session.next_cluster_id = max( + self.session.next_cluster_id, clustering.new_cluster_id() + ) + self.selection.restore(selection) + if propositions is not None: + self.merge_propositions.restore(propositions) + self._is_dirty = dirty + self.columns = columns + self._lineage_cache = None + if self.gui is not None: + self._reset_cluster_view() + self._restore_history_context( + selection, + { + 'mode': 'merge' if selection.is_merge_mode else 'post_merge', + 'tables': workflow, + }, + 'undo', + ) + self._refresh_bounded() + self._update_save_feedback() + raise + + def _require_capability(self, name): + if not self.session.manifest.get('capabilities', {}).get(name, False): + raise ValueError(f'This bounded package does not allow {name} operations.') + + def get_cluster_info(self, cluster_id, exclude=()): + info = super().get_cluster_info(cluster_id, exclude=exclude) + # Producer TSV labels may use metric names. They must not overwrite + # authoritative complete-sorting counts/rates with inherited values. + for name in ( + 'n_spikes', + 'fr', + 'full_count', + 'evidence_count', + 'waveform_count', + 'evidence', + ): + if name not in exclude and name in self.cluster_metrics: + info[name] = self.cluster_metrics[name](cluster_id) + return info + + def _select_after_merge(self, *args, **kwargs): + if self.gui is not None: + return super()._select_after_merge(*args, **kwargs) + + def _select_after_split(self, *args, **kwargs): + if self.gui is not None: + return super()._select_after_split(*args, **kwargs) + + def _restore_history_context(self, *args, **kwargs): + if self.gui is not None: + return super()._restore_history_context(*args, **kwargs) + + def _reset_cluster_view(self, *args, **kwargs): + if self.gui is not None: + return super()._reset_cluster_view(*args, **kwargs) + + def _workflow_context(self): + return super()._workflow_context() if self.gui is not None else {} + + def _validate_native(self, metadata=True): + if not np.array_equal(self.clustering.spike_clusters, self.session.assignments): + raise ValueError('Native assignments disagree with the bounded semantic command.') + if metadata: + active, expected, _ = interpret_events( + self.session.unit_ids, self.session.base_metadata, self.session.events + ) + for field in set(self.cluster_meta.fields) | set(expected): + for cluster in active: + actual = ( + self.cluster_meta.get(field, cluster) + if field in self.cluster_meta.fields + else None + ) + if actual != expected.get(field, {}).get(cluster): + raise ValueError( + 'Native metadata disagrees with the bounded semantic command.' + ) + + def _synchronize_metadata(self): + # Native ClusterMeta.undo does not retain nonstacked inherited merge + # labels. The semantic projection is authoritative in bounded mode. + data = {} + for field, values in self.session.metadata.items(): + if field not in self.cluster_meta.fields: + with self.policy.command(): + self.cluster_meta.add_field(field) + for cluster, value in values.items(): + data.setdefault(cluster, {})[field] = value + self.cluster_meta._data = data + + def _refresh_bounded(self): + if self.gui is not None: + data = self.cluster_info + self.cluster_view.change(data) + self.similarity_view.change(data) + if self.merge_view is not None: + self.merge_view.change(data) + emit('bounded_state_changed', self) + + def _capture(self): + self._validate_native() + self.session_history.capture() + self._global_history.add_to_current_action(self.session_history) + self._is_dirty = True + self._update_save_feedback() + self._synchronize_metadata() + self._refresh_bounded() + + @semantic_command + def merge(self, cluster_ids=None, to=None): + self._require_capability('merge') + parents = list(self.selected if cluster_ids is None else cluster_ids) + if len(parents) < 2: + return + target = max(self.clustering.new_cluster_id(), self.session.next_cluster_id) + target = target if to is None else int(to) + self.session.preflight_merge(parents, target) + with self.policy.command(): + out = super().merge(cluster_ids=cluster_ids, to=target) + if out is None: + return + resolved = { + field: self.cluster_meta.get(field, target) for field in self.cluster_meta.fields + } + self.session.record_merge(parents, target, metadata=resolved) + self._capture() + return out + + @semantic_command + def label(self, name, value, cluster_ids=None): + self._require_capability('label') + if name in ( + 'n_spikes', + 'fr', + 'full_count', + 'evidence_count', + 'waveform_count', + 'evidence', + ): + raise ValueError('Bounded count/rate/evidence columns are read-only metrics.') + ids = list(self.selected if cluster_ids is None else np.atleast_1d(cluster_ids)) + if not ids or self.selection.state.is_merge_mode: + return + # Preflight against a copy before touching native metadata. + before = self.session.snapshot() + self.session.record_label(name, ids, value) + self.session.restore(before) + with self.policy.command(): + super().label(name, value, cluster_ids=ids) + self.session.record_label(name, ids, value) + self._capture() + + def split(self, *args, **kwargs): + raise ValueError( + 'Lasso, feature and waveform splits cannot be replayed on unseen spikes. ' + 'Use Split by source amplitude for a full-source scalar predicate.' + ) + + @semantic_command + def split_source_amplitude(self, cluster_id, threshold): + if 'amplitudes.npy' not in self.session.manifest.get('capabilities', {}).get( + 'predicate_split_attributes', [] + ): + raise ValueError('This package does not support source-amplitude predicate splits.') + if self.selection.state.is_merge_mode: + raise ValueError('Leave Merge mode before splitting by source amplitude.') + if 'amplitudes.npy' not in self.session.attributes: + raise ValueError('The source has no scalar per-spike amplitudes.') + predicate = { + 'attribute': 'amplitudes.npy', + 'operator': '<', + 'threshold': float(threshold), + 'nonfinite': 'keep', + } + ids = self.clustering.spikes_per_cluster.get(int(cluster_id), np.array([], dtype=np.int64)) + selected = ids[ + evaluate_predicate(self.session.attributes['amplitudes.npy'][ids], predicate) + ] + if not len(selected) or len(selected) == len(ids): + raise ValueError('Both sides of the threshold must have evidence spikes.') + selected_id = max(self.clustering.new_cluster_id(), self.session.next_cluster_id) + remainder_id = selected_id + 1 + self.clustering._new_cluster_id = selected_id + metadata = { + field: { + child: self.cluster_meta.get(field, cluster_id) + for child in (selected_id, remainder_id) + } + for field in self.cluster_meta.fields + } + before = self.session.snapshot() + self.session.record_threshold_split( + int(cluster_id), [remainder_id, selected_id], predicate, metadata=metadata + ) + self.session.restore(before) + with self.policy.command(): + out = Supervisor.split(self, selected) + self.session.record_threshold_split( + int(cluster_id), [remainder_id, selected_id], predicate, metadata=metadata + ) + self._capture() + return out + + @semantic_command + def undo(self): + with self.policy.command(): + super().undo() + self._after_history() + + @semantic_command + def redo(self): + with self.policy.command(): + super().redo() + self._after_history() + + def _after_history(self): + self._validate_native(metadata=False) + # Never reuse allocated IDs after branching or restoring a snapshot. + self.session.next_cluster_id = max( + self.session.next_cluster_id, self.clustering.new_cluster_id() + ) + self._is_dirty = True + self._update_save_feedback() + self._synchronize_metadata() + self._refresh_bounded() + + def save(self): + self._validate_native() + documents = None + if self.merge_propositions is not None: + from .source import hash_file + + documents = { + 'curation_review.json': self.merge_propositions.catalog.review_mapping( + source_sha256=hash_file(self.session.path / 'curation.json') + ) + } + self.session.save( + assignments=self.clustering.spike_clusters, + metadata=self.session.metadata, + next_cluster_id=max(self.session.next_cluster_id, self.clustering.new_cluster_id()), + extra_documents=documents, + ) + if self.merge_propositions is not None: + self.merge_propositions.mark_saved() + self._is_dirty = False + self._update_save_feedback(saved=True) diff --git a/phy/apps/large_recording/tests/__init__.py b/phy/apps/large_recording/tests/__init__.py new file mode 100644 index 000000000..ba239ecd8 --- /dev/null +++ b/phy/apps/large_recording/tests/__init__.py @@ -0,0 +1 @@ +"""Long-recording workflow regression tests.""" diff --git a/phy/apps/large_recording/tests/test_adapter.py b/phy/apps/large_recording/tests/test_adapter.py new file mode 100644 index 000000000..76d4dda47 --- /dev/null +++ b/phy/apps/large_recording/tests/test_adapter.py @@ -0,0 +1,449 @@ +"""Bounded GUI correctness and ordinary-controller isolation.""" + +import json + +import numpy as np +import pytest +from click.testing import CliRunner +from phylib.utils import connect, unconnect + +from phy.apps import phycli +from phy.apps.large_recording.adapter import BoundedTemplateController +from phy.apps.large_recording.curation import BoundedSession +from phy.apps.large_recording.prepare import prepare_large +from phy.apps.large_recording.replay import apply_curation +from phy.apps.large_recording.tests.test_source import make_sorting +from phy.apps.template.gui import TemplateController + + +def make_controller(tmp_path, sparse=False): + params = make_gui_sorting(tmp_path / 'source') + output = prepare_large( + params, + tmp_path / 'bounded', + window_seconds=3 if sparse else 20, + stratum_seconds=20, + max_spikes=5 if sparse else 100, + bin_seconds=5, + batch_size=2, + ) + controller = TemplateController( + dir_path=output, + config_dir=tmp_path / 'config', + plugins=[], + enable_threading=False, + clear_cache=True, + ) + return params, controller + + +def make_gui_sorting(path): + params = make_sorting(path) + times = np.load(path / 'spike_times.npy') + np.save(path / 'spike_times.npy', times * 30000) + params.write_text( + params.read_text() + .replace('sample_rate = 1.0', 'sample_rate = 30000.0') + .replace('num_samples = 20', 'num_samples = 600000') + ) + return params + + +def test_manifest_factory_reserves_complete_ids_and_guards_mutations(tmp_path): + _, controller = make_controller(tmp_path, sparse=True) + try: + assert isinstance(controller, BoundedTemplateController) + with pytest.raises(ValueError): + controller.model.spike_clusters[0] = 123 + with pytest.raises(ValueError): + controller.model.spike_clusters = np.array([123]) + with pytest.raises(ValueError): + controller.model.sample_rate = 1.0 + with pytest.raises(TypeError): + controller.model.metadata['group'][0] = 'noise' + supervisor = controller.supervisor + assert 99 not in supervisor.clustering.cluster_ids + assert supervisor.clustering.new_cluster_id() == 100 + assert controller.get_full_count(0) == 4 + assert controller.get_mean_firing_rate(0) == 0.2 + for operation in ( + lambda: supervisor.clustering.merge([0, 1]), + lambda: supervisor.clustering.assign([0]), + lambda: supervisor.clustering.split([0]), + supervisor.clustering.undo, + supervisor.clustering.redo, + supervisor.clustering.reset, + lambda: supervisor.cluster_meta.set('group', [0], 'noise', add_to_stack=False), + lambda: supervisor.cluster_meta.from_dict({0: {'group': 'noise'}}), + lambda: controller.model.save_spike_clusters(supervisor.clustering.spike_clusters), + lambda: controller.model.save_metadata('group', {0: 'noise'}), + lambda: supervisor.split([0]), + ): + with pytest.raises(ValueError): + operation() + with pytest.raises(ValueError): + supervisor.clustering.spike_clusters[0] = 8 + assert not controller.bounded_session.events + with pytest.raises(ValueError): + supervisor.cluster_meta.add_field('group', 'noise') + with pytest.raises(TypeError): + supervisor.clustering.spikes_per_cluster[0] = np.array([0]) + with pytest.raises(ValueError): + supervisor.clustering.spikes_per_cluster[0][0] = 1 + assert 'CorrelogramView' not in controller.view_creator + assert 'FiringRateView' not in controller.view_creator + finally: + controller.close() + + +def test_loader_and_native_history_do_not_allocate_by_sparse_unit_id(qtbot, tmp_path): + params = make_gui_sorting(tmp_path / 'source') + labels = np.load(params.parent / 'spike_clusters.npy').astype(np.int64) + labels[labels == 99] = 2_000_000_000 + np.save(params.parent / 'spike_clusters.npy', labels) + (params.parent / 'cluster_group.tsv').write_text( + 'cluster_id\tgroup\n0\tgood\n1\tmua\n2\tgood\n2000000000\tnoise\n' + ) + output = prepare_large(params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20) + controller = TemplateController(dir_path=output, config_dir=tmp_path / 'config', plugins=[]) + try: + assert controller.model.sparse_clusters is None + first = controller.supervisor.merge([0, 1]) + assert first.added == [2_000_000_001] + controller.supervisor.merge([2, 2_000_000_001]) + controller.supervisor.undo() + assert 2_000_000_001 in controller.supervisor.clustering.cluster_ids + values = controller._window_correlograms([2, 2_000_000_001], 0.001, 0.004) + assert values.shape == (2, 2, 5) + controller.supervisor.save() + finally: + controller.close() + reopened = TemplateController(dir_path=output, config_dir=tmp_path / 'config', plugins=[]) + try: + assert 2_000_000_001 in reopened.model.cluster_ids + assert reopened.model.sparse_clusters is None + finally: + reopened.close() + + +def test_merge_label_history_save_reopen_replay_regenerate(tmp_path): + params, controller = make_controller(tmp_path) + path = controller.dir_path + try: + supervisor = controller.supervisor + original_file = (path / 'spike_clusters.npy').read_bytes() + up = supervisor.merge([0, 1]) + assert up.added == [100] + assert controller.get_full_count(100) == 7 + supervisor.label('group', 'noise', [100]) + supervisor.undo() + assert controller.bounded_session.metadata['group'][100] == 'good' + assert supervisor.cluster_meta.get('group', 100) == 'good' + supervisor.redo() + assert controller.bounded_session.metadata['group'][100] == 'noise' + supervisor.undo() + supervisor.undo() + np.testing.assert_array_equal( + supervisor.clustering.spike_clusters, controller.bounded_session.base_assignments + ) + supervisor.merge([1, 2]) + assert supervisor.clustering.new_cluster_id() == 102 + supervisor.label('quality', 'reviewed', [101]) + supervisor.save() + assert not supervisor.is_dirty() + assert (path / 'spike_clusters.npy').read_bytes() == original_file + expected = supervisor.clustering.spike_clusters.copy() + finally: + controller.close() + reopened = TemplateController( + dir_path=path, config_dir=tmp_path / 'config', plugins=[], enable_threading=False + ) + try: + np.testing.assert_array_equal(reopened.supervisor.clustering.spike_clusters, expected) + assert reopened.supervisor.clustering.new_cluster_id() == 102 + assert reopened.supervisor.cluster_meta.get('quality', 101) == 'reviewed' + assert reopened.supervisor._global_history.current_position == 0 + finally: + reopened.close() + output = apply_curation(params, path, tmp_path / 'curated', batch_size=2) + source = np.load(params.parent / 'spike_clusters.npy') + oracle = source.copy() + oracle[np.isin(source, [1, 2])] = 101 + np.testing.assert_array_equal(np.load(output / 'spike_clusters.npy'), oracle) + regenerated = prepare_large( + output / 'params.py', + tmp_path / 'regenerated', + window_seconds=20, + stratum_seconds=20, + bin_seconds=5, + ) + assert BoundedSession(regenerated).unit_ids == [0, 99, 101] + + +def test_source_threshold_replays_every_spike_and_full_metrics_pending(tmp_path): + params, controller = make_controller(tmp_path) + try: + supervisor = controller.supervisor + up = supervisor.split_source_amplitude(0, 4) + assert set(up.added) == {100, 101} + assert controller.get_full_count(100) is None + assert controller.get_mean_firing_rate(101) is None + supervisor.undo() + assert controller.get_full_count(0) == 4 + supervisor.redo() + supervisor.save() + output = apply_curation(params, controller.dir_path, tmp_path / 'curated', batch_size=2) + source = np.load(params.parent / 'spike_clusters.npy') + amplitudes = np.load(params.parent / 'amplitudes.npy') + oracle = source.copy() + oracle[(source == 0) & (amplitudes < 4)] = 100 + oracle[(source == 0) & ~(amplitudes < 4)] = 101 + np.testing.assert_array_equal(np.load(output / 'spike_clusters.npy'), oracle) + finally: + controller.close() + + +def test_failed_save_leaves_dirty_previous_checkpoint(tmp_path, monkeypatch): + _, controller = make_controller(tmp_path) + try: + supervisor = controller.supervisor + supervisor.merge([0, 1]) + supervisor.save() + supervisor.label('group', 'noise', [100]) + from phy.apps.large_recording import curation + + original = curation._write_json + + def fail(path, value): + if path.name == 'journal.json': + raise OSError('interrupted') + return original(path, value) + + monkeypatch.setattr(curation, '_write_json', fail) + with pytest.raises(OSError, match='interrupted'): + supervisor.save() + assert supervisor.is_dirty() + assert BoundedSession(controller.dir_path).metadata['group'][100] == 'good' + finally: + controller.close() + + +def test_standard_controller_retains_existing_path_and_save_format(qtbot, tmp_path): + params = make_gui_sorting(tmp_path / 'ordinary') + controller = TemplateController( + dir_path=params.parent, + config_dir=tmp_path / 'config', + plugins=[], + enable_threading=False, + sample_rate=30000, + n_channels_dat=3, + dtype='int16', + ) + gui = controller.create_gui(default_views=[], do_prompt_save=False) + try: + assert type(controller) is TemplateController + assert controller.gui_name == 'TemplateGUI' + assert controller.default_views[0] == 'WaveformView' + assert 'FiringRateView' in controller.view_creator + up = controller.supervisor.merge([0, 1]) + controller.supervisor.save() + assert up.added == [100] + assert not (params.parent / '.phy-curation').exists() + assert not (params.parent / '_phy_curation_manifest.json').exists() + assert 100 in np.load(params.parent / 'spike_clusters.npy') + finally: + gui.close() + controller.close() + + +def test_malformed_manifest_fails_before_model_loading(tmp_path): + params = make_sorting(tmp_path / 'ordinary') + (params.parent / '_phy_curation_manifest.json').write_text(json.dumps({'format_version': 99})) + with pytest.raises(ValueError, match='Unsupported'): + TemplateController(dir_path=params.parent, plugins=[]) + + +def test_cli_preparation_and_replay(tmp_path): + params = make_sorting(tmp_path / 'source') + runner = CliRunner() + output = tmp_path / 'bounded' + result = runner.invoke( + phycli, + [ + 'prepare-large', + str(params), + str(output), + '--window-seconds', + '20', + '--stratum-seconds', + '20', + '--bin-seconds', + '5', + '--batch-size', + '2', + ], + ) + assert result.exit_code == 0, result.output + session = BoundedSession(output) + session.record_merge([0, 1], 100, {'group': 'good'}) + session.save() + result = runner.invoke( + phycli, + [ + 'apply-curation', + str(params), + str(output), + str(tmp_path / 'curated'), + '--batch-size', + '2', + ], + ) + assert result.exit_code == 0, result.output + + +def test_stability_gui_and_exact_trends(qtbot, tmp_path): + _, controller = make_controller(tmp_path, sparse=True) + gui = controller.create_gui(default_views=['UnitStabilityView'], do_prompt_save=False) + try: + assert 'Bounded evidence' in gui.windowTitle() + view = gui.list_views('UnitStabilityView')[0] + data = view.get_data(0) + assert data.counts.sum() == 4 + np.testing.assert_allclose(data.trends['amplitude'], [1.5, np.nan, 6, 9], equal_nan=True) + view.on_select([99]) # Unit absent from evidence still has full-duration counts. + assert view.get_data(99).counts.sum() == 2 + controller.supervisor.merge([0, 1]) + controller.supervisor.block() + assert view.get_data(100).counts.sum() == 7 + finally: + gui.close() + controller.close() + + +def test_observer_mutation_rolls_back_native_and_semantic_state(tmp_path): + _, controller = make_controller(tmp_path) + supervisor = controller.supervisor + once = [] + + @connect(event='cluster', sender=supervisor.clustering) + def mutate(sender, up): + if not once: + once.append(True) + sender.assign([0]) + + try: + with pytest.raises(ValueError, match='Native assignments'): + supervisor.merge([0, 1]) + np.testing.assert_array_equal( + supervisor.clustering.spike_clusters, controller.bounded_session.base_assignments + ) + assert not controller.bounded_session.events + assert supervisor._global_history.current_position == 0 + assert not supervisor.is_dirty() + finally: + unconnect(mutate) + controller.close() + + +def test_evidence_navigation_never_reuses_old_feature_cache(tmp_path): + _, controller = make_controller(tmp_path, sparse=True) + try: + first = controller._get_features(0, channel_ids=[0, 1], load_all=True) + assert first.spike_ids.tolist() == [0, 3] + controller.inspect_evidence_window(1) + assert not controller._get_feature_view_spike_ids(0, load_all=True).size + assert not controller._get_features(0, channel_ids=[0, 1], load_all=True) + finally: + controller.close() + + +def test_late_clock_evidence_preserves_submillisecond_positions(tmp_path): + params = make_gui_sorting(tmp_path / 'source') + origin = 1000 * 3600 * 30000 + times = np.load(params.parent / 'spike_times.npy') + origin + times[1:4] = origin + np.array([30, 31, 60]) + np.save(params.parent / 'spike_times.npy', times) + params.write_text( + params.read_text().replace('num_samples = 600000', f'num_samples = {origin + 600000}') + + f'valid_intervals = [[{origin}, {origin + 600000}]]\n' + ) + output = prepare_large( + params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20, bin_seconds=3600 + ) + controller = TemplateController(dir_path=output, config_dir=tmp_path / 'config', plugins=[]) + try: + positions = controller._window_time(np.array([1, 2])) + assert positions[1] - positions[0] == pytest.approx(1 / 30000) + assert np.float32(positions[0]) != np.float32(positions[1]) + finally: + controller.close() + + +def test_empty_window_and_nonfinite_amplitude_measurements(qtbot, tmp_path): + _, controller = make_controller(tmp_path, sparse=True) + view = controller.create_evidence_amplitude_view() + try: + controller.inspect_evidence_window(1) + view.on_select([0]) + assert view.data_bounds == (0, 0, view.duration, 1) + assert all(not len(b.amplitudes) for b in view._displayed_bunchs) + finally: + view.canvas.close() + controller.close() + + +def test_proposition_reviews_checkpoint_and_failed_review_history(tmp_path): + from phy.cluster._propositions import MergePropositionController, decode_curation_mapping + + _, controller = make_controller(tmp_path) + supervisor = controller.supervisor + mapping = {'format_version': '2', 'unit_ids': [0, 1, 2, 99], 'merges': [{'unit_ids': [0, 1]}]} + producer = controller.dir_path / 'curation.json' + producer.write_text(json.dumps(mapping)) + supervisor.merge_propositions = MergePropositionController(decode_curation_mapping(mapping)) + key = supervisor.merge_propositions.catalog.propositions[0].key + try: + with pytest.raises(RuntimeError, match='observer'), supervisor._transaction(): + supervisor.merge_propositions.accept(key, [0, 1], 100) + raise RuntimeError('observer') + assert supervisor.merge_propositions._history.current_position == 0 + assert not supervisor.merge_propositions.catalog.reviews + supervisor.merge_propositions.reject(key) + supervisor.save() + controller.close() + reopened = TemplateController( + dir_path=producer.parent, config_dir=tmp_path / 'config', plugins=[] + ) + try: + assert ( + reopened.supervisor.merge_propositions.catalog.status_for(key).value == 'rejected' + ) + assert not reopened.supervisor.is_dirty() + finally: + reopened.close() + finally: + controller.close() + + +def test_nonfinite_source_amplitudes_are_omitted_from_plots_and_kept_in_replay(qtbot, tmp_path): + params = make_gui_sorting(tmp_path / 'source') + amplitudes = np.load(params.parent / 'amplitudes.npy') + amplitudes[3] = np.inf + amplitudes[6] = np.nan + np.save(params.parent / 'amplitudes.npy', amplitudes) + output = prepare_large(params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20) + controller = TemplateController(dir_path=output, config_dir=tmp_path / 'config', plugins=[]) + view = controller.create_evidence_amplitude_view() + try: + view.on_select([0]) + assert np.all(np.isfinite(view.data_bounds)) + assert all(np.all(np.isfinite(b.amplitudes)) for b in view._displayed_bunchs) + controller.supervisor.split_source_amplitude(0, 4) + controller.supervisor.save() + curated = apply_curation(params, output, tmp_path / 'curated') + labels = np.load(curated / 'spike_clusters.npy') + assert labels[0] == 100 + assert labels[3] == labels[6] == 101 + finally: + view.canvas.close() + controller.close() diff --git a/phy/apps/large_recording/tests/test_curation.py b/phy/apps/large_recording/tests/test_curation.py new file mode 100644 index 000000000..9997bc9cf --- /dev/null +++ b/phy/apps/large_recording/tests/test_curation.py @@ -0,0 +1,427 @@ +"""Exact bounded journal projection and transactional persistence regressions.""" + +import json +import shutil + +import numpy as np +import pytest + +from phy.apps.large_recording.curation import BASE_ASSIGNMENTS, BoundedSession +from phy.apps.large_recording.source import Source, hash_file +from phy.apps.tests.large_recording_fixtures import make_long_recording_fixture + + +def make_session(tmp_path): + fixture = make_long_recording_fixture(tmp_path) + # The highest complete-source unit is intentionally absent from evidence. + clusters = np.load(fixture.full / 'spike_clusters.npy') + clusters[-1] = 99 + np.save(fixture.full / 'spike_clusters.npy', clusters) + (fixture.full / 'cluster_group.tsv').write_text('cluster_id\tgroup\n0\tgood\n1\tnoise\n') + source = Source(fixture.full / 'params.py', batch_size=5) + rows = np.array([0, 2, 5, 8, 13, 21, 25], dtype=np.int64) + for name in ('spike_times.npy', 'spike_clusters.npy', 'spike_templates.npy', 'amplitudes.npy'): + np.save(fixture.bounded / name, source.arrays[name][rows]) + np.save(fixture.bounded / BASE_ASSIGNMENTS, clusters[rows]) + np.save(fixture.bounded / '_phy_spike_global_ids.npy', rows) + np.save(fixture.bounded / 'windows.npy', np.array([[0, source.num_samples]], dtype=np.int64)) + manifest = { + 'format': 'phy_bounded_curation', + 'format_version': 1, + 'source': { + 'identity': source.identity, + 'revision': source.revision, + 'path': str(source.params_path), + }, + 'sample_rate': source.sample_rate, + 'num_samples': source.num_samples, + 'valid_intervals': source.valid_intervals.tolist(), + 'num_spikes': source.n_spikes, + 'bounded_num_spikes': len(rows), + 'unit_ids': source.unit_ids.tolist(), + 'counts': source.counts.tolist(), + 'next_cluster_id': 100, + 'global_spike_ids_file': '_phy_spike_global_ids.npy', + 'windows_file': 'windows.npy', + 'base_metadata': source.metadata, + 'capabilities': { + 'merge': True, + 'label': True, + 'split': False, + 'predicate_split_attributes': ['amplitudes.npy'], + }, + 'bounded_files': { + name: hash_file(fixture.bounded / name) + for name in ( + BASE_ASSIGNMENTS, + '_phy_spike_global_ids.npy', + 'windows.npy', + 'spike_times.npy', + ) + }, + } + (fixture.bounded / '_phy_curation_manifest.json').write_text(json.dumps(manifest)) + return source, BoundedSession(fixture.bounded) + + +def test_merge_label_snapshot_branch_reopen(tmp_path): + source, session = make_session(tmp_path) + baseline = session.snapshot() + session.record_merge([0, 1], 100, {'group': 'good'}) + first = session.snapshot() + session.record_label('group', [100], 'mua') + second = session.snapshot() + session.restore(first) + assert session.metadata['group'][100] == 'good' + session.restore(second) + assert session.metadata['group'][100] == 'mua' + session.restore(baseline) + session.record_merge([1, 2], 100, {'group': 'noise'}) + session.record_label('custom', [99, 100], 'reviewed') + session.save(extra_documents={'curation_review.json': {'accepted': [100]}}) + reopened = BoundedSession(session.path) + assert reopened.next_cluster_id == 101 + assert reopened.events == session.events + assert reopened.metadata == session.metadata + assert np.array_equal(reopened.assignments, session.assignments) + assert reopened.extra_documents == {'curation_review.json': {'accepted': [100]}} + assert source.unit_ids.tolist() == [0, 1, 2, 99] + assert 99 not in session.base_assignments + + +@pytest.mark.parametrize('parents,target', [([0, 1], 99), ([0, 0], 100), ([0, 3], 100)]) +def test_merge_collision_and_invalid_parents_do_not_mutate(tmp_path, parents, target): + _, session = make_session(tmp_path) + before = session.snapshot() + with pytest.raises(ValueError): + session.record_merge(parents, target) + assert session.events == before['events'] + assert np.array_equal(session.assignments, before['assignments']) + + +def test_save_requires_exact_semantic_projection(tmp_path): + _, session = make_session(tmp_path) + session.record_merge([0, 1], 100, {'group': 'good'}) + with pytest.raises(ValueError, match='projection'): + session.save(assignments=session.base_assignments) + with pytest.raises(ValueError, match='metadata'): + session.save(metadata={}) + assert not (session.path / '.phy-curation' / 'current.json').exists() + + +@pytest.mark.parametrize('fail_at', ['journal.json', 'checkpoint.json', 'pointer']) +def test_interrupted_save_recovers_previous_generation(tmp_path, monkeypatch, fail_at): + from phy.apps.large_recording import curation + + _, session = make_session(tmp_path) + session.record_merge([0, 1], 100, {'group': 'good'}) + session.save() + previous = session.snapshot() + session.record_label('group', [100], 'mua') + original = curation._write_json + + def fail(path, value): + if path.name == fail_at or (fail_at == 'pointer' and path.name.startswith('.current-')): + raise OSError('injected interrupted write') + return original(path, value) + + monkeypatch.setattr(curation, '_write_json', fail) + with pytest.raises(OSError, match='interrupted'): + session.save() + reopened = BoundedSession(session.path) + assert reopened.events == previous['events'] + assert reopened.metadata == previous['metadata'] + assert session.metadata['group'][100] == 'mua' + + +def test_checkpoint_corruption_is_rejected(tmp_path): + _, session = make_session(tmp_path) + session.save() + (session.checkpoint_path / 'journal.json').write_text('[]\nchanged') + with pytest.raises(ValueError, match='corruption'): + BoundedSession(session.path) + + +def test_bound_base_corruption_is_rejected(tmp_path): + _, session = make_session(tmp_path) + session.close() + shutil.copyfile(session.path / 'spike_times.npy', session.path / BASE_ASSIGNMENTS) + with pytest.raises(ValueError, match='integrity'): + BoundedSession(session.path) + + +def test_threshold_projection_is_serialized_and_nonfinite_policy_explicit(tmp_path): + _, session = make_session(tmp_path) + predicate = { + 'attribute': 'amplitudes.npy', + 'operator': '<', + 'threshold': 15.0, + 'nonfinite': 'keep', + } + session.record_split(2, 100, 101, predicate, {'group': {100: 'noise', 101: 'good'}}) + expected = session.base_assignments.astype(np.int64).copy() + mask = expected == 2 + selected = np.load(session.path / 'amplitudes.npy') < 15 + expected[mask & selected] = 100 + expected[mask & ~selected] = 101 + assert np.array_equal(session.assignments, expected) + session.save() + assert BoundedSession(session.path).events == session.events + with pytest.raises(ValueError, match='Unsupported'): + session.record_split(0, 102, 103, dict(predicate, attribute='waveform_peak')) + + +@pytest.mark.parametrize('boundary', ['generation', 'pointer']) +def test_atomic_publication_failure_preserves_previous_commit(tmp_path, monkeypatch, boundary): + from phy.apps.large_recording import curation + + _, session = make_session(tmp_path) + session.save() + session.record_merge([0, 1], 100, {'group': 'good'}) + original = curation.os.replace + + def fail(source, destination): + if (boundary == 'generation' and destination.parent.name == 'generations') or ( + boundary == 'pointer' and destination.name == 'current.json' + ): + raise OSError('injected atomic publication failure') + return original(source, destination) + + monkeypatch.setattr(curation.os, 'replace', fail) + with pytest.raises(OSError, match='publication'): + session.save() + assert BoundedSession(session.path).events == [] + assert len(session.events) == 1 + + +def test_nonfinite_predicate_keep_and_reject_are_explicit(): + from phy.apps.large_recording.curation import evaluate_predicate + + values = np.array([1.0, np.nan, np.inf, -np.inf, 3.0]) + predicate = { + 'attribute': 'amplitudes.npy', + 'operator': '<', + 'threshold': 2.0, + 'nonfinite': 'keep', + } + assert evaluate_predicate(values, predicate).tolist() == [True, False, False, False, False] + with pytest.raises(ValueError, match='Nonfinite'): + evaluate_predicate(values, dict(predicate, nonfinite='reject')) + + +def test_semantically_corrupt_checkpoint_rejected_even_with_valid_file_hashes(tmp_path): + from phy.apps.large_recording.curation import _digest + + _, session = make_session(tmp_path) + session.record_merge([0, 1], 100, {'group': 'good'}) + session.save() + generation = session.checkpoint_path + journal = json.loads((generation / 'journal.json').read_text()) + journal[0]['target'] = 99 + (generation / 'journal.json').write_text(json.dumps(journal)) + descriptor = json.loads((generation / 'checkpoint.json').read_text()) + descriptor['files']['journal.json'] = hash_file(generation / 'journal.json') + (generation / 'checkpoint.json').write_text(json.dumps(descriptor)) + pointer = {'generation': generation.name, 'sha256': _digest(descriptor)} + (session.path / '.phy-curation' / 'current.json').write_text(json.dumps(pointer)) + with pytest.raises(ValueError, match='colliding'): + BoundedSession(session.path) + + +def test_compact_history_preserves_allocator_and_reconstructs_assignments(tmp_path): + _, session = make_session(tmp_path) + before = session.journal_state() + session.record_merge([0, 1], 100, {'group': 'good'}) + after = session.journal_state() + merged = session.assignments.copy() + assert set(after) == {'events', 'next_cluster_id'} + assert after['events'][0] is session.events[0] + session.restore_journal(before) + assert session.events == [] + assert session.next_cluster_id == 101 + assert np.array_equal(session.assignments, session.base_assignments) + session.restore_journal(after) + assert np.array_equal(session.assignments, merged) + session.restore_journal(before) + session.record_merge([1, 2], 101, {'group': 'noise'}) + assert session.next_cluster_id == 102 + + +def test_session_close_releases_own_file_mappings(tmp_path): + _, session = make_session(tmp_path) + mappings = [ + array._mmap + for array in ( + session.global_ids, + session.base_assignments, + *session.attributes.values(), + ) + ] + session.close() + assert all(mapping.closed for mapping in mappings) + session.close() + + +def test_stale_session_save_rejected_without_losing_first_commit(tmp_path): + _, first = make_session(tmp_path) + second = BoundedSession(first.path) + first.record_merge([0, 1], 100, {'group': 'good'}) + second.record_merge([1, 2], 100, {'group': 'noise'}) + first.save(extra_documents={'curation_review.json': {'accepted': [100]}}) + with pytest.raises(ValueError, match='reopen'): + second.save(extra_documents={'curation_review.json': {'rejected': [100]}}) + reopened = BoundedSession(first.path) + assert reopened.events == first.events + assert reopened.extra_documents == first.extra_documents + assert second.events != reopened.events + # The winning session can commit again without reopening. + first.record_label('group', [100], 'mua') + first.save() + assert BoundedSession(first.path).metadata['group'][100] == 'mua' + + +def test_failed_save_releases_commit_lock(tmp_path, monkeypatch): + from phy.apps.large_recording import curation + + _, session = make_session(tmp_path) + session.save() + other = BoundedSession(session.path) + original = curation._write_json + + def fail(path, value): + if path.name == 'journal.json': + raise OSError('injected locked save failure') + return original(path, value) + + with monkeypatch.context() as patch: + patch.setattr(curation, '_write_json', fail) + with pytest.raises(OSError, match='locked'): + session.save() + other.record_label('group', [0], 'mua') + other.save() + assert BoundedSession(session.path).metadata['group'][0] == 'mua' + + +def _concurrent_save_worker(path, barrier, results, label): + session = BoundedSession(path) + try: + session.record_label('group', [0], label) + barrier.wait(timeout=20) + try: + session.save() + except ValueError as error: + results.put(('conflict', str(error))) + else: + results.put(('saved', label)) + finally: + session.close() + + +def test_concurrent_process_saves_compare_generation_under_lock(tmp_path): + import multiprocessing + + _, session = make_session(tmp_path) + session.save() + context = multiprocessing.get_context('spawn') + barrier = context.Barrier(2) + results = context.Queue() + processes = [ + context.Process( + target=_concurrent_save_worker, + args=(str(session.path), barrier, results, label), + ) + for label in ('good', 'noise') + ] + for process in processes: + process.start() + try: + outcomes = [results.get(timeout=30) for _ in processes] + for process in processes: + process.join(timeout=10) + assert process.exitcode == 0 + assert sorted(status for status, _ in outcomes) == ['conflict', 'saved'] + winner = next(value for status, value in outcomes if status == 'saved') + assert 'reopen' in next(value for status, value in outcomes if status == 'conflict') + assert BoundedSession(session.path).metadata['group'][0] == winner + finally: + for process in processes: + if process.is_alive(): + process.terminate() + process.join(timeout=5) + results.close() + + +def _replace_capabilities(session, capabilities): + session.close() + path = session.path / '_phy_curation_manifest.json' + manifest = json.loads(path.read_text()) + if capabilities is None: + manifest.pop('capabilities', None) + else: + manifest['capabilities'] = capabilities + path.write_text(json.dumps(manifest)) + return BoundedSession(session.path) + + +@pytest.mark.parametrize( + 'capabilities', [None, {}, {'merge': False, 'label': False, 'split': False}] +) +def test_readonly_capabilities_reject_public_mutations_without_state_changes( + tmp_path, capabilities +): + _, session = make_session(tmp_path) + session = _replace_capabilities(session, capabilities) + before = session.snapshot() + predicate = { + 'attribute': 'amplitudes.npy', + 'operator': '<', + 'threshold': 15.0, + 'nonfinite': 'keep', + } + commands = [ + lambda: session.preflight_merge([0, 1], 100), + lambda: session.record_merge([0, 1], 100), + lambda: session.preflight_label('group', [0], 'good'), + lambda: session.record_label('group', [0], 'good'), + lambda: session.preflight_split(2, 100, 101, predicate), + lambda: session.record_split(2, 100, 101, predicate), + lambda: session.record_threshold_split(2, [101, 100], predicate), + ] + for command in commands: + with pytest.raises(ValueError, match='capability'): + command() + assert session.events == before['events'] + assert session.metadata == before['metadata'] + assert session.next_cluster_id == before['next_cluster_id'] + assert np.array_equal(session.assignments, before['assignments']) + + +@pytest.mark.parametrize( + 'capabilities', + [ + {'merge': 'true'}, + {'label': 'false'}, + {'stored_waveforms': 1}, + {'raw_traces': None}, + {'split': True}, + {'predicate_split_attributes': 'amplitudes.npy'}, + {'predicate_split_attributes': ['amplitudes.npy', 'amplitudes.npy']}, + {'predicate_split_attributes': ['spike_templates.npy']}, + {'predicate_split_attributes': ['waveform_peak']}, + {'predicate_split_attributes': [True]}, + {'full_attributes': 'amplitude'}, + ], +) +def test_malformed_or_unsupported_capabilities_fail_validation(tmp_path, capabilities): + _, session = make_session(tmp_path) + with pytest.raises(ValueError, match='capability|capabilities|unsupported|Unsupported'): + _replace_capabilities(session, capabilities) + + +def test_predicate_capability_requires_an_available_scalar_source_attribute(tmp_path): + _, session = make_session(tmp_path) + session.close() + values = np.load(session.path / 'amplitudes.npy') + np.save(session.path / 'amplitudes.npy', np.column_stack([values, values])) + with pytest.raises(ValueError, match='source scalar'): + BoundedSession(session.path) diff --git a/phy/apps/large_recording/tests/test_gui_lifecycle.py b/phy/apps/large_recording/tests/test_gui_lifecycle.py new file mode 100644 index 000000000..d2286b9eb --- /dev/null +++ b/phy/apps/large_recording/tests/test_gui_lifecycle.py @@ -0,0 +1,169 @@ +"""GUI transaction, proposition review, and pending-metric regressions.""" + +import json + +import numpy as np +import pytest +from phylib.utils import connect, unconnect + +from phy.apps.large_recording.prepare import prepare_large +from phy.apps.large_recording.tests.test_adapter import make_controller, make_gui_sorting +from phy.apps.template.gui import TemplateController + + +def test_window_time_labels_follow_navigation(qtbot, tmp_path): + _, controller = make_controller(tmp_path, sparse=True) + gui = controller.create_gui(default_views=['EvidenceAmplitudeView'], do_prompt_save=False) + try: + view = gui.get_view('EvidenceAmplitudeView') + assert view.time_axis_label.text() == 'Seconds from window start' + initial = view.window_label.text() + controller.inspect_evidence_window(1) + start, end = controller.current_evidence_window / controller.model.sample_rate + assert view.window_label.text() == ( + f'Evidence window 1 | Recording time: {start:.3f}–{end:.3f} s' + ) + assert view.window_label.text() != initial + assert view.duration == end - start + finally: + gui.close() + controller.close() + + +def test_gui_observer_mutation_rolls_back_views_selection_and_journal(qtbot, tmp_path): + _, controller = make_controller(tmp_path) + gui = controller.create_gui(default_views=[], do_prompt_save=False) + supervisor = controller.supervisor + once = [] + + try: + supervisor.select([0]) + supervisor.block() + assignments = supervisor.clustering.spike_clusters.copy() + events = list(controller.bounded_session.events) + selection = supervisor.selection.snapshot() + cluster_ids = supervisor.cluster_view.get_ids() + similarity_ids = supervisor.similarity_view.get_ids() + histories = tuple( + (history._index, len(history._history)) + for history in ( + supervisor.clustering._undo_stack, + supervisor.cluster_meta._undo_stack, + supervisor._global_history, + supervisor.session_history, + ) + ) + dirty = supervisor.is_dirty() + + @connect(event='cluster', sender=supervisor.clustering) + def mutate(sender, up): + if not once: + once.append(True) + sender.assign([0]) + + with pytest.raises((ValueError, KeyError)): + supervisor.merge([0, 1]) + supervisor.block() + + np.testing.assert_array_equal(supervisor.clustering.spike_clusters, assignments) + assert controller.bounded_session.events == events + assert supervisor.selection.snapshot() == selection + assert supervisor.cluster_view.get_ids() == cluster_ids + assert supervisor.similarity_view.get_ids() == similarity_ids + assert ( + tuple( + (history._index, len(history._history)) + for history in ( + supervisor.clustering._undo_stack, + supervisor.cluster_meta._undo_stack, + supervisor._global_history, + supervisor.session_history, + ) + ) + == histories + ) + assert supervisor.is_dirty() == dirty + finally: + if 'mutate' in locals(): + unconnect(mutate) + gui.close() + controller.close() + + +def test_gui_proposition_review_undo_redo_and_checkpoint_reopen(qtbot, tmp_path): + params = make_gui_sorting(tmp_path / 'source') + proposal = { + 'format_version': '2', + 'unit_ids': [0, 1, 2, 99], + 'merges': [{'unit_ids': [0, 1]}], + } + (params.parent / 'curation.json').write_text(json.dumps(proposal), encoding='utf-8') + bounded = prepare_large( + params, + tmp_path / 'bounded', + window_seconds=20, + stratum_seconds=20, + bin_seconds=5, + batch_size=2, + ) + config = tmp_path / 'config' + controller = TemplateController( + dir_path=bounded, config_dir=config, plugins=[], enable_threading=False + ) + gui = controller.create_gui(default_views=[], do_prompt_save=False) + supervisor = controller.supervisor + proposition = supervisor.merge_propositions.catalog.propositions[0] + key = proposition.key + + try: + supervisor._review_merge_proposition(supervisor.merge_propositions_view, key) + assert supervisor.selection.state.is_merge_mode + assert supervisor.selected_merge == [0, 1] + + merged = supervisor.merge() + supervisor.block() + assert merged.added == [100] + assert supervisor.merge_propositions.catalog.status_for(key).value == 'accepted' + assert controller.bounded_session.events[-1]['target'] == 100 + + supervisor.undo() + supervisor.block() + assert supervisor.merge_propositions.catalog.status_for(key).value == 'pending' + assert controller.bounded_session.events == [] + + supervisor.redo() + supervisor.block() + assert supervisor.merge_propositions.catalog.status_for(key).value == 'accepted' + assert controller.bounded_session.events[-1]['target'] == 100 + supervisor.save() + finally: + gui.close() + controller.close() + + reopened = TemplateController( + dir_path=bounded, config_dir=config, plugins=[], enable_threading=False + ) + try: + key = reopened.supervisor.merge_propositions.catalog.propositions[0].key + assert reopened.supervisor.merge_propositions.catalog.status_for(key).value == 'accepted' + assert reopened.bounded_session.events[-1]['target'] == 100 + finally: + reopened.close() + + +def test_default_bounded_views_handle_threshold_children_with_pending_metrics(qtbot, tmp_path): + _, controller = make_controller(tmp_path) + gui = controller.create_gui(do_prompt_save=False) + try: + supervisor = controller.supervisor + update = supervisor.split_source_amplitude(0, 4) + supervisor.block() + assert set(update.added) == {100, 101} + for cluster_id in update.added: + row = supervisor.cluster_view._model.row_by_id(cluster_id) + assert row['fr'] is None + assert row['full_count'] is None + assert row['evidence'] == 'pending replay' + finally: + gui.close() + controller.close() diff --git a/phy/apps/large_recording/tests/test_prepare.py b/phy/apps/large_recording/tests/test_prepare.py new file mode 100644 index 000000000..a3688fd8c --- /dev/null +++ b/phy/apps/large_recording/tests/test_prepare.py @@ -0,0 +1,218 @@ +"""Production bounded preparation: budgets, rows, waveforms, and publication.""" + +import json + +import numpy as np +import pytest +from phylib.io.model import load_model + +from phy.apps.large_recording.prepare import prepare_large +from phy.apps.large_recording.source import PER_SPIKE, hash_file +from phy.apps.large_recording.tests.test_source import make_sorting + + +@pytest.mark.parametrize('channel_map', [[0, 1, 2], [2, 0, 1]]) +def test_prepare_exact_arrays_waveforms_and_source_unchanged(tmp_path, channel_map): + params = make_sorting(tmp_path / 'full', columns=True) + np.save(params.parent / 'channel_map.npy', np.asarray(channel_map, dtype=np.int32)) + before = {p.name: hash_file(p) for p in params.parent.iterdir()} + output = prepare_large( + params, + tmp_path / 'bounded', + max_spikes=8, + window_seconds=3, + stratum_seconds=20, + seed=0, + bin_seconds=5, + batch_size=2, + ) + manifest = json.loads((output / '_phy_curation_manifest.json').read_text()) + windows = np.load(output / manifest['windows_file']) + samples = np.load(params.parent / 'spike_times.npy').reshape(-1) + expected_ids = np.array( + [ + index + for index, sample in enumerate(samples) + if any(start <= sample < end for start, end in windows) + ], + dtype=np.int64, + ) + global_ids = np.load(output / manifest['global_spike_ids_file']) + np.testing.assert_array_equal(global_ids, expected_ids) + assert len(global_ids) <= 8 + for name in PER_SPIKE: + if (params.parent / name).exists(): + np.testing.assert_array_equal( + np.load(output / name), np.load(params.parent / name)[expected_ids] + ) + source_waveform_ids = np.load(params.parent / '_phy_spikes_subset.spikes.npy') + expected_waveform_rows = [ + row for row, spike in enumerate(source_waveform_ids) if spike in expected_ids + ] + np.testing.assert_array_equal( + global_ids[np.load(output / '_phy_spikes_subset.spikes.npy')], + source_waveform_ids[expected_waveform_rows], + ) + for suffix in ('channels', 'waveforms'): + np.testing.assert_array_equal( + np.load(output / f'_phy_spikes_subset.{suffix}.npy'), + np.load(params.parent / f'_phy_spikes_subset.{suffix}.npy')[expected_waveform_rows], + ) + assert before == {p.name: hash_file(p) for p in params.parent.iterdir()} + assert 'dat_path = None' in (output / 'params.py').read_text() + assert manifest['next_cluster_id'] == 100 + assert sum(manifest['counts']) == len(samples) + assert sum(manifest['bounded_counts']) == len(global_ids) + model = load_model(output / 'params.py') + try: + ids = np.load(output / '_phy_spikes_subset.spikes.npy') + channels = np.load(output / '_phy_spikes_subset.channels.npy') + sparse = np.load(output / '_phy_spikes_subset.waveforms.npy') + expected = np.zeros((len(ids), 8, 3), dtype=np.float32) + for row, slots in enumerate(channels): + for slot, channel in enumerate(slots): + if channel >= 0: + expected[row, :, channel] = sparse[row, :, slot] + np.testing.assert_array_equal(model.get_waveforms(ids, channel_ids=[0, 1, 2]), expected) + finally: + model.close() + + +def test_deterministic_budgets_and_zero_evidence_catalog(tmp_path): + params = make_sorting(tmp_path / 'full') + outputs = [ + prepare_large( + params, + tmp_path / name, + max_spikes=5, + window_seconds=3, + stratum_seconds=20, + seed=0, + bin_seconds=5, + batch_size=2, + ) + for name in ('one', 'two') + ] + manifests = [ + json.loads((path / '_phy_curation_manifest.json').read_text()) for path in outputs + ] + assert manifests[0] == manifests[1] + assert manifests[0]['unit_ids'][-1] == 99 + assert manifests[0]['bounded_counts'][-1] == 0 + assert manifests[0]['omitted_windows'] + for path in outputs: + ids = np.load(path / '_phy_spike_global_ids.npy') + times = np.load(path / 'spike_times.npy') + assert len(ids) <= 5 + # All three tied samples at 1 retain their separate source identities. + assert np.count_nonzero(times == 1) == 2 + + +def test_rejected_budget_or_interrupted_publication_leaves_no_output(tmp_path, monkeypatch): + params = make_sorting(tmp_path / 'full') + with pytest.raises(ValueError, match='No bounded evidence'): + prepare_large( + params, tmp_path / 'empty', max_spikes=1, window_seconds=3, stratum_seconds=20, seed=0 + ) + with pytest.raises(ValueError, match='Output budget'): + prepare_large(params, tmp_path / 'too-big', max_bytes=10) + import phy.apps.large_recording.prepare as module + + def fail(*args, **kwargs): + raise OSError('simulated interruption') + + monkeypatch.setattr(module, 'build_summaries', fail) + with pytest.raises(OSError, match='interruption'): + prepare_large(params, tmp_path / 'interrupted', window_seconds=20, stratum_seconds=20) + assert not (tmp_path / 'interrupted').exists() + assert not list(tmp_path.glob('.interrupted.prepare-*')) + + +def test_prepare_missing_waveforms_is_explicit(tmp_path): + params = make_sorting(tmp_path / 'full') + for path in params.parent.glob('_phy_spikes_subset.*.npy'): + path.unlink() + output = prepare_large(params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20) + manifest = json.loads((output / '_phy_curation_manifest.json').read_text()) + assert not manifest['capabilities']['stored_waveforms'] + assert not (output / '_phy_spikes_subset.waveforms.npy').exists() + + +def test_spikeinterface_singular_sparse_template_indices(tmp_path): + params = make_sorting(tmp_path / 'full') + templates = np.load(params.parent / 'templates.npy')[:, :, :2] + templates[1, :, 1] = 0 + np.save(params.parent / 'templates.npy', templates) + indices = np.array([[0, 2], [1, -1], [2, 0]], dtype=np.int32) + np.save(params.parent / 'template_ind.npy', indices) + output = prepare_large(params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20) + model = load_model(output / 'params.py') + try: + np.testing.assert_array_equal(model.sparse_templates.cols, indices) + finally: + model.close() + + +def test_adjacent_acquisition_blocks_retain_observation_boundaries(tmp_path): + params = make_sorting(tmp_path / 'full') + output = prepare_large( + params, + tmp_path / 'bounded', + valid_intervals=[[0, 10], [10, 20]], + window_seconds=10, + stratum_seconds=10, + ) + manifest = json.loads((output / '_phy_curation_manifest.json').read_text()) + np.testing.assert_array_equal(np.load(output / manifest['windows_file']), [[0, 10], [10, 20]]) + assert manifest['windows'][0]['block'] != manifest['windows'][1]['block'] + + +def test_disk_preflight_rejects_before_writing_package(tmp_path, monkeypatch): + from types import SimpleNamespace + + from phy.apps.large_recording import prepare + + params = make_sorting(tmp_path / 'full') + monkeypatch.setattr(prepare.shutil, 'disk_usage', lambda path: SimpleNamespace(free=0)) + with pytest.raises(OSError, match='Insufficient disk'): + prepare_large(params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20) + assert not (tmp_path / 'bounded').exists() + assert not list(tmp_path.glob('.bounded.prepare-*')) + + +def test_final_windows_records_match_retained_samples_and_gaps(tmp_path): + params = make_sorting(tmp_path / 'full') + output = prepare_large( + params, + tmp_path / 'bounded', + valid_intervals=[[0, 7], [10, 20]], + window_seconds=5, + stratum_seconds=20, + bin_seconds=5, + ) + manifest = json.loads((output / '_phy_curation_manifest.json').read_text()) + intervals = np.load(output / manifest['windows_file']) + assert [[row['start'], row['end']] for row in manifest['windows']] == intervals.tolist() + for row in manifest['windows']: + start, end = row['start'], row['end'] + assert any(first <= start and end <= last for first, last in [[0, 7], [10, 20]]) + assert row['exposure_seconds'] == end - start + assert manifest['capabilities']['predicate_split_attributes'] == ['amplitudes.npy'] + + +def test_prepare_protects_loader_and_preserves_producer_documents(tmp_path): + params = make_sorting(tmp_path / 'full') + for name, value in ( + ('curation.json', {'merge': [[0, 1]]}), + ('curation_review.json', {'accepted': []}), + ): + (params.parent / name).write_text(json.dumps(value)) + output = prepare_large(params, tmp_path / 'bounded', window_seconds=20, stratum_seconds=20) + manifest = json.loads((output / '_phy_curation_manifest.json').read_text()) + assert manifest['bounded_files']['params.py'] == hash_file(output / 'params.py') + assert manifest['bounded_files']['_phy_cluster_metadata.json'] == hash_file( + output / '_phy_cluster_metadata.json' + ) + for name in ('curation.json', 'curation_review.json'): + assert (output / name).read_bytes() == (params.parent / name).read_bytes() + assert name not in manifest['excluded_files'] diff --git a/phy/apps/large_recording/tests/test_replay.py b/phy/apps/large_recording/tests/test_replay.py new file mode 100644 index 000000000..408f536fd --- /dev/null +++ b/phy/apps/large_recording/tests/test_replay.py @@ -0,0 +1,144 @@ +"""Full assignment oracle, revision binding and atomic publication regressions.""" + +import json + +import numpy as np +import pytest + +from phy.apps.large_recording.replay import apply_curation +from phy.apps.large_recording.source import Source, hash_file +from phy.apps.large_recording.tests.test_curation import make_session + + +def test_full_replay_chains_preserve_unseen_units_and_immutable_bytes(tmp_path): + source, session = make_session(tmp_path) + original_hashes = {name: hash_file(source.path / name) for name in source.arrays} + original = source.arrays['spike_clusters.npy'].copy() + session.record_merge([0, 1], 100, {'group': 'good', 'confidence': 0.8}) + session.record_label('group', [100], 'mua') + session.record_merge([100, 2], 101, {'group': 'good', 'confidence': None}) + session.record_label('custom', [99, 101], 'verified') + session.save() + output = apply_curation(source.params_path, session.path, tmp_path / 'curated', batch_size=3) + expected = np.where(np.isin(original, [0, 1, 2]), 101, original) + assert np.array_equal(np.load(output / 'spike_clusters.npy'), expected) + assert expected[-1] == 99 + curated = Source(output / 'params.py') + assert curated.metadata == session.metadata + assert curated.identity == source.identity + assert curated.revision != source.revision + lineage = json.loads((output / '_phy_curation_lineage.json').read_text()) + assert lineage['parent_revision'] == source.revision + assert lineage['revision'] == curated.revision + for name, digest in original_hashes.items(): + assert hash_file(source.path / name) == digest + if name != 'spike_clusters.npy': + assert hash_file(output / name) == digest + with pytest.raises(ValueError, match='revision'): + apply_curation(output / 'params.py', session.path, tmp_path / 'twice') + assert not (tmp_path / 'twice').exists() + second = apply_curation(source.params_path, session.path, tmp_path / 'second', batch_size=7) + assert hash_file(second / 'spike_clusters.npy') == hash_file(output / 'spike_clusters.npy') + + +def test_source_revision_and_feature_identity_are_independently_bound(tmp_path): + source, session = make_session(tmp_path) + session.save() + clusters = source.arrays['spike_clusters.npy'].copy() + amplitudes = source.arrays['amplitudes.npy'].copy() + source.close() + clusters[0] = 2 + np.save(source.path / 'spike_clusters.npy', clusters) + with pytest.raises(ValueError, match='revision'): + apply_curation(source.params_path, session.path, tmp_path / 'wrong_revision') + amplitudes[0] += 1 + np.save(source.path / 'amplitudes.npy', amplitudes) + with pytest.raises(ValueError, match='identity'): + apply_curation(source.params_path, session.path, tmp_path / 'wrong_feature') + assert not (tmp_path / 'wrong_revision').exists() + assert not (tmp_path / 'wrong_feature').exists() + + +def test_full_threshold_split_matches_scalar_oracle(tmp_path): + source, session = make_session(tmp_path) + original = source.arrays['spike_clusters.npy'].copy() + amplitudes = source.arrays['amplitudes.npy'].copy() + predicate = { + 'attribute': 'amplitudes.npy', + 'operator': '<', + 'threshold': 15.0, + 'nonfinite': 'keep', + } + session.record_split(2, 100, 101, predicate, {'group': {100: 'noise', 101: 'good'}}) + session.record_merge([0, 100], 102, {'group': 'mua'}) + session.save() + expected = original.astype(np.int64) + expected[(original == 2) & (amplitudes < 15)] = 100 + expected[(original == 2) & ~(amplitudes < 15)] = 101 + expected[np.isin(expected, [0, 100])] = 102 + output = apply_curation(source.params_path, session.path, tmp_path / 'split', batch_size=4) + assert np.array_equal(np.load(output / 'spike_clusters.npy'), expected) + assert Source(output / 'params.py').metadata == session.metadata + + +def test_failed_replay_does_not_publish_output(tmp_path, monkeypatch): + from phy.apps.large_recording import replay + + source, session = make_session(tmp_path) + session.save() + original = replay._write_json + + def fail(path, value): + if path.name == '_phy_curation_lineage.json': + raise OSError('injected replay interruption') + return original(path, value) + + monkeypatch.setattr(replay, '_write_json', fail) + with pytest.raises(OSError, match='interruption'): + apply_curation(source.params_path, session.path, tmp_path / 'interrupted') + assert not (tmp_path / 'interrupted').exists() + assert not list(tmp_path.glob('.interrupted-*')) + + +def test_replay_requires_committed_checkpoint_and_new_output(tmp_path): + source, session = make_session(tmp_path) + with pytest.raises(ValueError, match='Save'): + apply_curation(source.params_path, session.path, tmp_path / 'unsaved') + session.save() + with pytest.raises(FileExistsError): + apply_curation(source.params_path, session.path, source.path) + + +def test_replay_preserves_empty_split_children_and_allocator_lineage(tmp_path): + source, session = make_session(tmp_path) + session.record_split( + 2, + 100, + 101, + { + 'attribute': 'amplitudes.npy', + 'operator': '<', + 'threshold': 100.0, + 'nonfinite': 'keep', + }, + ) + session.save() + output = apply_curation(source.params_path, session.path, tmp_path / 'empty_child') + curated = Source(output / 'params.py') + assert 101 in curated.unit_ids + assert curated.counts[curated.unit_ids.tolist().index(101)] == 0 + assert curated.next_cluster_id == 102 + + +def test_bounded_threshold_attribute_must_match_exact_source_rows(tmp_path): + source, session = make_session(tmp_path) + session.save() + session.close() + # Even a manifest that does not hash this optional attribute must not allow + # evidence values different from those used by the full-data predicate. + values = np.load(session.path / 'amplitudes.npy') + values[0] += 1 + np.save(session.path / 'amplitudes.npy', values) + with pytest.raises(ValueError, match='does not match source rows'): + apply_curation(source.params_path, session.path, tmp_path / 'bad_evidence') + assert not (tmp_path / 'bad_evidence').exists() diff --git a/phy/apps/large_recording/tests/test_source.py b/phy/apps/large_recording/tests/test_source.py new file mode 100644 index 000000000..05578bf9d --- /dev/null +++ b/phy/apps/large_recording/tests/test_source.py @@ -0,0 +1,295 @@ +"""Exact source clock, alignment, and identity regressions.""" + +import json + +import numpy as np +import pytest + +from phy.apps.large_recording.source import Source, hash_file + + +def make_sorting(path, columns=False): + path.mkdir() + samples = np.array([0, 1, 1, 2, 5, 6, 10, 11, 14, 18, 19, 19], dtype=np.int64) + clusters = np.array([0, 1, 2, 0, 1, 2, 0, 1, 2, 0, 99, 99], dtype=np.int32) + templates = np.arange(len(samples), dtype=np.int32) % 3 + arrays = { + 'spike_times.npy': samples, + 'spike_clusters.npy': clusters, + 'spike_templates.npy': templates, + 'amplitudes.npy': np.arange(len(samples), dtype=np.float32), + 'spike_positions.npy': np.column_stack((np.zeros(len(samples)), samples * 2)), + 'spike_depths.npy': samples.astype(np.float64) * 2, + 'templates.npy': np.arange(3 * 8 * 3, dtype=np.float32).reshape(3, 8, 3), + 'channel_map.npy': np.array([2, 0, 1], dtype=np.int32), + 'channel_positions.npy': np.array([[0, 0], [0, 10], [0, 20]]), + 'pc_features.npy': np.arange(len(samples) * 4, dtype=np.float32).reshape(-1, 2, 2), + 'pc_feature_ind.npy': np.array([[0, 2], [2, 1], [1, 0]], dtype=np.int32), + } + for name, array in arrays.items(): + if columns and name in ('spike_times.npy', 'spike_clusters.npy', 'spike_templates.npy'): + array = array[:, None] + np.save(path / name, array) + ids = np.array([0, 2, 3, 8, 11], dtype=np.int64) + channels = np.array([[0, -1], [2, 1], [1, -1], [2, -1], [0, 1]], dtype=np.int32) + waveforms = np.zeros((len(ids), 8, 2), dtype=np.float32) + for row, spike in enumerate(ids): + for slot, channel in enumerate(channels[row]): + if channel != -1: + waveforms[row, :, slot] = spike + channel + np.arange(8) / 10 + for suffix, array in (('spikes', ids), ('channels', channels), ('waveforms', waveforms)): + np.save(path / f'_phy_spikes_subset.{suffix}.npy', array) + (path / 'params.py').write_text( + "sample_rate = 1.0\nnum_samples = 20\ndat_path = None\nn_channels_dat = 3\ndtype = 'int16'\n" + ) + (path / 'cluster_group.tsv').write_text( + 'cluster_id\tgroup\n0\tgood\n1\tmua\n2\tgood\n99\tnoise\n' + ) + return path / 'params.py' + + +@pytest.mark.parametrize('columns', [False, True]) +def test_source_batches_ties_counts_columns(tmp_path, columns): + source = Source(make_sorting(tmp_path / 'full', columns), batch_size=2) + assert source.n_spikes == 12 + assert source.unit_ids.tolist() == [0, 1, 2, 99] + assert source.counts.tolist() == [4, 3, 3, 2] + assert all(batch['spike_times.npy'].ndim == 1 for batch in source.iter_batches()) + assert all(isinstance(array, np.memmap) for array in source.arrays.values()) + assert source.metadata['group'][99] == 'noise' + + +def test_declared_waveform_source_identity_is_checked(tmp_path): + params = make_sorting(tmp_path / 'full') + manifest = { + 'format': 'phy_spikes_subset', + 'format_version': 1, + 'recording': {'num_samples': 20, 'sampling_frequency_hz': 1}, + 'source_sorting': {'num_spikes': 12, 'spike_times_sha256': '0' * 64}, + } + path = params.parent / '_phy_spikes_subset.manifest.json' + path.write_text(json.dumps(manifest)) + with pytest.raises(ValueError, match='fingerprint'): + Source(params) + manifest['source_sorting']['spike_times_sha256'] = hash_file(params.parent / 'spike_times.npy') + path.write_text(json.dumps(manifest)) + Source(params).close() + manifest['source_sorting']['num_spikes'] = 13 + path.write_text(json.dumps(manifest)) + with pytest.raises(ValueError, match='spike count'): + Source(params) + + +def test_failed_source_validation_closes_readers_and_rejects_unsupported_ids( + tmp_path, monkeypatch +): + params = make_sorting(tmp_path / 'full') + labels = np.load(params.parent / 'spike_clusters.npy').astype(np.int64) + labels[-1] = 2**31 + np.save(params.parent / 'spike_clusters.npy', labels) + original = np.load + readers = [] + + def read(*args, **kwargs): + array = original(*args, **kwargs) + readers.append(array) + return array + + monkeypatch.setattr(np, 'load', read) + with pytest.raises(ValueError, match='int32'): + Source(params) + assert all(array._mmap.closed for array in readers) + + +def test_identity_assignment_and_feature_layers(tmp_path): + params = make_sorting(tmp_path / 'full') + source = Source(params) + source.close() + clusters = np.load(params.parent / 'spike_clusters.npy') + clusters[clusters == 1] = 0 + np.save(params.parent / 'spike_clusters.npy', clusters) + merged = Source(params) + merged.close() + assert merged.identity == source.identity + assert merged.revision != source.revision + values = np.load(params.parent / 'amplitudes.npy') + values[0] = 99 + np.save(params.parent / 'amplitudes.npy', values) + altered = Source(params) + altered.close() + assert altered.identity != merged.identity + assert altered.revision == merged.revision + assert ( + hash_file(params.parent / 'amplitudes.npy') == altered.identity['arrays']['amplitudes.npy'] + ) + + +def test_exposure_gaps_and_missing_duration(tmp_path): + params = make_sorting(tmp_path / 'full') + source = Source(params, valid_intervals=[[0, 7], [10, 20]]) + assert source.valid_intervals.tolist() == [[0, 7], [10, 20]] + with pytest.raises(ValueError, match='outside valid'): + Source(params, valid_intervals=[[0, 5], [10, 20]]) + params.write_text(params.read_text().replace('num_samples = 20\n', '')) + with pytest.raises(ValueError, match='exposure is missing'): + Source(params) + assert Source(params, duration_seconds=20).num_samples == 20 + with pytest.raises(ValueError, match='integer sample count'): + Source(params, duration_seconds=20.5) + + +@pytest.mark.parametrize( + 'failure', ['order', 'alignment', 'unknown', 'segment', 'waveforms', 'ids'] +) +def test_reject_ambiguous_or_misaligned_sources(tmp_path, failure): + params = make_sorting(tmp_path / 'full') + path = params.parent + if failure == 'order': + np.save(path / 'spike_times.npy', np.arange(12)[::-1]) + elif failure == 'alignment': + np.save(path / 'pc_features.npy', np.zeros((11, 2, 2))) + elif failure == 'unknown': + np.save(path / 'mysterious.npy', np.ones(12)) + elif failure == 'segment': + params.write_text(params.read_text() + 'num_segments = 2\n') + elif failure == 'waveforms': + (path / '_phy_spikes_subset.channels.npy').unlink() + else: + np.save(path / '_phy_spikes_subset.spikes.npy', np.array([0, 2, 2, 8, 11], dtype=np.int64)) + with pytest.raises(ValueError): + Source(params, batch_size=2) + + +def test_metadata_json_preserves_typed_values(tmp_path): + params = make_sorting(tmp_path / 'full') + metadata = {'group': {'0': 'good'}, 'custom': {'99': None, '2': 4}} + (params.parent / '_phy_cluster_metadata.json').write_text(json.dumps(metadata)) + assert Source(params).metadata == {'group': {0: 'good'}, 'custom': {99: None, 2: 4}} + + +def test_replay_lineage_retains_zero_count_catalog_and_allocator(tmp_path): + params = make_sorting(tmp_path / 'full') + source = Source(params) + lineage = { + 'format': 'phy_curation_lineage', + 'format_version': 1, + 'coordinate_identity': source.identity, + 'revision': source.revision, + 'unit_ids': [0, 1, 2, 99, 103], + 'next_cluster_id': 105, + } + (params.parent / '_phy_curation_lineage.json').write_text(json.dumps(lineage)) + reopened = Source(params) + assert reopened.next_cluster_id == 105 + assert reopened.unit_ids[-1] == 103 + assert reopened.counts[-1] == 0 + assert reopened.identity == source.identity + lineage['next_cluster_id'] = 103 + (params.parent / '_phy_curation_lineage.json').write_text(json.dumps(lineage)) + with pytest.raises(ValueError, match='collides'): + Source(params) + + +def test_unknown_upstream_manifest_version_is_rejected(tmp_path): + params = make_sorting(tmp_path / 'full') + (params.parent / '_phy_spikes_subset.manifest.json').write_text( + json.dumps({'format': 'phy_spikes_subset', 'format_version': 2}) + ) + with pytest.raises(ValueError, match='Unsupported waveform manifest'): + Source(params) + + +def test_integer_clock_beyond_32_bits(tmp_path): + from phy.apps.tests.large_recording_fixtures import make_long_recording_fixture + + fixture = make_long_recording_fixture(tmp_path) + source = Source(fixture.full / 'params.py', batch_size=3) + assert source.num_samples == 1000 * 3600 * 30000 + assert source.arrays['spike_times.npy'][-1] > np.iinfo(np.int32).max + assert source.n_spikes == 32 + + +def test_source_heap_does_not_scale_with_spike_arrays(tmp_path): + import tracemalloc + + peaks = [] + for spikes in (1000, 250000): + path = tmp_path / str(spikes) + path.mkdir() + for name, dtype in ( + ('spike_times.npy', np.int64), + ('spike_clusters.npy', np.int32), + ('spike_templates.npy', np.int32), + ): + values = np.lib.format.open_memmap( + path / name, mode='w+', dtype=dtype, shape=(spikes,) + ) + for start in range(0, spikes, 1024): + stop = min(start + 1024, spikes) + values[start:stop] = np.arange(start, stop) if name == 'spike_times.npy' else 0 + values.flush() + del values + (path / 'params.py').write_text(f'sample_rate = 1\nnum_samples = {spikes}\n') + tracemalloc.start() + source = Source(path / 'params.py', batch_size=1024) + _, peak = tracemalloc.get_traced_memory() + tracemalloc.stop() + assert source.n_spikes == spikes + assert sum(source.counts) == spikes + peaks.append(peak) + # This checks allocated heap; file-backed OS pages are intentionally separate. + assert peaks[-1] < 4 * 1024 * 1024 + assert peaks[-1] - peaks[0] < 3 * 1024 * 1024 + + +def test_explicit_intervals_cannot_override_changed_producer_exposure(tmp_path): + params = make_sorting(tmp_path / 'full') + params.write_text(params.read_text() + 'valid_intervals = [[0, 7], [10, 20]]\n') + with pytest.raises(ValueError, match='disagree with acquisition metadata'): + Source(params, valid_intervals=[[0, 20]]) + with Source(params) as source: + arrays = list(source.arrays.values()) + assert source.valid_intervals.tolist() == [[0, 7], [10, 20]] + assert not source.arrays + assert all(array._mmap.closed for array in arrays) + + +def test_exporter_channel_sidecars_and_raw_map_validation(tmp_path): + params = make_sorting(tmp_path / 'full') + np.save( + params.parent / 'channel_map_si.npy', np.array(['hardware2', 'hardware0', 'hardware1']) + ) + np.save(params.parent / 'channel_groups.npy', np.array([1, 1, 0])) + with Source(params) as source: + assert 'channel_map_si.npy' in source.arrays + params.write_text(params.read_text().replace('n_channels_dat = 3', 'n_channels_dat = 2')) + with pytest.raises(ValueError, match='outside n_channels_dat'): + Source(params) + + +@pytest.mark.parametrize( + 'positions', [np.zeros((3, 1)), np.zeros((3, 3)), np.full((3, 2), np.nan)] +) +def test_channel_geometry_must_match_gui_contract(tmp_path, positions): + params = make_sorting(tmp_path / 'source') + np.save(params.parent / 'channel_positions.npy', positions) + with pytest.raises(ValueError, match='finite x/y'): + Source(params) + + +def test_exclusive_publish_preserves_racing_empty_destination(tmp_path): + from phy.apps.large_recording.source import publish_directory + + staged = tmp_path / 'stage' + staged.mkdir() + (staged / 'payload').write_text('complete') + destination = tmp_path / 'destination' + destination.mkdir() + with pytest.raises(FileExistsError): + publish_directory(staged, destination) + assert destination.is_dir() and not list(destination.iterdir()) + assert (staged / 'payload').read_text() == 'complete' + destination.rmdir() + publish_directory(staged, destination) + assert not staged.exists() + assert (destination / 'payload').read_text() == 'complete' diff --git a/phy/apps/large_recording/tests/test_statistics.py b/phy/apps/large_recording/tests/test_statistics.py new file mode 100644 index 000000000..1bee50516 --- /dev/null +++ b/phy/apps/large_recording/tests/test_statistics.py @@ -0,0 +1,43 @@ +"""Independent pair-count oracle and observation-boundary regressions.""" + +import numpy as np + +from phy.apps.large_recording.statistics import window_correlograms, window_isi + + +def test_window_isi_excludes_missing_observation_periods_and_preserves_ties(): + samples = np.array([0, 1, 1, 2, 1000, 1002], dtype=np.int64) + windows = np.array([[0, 3], [1000, 1003]], dtype=np.int64) + np.testing.assert_array_equal(window_isi(samples, windows, 1), [1, 0, 1, 2]) + + +def test_late_sample_ccg_matches_all_pairs_oracle_with_sparse_ids(): + origin = 1000 * 3600 * 30000 + samples = origin + np.array([0, 1, 30, 30, 31, 60, 300, 301, 330], dtype=np.int64) + clusters = np.array([2_000_000_000, 2, 2, 2_000_000_000, 2, 2, 2, 2_000_000_000, 2]) + windows = np.array([[origin, origin + 100], [origin + 300, origin + 400]], dtype=np.int64) + ids = [2_000_000_000, 2] + actual = window_correlograms(samples, clusters, windows, ids, 30000, 0.001, 0.004) + expected = np.zeros((2, 2, 5), dtype=np.int64) + lookup = {cluster: index for index, cluster in enumerate(ids)} + for start, end in windows: + rows = [i for i, tick in enumerate(samples) if start <= tick < end] + for index, left in enumerate(rows): + for right in rows[index + 1 :]: + delay = int(samples[right] - samples[left]) // 30 + if delay > 2: + continue + a, b = lookup[clusters[left]], lookup[clusters[right]] + expected[a, b, 2 + delay] += 1 + if a != b or delay != 0: + expected[b, a, 2 - delay] += 1 + np.testing.assert_array_equal(actual, expected) + + +def test_cross_unit_central_pairs_do_not_depend_on_tied_row_order(): + samples = np.array([0, 0, 10, 10], dtype=np.int64) + windows = np.array([[0, 20]], dtype=np.int64) + a = window_correlograms(samples, np.array([0, 1, 1, 0]), windows, [0, 1], 30000, 0.001, 0.004) + b = window_correlograms(samples, np.array([1, 0, 1, 0]), windows, [0, 1], 30000, 0.001, 0.004) + np.testing.assert_array_equal(a, b) + assert a[0, 1, 2] == 4 diff --git a/phy/apps/large_recording/tests/test_summaries.py b/phy/apps/large_recording/tests/test_summaries.py new file mode 100644 index 000000000..5eb4d03e2 --- /dev/null +++ b/phy/apps/large_recording/tests/test_summaries.py @@ -0,0 +1,67 @@ +"""Full counts and acquisition exposure compared with an independent tiny oracle.""" + +import numpy as np +import pytest + +from phy.apps.large_recording.source import Source +from phy.apps.large_recording.summaries import build_summaries +from phy.apps.large_recording.tests.test_source import make_sorting + + +def test_exact_counts_exposure_and_additive_measurements(tmp_path): + params = make_sorting(tmp_path / 'full', columns=True) + amplitude = np.load(params.parent / 'amplitudes.npy') + amplitude[1] = np.nan + np.save(params.parent / 'amplitudes.npy', amplitude[:, None]) + source = Source(params, valid_intervals=[[0, 7], [10, 20]], batch_size=2) + files = build_summaries(source, tmp_path / 'summary', bin_seconds=5) + counts = np.load(tmp_path / 'summary' / files['counts_file']) + expected = np.zeros((4, 4), dtype=np.int64) + expected_sum = np.zeros((4, 4)) + expected_n = np.zeros((4, 4), dtype=np.int64) + units = np.load(params.parent / 'spike_clusters.npy').reshape(-1) + times = np.load(params.parent / 'spike_times.npy').reshape(-1) + for unit, time, value in zip(units, times, amplitude, strict=True): + row = list(source.unit_ids).index(unit) + column = int(time // 5) + expected[row, column] += 1 + if np.isfinite(value): + expected_sum[row, column] += value + expected_n[row, column] += 1 + np.testing.assert_array_equal(counts, expected) + np.testing.assert_array_equal(counts.sum(axis=1), source.counts) + np.testing.assert_array_equal( + np.load(tmp_path / 'summary' / files['edges_file']), [0, 5, 10, 15, 20] + ) + np.testing.assert_array_equal( + np.load(tmp_path / 'summary' / files['exposure_file']), [5, 2, 5, 5] + ) + np.testing.assert_array_equal( + np.load(tmp_path / 'summary' / files['amplitude_sum_file']), expected_sum + ) + np.testing.assert_array_equal( + np.load(tmp_path / 'summary' / files['amplitude_count_file']), expected_n + ) + assert files['measurement_units']['amplitude'] == 'source amplitude units' + assert 'depth_sum_file' in files + + +def test_summary_budget_checked_before_output_allocation(tmp_path): + source = Source(make_sorting(tmp_path / 'full')) + with pytest.raises(ValueError, match='Summary budget'): + build_summaries(source, tmp_path / 'summary', bin_seconds=1, max_bytes=50) + assert not (tmp_path / 'summary').exists() + + +def test_position_coordinates_require_declared_depth_axis(tmp_path): + from phy.apps.large_recording.summaries import scalar_attributes + + params = make_sorting(tmp_path / 'full') + (params.parent / 'spike_depths.npy').unlink() + with Source(params) as source: + assert 'depth' not in scalar_attributes(source) + params.write_text(params.read_text() + 'depth_axis = 1\n') + with Source(params) as source: + np.testing.assert_array_equal( + scalar_attributes(source)['depth'], source.arrays['spike_positions.npy'][:, 1] + ) diff --git a/phy/apps/large_recording/views.py b/phy/apps/large_recording/views.py new file mode 100644 index 000000000..c96fabf09 --- /dev/null +++ b/phy/apps/large_recording/views.py @@ -0,0 +1,45 @@ +"""Explicitly sampled statistics and window-relative amplitude views.""" + +from phy.cluster.views import AmplitudeView, CorrelogramView, ISIView +from phy.gui.qt import QLabel, Qt + + +class SampledISIView(ISIView): + ex_status = 'Selected windows only; intervals crossing window boundaries excluded.' + + +class SampledCorrelogramView(CorrelogramView): + ex_status = 'Selected windows only; raw pair counts, no full-recording rate baseline.' + + +class EvidenceAmplitudeView(AmplitudeView): + """Raw source amplitude in an explicitly rebased evidence window.""" + + window_caption = 'Recording time unavailable' + + def set_window_context(self, index, start_seconds, end_seconds): + self.window_caption = ( + f'Evidence window {index} | Recording time: {start_seconds:.3f}–{end_seconds:.3f} s' + ) + if hasattr(self, 'window_label'): + self.window_label.setText(self.window_caption) + + def attach(self, gui): + super().attach(gui) + layout = self.dock.widget().layout() + self.window_label = QLabel(self.window_caption) + self.window_label.setWordWrap(True) + self.time_axis_label = QLabel('Seconds from window start') + for label in (self.window_label, self.time_axis_label): + label.setAlignment(Qt.AlignCenter) + label.setStyleSheet('color: white; background-color: black; padding: 4px;') + layout.insertWidget(0, self.window_label) + layout.insertWidget(layout.count() - 1, self.time_axis_label) + + def _get_data_bounds(self, bunchs): + if not any(len(bunch.amplitudes) for bunch in bunchs): + return (0, 0, self.duration, 1) + xmin, ymin, xmax, ymax = super()._get_data_bounds(bunchs) + if ymin == ymax: + ymin, ymax = ymin - 0.5, ymax + 0.5 + return xmin, ymin, xmax, ymax diff --git a/phy/apps/template/gui.py b/phy/apps/template/gui.py index bda606c42..d8c1db3ba 100644 --- a/phy/apps/template/gui.py +++ b/phy/apps/template/gui.py @@ -81,6 +81,14 @@ class TemplateController(WaveformMixin, FeatureMixin, TemplateMixin, TraceMixin, # Internal methods # ------------------------------------------------------------------------- + def __new__(cls, dir_path=None, *args, **kwargs): + if cls is TemplateController and dir_path is not None: + if (Path(dir_path) / '_phy_curation_manifest.json').exists(): + from ..large_recording.adapter import BoundedTemplateController + + return object.__new__(BoundedTemplateController) + return object.__new__(cls) + def _get_waveforms_dict(self): waveforms_dict = super()._get_waveforms_dict() # Remove waveforms and mean_waveforms if there is no raw data file. @@ -142,6 +150,7 @@ def _get_template_features(self, cluster_ids, load_all=False): def _set_view_creator(self): super()._set_view_creator() self.view_creator['TemplateFeatureView'] = self.create_template_feature_view + self.view_creator['UnitStabilityView'] = self.create_unit_stability_view # Public methods # ------------------------------------------------------------------------- @@ -255,6 +264,12 @@ def create_template_feature_view(self): return return TemplateFeatureView(coords=self._get_template_features) + def create_unit_stability_view(self): + """Optionally inspect selected-unit activity over the loaded recording.""" + from ..large_recording.stability import UnitStabilityView + + return UnitStabilityView(self) + # ------------------------------------------------------------------------------ # Template commands @@ -268,9 +283,14 @@ def template_gui(params_path, **kwargs): # pragma: no cover dir_path = p.parent _add_log_file(dir_path / 'phy.log') - model = load_model(params_path) + bounded = (dir_path / '_phy_curation_manifest.json').exists() + model = None if bounded else load_model(params_path) # Automatically export spike waveforms when using compressed raw ephys. - if model.spike_waveforms is None and isinstance(model.traces, MtscompEphysReader): + if ( + model is not None + and model.spike_waveforms is None + and isinstance(model.traces, MtscompEphysReader) + ): # TODO: customizable values below. model.save_spikes_subset_waveforms(max_n_spikes_per_template=500, max_n_channels=16) @@ -285,6 +305,22 @@ def template_gui(params_path, **kwargs): # pragma: no cover def template_describe(params_path): """Describe a template dataset.""" + if (Path(params_path).parent / '_phy_curation_manifest.json').exists(): + from ..large_recording.curation import BoundedSession + + session = BoundedSession(Path(params_path).parent) + try: + print(f'Bounded evidence: {len(session.assignments):,} spikes') + print( + f'Complete source: {session.manifest["num_spikes"]:,} spikes, {len(session.unit_ids):,} units' + ) + print( + f'Duration: {session.manifest["num_samples"] / session.manifest["sample_rate"]:.3f} s' + ) + print(f'Committed curation: {len(session.events)} semantic operations') + finally: + session.close() + return model = load_model(params_path) model.describe() model.close() diff --git a/phy/apps/tests/large_recording_fixtures.py b/phy/apps/tests/large_recording_fixtures.py new file mode 100644 index 000000000..7041213bd --- /dev/null +++ b/phy/apps/tests/large_recording_fixtures.py @@ -0,0 +1,153 @@ +"""Small, deterministic on-disk example of a 1000-hour two-stage Phy dataset. + +This is reference test data, not a prepare-large implementation. The fake +recording has 32 spikes spread across 1000 hours, not 1000 hours of raw traces. +The full and bounded folders use ordinary Phy filenames. +""" + +import hashlib +import json +import shutil +from dataclasses import dataclass +from pathlib import Path + +import numpy as np + + +@dataclass(frozen=True) +class LongRecordingFixture: + full: Path + bounded: Path + global_ids: np.ndarray + upstream_waveform_ids: np.ndarray + bounded_waveform_ids: np.ndarray + + +def _sha256(path): + return hashlib.sha256(path.read_bytes()).hexdigest() + + +def _reference_remap(global_ids, upstream_ids): + """Reference intersection for sorted, unique, nonnegative int64 row IDs.""" + global_ids = np.asarray(global_ids) + upstream_ids = np.asarray(upstream_ids) + for ids in (global_ids, upstream_ids): + if ids.ndim != 1 or ids.dtype != np.dtype('int64'): + raise ValueError('Row IDs must be one-dimensional int64 arrays.') + if np.any(ids < 0) or np.any(ids[1:] <= ids[:-1]): + raise ValueError('Row IDs must be nonnegative, sorted, and unique.') + if not len(global_ids) or not len(upstream_ids): + empty = np.array([], dtype=np.int64) + return empty, empty.copy() + positions = np.searchsorted(global_ids, upstream_ids) + keep = positions < len(global_ids) + keep[keep] &= global_ids[positions[keep]] == upstream_ids[keep] + return positions[keep].astype(np.int64), np.flatnonzero(keep).astype(np.int64) + + +def make_long_recording_fixture(root): + """Write one full Phy dataset and its expected bounded counterpart.""" + root = Path(root) + full, bounded = root / 'full', root / 'bounded' + full.mkdir(parents=True) + bounded.mkdir(parents=True) + + fs = 30_000 + num_samples = 1000 * 3600 * fs # Far beyond signed int32 sample indices. + samples = np.linspace(100, num_samples - 100, 32, dtype=np.int64) + samples[8] = samples[5] # Simultaneous spikes retain distinct source row identities. + samples[6:8] = samples[5] + units = (np.arange(32) % 3).astype(np.int32) + amplitudes = np.arange(1, 33, dtype=np.float32) + templates = np.zeros((3, 12, 4), dtype=np.float32) + for unit in range(3): + templates[unit, :, unit] = -np.arange(12, dtype=np.float32) + + for name, array in ( + ('spike_times.npy', samples), + ('spike_clusters.npy', units), + ('spike_templates.npy', units), + ('amplitudes.npy', amplitudes), + ('templates.npy', templates), + ('channel_map.npy', np.arange(4, dtype=np.int32)), + ('channel_positions.npy', np.array([[0, 0], [20, 0], [0, 20], [20, 20]])), + ): + np.save(full / name, array) + (full / 'params.py').write_text( + "sample_rate = 30000.0\ndat_path = None\nn_channels_dat = 4\ndtype = 'int16'\n", + encoding='utf-8', + ) + + upstream_ids = np.array([1, 2, 8, 13, 18, 21, 28, 31], dtype=np.int64) + channels = np.array( + [[0, 1, -1] if i % 2 == 0 else [2, 3, -1] for i in range(len(upstream_ids))], + dtype=np.int32, + ) + waveforms = np.zeros((len(upstream_ids), 12, 3), dtype=np.float32) + for row, spike_id in enumerate(upstream_ids): + waveforms[row, :, 0] = spike_id + np.arange(12) / 12 + waveforms[row, :, 1] = -(spike_id + np.arange(12) / 12) + stem = '_phy_spikes_subset' + np.save(full / f'{stem}.spikes.npy', upstream_ids) + np.save(full / f'{stem}.channels.npy', channels) + np.save(full / f'{stem}.waveforms.npy', waveforms) + + # Proposed, optional upstream window metadata: simple homogeneous arrays. + starts = np.minimum(samples[upstream_ids] - 10, num_samples - 30 * fs) + bounds = np.column_stack((starts, starts + 30 * fs)).astype(np.int64) + types = np.zeros(len(upstream_ids), dtype=np.uint8) + types[-1] = 1 # rescue + strata = (samples[upstream_ids] // (3600 * fs)).astype(np.int32) + strata[-1] = -1 # rescue has no survey stratum + np.save(full / f'{stem}.window_bounds.npy', bounds) + np.save(full / f'{stem}.window_types.npy', types) + np.save(full / f'{stem}.window_strata.npy', strata) + source_hash = _sha256(full / 'spike_times.npy') + (full / f'{stem}.manifest.json').write_text( + json.dumps( + { + 'format': 'phy_spikes_subset', + 'format_version': 1, + 'source_sorting': {'num_spikes': len(samples), 'spike_times_sha256': source_hash}, + 'recording': { + 'sampling_frequency_hz': fs, + 'num_samples': num_samples, + 'num_channels': 4, + 'num_segments': 1, + }, + 'selection': {'algorithm': 'synthetic_windows_v1', 'seed': 42}, + 'waveform_extraction': {'nbefore_samples': 4, 'nafter_samples': 8}, + 'window_types': {'0': 'survey', '1': 'rescue'}, + }, + indent=2, + ) + + '\n', + encoding='utf-8', + ) + + global_ids = np.array([0, 2, 5, 8, 13, 21, 25, 31], dtype=np.int64) + for name in ('spike_times.npy', 'spike_clusters.npy', 'spike_templates.npy', 'amplitudes.npy'): + np.save(bounded / name, np.load(full / name)[global_ids]) + for name in ('templates.npy', 'channel_map.npy', 'channel_positions.npy', 'params.py'): + shutil.copyfile(full / name, bounded / name) + np.save(bounded / '_phy_spike_global_ids.npy', global_ids) + + bounded_rows, waveform_rows = _reference_remap(global_ids, upstream_ids) + np.save(bounded / f'{stem}.spikes.npy', bounded_rows) + np.save(bounded / f'{stem}.channels.npy', channels[waveform_rows]) + np.save(bounded / f'{stem}.waveforms.npy', waveforms[waveform_rows]) + (bounded / '_phy_curation_manifest.json').write_text( + json.dumps( + { + 'format': 'phy_bounded_curation', + 'format_version': 1, + 'source': {'num_spikes': len(samples), 'spike_times_sha256': source_hash}, + 'global_spike_ids_file': '_phy_spike_global_ids.npy', + 'duration_seconds': 1000 * 3600, + }, + indent=2, + ) + + '\n', + encoding='utf-8', + ) + return LongRecordingFixture(full, bounded, global_ids, upstream_ids, bounded_rows) diff --git a/phy/apps/tests/test_large_recording_contract.py b/phy/apps/tests/test_large_recording_contract.py new file mode 100644 index 000000000..8e882afc1 --- /dev/null +++ b/phy/apps/tests/test_large_recording_contract.py @@ -0,0 +1,151 @@ +"""Reference contract tests for future prepare-large and waveform remapping. + +These tests deliberately do not exercise a production prepare-large command: +that command does not exist yet. Later implementations should reuse this +fixture and replace reference-only checks with command-level assertions. +""" + +import hashlib +import json + +import numpy as np +import pytest +from phylib.io.model import load_model + +from .large_recording_fixtures import _reference_remap, make_long_recording_fixture + + +@pytest.fixture +def long_recording(tmp_path): + return make_long_recording_fixture(tmp_path) + + +def test_long_recording_uses_exact_64_bit_sample_times(long_recording): + full = np.load(long_recording.full / 'spike_times.npy') + bounded = np.load(long_recording.bounded / 'spike_times.npy') + assert full.dtype == bounded.dtype == np.dtype('int64') + assert full[-1] > np.iinfo(np.int32).max + assert full[-1] / 30_000 > 999 * 3600 + np.testing.assert_array_equal(bounded, full[long_recording.global_ids]) + assert np.all(np.diff(long_recording.global_ids) > 0) + assert np.all(np.diff(bounded) >= 0) + assert bounded[2] == bounded[3] # equal times do not collapse distinct spikes + + +def test_two_stage_row_indices_and_waveform_data(long_recording): + full = long_recording.full + bounded = long_recording.bounded + global_ids = np.load(bounded / '_phy_spike_global_ids.npy') + full_waveform_ids = np.load(full / '_phy_spikes_subset.spikes.npy') + bounded_waveform_ids = np.load(bounded / '_phy_spikes_subset.spikes.npy') + assert global_ids.dtype == bounded_waveform_ids.dtype == np.dtype('int64') + np.testing.assert_array_equal( + global_ids[bounded_waveform_ids], np.intersect1d(full_waveform_ids, global_ids) + ) + upstream_rows = np.searchsorted(full_waveform_ids, global_ids[bounded_waveform_ids]) + for suffix in ('channels', 'waveforms'): + np.testing.assert_array_equal( + np.load(bounded / f'_phy_spikes_subset.{suffix}.npy'), + np.load(full / f'_phy_spikes_subset.{suffix}.npy')[upstream_rows], + ) + # A bounded row identifies one source row, and never a waveform tensor row. + assert bounded_waveform_ids.tolist() == [1, 3, 4, 5, 7] + assert global_ids[bounded_waveform_ids].tolist() == [2, 8, 13, 21, 31] + + +def test_ordinary_phy_waveform_files_are_aligned(long_recording): + bounded = long_recording.bounded + spike_ids = np.load(bounded / '_phy_spikes_subset.spikes.npy') + channels = np.load(bounded / '_phy_spikes_subset.channels.npy') + waveforms = np.load(bounded / '_phy_spikes_subset.waveforms.npy', mmap_mode='r') + assert spike_ids.ndim == 1 and np.all(np.diff(spike_ids) > 0) + assert channels.dtype == np.dtype('int32') + assert channels.shape == (len(spike_ids), waveforms.shape[2]) + assert waveforms.shape == (len(spike_ids), 12, 3) + assert np.all((channels == -1) | ((channels >= 0) & (channels < 4))) + assert np.all(waveforms[:, :, -1] == 0) # -1 means a zero-padded channel + assert spike_ids[-1] < len(np.load(bounded / 'spike_times.npy')) + + +def test_window_arrays_use_samples_and_simple_dtypes(long_recording): + full = long_recording.full + bounds = np.load(full / '_phy_spikes_subset.window_bounds.npy') + types = np.load(full / '_phy_spikes_subset.window_types.npy') + strata = np.load(full / '_phy_spikes_subset.window_strata.npy') + spike_samples = np.load(full / 'spike_times.npy') + assert bounds.shape == (len(types), 2) + assert len(types) == len(strata) + assert bounds.dtype == np.dtype('int64') + assert types.dtype == np.dtype('uint8') + assert strata.dtype == np.dtype('int32') + assert np.all(bounds[:, 1] > bounds[:, 0]) + assert np.all(bounds[:, 1] <= 1000 * 3600 * 30_000) + assert np.all((types == 0) | (types == 1)) + assert strata[-1] == -1 and types[-1] == 1 + source_waveform_times = spike_samples[long_recording.upstream_waveform_ids] + assert np.all((bounds[:, 0] <= source_waveform_times) & (source_waveform_times < bounds[:, 1])) + + +def test_reference_remap_handles_empty_and_disjoint_inputs(): + a = np.array([14, 19, 117], dtype=np.int64) + b = np.array([2, 4, 119], dtype=np.int64) + for left, right in ((a, b), (a, b[:0]), (a[:0], b), (a[:0], b[:0])): + bounded_rows, source_rows = _reference_remap(left, right) + assert bounded_rows.dtype == source_rows.dtype == np.dtype('int64') + assert bounded_rows.shape == source_rows.shape == (0,) + rows, upstream_rows = _reference_remap(a, np.array([14, 15, 117, 118])) + assert rows.tolist() == [0, 2] + assert upstream_rows.tolist() == [0, 2] + + +@pytest.mark.parametrize('ids', [[2, 1], [1, 1], [-1, 1]]) +def test_reference_remap_rejects_invalid_row_order(ids): + valid = np.array([1, 2], dtype=np.int64) + invalid = np.array(ids, dtype=np.int64) + for left, right in ((valid, invalid), (invalid, valid)): + with pytest.raises(ValueError, match='nonnegative, sorted, and unique'): + _reference_remap(left, right) + + +def test_manifests_bind_source_index_space_not_mutable_clusters(long_recording): + full = long_recording.full + bounded = long_recording.bounded + source_hash = hashlib.sha256((full / 'spike_times.npy').read_bytes()).hexdigest() + upstream = json.loads((full / '_phy_spikes_subset.manifest.json').read_text()) + downstream = json.loads((bounded / '_phy_curation_manifest.json').read_text()) + assert upstream['source_sorting']['spike_times_sha256'] == source_hash + assert downstream['source']['spike_times_sha256'] == source_hash + assert upstream['source_sorting']['num_spikes'] == len(np.load(full / 'spike_times.npy')) + assert downstream['global_spike_ids_file'] == '_phy_spike_global_ids.npy' + assert downstream['duration_seconds'] == 1000 * 3600 + assert hashlib.sha256((bounded / 'spike_times.npy').read_bytes()).hexdigest() != source_hash + + +def test_bounded_dataset_loads_with_existing_phylib(long_recording): + """The old Phy loader consumes bounded rows, not original global indices.""" + model = load_model(long_recording.bounded / 'params.py') + try: + np.testing.assert_array_equal( + model.spike_waveforms.spike_ids, long_recording.bounded_waveform_ids + ) + assert model.n_spikes == len(long_recording.global_ids) + assert model.n_samples_waveforms == 12 + assert model.spike_times[-1] > 999 * 3600 + np.testing.assert_array_equal( + model.get_waveforms(long_recording.bounded_waveform_ids, channel_ids=[0, 1, 2, 3]), + _expected_dense_waveforms(long_recording.bounded), + ) + finally: + model.close() + + +def _expected_dense_waveforms(path): + """Independent dense expansion of the fixture's sparse channel slots.""" + channels = np.load(path / '_phy_spikes_subset.channels.npy') + waveforms = np.load(path / '_phy_spikes_subset.waveforms.npy') + expected = np.zeros((len(channels), 12, 4), dtype=waveforms.dtype) + for row, channel_row in enumerate(channels): + for slot, channel in enumerate(channel_row): + if channel >= 0: + expected[row, :, channel] = waveforms[row, :, slot] + return expected diff --git a/phy/cluster/clustering.py b/phy/cluster/clustering.py index 662ca7439..32fa7480e 100644 --- a/phy/cluster/clustering.py +++ b/phy/cluster/clustering.py @@ -7,7 +7,8 @@ import logging import numpy as np -from phylib.io.array import _spikes_in_clusters, _spikes_per_cluster, _unique +from phylib.io.array import _spikes_in_clusters, _spikes_per_cluster +from phylib.io.array import _unique as _unique_dense from phylib.utils._types import _as_array, _is_array_like from phylib.utils.event import emit @@ -22,6 +23,23 @@ # ------------------------------------------------------------------------------ +def _unique(values): + """Preserve dense-ID performance without allocating by sparse absolute IDs.""" + if values is None or len(values) == 0: + return np.array([], dtype=np.int64) + values = np.asarray(values) + if int(values.max()) > max(1024, 4 * len(values)): + return np.unique(values[values >= 0]).astype(np.int64) + return _unique_dense(values) + + +def _largest_cluster(values): + if int(values.max()) > max(1024, 4 * len(values)): + ids, counts = np.unique(values, return_counts=True) + return int(ids[counts.argmax()]) + return int(np.bincount(values).argmax()) + + def _extend_spikes(spike_ids, spike_clusters): """Return all spikes belonging to the clusters containing the specified spikes.""" @@ -81,7 +99,7 @@ def _extend_assignment(spike_ids, old_spike_clusters, spike_clusters_rel, new_cl def _assign_update_info(spike_ids, old_spike_clusters, new_spike_clusters): old_clusters = _unique(old_spike_clusters) new_clusters = _unique(new_spike_clusters) - largest_old_cluster = np.bincount(old_spike_clusters).argmax() + largest_old_cluster = _largest_cluster(old_spike_clusters) descendants = list(set(zip(old_spike_clusters, new_spike_clusters))) update_info = UpdateInfo( description='assign', @@ -293,7 +311,9 @@ def _do_assign(self, spike_ids, new_spike_clusters): def _do_merge(self, spike_ids, cluster_ids, to): # Create the UpdateInfo instance here. descendants = [(cluster, to) for cluster in cluster_ids] - largest_old_cluster = np.bincount(self.spike_clusters[spike_ids]).argmax() + largest_old_cluster = max( + cluster_ids, key=lambda cluster: (len(self._spikes_per_cluster[cluster]), -cluster) + ) up = UpdateInfo( description='merge', spike_ids=list(spike_ids), diff --git a/phy/cluster/supervisor.py b/phy/cluster/supervisor.py index 54c5988d1..5b0b4a3f9 100644 --- a/phy/cluster/supervisor.py +++ b/phy/cluster/supervisor.py @@ -755,6 +755,11 @@ class Supervisor: skip_masked_clusters : bool Whether automatic navigation and similar-cluster selection skip noise and MUA clusters. The default is True. + clustering_factory : callable, optional + Construct the clustering model. Defaults to Clustering; bounded curation + supplies a guarded implementation with the same constructor interface. + cluster_meta_factory : callable, optional + Construct metadata from groups and labels. Defaults to create_cluster_meta. Events ------ @@ -790,6 +795,8 @@ def __init__( n_similar_clusters_to_select=None, skip_masked_clusters=True, merge_propositions=None, + clustering_factory=None, + cluster_meta_factory=None, ): super().__init__() self.context = context @@ -833,7 +840,7 @@ def __init__( # Create Clustering and ClusterMeta. # Load the cached spikes_per_cluster array. spc = context.load('spikes_per_cluster') if context else None - self.clustering = Clustering( + self.clustering = (clustering_factory or Clustering)( spike_clusters, spikes_per_cluster=spc, new_cluster_id=new_cluster_id ) @@ -841,7 +848,9 @@ def __init__( self._save_spikes_per_cluster() # Create the ClusterMeta instance. - self.cluster_meta = create_cluster_meta(cluster_groups or {}, self.cluster_labels) + self.cluster_meta = (cluster_meta_factory or create_cluster_meta)( + cluster_groups or {}, self.cluster_labels + ) # Create the GlobalHistory instance. self._global_history = GlobalHistory( diff --git a/phy/cluster/tests/test_clustering.py b/phy/cluster/tests/test_clustering.py index 9557221d7..ef8d03beb 100644 --- a/phy/cluster/tests/test_clustering.py +++ b/phy/cluster/tests/test_clustering.py @@ -26,6 +26,28 @@ # ------------------------------------------------------------------------------ +def test_sparse_absolute_cluster_ids_do_not_allocate_dense_histograms(monkeypatch): + original = np.bincount + + def bounded_bincount(values, *args, **kwargs): + assert not len(values) or np.max(values) < 1024 + return original(values, *args, **kwargs) + + monkeypatch.setattr(np, 'bincount', bounded_bincount) + high = 2_000_000_000 + clustering = Clustering(np.array([0, 1, 1, high], dtype=np.int64)) + first = clustering.merge([0, 1]) + assert first.largest_old_cluster == 1 + assert first.added == [high + 1] + second = clustering.merge([high, high + 1]) + assert second.largest_old_cluster == high + 1 + clustering.undo() + ae(clustering.spike_clusters, [high + 1, high + 1, high + 1, high]) + clustering.split([0]) + clustering.undo() + ae(clustering.spike_clusters, [high + 1, high + 1, high + 1, high]) + + def test_extend_spikes_simple(): spike_clusters = np.array([3, 5, 2, 9, 5, 5, 2]) spike_ids = np.array([2, 4, 0]) diff --git a/tools/benchmark_long_recording.py b/tools/benchmark_long_recording.py new file mode 100644 index 000000000..b6bfabeef --- /dev/null +++ b/tools/benchmark_long_recording.py @@ -0,0 +1,398 @@ +"""Reproducible, file-backed capacity benchmark for long-recording services. + +Run from a development environment, for example:: + + uv run python tools/benchmark_long_recording.py --spikes 1000000 + +The benchmark creates synthetic Phy arrays without raw traces, times preparation +and replay in separate child processes, and checks replay against an independent +streamed oracle. It is intentionally capped at two million source spikes. +""" + +import argparse +import json +import os +import platform +import subprocess +import sys +import tempfile +import threading +import time +from pathlib import Path + +import numpy as np + +REPOSITORY = Path(__file__).resolve().parents[1] +if str(REPOSITORY) not in sys.path: + sys.path.insert(0, str(REPOSITORY)) + + +def _file_bytes(path): + path = Path(path) + return sum(item.stat().st_size for item in path.rglob('*') if item.is_file()) + + +def _make_source(path, spikes, units, duration_hours, batch_size): + """Write canonical arrays in bounded chunks, without allocating full vectors.""" + path.mkdir(parents=True) + duration_seconds = duration_hours * 3600 + sample_rate = 30_000 + num_samples = round(duration_seconds * sample_rate) + if num_samples <= 0 or num_samples > np.iinfo(np.int64).max: + raise ValueError('duration-hours produces an invalid int64 sample axis.') + + arrays = { + 'spike_times.npy': np.lib.format.open_memmap( + path / 'spike_times.npy', mode='w+', dtype=np.int64, shape=(spikes,) + ), + 'spike_clusters.npy': np.lib.format.open_memmap( + path / 'spike_clusters.npy', mode='w+', dtype=np.int32, shape=(spikes,) + ), + 'spike_templates.npy': np.lib.format.open_memmap( + path / 'spike_templates.npy', mode='w+', dtype=np.int32, shape=(spikes,) + ), + } + for start in range(0, spikes, batch_size): + stop = min(start + batch_size, spikes) + rows = np.arange(start, stop, dtype=np.int64) + arrays['spike_times.npy'][start:stop] = rows * (num_samples - 1) // spikes + assignments = (rows % units).astype(np.int32) + arrays['spike_clusters.npy'][start:stop] = assignments + arrays['spike_templates.npy'][start:stop] = assignments + for array in arrays.values(): + array.flush() + array._mmap.close() + + templates = np.lib.format.open_memmap( + path / 'templates.npy', mode='w+', dtype=np.float32, shape=(units, 4, 2) + ) + templates[:] = 0 + templates.flush() + templates._mmap.close() + np.save(path / 'channel_map.npy', np.array([0, 1], dtype=np.int32)) + np.save(path / 'channel_positions.npy', np.array([[0.0, 0.0], [20.0, 0.0]], dtype=np.float32)) + (path / 'params.py').write_text( + f'sample_rate = {sample_rate}.0\n' + f'num_samples = {num_samples}\n' + 'n_channels_dat = 2\n' + "dtype = 'int16'\n" + 'dat_path = None\n', + encoding='utf-8', + ) + return duration_seconds + + +def _peak_rss_bytes(): + """Return process peak RSS, or an explicit unavailable result.""" + if os.name != 'nt': + try: + import resource + + value = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + # Linux and the BSDs report KiB; macOS reports bytes. + return int(value if sys.platform == 'darwin' else value * 1024), 'resource.ru_maxrss' + except (ImportError, AttributeError, OSError): + pass + + try: + import psutil + except ImportError: + return None, 'unavailable (install psutil for sampled peak RSS on Windows)' + + process = psutil.Process() + peak = [process.memory_info().rss] + stopped = threading.Event() + + def sample(): + while not stopped.wait(0.02): + try: + peak[0] = max(peak[0], process.memory_info().rss) + except psutil.Error: + return + + sampler = threading.Thread(target=sample, daemon=True) + sampler.start() + return _SampledRSS(stopped, sampler, peak) + + +class _SampledRSS: + def __init__(self, stopped, sampler, peak): + self.stopped = stopped + self.sampler = sampler + self.peak = peak + + def finish(self): + self.stopped.set() + self.sampler.join() + return int(self.peak[0]), 'psutil sampled every 20 ms' + + +def _measure_rss_start(): + value = _peak_rss_bytes() + return value if isinstance(value, _SampledRSS) else value + + +def _rss_result(handle): + return handle.finish() if isinstance(handle, _SampledRSS) else handle + + +def _phase_prepare(args): + from phy.apps.large_recording.prepare import prepare_large + + rss = _measure_rss_start() + started = time.perf_counter() + bounded = prepare_large( + args.params, + args.bounded_path, + duration_seconds=args.duration_seconds, + max_spikes=args.max_spikes, + max_bytes=args.max_bytes, + window_seconds=args.window_seconds, + stratum_seconds=args.stratum_seconds, + seed=args.seed, + batch_size=args.batch_size, + bin_seconds=args.bin_seconds, + ) + elapsed = time.perf_counter() - started + peak, method = _rss_result(rss) + print( + json.dumps( + { + 'phase': 'prepare', + 'wall_seconds': elapsed, + 'peak_rss_bytes': peak, + 'peak_rss_method': method, + 'bounded_path': str(bounded), + } + ) + ) + + +def _phase_curate(args): + from phy.apps.large_recording.curation import BoundedSession + + session = BoundedSession(args.bounded_path) + try: + if 0 not in session.unit_ids or 1 not in session.unit_ids: + raise ValueError('Synthetic source must contain units 0 and 1.') + session.record_merge([0, 1], session.next_cluster_id) + session.save() + finally: + session.close() + print(json.dumps({'phase': 'curate', 'saved': True})) + + +def _phase_replay(args): + from phy.apps.large_recording.replay import apply_curation + + rss = _measure_rss_start() + started = time.perf_counter() + output = apply_curation( + args.params, args.bounded_path, args.curated_path, batch_size=args.batch_size + ) + elapsed = time.perf_counter() - started + peak, method = _rss_result(rss) + print( + json.dumps( + { + 'phase': 'replay', + 'wall_seconds': elapsed, + 'peak_rss_bytes': peak, + 'peak_rss_method': method, + 'curated_path': str(output), + } + ) + ) + + +def _oracle(source_path, curated_path, spikes, units, batch_size): + """Independent streamed oracle for the benchmark's one known merge.""" + source = np.load(source_path / 'spike_clusters.npy', mmap_mode='r', allow_pickle=False) + output = np.load(curated_path / 'spike_clusters.npy', mmap_mode='r', allow_pickle=False) + expected_target = units + checked = 0 + matches = source.shape == output.shape == (spikes,) + for start in range(0, spikes, batch_size): + stop = min(start + batch_size, spikes) + expected = np.array(source[start:stop], copy=True) + expected[(expected == 0) | (expected == 1)] = expected_target + matches &= np.array_equal(expected, output[start:stop]) + checked += stop - start + source._mmap.close() + output._mmap.close() + return {'matches': bool(matches), 'rows_checked': checked, 'merge_target': expected_target} + + +def _run_phase( + script, + phase, + params, + bounded, + curated, + duration_seconds, + max_spikes, + max_bytes, + window_seconds, + stratum_seconds, + seed, + batch_size, + bin_seconds, +): + command = [ + sys.executable, + str(script), + '--_phase', + phase, + '--params', + str(params), + '--bounded-path', + str(bounded), + '--curated-path', + str(curated), + '--duration-seconds', + str(duration_seconds), + '--max-spikes', + str(max_spikes), + '--max-bytes', + str(max_bytes), + '--window-seconds', + str(window_seconds), + '--stratum-seconds', + str(stratum_seconds), + '--seed', + str(seed), + '--batch-size', + str(batch_size), + '--bin-seconds', + str(bin_seconds), + ] + try: + completed = subprocess.run( + command, cwd=REPOSITORY, check=True, capture_output=True, text=True + ) + except subprocess.CalledProcessError as error: + raise RuntimeError(error.stderr.strip() or error.stdout.strip()) from error + return json.loads(completed.stdout.strip().splitlines()[-1]) + + +def _run(args): + if args.spikes <= 0 or args.spikes > 2_000_000: + raise ValueError('--spikes must be between 1 and 2,000,000.') + if args.units < 2 or args.units > args.spikes: + raise ValueError('--units must be at least 2 and no greater than --spikes.') + if args.batch_size <= 0 or args.duration_hours <= 0: + raise ValueError('--batch-size and --duration-hours must be positive.') + + temporary = None + if args.output_dir: + root = Path(args.output_dir).resolve() + root.mkdir(parents=True, exist_ok=True) + else: + temporary = tempfile.TemporaryDirectory(prefix='phy-long-recording-benchmark-') + root = Path(temporary.name) + source_dir, bounded_dir, curated_dir = ( + root / name for name in ('source', 'bounded', 'curated') + ) + duration_seconds = _make_source( + source_dir, args.spikes, args.units, args.duration_hours, args.batch_size + ) + source_bytes = _file_bytes(source_dir) + common = { + 'params': source_dir / 'params.py', + 'bounded': bounded_dir, + 'curated': curated_dir, + 'duration_seconds': duration_seconds, + 'max_spikes': args.max_spikes, + 'max_bytes': args.max_bytes, + 'window_seconds': args.window_seconds, + 'stratum_seconds': args.stratum_seconds, + 'seed': args.seed, + 'batch_size': args.batch_size, + 'bin_seconds': args.bin_seconds, + } + prepare = _run_phase(Path(__file__).resolve(), 'prepare', **common) + _run_phase(Path(__file__).resolve(), 'curate', **common) + replay = _run_phase(Path(__file__).resolve(), 'replay', **common) + oracle = _oracle(source_dir, curated_dir, args.spikes, args.units, args.batch_size) + if not oracle['matches']: + raise RuntimeError('Independent streamed replay oracle did not match curated assignments.') + + result = { + 'benchmark': 'phy_long_recording_v1', + 'machine': { + 'platform': platform.platform(), + 'machine': platform.machine(), + 'processor': platform.processor(), + 'python': platform.python_version(), + 'numpy': np.__version__, + 'cpu_count': os.cpu_count(), + }, + 'parameters': { + 'spikes': args.spikes, + 'units': args.units, + 'duration_hours': args.duration_hours, + 'sample_rate_hz': 30000, + 'batch_size': args.batch_size, + 'max_spikes': args.max_spikes, + 'max_bytes': args.max_bytes, + 'window_seconds': args.window_seconds, + 'stratum_seconds': args.stratum_seconds, + 'seed': args.seed, + 'bin_seconds': args.bin_seconds, + }, + 'source_file_bytes': source_bytes, + 'bounded_file_bytes': _file_bytes(bounded_dir), + 'curated_file_bytes': _file_bytes(curated_dir), + 'bounded_spikes': int( + json.loads((bounded_dir / '_phy_curation_manifest.json').read_text(encoding='utf-8'))[ + 'bounded_num_spikes' + ] + ), + 'preparation': prepare, + 'replay': replay, + 'independent_streamed_replay_oracle': oracle, + 'output_dir': str(root), + 'output_dir_retained': temporary is None, + } + print(json.dumps(result, indent=2, sort_keys=True)) + if temporary is not None: + temporary.cleanup() + + +def _parser(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--spikes', type=int, default=1_000_000) + parser.add_argument('--units', type=int, default=1000) + parser.add_argument('--duration-hours', type=float, default=1000) + parser.add_argument('--batch-size', type=int, default=100_000) + parser.add_argument('--max-spikes', type=int, default=100_000) + parser.add_argument('--max-bytes', type=int, default=536_870_912) + parser.add_argument('--window-seconds', type=float, default=30) + parser.add_argument('--stratum-seconds', type=float, default=3600) + parser.add_argument('--seed', type=int, default=0) + parser.add_argument('--bin-seconds', type=float, default=3600) + parser.add_argument('--output-dir', type=Path) + parser.add_argument( + '--_phase', choices=('prepare', 'curate', 'replay'), help=argparse.SUPPRESS + ) + parser.add_argument('--params', type=Path, help=argparse.SUPPRESS) + parser.add_argument('--bounded-path', type=Path, help=argparse.SUPPRESS) + parser.add_argument('--curated-path', type=Path, help=argparse.SUPPRESS) + parser.add_argument('--duration-seconds', type=float, help=argparse.SUPPRESS) + return parser + + +def main(): + args = _parser().parse_args() + if args._phase == 'prepare': + _phase_prepare(args) + elif args._phase == 'curate': + _phase_curate(args) + elif args._phase == 'replay': + _phase_replay(args) + else: + _run(args) + + +if __name__ == '__main__': + main() diff --git a/tools/capture_long_recording_views.py b/tools/capture_long_recording_views.py new file mode 100644 index 000000000..1c40089c8 --- /dev/null +++ b/tools/capture_long_recording_views.py @@ -0,0 +1,254 @@ +"""Capture reproducible screenshots of the bounded-evidence Qt views. + +The script builds a temporary synthetic Phy sorting, prepares it with the real +long-recording workflow, opens the real GUI, and saves two screenshots. It does +not download or modify a recording. + +Run with a display, or under a virtual X server:: + + QT_QPA_PLATFORM=xcb LIBGL_ALWAYS_SOFTWARE=1 \ + xvfb-run -a -s '-screen 0 2400x1600x24 +extension GLX' \ + uv run python tools/capture_long_recording_views.py +""" + +import argparse +import sys +import tempfile +from pathlib import Path + +import numpy as np +from PyQt5.QtTest import QTest + +from phy.apps.large_recording.prepare import prepare_large +from phy.apps.template.gui import TemplateController +from phy.gui.qt import QApplication, Qt + +REPOSITORY = Path(__file__).resolve().parents[1] +DEFAULT_OUTPUT = REPOSITORY / 'docs' / 'images' +SAMPLE_RATE = 30_000 +HOURS = 8 +CLUSTER_IDS = np.array([0, 2, 7, 99], dtype=np.int32) + + +def _make_source(path): + """Write a varied but deterministic synthetic source sorting.""" + path.mkdir(parents=True) + rng = np.random.default_rng(20261013) + samples = SAMPLE_RATE * 3600 * HOURS + intervals = np.array( + [[0, 3 * 3600 * SAMPLE_RATE], [int(3.5 * 3600 * SAMPLE_RATE), samples]], + dtype=np.int64, + ) + valid = intervals[:, 1] - intervals[:, 0] + cumulative = np.cumsum(valid) + offsets = rng.integers(0, cumulative[-1], size=120_000, dtype=np.int64) + interval_index = np.searchsorted(cumulative, offsets, side='right') + previous = np.where(interval_index == 0, 0, cumulative[np.maximum(interval_index - 1, 0)]) + spike_times = intervals[interval_index, 0] + offsets - previous + order = np.argsort(spike_times, kind='stable') + spike_times = spike_times[order] + + unit_index = np.arange(len(spike_times), dtype=np.int64) % len(CLUSTER_IDS) + unit_index = unit_index[order] + elapsed = spike_times / SAMPLE_RATE + hour_phase = elapsed / 3600 + amplitudes = ( + 35 + + unit_index * 9 + + 5 * np.sin(hour_phase * np.pi / 3) + + 2 * hour_phase + + rng.normal(0, 2.5, len(spike_times)) + ).astype(np.float32) + depths = ( + 250 + + unit_index * 180 + + 20 * np.sin(hour_phase * np.pi / 2) + + rng.normal(0, 8, len(spike_times)) + ).astype(np.float32) + templates_for_spike = unit_index.astype(np.int32) + + np.save(path / 'spike_times.npy', spike_times.astype(np.int64)) + np.save(path / 'spike_clusters.npy', CLUSTER_IDS[unit_index]) + np.save(path / 'spike_templates.npy', templates_for_spike) + np.save(path / 'amplitudes.npy', amplitudes) + np.save(path / 'spike_depths.npy', depths) + + channel_map = np.arange(4, dtype=np.int32) + channel_positions = np.array([[0, 0], [20, 100], [0, 200], [20, 300]], dtype=np.float32) + np.save(path / 'channel_map.npy', channel_map) + np.save(path / 'channel_positions.npy', channel_positions) + + templates = np.zeros((len(CLUSTER_IDS), 40, len(channel_map)), dtype=np.float32) + time = np.arange(40, dtype=np.float32) + for unit in range(len(CLUSTER_IDS)): + center = 13 + unit * 3 + pulse = -np.exp(-0.5 * ((time - center) / 3.5) ** 2) * (35 + 8 * unit) + for channel in range(len(channel_map)): + templates[unit, :, channel] = pulse * np.exp(-0.7 * abs(channel - unit)) + np.save(path / 'templates.npy', templates) + + pc_feature_ind = np.array( + [ + [unit, (unit + 1) % len(channel_map), (unit + 2) % len(channel_map)] + for unit in range(len(CLUSTER_IDS)) + ], + dtype=np.int32, + ) + np.save(path / 'pc_feature_ind.npy', pc_feature_ind) + features = np.empty((len(spike_times), 3, 3), dtype=np.float32) + features[:, 0, :] = unit_index[:, None] * 3.5 + rng.normal(0, 0.65, (len(spike_times), 3)) + features[:, 1, :] = np.sin(hour_phase[:, None] * np.pi / 2) + rng.normal( + 0, 0.7, (len(spike_times), 3) + ) + features[:, 2, :] = np.cos(hour_phase[:, None] * np.pi / 3) + rng.normal( + 0, 0.6, (len(spike_times), 3) + ) + np.save(path / 'pc_features.npy', features) + + # Stored waveforms are a display reservoir, aligned to sorted source rows. + waveform_ids = np.arange(0, len(spike_times), 8, dtype=np.int64) + waveform_channels = np.tile( + np.arange(len(channel_map), dtype=np.int32), (len(waveform_ids), 1) + ) + waveforms = templates[unit_index[waveform_ids]].copy() + waveforms += rng.normal(0, 2.2, waveforms.shape).astype(np.float32) + stem = '_phy_spikes_subset' + np.save(path / f'{stem}.spikes.npy', waveform_ids) + np.save(path / f'{stem}.channels.npy', waveform_channels) + np.save(path / f'{stem}.waveforms.npy', waveforms) + + with (path / 'cluster_group.tsv').open('w', encoding='utf-8') as stream: + stream.write('cluster_id\tgroup\n0\tgood\n2\tgood\n7\tmua\n99\tgood\n') + (path / 'params.py').write_text( + f'sample_rate = {SAMPLE_RATE}.0\n' + f'num_samples = {samples}\n' + 'n_channels_dat = 4\n' + "dtype = 'int16'\n" + 'dat_path = None\n' + f'valid_intervals = {intervals.tolist()!r}\n', + encoding='utf-8', + ) + return path / 'params.py', intervals + + +def _wait_for_render(): + QApplication.processEvents() + QTest.qWait(1200) + QApplication.processEvents() + + +def _capture(gui, destination): + _wait_for_render() + # Phy's OpenGL canvases are native child windows. QWidget.grab() omits + # those children on some Qt/X11 combinations, so grab the rendered screen + # and crop to the application window instead. + screen = QApplication.primaryScreen() + frame = gui.frameGeometry() + pixmap = screen.grabWindow(0).copy(frame.x(), frame.y(), frame.width(), frame.height()) + if not pixmap.save(str(destination), 'PNG'): + raise OSError(f'Could not save screenshot: {destination}') + + +def _capture_views(bounded, config_dir, output_dir, stability): + views = ( + ['UnitStabilityView'] + if stability + else [ + 'EvidenceAmplitudeView', + 'FeatureView', + 'SampledISIView', + 'SampledCorrelogramView', + ] + ) + controller = TemplateController( + dir_path=bounded, + config_dir=config_dir, + plugins=[], + enable_threading=False, + ) + gui = controller.create_gui(default_views=views, do_prompt_save=False) + try: + gui.resize(2100, 1200) + gui.show() + _wait_for_render() + # Focus the captures on the plotted evidence rather than the GUI tables. + similarity = gui.get_view('SimilarityView') + if similarity is not None: + similarity.dock.hide() + cluster_table = gui.get_view('ClusterView') + if cluster_table is not None: + cluster_table.dock.hide() + if not stability: + amplitude, feature, isi, correlogram = ( + gui.get_view(name) + for name in ( + 'EvidenceAmplitudeView', + 'FeatureView', + 'SampledISIView', + 'SampledCorrelogramView', + ) + ) + gui.splitDockWidget(amplitude.dock, feature.dock, Qt.Vertical) + gui.splitDockWidget(amplitude.dock, isi.dock, Qt.Horizontal) + gui.splitDockWidget(isi.dock, correlogram.dock, Qt.Vertical) + gui.resizeDocks([amplitude.dock, isi.dock], [1000, 1000], Qt.Horizontal) + gui.resizeDocks([amplitude.dock, feature.dock], [550, 550], Qt.Vertical) + gui.resizeDocks([isi.dock, correlogram.dock], [550, 550], Qt.Vertical) + _wait_for_render() + windows = controller.evidence_windows + times = controller.model.spike_samples + counts = [ + int(np.searchsorted(times, end) - np.searchsorted(times, start)) + for start, end in windows + ] + best_window = int(np.argmax(counts)) + controller.inspect_evidence_window(best_window) + selected = [0] if stability else [0, 2] + controller.supervisor.select(selected) + controller.supervisor.block() + _wait_for_render() + if stability: + filename = 'long-recording-overview.png' + else: + filename = 'long-recording-evidence.png' + destination = output_dir / filename + _capture(gui, destination) + print( + f'{destination} ({gui.width()}x{gui.height()}, window {best_window}, ' + f'{counts[best_window]} evidence spikes)' + ) + finally: + gui.close() + controller.close() + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--output-dir', type=Path, default=DEFAULT_OUTPUT) + args = parser.parse_args() + args.output_dir.mkdir(parents=True, exist_ok=True) + app = QApplication.instance() or QApplication(sys.argv[:1]) + + with tempfile.TemporaryDirectory(prefix='phy-long-recording-screenshots-') as temporary: + root = Path(temporary) + params, intervals = _make_source(root / 'source') + bounded = prepare_large( + params, + root / 'bounded', + duration_seconds=SAMPLE_RATE * 3600 * HOURS / SAMPLE_RATE, + valid_intervals=intervals, + max_spikes=20_000, + max_bytes=268_435_456, + window_seconds=30, + stratum_seconds=3600, + seed=0, + batch_size=50_000, + bin_seconds=1800, + ) + _capture_views(bounded, root / 'config-stability', args.output_dir, stability=True) + _capture_views(bounded, root / 'config-evidence', args.output_dir, stability=False) + app.quit() + + +if __name__ == '__main__': + main()