Add adata convert: matrix dtype, layout and density on disk - #15
Conversation
Closes #13, which asked to change the dtype of X/data and X/indices on disk because count matrices stored as float64 waste disk and read time. Generalised to any matrix, plus the CSR<->CSC conversion that concat has been telling people to do since it was written, with nothing to do it. adata convert data.h5ad X -o out.h5ad --dtype float32 adata convert data.h5ad X -o out.h5ad --layout csc adata convert data.h5ad --all -o out.h5ad --dtype float32 Two refusals, both before the destination store exists. A cast is checked by casting every value and casting it back, not by comparing dtypes: whether float64 counts survive float32 depends on the counts, which is the actual question #13 is asking, and the answer is worth printing even when it is yes. A densification is checked against the projected size and refused above 4x growth. --force overrides either and both messages say what was measured. Correctness is held to scipy: every path is compared against what scipy would have produced from the same matrix -- indptr, indices and data each compared exactly, not approximately -- for both transpose directions, both implementations, and across bucket counts from one to hundreds. CSR->CSC->CSR is asserted to be the identity. Transposing streams by default. It cannot be done in one pass, so: count nonzeros per output major to build indptr, scatter the input into buckets by where it lands, then sort and write each bucket in order. Two extra passes over nnz buys a peak set by the chunk rather than by the matrix. --in-memory does it with one lexsort for when the matrix fits. Three things this got wrong on the way, all caught by measuring rather than by reading the code: The first version made files eight times *larger* on a conversion asked for to make them smaller. _growable creates with a fixed 65,536-element chunk, so a 2,400-nonzero output allocated 786 KB of mostly empty chunk; and the index dtype defaulted to int64, silently widening every int32 store that passed through. Outputs are now sized up front -- nnz is known, and is invariant under a transpose -- with chunking and compression forwarded from the source, and the index dtype is preserved unless asked. A test asserts the file shrinks, because that is the feature. The docstring claimed nothing is written before the checks and that was false: the destination was created, and obs and var copied into it, before convert_matrix ran its first check. plan_conversion now resolves and checks every target before the output store is opened, and two tests assert no output file exists after a refusal. The streaming transpose held the whole matrix below 4M nonzeros, because the bucket size was a constant rather than derived from the chunk, so everything smaller went into a single bucket. The performance guard caught it: peak grew 15.9x for 32x the nonzeros on a path whose whole justification is that it does not. Buckets now follow `chunk`. concat's check is unchanged in behaviour but now names the command to run, and is tested for the first time -- CSR+CSC, CSR+dense, CSC+dense, X and layers, including a test that follows the suggested fix and concatenates successfully afterwards. The unreachable duplicate raise is kept and labelled a backstop. Guards: reads linear in nnz for a cast and in n_obs for a layout change; streaming transpose peak bounded while the in-memory one grows, compared in the same run; and an assertion that the default really is the streaming path, since the two produce identical output and nothing else can tell them apart. Benchmark cases convert-dtype and convert-layout, the latter running both implementations so the trade is visible. 952 tests pass, coverage 92.65%, 204 compatibility tests across six anndata releases. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The command list and the worked examples both stopped at concat. The new section is its own, not folded into "Filtering without a name list", because shrinking a store and making encodings agree is not filtering. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Codex Review SummaryThis comment shows the latest Codex review activity on this pull request.
ℹ️ About Codex in GitHubYour team has set up Codex to review pull requests in this repo. Reviews are triggered when you
Codex reacts with 👀 while any review is running, comments if it has suggestions, and reacts with 👍 once all reviews finish with no findings. |
![chatgpt-codex-connector[bot]](https://avatars.githubusercontent.com/in/1144995?s=60&v=4){kind=small}
There was a problem hiding this comment.
💡 Codex Review
Here are some automated review suggestions for this pull request.
Reviewed commit: 066fec844c
ℹ️ About Codex in GitHub
Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you
- Open a pull request for review
- Mark a draft as ready
- Comment "@codex review".
If Codex has suggestions, it will comment; otherwise it will react with 👍.
Codex can also answer questions or update the PR. Try commenting "@codex address that feedback".
| + ", ".join(targets) | ||
| ) | ||
|
|
||
| with open_store(dst_path, "w", zarr_format=zarr_format) as dst_store: |
There was a problem hiding this comment.
When -o names the input .zarr directory (or an alias resolving to it), opening dst_path with mode "w" clears the hierarchy while src still refers to that same store. The subsequent copy therefore fails or produces an empty/invalid store after destroying the original data; validate that source and destination are distinct before opening the destination.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Confirmed and fixed in 972fed5. Worse than described for Zarr: HDF5 refuses the second open, but Zarr completed successfully and left the store with zero nonzeros where the data had been. Resolved paths are now compared (symlinks and .. included) and the error points at --inplace, which writes to a temp file and swaps. Two tests, including the relative-path alias.
| if key in ("layers", "raw") and any( | ||
| p == key or p.startswith(f"{key}/") for p in by_parent | ||
| ): | ||
| _write_group(src[key], dst, key, by_parent, **options) |
There was a problem hiding this comment.
For an advertised explicit path such as obsm/X_pca, planning resolves and validates the dataset, but this branch only rebuilds layers and raw; obsm is copied wholesale instead. The command exits successfully while leaving the requested matrix unchanged, so _write needs to recurse through every parent represented in by_parent, not only these two names.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Confirmed and fixed in 972fed5. obsm/X_pca was planned, validated, then copied over unchanged with a success message — the docs advertise that path, so this was a silent no-op on an advertised feature. _write now descends into every parent that has a target, not just layers and raw. Test asserts the obsm matrix converts and that X is left alone.
| if csr: | ||
| block[major, minor] = data | ||
| else: | ||
| block[minor, major] = data |
There was a problem hiding this comment.
If a valid CSR/CSC input contains duplicate indices within a major axis, NumPy advanced assignment keeps only one of the repeated values here, whereas the sparse matrix represents their sum. Consequently --layout dense silently changes matrix values for non-canonical sparse inputs; the writes must accumulate repeated coordinates, for example with indexed addition.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Confirmed and fixed in 972fed5, using np.add.at.
Worth recording how nearly this was dismissed: my first reproduction went through ad.AnnData(X=dup), which canonicalises on write, so the duplicates never reached disk and the output looked correct. It only shows with a hand-built file. Verified: scipy's toarray()[0,1] is 3.0, the old code gave 2.0, the new code gives 3.0. The CSR→CSC path was already correct, since it moves duplicates rather than resolving them.
| per_bucket = max(MIN_BUCKET_ENTRIES, int(bucket_entries or chunk)) | ||
| n_buckets = max(1, int(np.ceil(nnz / per_bucket))) if nnz else 1 | ||
| n_buckets = min(n_buckets, n_major_out) or 1 | ||
| bounds = np.linspace(0, n_major_out, n_buckets + 1).astype(np.int64) |
There was a problem hiding this comment.
For a skewed matrix where most nonzeros land in a small range of output majors, these equal-width coordinate bounds place nearly all entries in one bucket. Pass 3 then loads that entire bucket and its sort arrays into memory, so the default “streaming” transpose can still OOM despite a small chunk; use the already-computed per-major counts to choose approximately equal-nnz bucket boundaries.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Confirmed and fixed in 972fed5. Measured on a matrix with 95% of entries in the first five columns: the largest of three equal-width buckets held 88% of the nonzeros, so the bucket rather than the chunk set the peak — the streaming path could still exhaust memory on exactly the shape it exists for, since single-cell matrices are skewed.
Bounds now come from the cumulative counts already computed in pass 1, so it costs nothing. An explicit bucket_entries is also honoured as given rather than floored, which is both more correct and what let the regression test exercise the bucketing at all.
| elif src.sparse: | ||
| # Keep what the source used. Defaulting to int64 silently doubled | ||
| # the index arrays of every int32 store that passed through. | ||
| idx_dtype = np.dtype(src.obj["indices"].dtype) |
There was a problem hiding this comment.
For a large sparse store with indices as int32 but indptr as int64 (for example, a narrow matrix with more than 2^31 nonzeros), the default inference selects int32 solely from indices; all writers then cast indptr to that dtype. Because the range check only runs when the user explicitly supplies --indices-dtype, an ordinary dtype or layout conversion silently overflows row/column offsets and corrupts the matrix.
Useful? React with 👍 / 👎.
There was a problem hiding this comment.
Confirmed and fixed in 972fed5. Reproduced with an int32-indices/int64-indptr store: the output narrowed indptr to int32.
The two widths are now tracked separately and each preserved from the source. I also took the second half of your point — the range check previously ran only when --indices-dtype was passed, but an inferred dtype can be too narrow just as easily, so it now runs always, and checks indices against the larger dimension and indptr against nnz rather than one bound for both.
All five reproduced before fixing and re-measured after. Two would have corrupted or destroyed data while reporting success, which is the worst shape a bug can take in something that rewrites people's matrices. An output path naming the input destroyed the store. HDF5 happens to refuse the second open; Zarr does not, so `-o` pointing at the input completed "successfully" and left a store with zero nonzeros where the data had been. Resolved paths are now compared, symlinks and `..` included, and the error points at --inplace, which writes to a temporary file and swaps. An explicitly named matrix outside layers and raw was silently not converted. `_write` recursed into those two groups only, so `obsm/X_pca` -- a path the docs advertise -- was planned, validated, then copied over unchanged with a success message. It now descends into every parent that has a target. Densifying dropped duplicate coordinates instead of summing them. Repeated indices within a major axis are legal in a CSR store and mean their sum, which is what scipy's own `toarray` gives; plain advanced assignment keeps whichever came last. np.add.at accumulates. Reaching this takes a hand-built file, since anndata canonicalises on write -- my first attempt to reproduce it failed for exactly that reason and looked like a false positive. indptr inherited the dtype of indices. A narrow matrix with more than 2^31 nonzeros needs int64 offsets over int32 coordinates; inferring one from the other narrowed the offsets and corrupted the matrix. The two widths are now tracked separately, each preserved from the source, and both range-checked -- always, not only when --indices-dtype is passed, since an inferred dtype can be too narrow just as easily. Transpose buckets were split by coordinate range, not by nonzero count. A skewed matrix therefore landed in one bucket, so the bucket rather than the chunk set the peak and the streaming path could still exhaust memory. Measured on a matrix with 95% of entries in the first five columns, the largest of three equal-width buckets held 88% of them. Bounds now come from the cumulative counts already computed in pass 1, and an explicit `bucket_entries` is honoured as given rather than floored, which is both more correct and what let the test exercise it. Each finding has a test that fails against the previous code. 961 tests pass, coverage 92.98%. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Closes #13.
Why
The issue asks to change the dtype of
X/dataandX/indiceson disk — count matrices stored as float64 waste disk and read time. Generalised here to any matrix, plus the CSR↔CSC conversion thatconcathas been telling people to do since it was written, with nothing to do it with.Works on
X, any layer,raw/X, or any 2-D array by path.Two refusals, both before the destination exists
--forceoverrides either, and both messages say what was measured.Transposing streams by default
It can't be done in one pass, so: count nonzeros per output major to build
indptr, scatter the input into buckets by where it lands, sort and write each bucket in order. Two extra passes over nnz buys a peak set by the chunk rather than by the matrix.--in-memorydoes it with onelexsortfor when the matrix fits.Correctness is held to scipy
Every path is compared against what scipy would have produced from the same matrix —
indptr,indicesanddataeach compared exactly, for both directions, both implementations, and across bucket counts from one to hundreds. CSR→CSC→CSR is asserted to be the identity.Three things this got wrong, all caught by measuring
Worth reading, because they're the interesting part.
The first version made files 8× larger on a conversion asked for to make them smaller.
_growablecreates with a fixed 65,536-element chunk, so a 2,400-nonzero output allocated 786 KB of mostly empty chunk; and the index dtype defaulted to int64, silently widening every int32 store. Outputs are now sized up front (nnz is known, and invariant under a transpose) with chunking and compression forwarded, and the index dtype is preserved unless asked. A test asserts the file shrinks, because that is the feature.The docstring claimed nothing is written before the checks, and that was false — the destination was created and obs/var copied into it before the first check ran.
plan_conversionnow resolves and checks every target before the output store is opened; two tests assert no output exists after a refusal.The streaming transpose held the whole matrix below 4M nonzeros, because the bucket size was a constant rather than derived from the chunk. The performance guard caught it: peak grew 15.9× for 32× the nonzeros, on the one path whose entire justification is that it doesn't. Buckets now follow
chunk.concat
Behaviour unchanged, but it now names the command to run — and is tested for the first time: CSR+CSC, CSR+dense, CSC+dense,
Xand layers, including a test that follows the suggested fix and concatenates successfully afterwards. The unreachable duplicate raise is kept and labelled a backstop.Guards and benchmark
Reads linear in nnz for a cast and in n_obs for a layout change; streaming transpose peak bounded while the in-memory one grows, compared in the same run; and an assertion that the default really is the streaming path — the two produce identical output, so nothing else can tell them apart.
convert-dtypeandconvert-layoutbenchmark cases, the latter running both implementations so the trade is visible.Verification
Review notes
--allmeansX, layers andraw/X— not obsm/varm, whose dtype is rarely what anyone is shrinking. Those still convert by explicit path.--layout sparsekeeps the source's major axis, so dense→sparse→dense is not a silent transpose.🤖 Generated with Claude Code