[tien-tong]

Summary

Rewrites the modelarrayio to-modelarray conversion path for speed while keeping the output layout that ModelArray reads (dataset paths, shapes, source names and metadata are unchanged; see Breaking changes for behaviour that did change). On real RBC freesurferpost data on CUBIC, one full 8,693-subject CIFTI scalar drops from about 3.3 h to about 7 min at one worker (27.5×), and to 69–91 s at eight workers. All 70 production scalars now convert in ≈22 min of job-array wall-clock, where the original code took 112 h for its first scalar and stalled at 15/70. Every output is bitwise-equal to the deployed reference.

The design, measurements and rejected alternatives are in design/TO_MODELARRAY_RUNTIME_OPTIMIZATION_PLAN.md.

Breaking changes

• --workers has a new meaning. It is now the maximum number of input sources loaded concurrently, for every modality and backend; it previously controlled concurrent TileDB scalar writes. Zero or negative values are rejected as a usage error.
• TileDB compresses with zstd level 5 by default (previously gzip level 4). HDF5 stays gzip level 4. Pass --compression gzip for the old TileDB filters.
• Invalid compression settings raise instead of being adjusted. Codec/backend pairs that do not exist (for example lzf with TileDB, which used to warn and continue) and out-of-range gzip levels (which used to be clamped) now raise ValueError before any data is read.
• Stricter CIFTI validation. Every subject's spatial axes must match the first subject's exactly. main silently accepted, and so misaligned, cohorts where subjects had the same structure names but different vertices, a pconn with permuted column parcels, or a second scalar on a different volume grid; these are now rejected.
• --s3-workers is deprecated (FutureWarning) and will be removed in the next breaking release. A value matching --workers is accepted; a conflicting value is rejected.
• h5_to_cifti writes float32. Results are no longer re-encoded in the template file's on-disk dtype (an int16 template with scaling used to quantize p-values).

Converting into an existing TileDB directory still merges into it, as on main: the scalars and parcel arrays written by a run replace their namesakes and every other scalar is kept.

What was slow

• HDF5 writes were capped by single-threaded chunk compression (~70 MB/s), not disk.
• TileDB writes reopened the array once per spatial stripe: one scalar became ~92 opens and ~92 fragments across hundreds of files.
• CIFTI loading parsed every file's XML header up front, GIL-bound, so extra workers made it slower.
• Every loader built the whole cohort across all scalars in memory before writing anything.
• NIfTI row assembly decoded full-volume float64 temporaries for the scalar and its mask on every row.
• S3 loading created a new client per file.

What changed

Loader (utils/loading.py, utils/{cifti,nifti,mif}.py, utils/s3_utils.py)

• One scalar at a time is staged into a typed, C-contiguous matrix through a bounded ThreadPoolExecutor.
• CIFTI values are read straight from each file's data block when its on-disk NIfTI-2 header (dims, dtype, no scaling) and CIFTI XML layout (every non-scalar MatrixIndicesMap, byte for byte) match the scalar's first file, skipping the per-file XML parse. dscalar, pscalar and pconn are supported; pconn blocks are reordered from NIfTI's column-major layout. Scaled data and any file that does not match is decoded by nibabel instead, and both paths produce identical values. Rows served this way are counted as fast_path_rows; MODELARRAYIO_CIFTI_RAW=0 disables the fast path.
• CIFTI axis validation compares a cheap spatial signature (all spatial-axis properties) instead of full label arrays.
• NIfTI masks are cached and selected with flat indices; decoding happens in the storage dtype.
• One shared S3 connection pool sized to the loader concurrency; MODELARRAYIO_S3_ANON=1 for public buckets. Broken git-annex symlinks report an actionable error instead of being read.

Storage (storage/{h5_storage,tiledb_storage,utils}.py)

• HDF5 gzip chunks are compressed on several threads and written with direct chunk writes. The thread count defaults to every CPU available to the process; MODELARRAYIO_HDF5_COMPRESSION_WORKERS caps it. Other HDF5 filters (lzf, none) use ordinary h5py writes.
• Each TileDB scalar matrix is written in a single dense submission: 1 fragment instead of one per stripe (350 → 8 files for MIF, 596 → 38 for CIFTI on the CUBIC comparison).
• TileDB's internal thread pools are recorded in telemetry but never capped.
• Compression codecs and levels are resolved and validated in one place for both backends; byte shuffle is on for compressed outputs and --compression none disables the filter pipeline.

CLI (cli/runtime.py, cli/to_modelarray.py, cli/{cifti,nifti,mif}_to_h5.py, cli/parser_utils.py, cli/h5_to_cifti.py)

• --workers N as described above; the resolved worker count is logged when it is above one.
• --dtype float32 stays the default; --dtype float64 is the documented opt-in that keeps source precision end to end.
• New --timing-jsonl PATH appends structured per-scalar and total runtime records (load/write phases, peak RSS, fragment counts, fast-path rows, TileDB pool settings).
• Outputs are written to a same-filesystem temporary path and validated before publication. An HDF5 file is replaced atomically; an existing TileDB directory is merged into with per-entry rollback.
• The group mask, index and directions inputs accept s3:// URIs. The cohort file must be local.
• tqdm dependency removed.

Docs

• README gains a "Fast conversion and input sources" section covering the worker semantics, dtype trade-off, S3/annex behaviour, TileDB merging and timing output.
• Sphinx: autosectionlabel prefixing fixes duplicate-label warnings; the sphinx-argparse patch tolerates versions without ArgParseDomain.

Benchmarks (benchmarks/)

New self-contained suite under benchmarks/to_modelarray/: a harness that runs the full CLI across NIfTI/CIFTI/MIF × HDF5/TileDB, cohort generators (public ABIDE from s3://fcp-indi, the MIF walkthrough demo, a reproducible CIFTI synthesizer), an old-commit launcher pinned to 9ab13dd, output-equality and R read gates, a comparison-report generator, and the CUBIC SLURM ladder (container env, sbatch templates, scalar-sharded job array). Site-specific locations are read from environment variables (see cubic/README.md). Pinned results live in baselines/; report PDFs are regenerated from those CSVs and are not committed.

Cross-modality, 100 subjects, original 9ab13dd vs this branch, one worker, CUBIC (baselines/cubic_cross_modality_9ab13dd_vs_optimized_20260922.csv):

Modality	Backend	Original (s)	Optimized (s)	Speedup	Output files	Values equal
NIfTI	HDF5	3.02	2.17	1.40×	1 → 1	yes
NIfTI	TileDB	3.41	2.33	1.46×	44 → 8	yes
MIF	HDF5	10.76	2.91	3.70×	1 → 1	yes
MIF	TileDB	14.52	3.53	4.11×	350 → 8	yes
CIFTI (3 scalars)	HDF5	385.78	12.12	31.8×	1 → 1	yes
CIFTI (3 scalars)	TileDB	396.47	12.56	31.6×	596 → 38	yes

Production scale, one full CIFTI scalar (area, 8,693 subjects), CUBIC (baselines/cubic_production_scale_9ab13dd_vs_optimized_20260923.csv):

Backend	Original, 1 worker (s)	Optimized, 1 worker (s)	4 workers (s)	8 workers (s)	Peak RSS
HDF5	11,854.8	429.9	120.9	90.6	~12 GB
TileDB	11,501.6	418.1	83.2	68.6	~20 GB

Peak memory is flat across workers because the staging matrix dominates. Other modalities were not swept at production scale, so the CLI keeps the conservative default of one worker.

The CUBIC numbers were measured on this PR's first commit. The review follow-up keeps the dscalar fast path (still one nibabel decode per scalar, with a slightly lower per-file cost on 91,282-greyordinate files), adds the pconn fast path, and was checked against nibabel and main on dscalar, pscalar and pconn cohorts (NaNs, float64, big-endian, scaled int16 and mixed encodings), on both backends and both dtypes, with the fast path on and off.

Deferred

• Overlapping load and write: loading dominates (≈1.7 s load vs ≈0.1 s write per example scalar), so hiding the write buys about five percent; its real value would be hiding S3 latency.
• --only-scalar: per-scalar cohort CSVs select the scalar for the sharded HPC recipe instead.
• A CLI option for the HDF5 compression thread count (for example defaulting to SLURM_CPUS_PER_TASK); for now only MODELARRAYIO_HDF5_COMPRESSION_WORKERS caps it.

Tests

New: test_runtime.py, test_runtime_integration_matrix.py, test_loading.py, test_storage_utils.py, test_s3_utils.py, test_benchmark_harness.py, test_comparison_report.py, test_cubic_tools.py, test_make_cifti_cohort.py. Extended: CIFTI/NIfTI utils, HDF5/TileDB storage, CLI utils and parser tests, including regression tests for the CIFTI fast path (pconn order and every map's digest, rejection branches, no per-file nibabel decoding), sequential TileDB conversions into one root, and CIFTI export dtype and round trip.

🤖 Generated with Claude Code

[codecov-commenter]

Codecov Report

❌ Patch coverage is 87.12999% with 200 lines in your changes missing coverage. Please review.
✅ Project coverage is 89.09%. Comparing base (9ab13dd) to head (3b5b43f).

Files with missing lines	Patch %	Lines
src/modelarrayio/cli/runtime.py	84.22%	32 Missing and 24 partials ⚠️
src/modelarrayio/utils/cifti.py	89.20%	22 Missing and 12 partials ⚠️
src/modelarrayio/utils/loading.py	88.31%	17 Missing and 10 partials ⚠️
src/modelarrayio/storage/h5_storage.py	84.90%	15 Missing and 9 partials ⚠️
src/modelarrayio/storage/utils.py	75.00%	9 Missing and 9 partials ⚠️
src/modelarrayio/utils/nifti.py	86.29%	11 Missing and 6 partials ⚠️
src/modelarrayio/utils/s3_utils.py	84.00%	7 Missing and 5 partials ⚠️
src/modelarrayio/storage/tiledb_storage.py	90.16%	4 Missing and 2 partials ⚠️
src/modelarrayio/utils/mif.py	79.31%	4 Missing and 2 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #110      +/-   ##
==========================================
- Coverage   90.81%   89.09%   -1.73%     
==========================================
  Files          21       23       +2     
  Lines        1470     2649    +1179     
  Branches      305      538     +233     
==========================================
+ Hits         1335     2360    +1025     
- Misses         76      177     +101     
- Partials       59      112      +53     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.