From e9e60973b45f6d14e0439a39f5002f46204dc167 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 26 Mar 2026 14:45:30 -0400 Subject: [PATCH 01/20] Support scalar-columns in other converters. --- src/modelarrayio/cli/cifti_to_h5.py | 27 ++++----- src/modelarrayio/cli/mif_to_h5.py | 77 ++++++++++++++++--------- src/modelarrayio/cli/nifti_to_h5.py | 75 +++++++++++++++++------- src/modelarrayio/cli/parser_utils.py | 22 +++---- src/modelarrayio/utils/cifti.py | 49 +--------------- src/modelarrayio/utils/fixels.py | 66 +++++++++++++++++++++ src/modelarrayio/utils/misc.py | 86 ++++++++++++++++++++++++++++ src/modelarrayio/utils/voxels.py | 14 ++++- 8 files changed, 289 insertions(+), 127 deletions(-) create mode 100644 src/modelarrayio/utils/misc.py diff --git a/src/modelarrayio/cli/cifti_to_h5.py b/src/modelarrayio/cli/cifti_to_h5.py index 6fec34b..d316e08 100644 --- a/src/modelarrayio/cli/cifti_to_h5.py +++ b/src/modelarrayio/cli/cifti_to_h5.py @@ -13,14 +13,13 @@ from tqdm import tqdm from modelarrayio.cli import utils as cli_utils -from modelarrayio.cli.parser_utils import add_scalar_columns_arg, add_to_modelarray_args +from modelarrayio.cli.parser_utils import add_to_modelarray_args from modelarrayio.utils.cifti import ( - _build_scalar_sources, - _cohort_to_long_dataframe, - _load_cohort_cifti, brain_names_to_dataframe, extract_cifti_scalar_data, + load_cohort_cifti, ) +from modelarrayio.utils.misc import build_scalar_sources, cohort_to_long_dataframe logger = logging.getLogger(__name__) @@ -47,7 +46,7 @@ def cifti_to_h5( Path to a csv with demographic info and paths to data backend : :obj:`str` Backend to use for storage (``'hdf5'`` or ``'tiledb'``) - output : :obj:`str` + output : :obj:`pathlib.Path` Output path. For the hdf5 backend, path to an .h5 file; for the tiledb backend, path to a .tdb directory. storage_dtype : :obj:`str` @@ -77,19 +76,18 @@ def cifti_to_h5( 0 if successful, 1 if failed. """ cohort_df = pd.read_csv(cohort_file) - cohort_long = _cohort_to_long_dataframe(cohort_df, scalar_columns=scalar_columns) - output_path = Path(output) + cohort_long = cohort_to_long_dataframe(cohort_df, scalar_columns=scalar_columns) if cohort_long.empty: raise ValueError('Cohort file does not contain any scalar entries after normalization.') - scalar_sources = _build_scalar_sources(cohort_long) + scalar_sources = build_scalar_sources(cohort_long) if not scalar_sources: raise ValueError('Unable to derive scalar sources from cohort file.') if backend == 'hdf5': - scalars, last_brain_names = _load_cohort_cifti(cohort_long, s3_workers) + scalars, last_brain_names = load_cohort_cifti(cohort_long, s3_workers) greyordinate_table, structure_names = brain_names_to_dataframe(last_brain_names) - output_path = cli_utils.prepare_output_parent(output_path) - with h5py.File(output_path, 'w') as h5_file: + output = cli_utils.prepare_output_parent(output) + with h5py.File(output, 'w') as h5_file: cli_utils.write_table_dataset( h5_file, 'greyordinates', @@ -107,9 +105,9 @@ def cifti_to_h5( chunk_voxels=chunk_voxels, target_chunk_mb=target_chunk_mb, ) - return int(not output_path.exists()) + return int(not output.exists()) - output_path.mkdir(parents=True, exist_ok=True) + output.mkdir(parents=True, exist_ok=True) if not scalar_sources: return 0 @@ -127,7 +125,7 @@ def _process_scalar_job(scalar_name, source_files): if rows: cli_utils.write_tiledb_scalar_matrices( - output_path, + output, {scalar_name: rows}, {scalar_name: source_files}, storage_dtype=storage_dtype, @@ -178,5 +176,4 @@ def _parse_cifti_to_h5(): formatter_class=argparse.ArgumentDefaultsHelpFormatter, ) add_to_modelarray_args(parser, default_output='greyordinatearray.h5') - add_scalar_columns_arg(parser) return parser diff --git a/src/modelarrayio/cli/mif_to_h5.py b/src/modelarrayio/cli/mif_to_h5.py index 38ed1cc..b651ddf 100644 --- a/src/modelarrayio/cli/mif_to_h5.py +++ b/src/modelarrayio/cli/mif_to_h5.py @@ -4,7 +4,8 @@ import argparse import logging -from collections import defaultdict +import os +from concurrent.futures import ThreadPoolExecutor, as_completed from functools import partial from pathlib import Path @@ -14,7 +15,8 @@ from modelarrayio.cli import utils as cli_utils from modelarrayio.cli.parser_utils import _is_file, add_to_modelarray_args -from modelarrayio.utils.fixels import gather_fixels, mif_to_nifti2 +from modelarrayio.utils.fixels import gather_fixels, load_cohort_mif +from modelarrayio.utils.misc import cohort_to_long_dataframe logger = logging.getLogger(__name__) @@ -33,6 +35,7 @@ def mif_to_h5( target_chunk_mb=2.0, workers=None, s3_workers=1, + scalar_columns=None, ): """Load all fixeldb data and write to an HDF5 or TileDB file. @@ -75,25 +78,20 @@ def mif_to_h5( """ # gather fixel data fixel_table, voxel_table = gather_fixels(index_file, directions_file) - output_path = Path(output) - # gather cohort data cohort_df = pd.read_csv(cohort_file) + cohort_long = cohort_to_long_dataframe(cohort_df, scalar_columns=scalar_columns) + if cohort_long.empty: + raise ValueError('Cohort file does not contain any scalar entries after normalization.') - # upload each cohort's data - scalars = defaultdict(list) - sources_lists = defaultdict(list) logger.info('Extracting .mif data...') - for row in tqdm(cohort_df.itertuples(index=False), total=cohort_df.shape[0]): - scalar_file = row.source_file - _scalar_img, scalar_data = mif_to_nifti2(scalar_file) - scalars[row.scalar_name].append(scalar_data) - sources_lists[row.scalar_name].append(row.source_file) + scalars, sources_lists = load_cohort_mif(cohort_long, s3_workers) + if not sources_lists: + raise ValueError('Unable to derive scalar sources from cohort file.') - # Write the output if backend == 'hdf5': - output_path = cli_utils.prepare_output_parent(output_path) - with h5py.File(output_path, 'w') as h5_file: + output = cli_utils.prepare_output_parent(output) + with h5py.File(output, 'w') as h5_file: cli_utils.write_table_dataset(h5_file, 'fixels', fixel_table) cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table) cli_utils.write_hdf5_scalar_matrices( @@ -107,19 +105,42 @@ def mif_to_h5( chunk_voxels=chunk_voxels, target_chunk_mb=target_chunk_mb, ) - return int(not output_path.exists()) - - cli_utils.write_tiledb_scalar_matrices( - output_path, - scalars, - sources_lists, - storage_dtype=storage_dtype, - compression=compression, - compression_level=compression_level, - shuffle=shuffle, - chunk_voxels=chunk_voxels, - target_chunk_mb=target_chunk_mb, - ) + return int(not output.exists()) + + output.mkdir(parents=True, exist_ok=True) + + scalar_names = list(sources_lists.keys()) + worker_count = workers if isinstance(workers, int) and workers > 0 else None + if worker_count is None: + cpu_count = os.cpu_count() or 1 + worker_count = min(len(scalar_names), max(1, cpu_count)) + else: + worker_count = min(len(scalar_names), worker_count) + + def _write_scalar_job(scalar_name): + cli_utils.write_tiledb_scalar_matrices( + output, + {scalar_name: scalars[scalar_name]}, + {scalar_name: sources_lists[scalar_name]}, + storage_dtype=storage_dtype, + compression=compression, + compression_level=compression_level, + shuffle=shuffle, + chunk_voxels=chunk_voxels, + target_chunk_mb=target_chunk_mb, + ) + + if worker_count <= 1: + for scalar_name in scalar_names: + _write_scalar_job(scalar_name) + else: + with ThreadPoolExecutor(max_workers=worker_count) as executor: + futures = { + executor.submit(_write_scalar_job, scalar_name): scalar_name + for scalar_name in scalar_names + } + for future in tqdm(as_completed(futures), total=len(futures), desc='TileDB scalars'): + future.result() return 0 diff --git a/src/modelarrayio/cli/nifti_to_h5.py b/src/modelarrayio/cli/nifti_to_h5.py index 28e2422..38cee30 100644 --- a/src/modelarrayio/cli/nifti_to_h5.py +++ b/src/modelarrayio/cli/nifti_to_h5.py @@ -4,6 +4,8 @@ import argparse import logging +import os +from concurrent.futures import ThreadPoolExecutor, as_completed from functools import partial from pathlib import Path @@ -11,10 +13,12 @@ import nibabel as nb import numpy as np import pandas as pd +from tqdm import tqdm from modelarrayio.cli import utils as cli_utils -from modelarrayio.cli.parser_utils import _is_file, add_to_modelarray_args -from modelarrayio.utils.voxels import _load_cohort_voxels +from modelarrayio.cli.parser_utils import _is_file, add_scalar_columns_arg, add_to_modelarray_args +from modelarrayio.utils.misc import cohort_to_long_dataframe +from modelarrayio.utils.voxels import load_cohort_voxels logger = logging.getLogger(__name__) @@ -32,6 +36,7 @@ def nifti_to_h5( target_chunk_mb=2.0, workers=None, s3_workers=1, + scalar_columns=None, ): """Load all volume data and write to an HDF5 or TileDB file. @@ -43,7 +48,7 @@ def nifti_to_h5( Path to a CSV with demographic info and paths to data. backend : :obj:`str` Storage backend (``'hdf5'`` or ``'tiledb'``). - output : :obj:`str` + output : :obj:`pathlib.Path` Output path. For the hdf5 backend, path to an .h5 file; for the tiledb backend, path to a .tdb directory. storage_dtype : :obj:`str` @@ -65,13 +70,14 @@ def nifti_to_h5( s3_workers : :obj:`int` Number of parallel workers for S3 downloads. Default 1. """ - cohort_df = pd.read_csv(cohort_file) - output_path = Path(output) - group_mask_img = nb.load(group_mask_file) group_mask_matrix = group_mask_img.get_fdata() > 0 voxel_coords = np.column_stack(np.nonzero(group_mask_matrix)) + cohort_df = pd.read_csv(cohort_file) + cohort_long = cohort_to_long_dataframe(cohort_df, scalar_columns=scalar_columns) + if cohort_long.empty: + raise ValueError('Cohort file does not contain any scalar entries after normalization.') voxel_table = pd.DataFrame( { 'voxel_id': np.arange(voxel_coords.shape[0]), @@ -82,11 +88,13 @@ def nifti_to_h5( ) logger.info('Extracting NIfTI data...') - scalars, sources_lists = _load_cohort_voxels(cohort_df, group_mask_matrix, s3_workers) + scalars, sources_lists = load_cohort_voxels(cohort_long, group_mask_matrix, s3_workers) + if not sources_lists: + raise ValueError('Unable to derive scalar sources from cohort file.') if backend == 'hdf5': - output_path = cli_utils.prepare_output_parent(output_path) - with h5py.File(output_path, 'w') as h5_file: + output = cli_utils.prepare_output_parent(output) + with h5py.File(output, 'w') as h5_file: cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table) cli_utils.write_hdf5_scalar_matrices( h5_file, @@ -99,19 +107,42 @@ def nifti_to_h5( chunk_voxels=chunk_voxels, target_chunk_mb=target_chunk_mb, ) - return int(not output_path.exists()) - - cli_utils.write_tiledb_scalar_matrices( - output_path, - scalars, - sources_lists, - storage_dtype=storage_dtype, - compression=compression, - compression_level=compression_level, - shuffle=shuffle, - chunk_voxels=chunk_voxels, - target_chunk_mb=target_chunk_mb, - ) + return int(not output.exists()) + + output.mkdir(parents=True, exist_ok=True) + + scalar_names = list(sources_lists.keys()) + worker_count = workers if isinstance(workers, int) and workers > 0 else None + if worker_count is None: + cpu_count = os.cpu_count() or 1 + worker_count = min(len(scalar_names), max(1, cpu_count)) + else: + worker_count = min(len(scalar_names), worker_count) + + def _write_scalar_job(scalar_name): + cli_utils.write_tiledb_scalar_matrices( + output, + {scalar_name: scalars[scalar_name]}, + {scalar_name: sources_lists[scalar_name]}, + storage_dtype=storage_dtype, + compression=compression, + compression_level=compression_level, + shuffle=shuffle, + chunk_voxels=chunk_voxels, + target_chunk_mb=target_chunk_mb, + ) + + if worker_count <= 1: + for scalar_name in scalar_names: + _write_scalar_job(scalar_name) + else: + with ThreadPoolExecutor(max_workers=worker_count) as executor: + futures = { + executor.submit(_write_scalar_job, scalar_name): scalar_name + for scalar_name in scalar_names + } + for future in tqdm(as_completed(futures), total=len(futures), desc='TileDB scalars'): + future.result() return 0 diff --git a/src/modelarrayio/cli/parser_utils.py b/src/modelarrayio/cli/parser_utils.py index 55c2c4a..05865fe 100644 --- a/src/modelarrayio/cli/parser_utils.py +++ b/src/modelarrayio/cli/parser_utils.py @@ -21,6 +21,15 @@ def add_to_modelarray_args(parser, default_output='output.h5'): ), default=default_output, ) + parser.add_argument( + '--scalar-columns', + '--scalar_columns', + nargs='+', + help=( + 'Column names containing scalar file paths when the cohort table is in wide format. ' + 'If omitted, the cohort file must include "scalar_name" and "source_file" columns.' + ), + ) parser.add_argument( '--backend', help='Storage backend for subject-by-element matrix', @@ -110,19 +119,6 @@ def add_to_modelarray_args(parser, default_output='output.h5'): return parser -def add_scalar_columns_arg(parser): - parser.add_argument( - '--scalar-columns', - '--scalar_columns', - nargs='+', - help=( - 'Column names containing scalar file paths when the cohort table is in wide format. ' - "If omitted, the cohort file must include 'scalar_name' and 'source_file' columns." - ), - ) - return parser - - def add_log_level_arg(parser): parser.add_argument( '--log-level', diff --git a/src/modelarrayio/utils/cifti.py b/src/modelarrayio/utils/cifti.py index 6dfdfe8..b3b9a66 100644 --- a/src/modelarrayio/utils/cifti.py +++ b/src/modelarrayio/utils/cifti.py @@ -2,7 +2,7 @@ from __future__ import annotations -from collections import OrderedDict, defaultdict +from collections import defaultdict from concurrent.futures import ThreadPoolExecutor, as_completed from pathlib import Path @@ -14,51 +14,6 @@ from modelarrayio.utils.s3_utils import load_nibabel -def _cohort_to_long_dataframe(cohort_df, scalar_columns=None): - scalar_columns = [col for col in (scalar_columns or []) if col] - if scalar_columns: - missing = [col for col in scalar_columns if col not in cohort_df.columns] - if missing: - raise ValueError(f'Wide-format cohort is missing scalar columns: {missing}') - records = [] - selected_columns = cohort_df[scalar_columns] - for row_values in selected_columns.itertuples(index=False, name=None): - for scalar_col, source_val in zip(scalar_columns, row_values, strict=True): - if pd.isna(source_val) or source_val is None: - continue - source_str = str(source_val).strip() - if not source_str: - continue - records.append({'scalar_name': scalar_col, 'source_file': source_str}) - return pd.DataFrame.from_records(records, columns=['scalar_name', 'source_file']) - - required = {'scalar_name', 'source_file'} - missing = required - set(cohort_df.columns) - if missing: - raise ValueError( - f'Cohort file must contain columns {sorted(required)} when ' - '--scalar-columns is not used.' - ) - - long_df = cohort_df[list(required)].copy() - long_df = long_df.dropna(subset=['scalar_name', 'source_file']) - long_df['scalar_name'] = long_df['scalar_name'].astype(str).str.strip() - long_df['source_file'] = long_df['source_file'].astype(str).str.strip() - long_df = long_df[(long_df['scalar_name'] != '') & (long_df['source_file'] != '')] - return long_df.reset_index(drop=True) - - -def _build_scalar_sources(long_df): - scalar_sources = OrderedDict() - for row in long_df.itertuples(index=False): - scalar = str(row.scalar_name) - source = str(row.source_file) - if not scalar or not source: - continue - scalar_sources.setdefault(scalar, []).append(source) - return scalar_sources - - def extract_cifti_scalar_data(cifti_file, reference_brain_names=None): """Load a scalar cifti file and get its data and mapping @@ -137,7 +92,7 @@ def brain_names_to_dataframe(brain_names): return greyordinate_df, structure_name_strings -def _load_cohort_cifti(cohort_long, s3_workers): +def load_cohort_cifti(cohort_long, s3_workers): """Load all CIFTI scalar rows from the cohort, optionally in parallel. The first file is always loaded serially to obtain the reference brain diff --git a/src/modelarrayio/utils/fixels.py b/src/modelarrayio/utils/fixels.py index 38e650e..43755b3 100644 --- a/src/modelarrayio/utils/fixels.py +++ b/src/modelarrayio/utils/fixels.py @@ -3,11 +3,14 @@ import shutil import subprocess import tempfile +from collections import defaultdict +from concurrent.futures import ThreadPoolExecutor, as_completed from pathlib import Path import nibabel as nb import numpy as np import pandas as pd +from tqdm import tqdm def find_mrconvert(): @@ -91,6 +94,69 @@ def mif_to_nifti2(mif_file): return nifti2_img, data +def load_cohort_mif(cohort_long, s3_workers): + """Load all MIF scalar rows from the cohort, optionally in parallel. + + When s3_workers > 1, a ThreadPoolExecutor is used to run mrconvert + calls concurrently (subprocess calls release the GIL). Results arrive + via as_completed and are indexed by (scalar_name, subj_idx) so the + final ordered lists are reconstructed correctly regardless of completion + order. + + Parameters + ---------- + cohort_long : :obj:`pandas.DataFrame` + Long-format cohort dataframe with columns 'scalar_name' and 'source_file'. + s3_workers : :obj:`int` + Number of parallel workers for loading. + + Returns + ------- + scalars : dict[str, list[np.ndarray]] + Per-scalar ordered list of 1-D subject arrays, ready for stripe-write. + sources_lists : dict[str, list[str]] + Per-scalar ordered list of source file paths (for HDF5 metadata). + """ + scalar_subj_counter = defaultdict(int) + jobs = [] + sources_lists = defaultdict(list) + + for row in cohort_long.itertuples(index=False): + sn = row.scalar_name + subj_idx = scalar_subj_counter[sn] + scalar_subj_counter[sn] += 1 + src = row.source_file + jobs.append((sn, subj_idx, src)) + sources_lists[sn].append(src) + + def _worker(job): + sn, subj_idx, src = job + _img, data = mif_to_nifti2(src) + return sn, subj_idx, data + + if s3_workers > 1: + results = defaultdict(dict) + with ThreadPoolExecutor(max_workers=s3_workers) as pool: + futures = {pool.submit(_worker, job): job for job in jobs} + for future in tqdm( + as_completed(futures), + total=len(futures), + desc='Loading MIF data', + ): + sn, subj_idx, data = future.result() + results[sn][subj_idx] = data + scalars = { + sn: [results[sn][i] for i in range(cnt)] for sn, cnt in scalar_subj_counter.items() + } + else: + scalars = defaultdict(list) + for job in tqdm(jobs, desc='Loading MIF data'): + sn, subj_idx, data = _worker(job) + scalars[sn].append(data) + + return scalars, sources_lists + + def gather_fixels(index_file, directions_file): """Load the index and directions files to get lookup tables. diff --git a/src/modelarrayio/utils/misc.py b/src/modelarrayio/utils/misc.py new file mode 100644 index 0000000..f7fcd61 --- /dev/null +++ b/src/modelarrayio/utils/misc.py @@ -0,0 +1,86 @@ +"""Miscellaneous utility functions.""" + +from __future__ import annotations + +from collections import OrderedDict + +import pandas as pd + + +def cohort_to_long_dataframe(cohort_df, scalar_columns=None): + """Convert a wide-format cohort dataframe to a long-format dataframe. + + Parameters + ---------- + cohort_df : :obj:`pandas.DataFrame` + Wide-format cohort dataframe + scalar_columns : :obj:`list` + List of scalar columns to use. If provided, these columns are treated as + file-path columns and melted into 'scalar_name'/'source_file' rows. All + remaining columns (e.g. 'source_mask_file') are broadcast to every output + row. If not provided, the dataframe is treated as already long-format. + + Returns + ------- + long_df : :obj:`pandas.DataFrame` + Long-format cohort dataframe with columns 'scalar_name', 'source_file', + and any non-scalar columns from the input. + """ + scalar_columns = [col for col in (scalar_columns or []) if col] + if scalar_columns: + missing = [col for col in scalar_columns if col not in cohort_df.columns] + if missing: + raise ValueError(f'Wide-format cohort is missing scalar columns: {missing}') + extra_columns = [col for col in cohort_df.columns if col not in scalar_columns] + records = [] + for _, row in cohort_df.iterrows(): + extra = {col: row[col] for col in extra_columns} + for scalar_col in scalar_columns: + source_val = row[scalar_col] + if pd.isna(source_val) or source_val is None: + continue + source_str = str(source_val).strip() + if not source_str: + continue + records.append({'scalar_name': scalar_col, 'source_file': source_str, **extra}) + output_columns = ['scalar_name', 'source_file'] + extra_columns + return pd.DataFrame.from_records(records, columns=output_columns) + + required = {'scalar_name', 'source_file'} + missing = required - set(cohort_df.columns) + if missing: + raise ValueError( + f'Cohort file must contain columns {sorted(required)} when ' + '--scalar-columns is not used.' + ) + + long_df = cohort_df.copy() + long_df = long_df.dropna(subset=['scalar_name', 'source_file']) + long_df['scalar_name'] = long_df['scalar_name'].astype(str).str.strip() + long_df['source_file'] = long_df['source_file'].astype(str).str.strip() + long_df = long_df[(long_df['scalar_name'] != '') & (long_df['source_file'] != '')] + return long_df.reset_index(drop=True) + + +def build_scalar_sources(long_df): + """Build a dictionary of scalar sources from a long dataframe. + + Parameters + ---------- + long_df : :obj:`pandas.DataFrame` + Long-format cohort dataframe with columns 'scalar_name' and 'source_file'. + + Returns + ------- + scalar_sources : :obj:`OrderedDict` + Dictionary of scalar sources. + Keys are scalar names, values are lists of source files. + """ + scalar_sources = OrderedDict() + for row in long_df.itertuples(index=False): + scalar = str(row.scalar_name) + source = str(row.source_file) + if not scalar or not source: + continue + scalar_sources.setdefault(scalar, []).append(source) + return scalar_sources diff --git a/src/modelarrayio/utils/voxels.py b/src/modelarrayio/utils/voxels.py index 0912b84..663eaaa 100644 --- a/src/modelarrayio/utils/voxels.py +++ b/src/modelarrayio/utils/voxels.py @@ -13,7 +13,7 @@ from modelarrayio.utils.s3_utils import load_nibabel -def _load_cohort_voxels(cohort_df, group_mask_matrix, s3_workers): +def load_cohort_voxels(cohort_long, group_mask_matrix, s3_workers): """Load all voxel rows from the cohort, optionally in parallel. When s3_workers > 1, a ThreadPoolExecutor is used. Threads share memory so @@ -21,6 +21,16 @@ def _load_cohort_voxels(cohort_df, group_mask_matrix, s3_workers): arrive via as_completed and are indexed by (scalar_name, subj_idx) so the final ordered lists are reconstructed correctly regardless of completion order. + Parameters + ---------- + cohort_long : :obj:`pandas.DataFrame` + Long-format cohort dataframe with columns 'scalar_name', 'source_file', + and 'source_mask_file'. + group_mask_matrix : :obj:`numpy.ndarray` + Boolean group mask array. + s3_workers : :obj:`int` + Number of parallel workers for loading. + Returns ------- scalars : dict[str, list[np.ndarray]] @@ -32,7 +42,7 @@ def _load_cohort_voxels(cohort_df, group_mask_matrix, s3_workers): jobs = [] sources_lists = defaultdict(list) - for row in cohort_df.itertuples(index=False): + for row in cohort_long.itertuples(index=False): sn = row.scalar_name subj_idx = scalar_subj_counter[sn] scalar_subj_counter[sn] += 1 From 66128ad5eec05b494440ae2ad77577175062c36a Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 26 Mar 2026 14:45:41 -0400 Subject: [PATCH 02/20] Run ruff. --- src/modelarrayio/cli/nifti_to_h5.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/modelarrayio/cli/nifti_to_h5.py b/src/modelarrayio/cli/nifti_to_h5.py index 38cee30..e472f7c 100644 --- a/src/modelarrayio/cli/nifti_to_h5.py +++ b/src/modelarrayio/cli/nifti_to_h5.py @@ -16,7 +16,7 @@ from tqdm import tqdm from modelarrayio.cli import utils as cli_utils -from modelarrayio.cli.parser_utils import _is_file, add_scalar_columns_arg, add_to_modelarray_args +from modelarrayio.cli.parser_utils import _is_file, add_to_modelarray_args from modelarrayio.utils.misc import cohort_to_long_dataframe from modelarrayio.utils.voxels import load_cohort_voxels From 962a807f742bac4903b67e1e91678dac4396d9f1 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 26 Mar 2026 14:48:05 -0400 Subject: [PATCH 03/20] Update test_cifti_cohort.py --- test/test_cifti_cohort.py | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/test/test_cifti_cohort.py b/test/test_cifti_cohort.py index 8117387..a5c3b9b 100644 --- a/test/test_cifti_cohort.py +++ b/test/test_cifti_cohort.py @@ -6,11 +6,8 @@ import pandas as pd import pytest -from modelarrayio.utils.cifti import ( - _build_scalar_sources, - _cohort_to_long_dataframe, - brain_names_to_dataframe, -) +from modelarrayio.utils.cifti import brain_names_to_dataframe +from modelarrayio.utils.misc import build_scalar_sources, cohort_to_long_dataframe def test_cohort_long_format_preserves_rows() -> None: @@ -21,7 +18,7 @@ def test_cohort_long_format_preserves_rows() -> None: 'extra_col': [1, 2], } ) - long_df = _cohort_to_long_dataframe(df) + long_df = cohort_to_long_dataframe(df) assert len(long_df) == 2 assert set(long_df.columns) == {'scalar_name', 'source_file'} assert long_df.iloc[0]['scalar_name'] == 'THICK' @@ -34,7 +31,7 @@ def test_cohort_long_format_strips_and_drops_empty() -> None: 'source_file': [' a.nii ', ' b.nii '], } ) - long_df = _cohort_to_long_dataframe(df) + long_df = cohort_to_long_dataframe(df) assert len(long_df) == 1 assert long_df.iloc[0]['scalar_name'] == 'THICK' @@ -47,7 +44,7 @@ def test_cohort_wide_format_expands_columns() -> None: 'FA': ['f1.nii', ''], } ) - long_df = _cohort_to_long_dataframe(df, scalar_columns=['THICK', 'FA']) + long_df = cohort_to_long_dataframe(df, scalar_columns=['THICK', 'FA']) # Row 2 has empty FA — skipped assert len(long_df) == 3 scalars = set(long_df['scalar_name']) @@ -57,13 +54,13 @@ def test_cohort_wide_format_expands_columns() -> None: def test_cohort_wide_format_missing_scalar_column_raises() -> None: df = pd.DataFrame({'THICK': ['a.nii']}) with pytest.raises(ValueError, match='missing scalar columns'): - _cohort_to_long_dataframe(df, scalar_columns=['THICK', 'MISSING']) + cohort_to_long_dataframe(df, scalar_columns=['THICK', 'MISSING']) def test_cohort_long_missing_required_raises() -> None: df = pd.DataFrame({'only_this': [1]}) with pytest.raises(ValueError, match='scalar_name'): - _cohort_to_long_dataframe(df) + cohort_to_long_dataframe(df) def test_build_scalar_sources_ordering() -> None: @@ -73,7 +70,7 @@ def test_build_scalar_sources_ordering() -> None: 'source_file': ['x1', 'x2', 'y1'], } ) - src = _build_scalar_sources(long_df) + src = build_scalar_sources(long_df) assert list(src.keys()) == ['A', 'B'] assert src['A'] == ['x1', 'x2'] assert src['B'] == ['y1'] From 367541c21035e52d308303216585e7d848b2a866 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Thu, 26 Mar 2026 14:50:27 -0400 Subject: [PATCH 04/20] Update test_cifti_cohort.py --- test/test_cifti_cohort.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_cifti_cohort.py b/test/test_cifti_cohort.py index a5c3b9b..eee8e00 100644 --- a/test/test_cifti_cohort.py +++ b/test/test_cifti_cohort.py @@ -20,7 +20,7 @@ def test_cohort_long_format_preserves_rows() -> None: ) long_df = cohort_to_long_dataframe(df) assert len(long_df) == 2 - assert set(long_df.columns) == {'scalar_name', 'source_file'} + assert set(long_df.columns) == {'extra_col', 'scalar_name', 'source_file'} assert long_df.iloc[0]['scalar_name'] == 'THICK' From 7101570613decac090bb52c48172b7631399db7b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 09:52:37 -0400 Subject: [PATCH 05/20] Rename util submodules. --- src/modelarrayio/utils/{fixels.py => mif.py} | 0 src/modelarrayio/utils/{voxels.py => nifti.py} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename src/modelarrayio/utils/{fixels.py => mif.py} (100%) rename src/modelarrayio/utils/{voxels.py => nifti.py} (100%) diff --git a/src/modelarrayio/utils/fixels.py b/src/modelarrayio/utils/mif.py similarity index 100% rename from src/modelarrayio/utils/fixels.py rename to src/modelarrayio/utils/mif.py diff --git a/src/modelarrayio/utils/voxels.py b/src/modelarrayio/utils/nifti.py similarity index 100% rename from src/modelarrayio/utils/voxels.py rename to src/modelarrayio/utils/nifti.py From c15a2d5de068992b4930ff8986852ff6f9894c9b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 09:53:19 -0400 Subject: [PATCH 06/20] Fix imports. --- src/modelarrayio/cli/h5_to_mif.py | 2 +- src/modelarrayio/cli/mif_to_h5.py | 2 +- src/modelarrayio/cli/nifti_to_h5.py | 2 +- test/test_voxels_utils.py | 2 +- 4 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/modelarrayio/cli/h5_to_mif.py b/src/modelarrayio/cli/h5_to_mif.py index bb6bfed..99e0d3b 100644 --- a/src/modelarrayio/cli/h5_to_mif.py +++ b/src/modelarrayio/cli/h5_to_mif.py @@ -14,7 +14,7 @@ from modelarrayio.cli import utils as cli_utils from modelarrayio.cli.parser_utils import _is_file, add_from_modelarray_args, add_log_level_arg -from modelarrayio.utils.fixels import mif_to_nifti2, nifti2_to_mif +from modelarrayio.utils.mif import mif_to_nifti2, nifti2_to_mif logger = logging.getLogger(__name__) diff --git a/src/modelarrayio/cli/mif_to_h5.py b/src/modelarrayio/cli/mif_to_h5.py index b651ddf..326cb87 100644 --- a/src/modelarrayio/cli/mif_to_h5.py +++ b/src/modelarrayio/cli/mif_to_h5.py @@ -15,7 +15,7 @@ from modelarrayio.cli import utils as cli_utils from modelarrayio.cli.parser_utils import _is_file, add_to_modelarray_args -from modelarrayio.utils.fixels import gather_fixels, load_cohort_mif +from modelarrayio.utils.mif import gather_fixels, load_cohort_mif from modelarrayio.utils.misc import cohort_to_long_dataframe logger = logging.getLogger(__name__) diff --git a/src/modelarrayio/cli/nifti_to_h5.py b/src/modelarrayio/cli/nifti_to_h5.py index e472f7c..34a12f5 100644 --- a/src/modelarrayio/cli/nifti_to_h5.py +++ b/src/modelarrayio/cli/nifti_to_h5.py @@ -18,7 +18,7 @@ from modelarrayio.cli import utils as cli_utils from modelarrayio.cli.parser_utils import _is_file, add_to_modelarray_args from modelarrayio.utils.misc import cohort_to_long_dataframe -from modelarrayio.utils.voxels import load_cohort_voxels +from modelarrayio.utils.nifti import load_cohort_voxels logger = logging.getLogger(__name__) diff --git a/test/test_voxels_utils.py b/test/test_voxels_utils.py index 82d7c92..dc56625 100644 --- a/test/test_voxels_utils.py +++ b/test/test_voxels_utils.py @@ -6,7 +6,7 @@ import numpy as np import pytest -from modelarrayio.utils.voxels import flattened_image +from modelarrayio.utils.nifti import flattened_image def _eye_affine(): From 9d3cc3bccdd197ef1ab8efecca83f93d72f656b2 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 09:55:21 -0400 Subject: [PATCH 07/20] Update test_fixels_utils.py --- test/test_fixels_utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/test/test_fixels_utils.py b/test/test_fixels_utils.py index 179f5bf..91bb270 100644 --- a/test/test_fixels_utils.py +++ b/test/test_fixels_utils.py @@ -6,7 +6,7 @@ import numpy as np import pytest -from modelarrayio.utils import fixels +from modelarrayio.utils import mif def _make_nifti2(shape=(2, 1, 1)) -> nb.Nifti2Image: @@ -15,14 +15,14 @@ def _make_nifti2(shape=(2, 1, 1)) -> nb.Nifti2Image: def test_nifti2_to_mif_raises_when_mrconvert_missing(tmp_path, monkeypatch) -> None: - monkeypatch.setattr(fixels, 'find_mrconvert', lambda: None) + monkeypatch.setattr(mif, 'find_mrconvert', lambda: None) with pytest.raises(FileNotFoundError, match='mrconvert'): - fixels.nifti2_to_mif(_make_nifti2(), tmp_path / 'out.mif') + mif.nifti2_to_mif(_make_nifti2(), tmp_path / 'out.mif') def test_mif_to_nifti2_raises_when_mrconvert_missing(monkeypatch) -> None: - monkeypatch.setattr(fixels, 'find_mrconvert', lambda: None) + monkeypatch.setattr(mif, 'find_mrconvert', lambda: None) with pytest.raises(FileNotFoundError, match='mrconvert'): - fixels.mif_to_nifti2('missing_input.mif') + mif.mif_to_nifti2('missing_input.mif') From 4c8e19d7ec5a8edbaa9f82b389e8628294a9f934 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 10:00:12 -0400 Subject: [PATCH 08/20] Update docstrings. --- src/modelarrayio/storage/h5_storage.py | 121 +++++++++++++++++++++ src/modelarrayio/storage/tiledb_storage.py | 110 ++++++++++++++++++- src/modelarrayio/utils/mif.py | 9 +- src/modelarrayio/utils/nifti.py | 18 ++- 4 files changed, 251 insertions(+), 7 deletions(-) diff --git a/src/modelarrayio/storage/h5_storage.py b/src/modelarrayio/storage/h5_storage.py index fedf62a..7c08e1c 100644 --- a/src/modelarrayio/storage/h5_storage.py +++ b/src/modelarrayio/storage/h5_storage.py @@ -16,10 +16,38 @@ def resolve_dtype(storage_dtype): + """Resolve a storage dtype to a supported NumPy floating type. + + Parameters + ---------- + storage_dtype : :obj:`str` + Storage dtype. + + Returns + ------- + :obj:`numpy.dtype` + Supported NumPy floating type. + """ return storage_utils.resolve_dtype(storage_dtype) def resolve_compression(compression, compression_level, shuffle): + """Resolve a compression method to a supported compression method. + + Parameters + ---------- + compression : :obj:`str` + Compression method. + compression_level : :obj:`int` + Compression level. + shuffle : :obj:`bool` + Whether to shuffle the data. + + Returns + ------- + :obj:`tuple` + Compression method, compression level, and whether to shuffle the data. + """ comp = ( None if compression is None or str(compression).lower() == 'none' @@ -39,6 +67,26 @@ def resolve_compression(compression, compression_level, shuffle): def compute_chunk_shape_full_subjects( num_subjects, num_items, item_chunk, target_chunk_mb, storage_np_dtype ): + """Compute a chunk shape for a full subject. + + Parameters + ---------- + num_subjects : :obj:`int` + Number of subjects. + num_items : :obj:`int` + Number of items. + item_chunk : :obj:`int` + Item chunk. + target_chunk_mb : :obj:`float` + Target chunk size in MB. + storage_np_dtype : :obj:`numpy.dtype` + Storage numpy dtype. + + Returns + ------- + :obj:`tuple` + Chunk shape. + """ chunk = storage_utils.compute_full_subject_chunk_shape( num_subjects=num_subjects, num_items=num_items, @@ -69,6 +117,36 @@ def create_scalar_matrix_dataset( chunk_voxels=0, target_chunk_mb=2.0, ): + """Create a scalar matrix dataset in an HDF5 file. + + Parameters + ---------- + h5file : :obj:`h5py.File` + HDF5 file. + dataset_path : :obj:`str` + Dataset path. + stacked_values : :obj:`numpy.ndarray` + Stacked values. + sources_list : :obj:`list` + Sources list. + storage_dtype : :obj:`str` + Storage dtype. + compression : :obj:`str` + Compression method. + compression_level : :obj:`int` + Compression level. + shuffle : :obj:`bool` + Whether to shuffle the data. + chunk_voxels : :obj:`int` + Chunk voxels. + target_chunk_mb : :obj:`float` + Target chunk size in MB. + + Returns + ------- + :obj:`h5py.Dataset` + Scalar matrix dataset. + """ storage_np_dtype = resolve_dtype(storage_dtype) comp, comp_opts, use_shuffle = resolve_compression(compression, compression_level, shuffle) @@ -118,6 +196,38 @@ def create_empty_scalar_matrix_dataset( target_chunk_mb=2.0, sources_list=None | pd.Series | list, ): + """Create an empty scalar matrix dataset in an HDF5 file. + + Parameters + ---------- + h5file : :obj:`h5py.File` + HDF5 file. + dataset_path : :obj:`str` + Dataset path. + num_subjects : :obj:`int` + Number of subjects. + num_items : :obj:`int` + Number of items. + storage_dtype : :obj:`str` + Storage dtype. + compression : :obj:`str` + Compression method. + compression_level : :obj:`int` + Compression level. + shuffle : :obj:`bool` + Whether to shuffle the data. + chunk_voxels : :obj:`int` + Chunk voxels. + target_chunk_mb : :obj:`float` + Target chunk size in MB. + sources_list : :obj:`list` + Sources list. + + Returns + ------- + :obj:`h5py.Dataset` + Empty scalar matrix dataset. + """ storage_np_dtype = resolve_dtype(storage_dtype) comp, comp_opts, use_shuffle = resolve_compression(compression, compression_level, shuffle) @@ -150,6 +260,17 @@ def create_empty_scalar_matrix_dataset( def write_column_names(h5_file: h5py.File, scalar: str, sources: pd.Series | list): + """Write column names to an HDF5 file. + + Parameters + ---------- + h5_file : :obj:`h5py.File` + HDF5 file. + scalar : :obj:`str` + Scalar name. + sources : :obj:`list` + Sources list. + """ values = np.array(storage_utils.normalize_column_names(sources), dtype=object) grp = h5_file.require_group(f'scalars/{scalar}') diff --git a/src/modelarrayio/storage/tiledb_storage.py b/src/modelarrayio/storage/tiledb_storage.py index df39374..e8d1f7f 100644 --- a/src/modelarrayio/storage/tiledb_storage.py +++ b/src/modelarrayio/storage/tiledb_storage.py @@ -16,6 +16,18 @@ def resolve_dtype(storage_dtype): + """Resolve a storage dtype to a supported NumPy floating type. + + Parameters + ---------- + storage_dtype : :obj:`str` + Storage dtype. + + Returns + ------- + :obj:`numpy.dtype` + Supported NumPy floating type. + """ return storage_utils.resolve_dtype(storage_dtype) @@ -46,6 +58,26 @@ def _build_filter_list(compression: str | None, compression_level: int | None, s def compute_tile_shape_full_subjects( num_subjects, num_items, item_tile, target_tile_mb, storage_np_dtype ): + """Compute a tile shape for a full subject. + + Parameters + ---------- + num_subjects : :obj:`int` + Number of subjects. + num_items : :obj:`int` + Number of items. + item_tile : :obj:`int` + Item tile. + target_tile_mb : :obj:`float` + Target tile size in MB. + storage_np_dtype : :obj:`numpy.dtype` + Storage numpy dtype. + + Returns + ------- + :obj:`tuple` + Tile shape. + """ tile = storage_utils.compute_full_subject_chunk_shape( num_subjects=num_subjects, num_items=num_items, @@ -82,6 +114,36 @@ def create_scalar_matrix_array( tile_voxels=0, target_tile_mb=2.0, ): + """Create a scalar matrix array in a TileDB directory. + + Parameters + ---------- + base_uri : :obj:`str` + Base URI. + dataset_path : :obj:`str` + Dataset path. + stacked_values : :obj:`numpy.ndarray` + Stacked values. + sources_list : :obj:`list` + Sources list. + storage_dtype : :obj:`str` + Storage dtype. + compression : :obj:`str` + Compression method. + compression_level : :obj:`int` + Compression level. + shuffle : :obj:`bool` + Whether to shuffle the data. + tile_voxels : :obj:`int` + Tile voxels. + target_tile_mb : :obj:`float` + Target tile size in MB. + + Returns + ------- + :obj:`str` + URI of the created array. + """ storage_np_dtype = resolve_dtype(storage_dtype) if stacked_values.dtype != storage_np_dtype: stacked_values = stacked_values.astype(storage_np_dtype) @@ -144,6 +206,38 @@ def create_empty_scalar_matrix_array( target_tile_mb=2.0, sources_list: Sequence[str] | None = None, ): + """Create an empty scalar matrix array in a TileDB directory. + + Parameters + ---------- + base_uri : :obj:`str` + Base URI. + dataset_path : :obj:`str` + Dataset path. + num_subjects : :obj:`int` + Number of subjects. + num_items : :obj:`int` + Number of items. + storage_dtype : :obj:`str` + Storage dtype. + compression : :obj:`str` + Compression method. + compression_level : :obj:`int` + Compression level. + shuffle : :obj:`bool` + Whether to shuffle the data. + tile_voxels : :obj:`int` + Tile voxels. + target_tile_mb : :obj:`float` + Target tile size in MB. + sources_list : :obj:`list` + Sources list. + + Returns + ------- + :obj:`str` + URI of the created array. + """ storage_np_dtype = resolve_dtype(storage_dtype) tile_shape = compute_tile_shape_full_subjects( num_subjects, num_items, tile_voxels, target_tile_mb, storage_np_dtype @@ -185,8 +279,7 @@ def create_empty_scalar_matrix_array( def write_rows_in_column_stripes(uri: str, rows: Sequence[np.ndarray]): - """ - Fill a 2D TileDB dense array by buffering column-aligned stripes to minimize + """Fill a 2D TileDB dense array by buffering column-aligned stripes to minimize tile writes, using about one tile's worth of memory. Parameters @@ -226,9 +319,16 @@ def write_rows_in_column_stripes(uri: str, rows: Sequence[np.ndarray]): def write_column_names(base_uri: str, scalar: str, sources: Sequence[str]): - """ - Store column names as a 1D dense TileDB array for the given scalar. - This mirrors the HDF5 dataset approach and scales to large cohorts. + """Store column names as a 1D dense TileDB array for the given scalar. + + Parameters + ---------- + base_uri : :obj:`str` + Base URI. + scalar : :obj:`str` + Scalar name. + sources : :obj:`list` + Sources list. """ sources = storage_utils.normalize_column_names(sources) uri = os.path.join(base_uri, 'scalars', scalar, 'column_names') diff --git a/src/modelarrayio/utils/mif.py b/src/modelarrayio/utils/mif.py index 43755b3..aa9d913 100644 --- a/src/modelarrayio/utils/mif.py +++ b/src/modelarrayio/utils/mif.py @@ -1,4 +1,4 @@ -"""Utility functions for fixel-wise data.""" +"""Utility functions for MIF data.""" import shutil import subprocess @@ -14,6 +14,13 @@ def find_mrconvert(): + """Find the mrconvert executable on the system. + + Returns + ------- + :obj:`str` + Path to the mrconvert executable. + """ return shutil.which('mrconvert') diff --git a/src/modelarrayio/utils/nifti.py b/src/modelarrayio/utils/nifti.py index 663eaaa..9e18ff7 100644 --- a/src/modelarrayio/utils/nifti.py +++ b/src/modelarrayio/utils/nifti.py @@ -1,4 +1,4 @@ -"""Utility functions for voxel-wise data.""" +"""Utility functions for NIfTI data.""" from __future__ import annotations @@ -82,6 +82,22 @@ def _worker(job): def flattened_image(scalar_image, scalar_mask, group_mask_matrix): + """Flatten a scalar image to a 1-D array. + + Parameters + ---------- + scalar_image : :obj:`nibabel.Nifti1Image` + Scalar image. + scalar_mask : :obj:`nibabel.Nifti1Image` + Scalar mask. + group_mask_matrix : :obj:`numpy.ndarray` + Group mask matrix. + + Returns + ------- + :obj:`numpy.ndarray` + Flattened scalar image. + """ scalar_mask_img = ( scalar_mask if hasattr(scalar_mask, 'get_fdata') else nb.load(Path(scalar_mask)) ) From 275a2b05371470c7affee2037bcd52284f0a85e5 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 12:14:42 -0400 Subject: [PATCH 09/20] Fix bugs. np.unicode_ isn't available in numpy 2. --- src/modelarrayio/cli/parser_utils.py | 1 + src/modelarrayio/storage/tiledb_storage.py | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/src/modelarrayio/cli/parser_utils.py b/src/modelarrayio/cli/parser_utils.py index 05865fe..a083699 100644 --- a/src/modelarrayio/cli/parser_utils.py +++ b/src/modelarrayio/cli/parser_utils.py @@ -20,6 +20,7 @@ def add_to_modelarray_args(parser, default_output='output.h5'): 'for the tiledb backend, path to a .tdb directory.' ), default=default_output, + type=Path, ) parser.add_argument( '--scalar-columns', diff --git a/src/modelarrayio/storage/tiledb_storage.py b/src/modelarrayio/storage/tiledb_storage.py index e8d1f7f..33be392 100644 --- a/src/modelarrayio/storage/tiledb_storage.py +++ b/src/modelarrayio/storage/tiledb_storage.py @@ -339,7 +339,7 @@ def write_column_names(base_uri: str, scalar: str, sources: Sequence[str]): name='idx', domain=(0, max(n - 1, 0)), tile=max(1, min(n, 1024)), dtype=np.int64 ) dom = tiledb.Domain(dim_idx) - attr_values = tiledb.Attr(name='values', dtype=np.unicode_) + attr_values = tiledb.Attr(name='values', dtype=np.str_) schema = tiledb.ArraySchema(domain=dom, attrs=[attr_values], sparse=False) if tiledb.object_type(uri): From bb99a7b7bd682741c24c6e48bedfeb5f22fa81b5 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 12:17:14 -0400 Subject: [PATCH 10/20] Update test_parser_utils.py --- test/test_parser_utils.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_parser_utils.py b/test/test_parser_utils.py index 4aed56a..278667a 100644 --- a/test/test_parser_utils.py +++ b/test/test_parser_utils.py @@ -3,6 +3,7 @@ from __future__ import annotations import argparse +from pathlib import Path from modelarrayio.cli import parser_utils @@ -66,7 +67,7 @@ def test_output_hdf5_default_name_override(tmp_path_factory) -> None: p = argparse.ArgumentParser() parser_utils.add_to_modelarray_args(p, default_output='custom.h5') args = p.parse_args(['--cohort-file', str(cohort_file)]) - assert args.output == 'custom.h5' + assert args.output == Path('custom.h5') def test_tiledb_args_group(tmp_path_factory) -> None: @@ -77,7 +78,7 @@ def test_tiledb_args_group(tmp_path_factory) -> None: parser_utils.add_to_modelarray_args(p, default_output='arrays.tdb') args = p.parse_args(['--cohort-file', str(cohort_file), '--backend', 'tiledb']) assert args.output == 'arrays.tdb' - assert args.backend == 'tiledb' + assert args.backend == Path('tiledb') assert args.workers == 0 assert args.s3_workers == 1 assert args.log_level == 'INFO' From cecfb41d6e00273cffbf51aa485c98be436635d3 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 12:18:39 -0400 Subject: [PATCH 11/20] Update test_parser_utils.py --- test/test_parser_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parser_utils.py b/test/test_parser_utils.py index 278667a..7b22509 100644 --- a/test/test_parser_utils.py +++ b/test/test_parser_utils.py @@ -77,7 +77,7 @@ def test_tiledb_args_group(tmp_path_factory) -> None: p = argparse.ArgumentParser() parser_utils.add_to_modelarray_args(p, default_output='arrays.tdb') args = p.parse_args(['--cohort-file', str(cohort_file), '--backend', 'tiledb']) - assert args.output == 'arrays.tdb' + assert args.output == Path('arrays.tdb') assert args.backend == Path('tiledb') assert args.workers == 0 assert args.s3_workers == 1 From f1cdfc07d14b7eb980317d08ca8f75e4842dd257 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Fri, 27 Mar 2026 12:19:59 -0400 Subject: [PATCH 12/20] Update test_parser_utils.py --- test/test_parser_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parser_utils.py b/test/test_parser_utils.py index 7b22509..6fa13d4 100644 --- a/test/test_parser_utils.py +++ b/test/test_parser_utils.py @@ -78,7 +78,7 @@ def test_tiledb_args_group(tmp_path_factory) -> None: parser_utils.add_to_modelarray_args(p, default_output='arrays.tdb') args = p.parse_args(['--cohort-file', str(cohort_file), '--backend', 'tiledb']) assert args.output == Path('arrays.tdb') - assert args.backend == Path('tiledb') + assert args.backend == 'tiledb' assert args.workers == 0 assert args.s3_workers == 1 assert args.log_level == 'INFO' From 297a224033a5e22694158bebd6617fe7239fbc7f Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 11:07:06 -0400 Subject: [PATCH 13/20] Split scalar outputs. --- src/modelarrayio/cli/cifti_to_h5.py | 40 +++++++++++++-- src/modelarrayio/cli/mif_to_h5.py | 34 +++++++++++-- src/modelarrayio/cli/nifti_to_h5.py | 33 ++++++++++-- src/modelarrayio/cli/utils.py | 7 +++ test/test_cifti_cli.py | 64 +++++++++++++++++++++++ test/test_mif_cli.py | 77 ++++++++++++++++++++++++++++ test/test_voxels_cli.py | 79 +++++++++++++++++++++++++++++ 7 files changed, 323 insertions(+), 11 deletions(-) create mode 100644 test/test_mif_cli.py diff --git a/src/modelarrayio/cli/cifti_to_h5.py b/src/modelarrayio/cli/cifti_to_h5.py index d316e08..5b8ccd0 100644 --- a/src/modelarrayio/cli/cifti_to_h5.py +++ b/src/modelarrayio/cli/cifti_to_h5.py @@ -82,8 +82,39 @@ def cifti_to_h5( scalar_sources = build_scalar_sources(cohort_long) if not scalar_sources: raise ValueError('Unable to derive scalar sources from cohort file.') + scalar_names = list(scalar_sources.keys()) + split_scalar_outputs = bool(scalar_columns) if backend == 'hdf5': + if split_scalar_outputs: + scalars, last_brain_names = load_cohort_cifti(cohort_long, s3_workers) + greyordinate_table, structure_names = brain_names_to_dataframe(last_brain_names) + outputs: list[Path] = [] + for scalar_name in scalar_names: + scalar_output = cli_utils.prepare_output_parent( + cli_utils.prefixed_output_path(output, scalar_name) + ) + with h5py.File(scalar_output, 'w') as h5_file: + cli_utils.write_table_dataset( + h5_file, + 'greyordinates', + greyordinate_table, + extra_attrs={'structure_names': structure_names}, + ) + cli_utils.write_hdf5_scalar_matrices( + h5_file, + {scalar_name: scalars[scalar_name]}, + {scalar_name: scalar_sources[scalar_name]}, + storage_dtype=storage_dtype, + compression=compression, + compression_level=compression_level, + shuffle=shuffle, + chunk_voxels=chunk_voxels, + target_chunk_mb=target_chunk_mb, + ) + outputs.append(scalar_output) + return int(not all(path.exists() for path in outputs)) + scalars, last_brain_names = load_cohort_cifti(cohort_long, s3_workers) greyordinate_table, structure_names = brain_names_to_dataframe(last_brain_names) output = cli_utils.prepare_output_parent(output) @@ -107,7 +138,6 @@ def cifti_to_h5( ) return int(not output.exists()) - output.mkdir(parents=True, exist_ok=True) if not scalar_sources: return 0 @@ -124,8 +154,13 @@ def _process_scalar_job(scalar_name, source_files): rows.append(cifti_data) if rows: + scalar_output = ( + cli_utils.prefixed_output_path(output, scalar_name) + if split_scalar_outputs + else output + ) cli_utils.write_tiledb_scalar_matrices( - output, + scalar_output, {scalar_name: rows}, {scalar_name: source_files}, storage_dtype=storage_dtype, @@ -138,7 +173,6 @@ def _process_scalar_job(scalar_name, source_files): ) return scalar_name - scalar_names = list(scalar_sources.keys()) worker_count = workers if isinstance(workers, int) and workers > 0 else None if worker_count is None: cpu_count = os.cpu_count() or 1 diff --git a/src/modelarrayio/cli/mif_to_h5.py b/src/modelarrayio/cli/mif_to_h5.py index 326cb87..b1bff4d 100644 --- a/src/modelarrayio/cli/mif_to_h5.py +++ b/src/modelarrayio/cli/mif_to_h5.py @@ -89,7 +89,33 @@ def mif_to_h5( if not sources_lists: raise ValueError('Unable to derive scalar sources from cohort file.') + scalar_names = list(sources_lists.keys()) + split_scalar_outputs = bool(scalar_columns) + if backend == 'hdf5': + if split_scalar_outputs: + outputs: list[Path] = [] + for scalar_name in scalar_names: + scalar_output = cli_utils.prepare_output_parent( + cli_utils.prefixed_output_path(output, scalar_name) + ) + with h5py.File(scalar_output, 'w') as h5_file: + cli_utils.write_table_dataset(h5_file, 'fixels', fixel_table) + cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table) + cli_utils.write_hdf5_scalar_matrices( + h5_file, + {scalar_name: scalars[scalar_name]}, + {scalar_name: sources_lists[scalar_name]}, + storage_dtype=storage_dtype, + compression=compression, + compression_level=compression_level, + shuffle=shuffle, + chunk_voxels=chunk_voxels, + target_chunk_mb=target_chunk_mb, + ) + outputs.append(scalar_output) + return int(not all(path.exists() for path in outputs)) + output = cli_utils.prepare_output_parent(output) with h5py.File(output, 'w') as h5_file: cli_utils.write_table_dataset(h5_file, 'fixels', fixel_table) @@ -107,9 +133,6 @@ def mif_to_h5( ) return int(not output.exists()) - output.mkdir(parents=True, exist_ok=True) - - scalar_names = list(sources_lists.keys()) worker_count = workers if isinstance(workers, int) and workers > 0 else None if worker_count is None: cpu_count = os.cpu_count() or 1 @@ -118,8 +141,11 @@ def mif_to_h5( worker_count = min(len(scalar_names), worker_count) def _write_scalar_job(scalar_name): + scalar_output = ( + cli_utils.prefixed_output_path(output, scalar_name) if split_scalar_outputs else output + ) cli_utils.write_tiledb_scalar_matrices( - output, + scalar_output, {scalar_name: scalars[scalar_name]}, {scalar_name: sources_lists[scalar_name]}, storage_dtype=storage_dtype, diff --git a/src/modelarrayio/cli/nifti_to_h5.py b/src/modelarrayio/cli/nifti_to_h5.py index 34a12f5..00c4f4d 100644 --- a/src/modelarrayio/cli/nifti_to_h5.py +++ b/src/modelarrayio/cli/nifti_to_h5.py @@ -92,7 +92,32 @@ def nifti_to_h5( if not sources_lists: raise ValueError('Unable to derive scalar sources from cohort file.') + scalar_names = list(sources_lists.keys()) + split_scalar_outputs = bool(scalar_columns) + if backend == 'hdf5': + if split_scalar_outputs: + outputs: list[Path] = [] + for scalar_name in scalar_names: + scalar_output = cli_utils.prepare_output_parent( + cli_utils.prefixed_output_path(output, scalar_name) + ) + with h5py.File(scalar_output, 'w') as h5_file: + cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table) + cli_utils.write_hdf5_scalar_matrices( + h5_file, + {scalar_name: scalars[scalar_name]}, + {scalar_name: sources_lists[scalar_name]}, + storage_dtype=storage_dtype, + compression=compression, + compression_level=compression_level, + shuffle=shuffle, + chunk_voxels=chunk_voxels, + target_chunk_mb=target_chunk_mb, + ) + outputs.append(scalar_output) + return int(not all(path.exists() for path in outputs)) + output = cli_utils.prepare_output_parent(output) with h5py.File(output, 'w') as h5_file: cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table) @@ -109,9 +134,6 @@ def nifti_to_h5( ) return int(not output.exists()) - output.mkdir(parents=True, exist_ok=True) - - scalar_names = list(sources_lists.keys()) worker_count = workers if isinstance(workers, int) and workers > 0 else None if worker_count is None: cpu_count = os.cpu_count() or 1 @@ -120,8 +142,11 @@ def nifti_to_h5( worker_count = min(len(scalar_names), worker_count) def _write_scalar_job(scalar_name): + scalar_output = ( + cli_utils.prefixed_output_path(output, scalar_name) if split_scalar_outputs else output + ) cli_utils.write_tiledb_scalar_matrices( - output, + scalar_output, {scalar_name: scalars[scalar_name]}, {scalar_name: sources_lists[scalar_name]}, storage_dtype=storage_dtype, diff --git a/src/modelarrayio/cli/utils.py b/src/modelarrayio/cli/utils.py index 36918b4..e1614b8 100644 --- a/src/modelarrayio/cli/utils.py +++ b/src/modelarrayio/cli/utils.py @@ -37,6 +37,13 @@ def prepare_output_parent(output_file: str | Path) -> Path: return output_path +def prefixed_output_path(output_path: str | Path, prefix: str) -> Path: + """Return output path with a sanitized prefix added to its filename.""" + path = Path(output_path) + safe_prefix = sanitize_result_name(prefix) + return path.with_name(f'{safe_prefix}_{path.name}') + + def write_table_dataset( h5_file: h5py.File, dataset_name: str, diff --git a/test/test_cifti_cli.py b/test/test_cifti_cli.py index 2fb3fba..557b2d3 100644 --- a/test/test_cifti_cli.py +++ b/test/test_cifti_cli.py @@ -107,3 +107,67 @@ def test_concifti_cli_creates_expected_hdf5(tmp_path, monkeypatch): # Spot-check a couple values assert np.isclose(float(dset[0, 0]), 0.0) assert np.isclose(float(dset[1, 0]), 1.0) + + +def test_cifti_to_h5_scalar_columns_writes_prefixed_outputs(tmp_path, monkeypatch): + vol_shape = (2, 2, 2) + mask = np.zeros(vol_shape, dtype=bool) + mask[(0, 0, 0)] = True + mask[(1, 1, 1)] = True + n_grayordinates = int(mask.sum()) + + rows = [] + for sidx in range(2): + alpha_vals = np.arange(n_grayordinates, dtype=np.float32) + sidx + beta_vals = np.arange(n_grayordinates, dtype=np.float32) + 10 + sidx + + alpha_img = _make_synthetic_cifti_dscalar(mask, alpha_vals) + beta_img = _make_synthetic_cifti_dscalar(mask, beta_vals) + + alpha_path = tmp_path / f'sub-{sidx + 1}_alpha.dscalar.nii' + beta_path = tmp_path / f'sub-{sidx + 1}_beta.dscalar.nii' + alpha_img.to_filename(alpha_path) + beta_img.to_filename(beta_path) + + rows.append( + {'subject_id': f'sub-{sidx + 1}', 'alpha': alpha_path.name, 'beta': beta_path.name} + ) + + cohort_csv = tmp_path / 'cohort_cifti_wide.csv' + with cohort_csv.open('w', newline='') as f: + writer = csv.DictWriter(f, fieldnames=['subject_id', 'alpha', 'beta']) + writer.writeheader() + writer.writerows(rows) + + out_h5 = tmp_path / 'greyordinatearray.h5' + alpha_out = tmp_path / 'alpha_greyordinatearray.h5' + beta_out = tmp_path / 'beta_greyordinatearray.h5' + + monkeypatch.chdir(tmp_path) + assert ( + modelarrayio_main( + [ + 'cifti-to-h5', + '--cohort-file', + str(cohort_csv), + '--scalar-columns', + 'alpha', + 'beta', + '--output', + str(out_h5), + ] + ) + == 0 + ) + + assert alpha_out.exists() + assert beta_out.exists() + assert not out_h5.exists() + + with h5py.File(alpha_out, 'r') as h5: + assert 'greyordinates' in h5 + assert sorted(h5['scalars'].keys()) == ['alpha'] + + with h5py.File(beta_out, 'r') as h5: + assert 'greyordinates' in h5 + assert sorted(h5['scalars'].keys()) == ['beta'] diff --git a/test/test_mif_cli.py b/test/test_mif_cli.py new file mode 100644 index 0000000..1bf159a --- /dev/null +++ b/test/test_mif_cli.py @@ -0,0 +1,77 @@ +import csv + +import h5py +import numpy as np +import pandas as pd + +import modelarrayio.cli.mif_to_h5 as mif_cli +from modelarrayio.cli.main import main as modelarrayio_main + + +def test_mif_to_h5_scalar_columns_writes_prefixed_outputs(tmp_path, monkeypatch): + index_file = tmp_path / 'index.nii.gz' + directions_file = tmp_path / 'directions.nii.gz' + index_file.write_bytes(b'index') + directions_file.write_bytes(b'directions') + + cohort_csv = tmp_path / 'cohort_mif_wide.csv' + with cohort_csv.open('w', newline='') as f: + writer = csv.DictWriter(f, fieldnames=['subject_id', 'alpha', 'beta']) + writer.writeheader() + writer.writerow({'subject_id': 'sub-1', 'alpha': 'a1.mif', 'beta': 'b1.mif'}) + writer.writerow({'subject_id': 'sub-2', 'alpha': 'a2.mif', 'beta': 'b2.mif'}) + + fixel_table = pd.DataFrame({'fixel_id': [0, 1], 'x': [0.0, 1.0], 'y': [0.0, 1.0], 'z': [0.0, 1.0]}) + voxel_table = pd.DataFrame({'voxel_id': [0], 'i': [0], 'j': [0], 'k': [0]}) + + def fake_gather_fixels(_index_file, _directions_file): + return fixel_table, voxel_table + + def fake_load_cohort_mif(cohort_long, _s3_workers): + scalars = {} + sources = {} + for scalar_name, group in cohort_long.groupby('scalar_name'): + n_rows = group.shape[0] + scalars[scalar_name] = [np.array([1.0, 2.0], dtype=np.float32) for _ in range(n_rows)] + sources[scalar_name] = group['source_file'].tolist() + return scalars, sources + + monkeypatch.setattr(mif_cli, 'gather_fixels', fake_gather_fixels) + monkeypatch.setattr(mif_cli, 'load_cohort_mif', fake_load_cohort_mif) + + out_h5 = tmp_path / 'fixelarray.h5' + alpha_out = tmp_path / 'alpha_fixelarray.h5' + beta_out = tmp_path / 'beta_fixelarray.h5' + assert ( + modelarrayio_main( + [ + 'mif-to-h5', + '--index-file', + str(index_file), + '--directions-file', + str(directions_file), + '--cohort-file', + str(cohort_csv), + '--scalar-columns', + 'alpha', + 'beta', + '--output', + str(out_h5), + ] + ) + == 0 + ) + + assert alpha_out.exists() + assert beta_out.exists() + assert not out_h5.exists() + + with h5py.File(alpha_out, 'r') as h5: + assert 'fixels' in h5 + assert 'voxels' in h5 + assert sorted(h5['scalars'].keys()) == ['alpha'] + + with h5py.File(beta_out, 'r') as h5: + assert 'fixels' in h5 + assert 'voxels' in h5 + assert sorted(h5['scalars'].keys()) == ['beta'] diff --git a/test/test_voxels_cli.py b/test/test_voxels_cli.py index c75acf7..1357fc6 100644 --- a/test/test_voxels_cli.py +++ b/test/test_voxels_cli.py @@ -204,3 +204,82 @@ def test_h5_to_nifti_cli_writes_results_with_dataset_column_names(tmp_path): assert effect_data[coord] == pytest.approx([0.1, 0.2, 0.3][idx]) assert pvalue_data[coord] == pytest.approx([0.9, 0.8, 0.7][idx]) assert inv_pvalue_data[coord] == pytest.approx([0.1, 0.2, 0.3][idx]) + + +def test_nifti_to_h5_scalar_columns_writes_prefixed_outputs(tmp_path, monkeypatch): + shape = (3, 3, 3) + group_mask = np.zeros(shape, dtype=bool) + true_coords = [(0, 0, 1), (1, 1, 1), (2, 2, 0)] + for i, j, k in true_coords: + group_mask[i, j, k] = True + + group_mask_file = tmp_path / 'group_mask.nii.gz' + _make_nifti(group_mask.astype(np.uint8)).to_filename(group_mask_file) + + rows = [] + for sidx in range(2): + subj_mask_file = tmp_path / f'sub-{sidx + 1}_mask.nii.gz' + _make_nifti(group_mask.astype(np.uint8)).to_filename(subj_mask_file) + + alpha_data = np.zeros(shape, dtype=np.float32) + beta_data = np.zeros(shape, dtype=np.float32) + for i, j, k in true_coords: + alpha_data[i, j, k] = 10.0 + sidx + beta_data[i, j, k] = 20.0 + sidx + + alpha_file = tmp_path / f'sub-{sidx + 1}_alpha.nii.gz' + beta_file = tmp_path / f'sub-{sidx + 1}_beta.nii.gz' + _make_nifti(alpha_data).to_filename(alpha_file) + _make_nifti(beta_data).to_filename(beta_file) + + rows.append( + { + 'subject_id': f'sub-{sidx + 1}', + 'alpha': alpha_file.name, + 'beta': beta_file.name, + 'source_mask_file': subj_mask_file.name, + } + ) + + cohort_csv = tmp_path / 'cohort_wide.csv' + with cohort_csv.open('w', newline='') as f: + writer = csv.DictWriter( + f, fieldnames=['subject_id', 'alpha', 'beta', 'source_mask_file'] + ) + writer.writeheader() + writer.writerows(rows) + + out_h5 = tmp_path / 'voxelarray.h5' + alpha_out = tmp_path / 'alpha_voxelarray.h5' + beta_out = tmp_path / 'beta_voxelarray.h5' + + monkeypatch.chdir(tmp_path) + assert ( + modelarrayio_main( + [ + 'nifti-to-h5', + '--group-mask-file', + str(group_mask_file), + '--cohort-file', + str(cohort_csv), + '--scalar-columns', + 'alpha', + 'beta', + '--output', + str(out_h5), + ] + ) + == 0 + ) + + assert alpha_out.exists() + assert beta_out.exists() + assert not out_h5.exists() + + with h5py.File(alpha_out, 'r') as h5: + assert 'voxels' in h5 + assert sorted(h5['scalars'].keys()) == ['alpha'] + + with h5py.File(beta_out, 'r') as h5: + assert 'voxels' in h5 + assert sorted(h5['scalars'].keys()) == ['beta'] From c7f862bf9e7e3bee62a4d8b161702d1c27317d84 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 11:09:23 -0400 Subject: [PATCH 14/20] Run ruff. --- test/test_mif_cli.py | 4 +++- test/test_voxels_cli.py | 4 +--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/test/test_mif_cli.py b/test/test_mif_cli.py index 1bf159a..59d8603 100644 --- a/test/test_mif_cli.py +++ b/test/test_mif_cli.py @@ -21,7 +21,9 @@ def test_mif_to_h5_scalar_columns_writes_prefixed_outputs(tmp_path, monkeypatch) writer.writerow({'subject_id': 'sub-1', 'alpha': 'a1.mif', 'beta': 'b1.mif'}) writer.writerow({'subject_id': 'sub-2', 'alpha': 'a2.mif', 'beta': 'b2.mif'}) - fixel_table = pd.DataFrame({'fixel_id': [0, 1], 'x': [0.0, 1.0], 'y': [0.0, 1.0], 'z': [0.0, 1.0]}) + fixel_table = pd.DataFrame( + {'fixel_id': [0, 1], 'x': [0.0, 1.0], 'y': [0.0, 1.0], 'z': [0.0, 1.0]} + ) voxel_table = pd.DataFrame({'voxel_id': [0], 'i': [0], 'j': [0], 'k': [0]}) def fake_gather_fixels(_index_file, _directions_file): diff --git a/test/test_voxels_cli.py b/test/test_voxels_cli.py index 1357fc6..729d82f 100644 --- a/test/test_voxels_cli.py +++ b/test/test_voxels_cli.py @@ -243,9 +243,7 @@ def test_nifti_to_h5_scalar_columns_writes_prefixed_outputs(tmp_path, monkeypatc cohort_csv = tmp_path / 'cohort_wide.csv' with cohort_csv.open('w', newline='') as f: - writer = csv.DictWriter( - f, fieldnames=['subject_id', 'alpha', 'beta', 'source_mask_file'] - ) + writer = csv.DictWriter(f, fieldnames=['subject_id', 'alpha', 'beta', 'source_mask_file']) writer.writeheader() writer.writerows(rows) From c541c2bfb4ed578bfcfad8c586cabe2396ee319b Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 11:20:13 -0400 Subject: [PATCH 15/20] Add outputs page. --- docs/index.rst | 1 + docs/outputs.rst | 175 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 176 insertions(+) create mode 100644 docs/outputs.rst diff --git a/docs/index.rst b/docs/index.rst index 1dfedb2..2d62c0e 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -14,4 +14,5 @@ installation auto_examples/index usage + outputs api diff --git a/docs/outputs.rst b/docs/outputs.rst new file mode 100644 index 0000000..63324c8 --- /dev/null +++ b/docs/outputs.rst @@ -0,0 +1,175 @@ +####### +Outputs +####### + +This page describes what each CLI command writes, how files are named, and what data +is stored inside each output artifact. + + +***************** +Commands Overview +***************** + +The commands fall into two groups: + +- ``*-to-h5`` commands: convert input neuroimaging data into either: + - one or more HDF5 files (``--backend hdf5``), or + - one or more TileDB directories (``--backend tiledb``). +- ``h5-to-*`` commands: convert analysis results stored in an HDF5 file into image files. + + +********************* +nifti-to-h5 (volumes) +********************* + +Default output name (HDF5 backend): ``voxelarray.h5``. + +HDF5 output contents: + +- ``voxels`` dataset: + - transposed voxel table with rows for ``voxel_id``, ``i``, ``j``, ``k``. + - attribute ``column_names = ['voxel_id', 'i', 'j', 'k']``. +- Per scalar: + - ``scalars//values`` with shape ``(n_subjects, n_voxels)``. + - ``scalars//column_names`` listing source file names. + +TileDB output contents: + +- Per scalar dense array at ``scalars//values`` with shape + ``(n_subjects, n_voxels)``. +- Column names are stored in array metadata (``column_names``). + +When ``--scalar-columns`` is provided: + +- Output is split by scalar column name. +- Example: ``--scalar-columns alpha beta --output voxelarray.h5`` writes: + - ``alpha_voxelarray.h5`` + - ``beta_voxelarray.h5`` +- The same prefix rule also applies to TileDB output paths. + + +******************** +cifti-to-h5 (CIFTI) +******************** + +Default output name (HDF5 backend): ``greyordinatearray.h5``. + +HDF5 output contents: + +- ``greyordinates`` dataset: + - transposed table with rows for ``vertex_id`` and ``structure_id``. + - attribute ``column_names = ['vertex_id', 'structure_id']``. + - attribute ``structure_names`` listing CIFTI brain structures. +- Per scalar: + - ``scalars//values`` with shape ``(n_subjects, n_greyordinates)``. + - ``scalars//column_names`` listing source file names. + +TileDB output contents: + +- Per scalar dense array at ``scalars//values`` with shape + ``(n_subjects, n_greyordinates)``. +- Column names metadata is written on each scalar matrix. +- An explicit TileDB array is also written at ``scalars//column_names``. + +When ``--scalar-columns`` is provided: + +- Output is split by scalar column name. +- Example: ``--scalar-columns alpha beta --output greyordinatearray.h5`` writes: + - ``alpha_greyordinatearray.h5`` + - ``beta_greyordinatearray.h5`` +- The same prefix rule also applies to TileDB output paths. + + +****************** +mif-to-h5 (fixels) +****************** + +Default output name (HDF5 backend): ``fixelarray.h5``. + +HDF5 output contents: + +- ``fixels`` dataset: + - transposed fixel table (``fixel_id``, coordinates/directions metadata from input fixel DB). + - attribute ``column_names`` containing table column names. +- ``voxels`` dataset: + - transposed voxel table with ``voxel_id``, ``i``, ``j``, ``k``. + - attribute ``column_names`` containing table column names. +- Per scalar: + - ``scalars//values`` with shape ``(n_subjects, n_fixels)``. + - ``scalars//column_names`` listing source file names. + +TileDB output contents: + +- Per scalar dense array at ``scalars//values`` with shape + ``(n_subjects, n_fixels)``. +- Column names are stored in array metadata (``column_names``). + +When ``--scalar-columns`` is provided: + +- Output is split by scalar column name. +- Example: ``--scalar-columns alpha beta --output fixelarray.h5`` writes: + - ``alpha_fixelarray.h5`` + - ``beta_fixelarray.h5`` +- The same prefix rule also applies to TileDB output paths. + + +*********************************** +h5-to-* commands (result exporters) +*********************************** + +These commands read statistical results from: + +- ``results//results_matrix`` (shape: ``(n_results, n_elements)``). + +Result names are read in this order: + +- ``results_matrix.attrs['colnames']`` (if present), +- ``results//column_names`` dataset, +- ``results//results_matrix/column_names`` dataset, +- fallback names: ``component001``, ``component002``, ... + +Any spaces or ``/`` in result names are replaced with ``_`` in filenames. + + +h5-to-nifti +=========== + +Writes one file per result to ``--output-dir``: + +- ``_`` (default extension ``.nii.gz``). +- If a result name contains ``p.value``, an additional file is written: + ``_``, + containing ``1 - p.value``. + +Each output volume uses ``--group-mask-file`` to map vectorized results back into 3D space. + + +h5-to-cifti +=========== + +Writes one CIFTI dscalar file per result to ``--output-dir``: + +- ``_.dscalar.nii``. +- If a result name contains ``p.value``, also writes the ``1 - p.value`` companion file + with ``1m.p.value`` in its name. + +The header is taken from ``--example-cifti`` (or from the first cohort ``source_file`` if +``--cohort-file`` is used instead). + + +h5-to-mif +========= + +Writes one MIF file per result to ``--output-dir``: + +- ``_.mif``. +- If a result name contains ``p.value``, also writes the ``1 - p.value`` companion file + with ``1m.p.value`` in its name. + +Also copies these files into ``--output-dir``: + +- ``--index-file`` +- ``--directions-file`` + +The output MIF geometry/header template is taken from ``--example-mif`` (or from the first +cohort ``source_file`` if ``--cohort-file`` is used instead). From 6335c11c9a940f49f1aa255749600a2feb89a1b0 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 12:39:50 -0400 Subject: [PATCH 16/20] Update docs/outputs.rst Co-authored-by: Matt Cieslak --- docs/outputs.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/outputs.rst b/docs/outputs.rst index 63324c8..092bdd8 100644 --- a/docs/outputs.rst +++ b/docs/outputs.rst @@ -48,9 +48,9 @@ When ``--scalar-columns`` is provided: - The same prefix rule also applies to TileDB output paths. -******************** +******************* cifti-to-h5 (CIFTI) -******************** +******************* Default output name (HDF5 backend): ``greyordinatearray.h5``. From 571ea88e78fd51bb591c6fff25ae349afc1550a7 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 12:51:13 -0400 Subject: [PATCH 17/20] Rename variables. num_subjects --> n_files num_items --> n_elements --- src/modelarrayio/cli/utils.py | 20 ++++----- src/modelarrayio/storage/h5_storage.py | 48 ++++++++++---------- src/modelarrayio/storage/tiledb_storage.py | 52 +++++++++++----------- src/modelarrayio/storage/utils.py | 18 ++++---- test/test_cifti_cli.py | 6 +-- test/test_h5_storage.py | 28 +++++------- test/test_storage_utils.py | 4 +- test/test_voxels_cli.py | 4 +- test/test_voxels_s3.py | 6 +-- 9 files changed, 91 insertions(+), 95 deletions(-) diff --git a/src/modelarrayio/cli/utils.py b/src/modelarrayio/cli/utils.py index e1614b8..85ee3a8 100644 --- a/src/modelarrayio/cli/utils.py +++ b/src/modelarrayio/cli/utils.py @@ -73,15 +73,15 @@ def write_hdf5_scalar_matrices( ) -> None: """Write per-scalar matrices into an open HDF5 file.""" for scalar_name, rows in scalars.items(): - num_subjects = len(rows) - if num_subjects == 0: + n_files = len(rows) + if n_files == 0: continue - num_items = rows[0].shape[0] + n_elements = rows[0].shape[0] dataset = h5_storage.create_empty_scalar_matrix_dataset( h5_file, f'scalars/{scalar_name}/values', - num_subjects, - num_items, + n_files, + n_elements, storage_dtype=storage_dtype, compression=compression, compression_level=compression_level, @@ -111,16 +111,16 @@ def write_tiledb_scalar_matrices( output_path.mkdir(parents=True, exist_ok=True) for scalar_name, rows in scalars.items(): - num_subjects = len(rows) - if num_subjects == 0: + n_files = len(rows) + if n_files == 0: continue - num_items = rows[0].shape[0] + n_elements = rows[0].shape[0] dataset_path = f'scalars/{scalar_name}/values' tiledb_storage.create_empty_scalar_matrix_array( str(output_path), dataset_path, - num_subjects, - num_items, + n_files, + n_elements, storage_dtype=storage_dtype, compression=compression, compression_level=compression_level, diff --git a/src/modelarrayio/storage/h5_storage.py b/src/modelarrayio/storage/h5_storage.py index 7c08e1c..c9bf260 100644 --- a/src/modelarrayio/storage/h5_storage.py +++ b/src/modelarrayio/storage/h5_storage.py @@ -65,15 +65,15 @@ def resolve_compression(compression, compression_level, shuffle): def compute_chunk_shape_full_subjects( - num_subjects, num_items, item_chunk, target_chunk_mb, storage_np_dtype + n_files, n_elements, item_chunk, target_chunk_mb, storage_np_dtype ): """Compute a chunk shape for a full subject. Parameters ---------- - num_subjects : :obj:`int` + n_files : :obj:`int` Number of subjects. - num_items : :obj:`int` + n_elements : :obj:`int` Number of items. item_chunk : :obj:`int` Item chunk. @@ -88,8 +88,8 @@ def compute_chunk_shape_full_subjects( Chunk shape. """ chunk = storage_utils.compute_full_subject_chunk_shape( - num_subjects=num_subjects, - num_items=num_items, + n_files=n_files, + n_elements=n_elements, item_chunk=item_chunk, target_chunk_mb=target_chunk_mb, storage_np_dtype=storage_np_dtype, @@ -97,8 +97,8 @@ def compute_chunk_shape_full_subjects( logger.debug( 'Computed chunk shape: %s (subjects=%d, items=%d, item_chunk=%s, target_chunk_mb=%.2f)', chunk, - num_subjects, - num_items, + n_files, + n_elements, str(item_chunk), float(target_chunk_mb), ) @@ -153,22 +153,22 @@ def create_scalar_matrix_dataset( if stacked_values.dtype != storage_np_dtype: stacked_values = stacked_values.astype(storage_np_dtype) - num_subjects, num_items = stacked_values.shape + n_files, n_elements = stacked_values.shape chunk_shape = compute_chunk_shape_full_subjects( - num_subjects, num_items, chunk_voxels, target_chunk_mb, storage_np_dtype + n_files, n_elements, chunk_voxels, target_chunk_mb, storage_np_dtype ) logger.info( 'Creating dataset %s with shape (%d, %d), dtype=%s, chunks=%s, compression=%s', dataset_path, - num_subjects, - num_items, + n_files, + n_elements, storage_np_dtype, chunk_shape, str(comp), ) dset = h5file.create_dataset( dataset_path, - shape=(num_subjects, num_items), + shape=(n_files, n_elements), dtype=storage_np_dtype, chunks=chunk_shape, compression=comp, @@ -186,8 +186,8 @@ def create_scalar_matrix_dataset( def create_empty_scalar_matrix_dataset( h5file, dataset_path, - num_subjects, - num_items, + n_files, + n_elements, storage_dtype='float32', compression='gzip', compression_level=4, @@ -204,9 +204,9 @@ def create_empty_scalar_matrix_dataset( HDF5 file. dataset_path : :obj:`str` Dataset path. - num_subjects : :obj:`int` + n_files : :obj:`int` Number of subjects. - num_items : :obj:`int` + n_elements : :obj:`int` Number of items. storage_dtype : :obj:`str` Storage dtype. @@ -232,20 +232,20 @@ def create_empty_scalar_matrix_dataset( comp, comp_opts, use_shuffle = resolve_compression(compression, compression_level, shuffle) chunk_shape = compute_chunk_shape_full_subjects( - num_subjects, num_items, chunk_voxels, target_chunk_mb, storage_np_dtype + n_files, n_elements, chunk_voxels, target_chunk_mb, storage_np_dtype ) logger.info( 'Creating empty dataset %s with shape (%d, %d), dtype=%s, chunks=%s, compression=%s', dataset_path, - num_subjects, - num_items, + n_files, + n_elements, storage_np_dtype, chunk_shape, str(comp), ) dset = h5file.create_dataset( dataset_path, - shape=(num_subjects, num_items), + shape=(n_files, n_elements), dtype=storage_np_dtype, chunks=chunk_shape, compression=comp, @@ -293,13 +293,13 @@ def write_rows_in_column_stripes(dset, rows): Parameters ---------- dset : h5py.Dataset - Target dataset with shape (num_subjects, num_elements) and chunking set. + Target dataset with shape (n_files, num_elements) and chunking set. rows : Sequence[np.ndarray] List/sequence of 1D arrays, one per subject, length == num_elements. Each will be cast on write to dset.dtype if needed. """ - num_subjects, num_elements = dset.shape - if len(rows) != num_subjects: + n_files, num_elements = dset.shape + if len(rows) != n_files: raise ValueError('rows length does not match dataset subjects dimension') stripe_width = dset.chunks[1] if dset.chunks is not None else max(1, num_elements // 8) logger.info( @@ -309,7 +309,7 @@ def write_rows_in_column_stripes(dset, rows): str(dset.chunks), ) - buf = np.empty((num_subjects, stripe_width), dtype=dset.dtype) + buf = np.empty((n_files, stripe_width), dtype=dset.dtype) with logging_redirect_tqdm(): for start in tqdm( range(0, num_elements, stripe_width), diff --git a/src/modelarrayio/storage/tiledb_storage.py b/src/modelarrayio/storage/tiledb_storage.py index 33be392..67ef6a5 100644 --- a/src/modelarrayio/storage/tiledb_storage.py +++ b/src/modelarrayio/storage/tiledb_storage.py @@ -56,15 +56,15 @@ def _build_filter_list(compression: str | None, compression_level: int | None, s def compute_tile_shape_full_subjects( - num_subjects, num_items, item_tile, target_tile_mb, storage_np_dtype + n_files, n_elements, item_tile, target_tile_mb, storage_np_dtype ): """Compute a tile shape for a full subject. Parameters ---------- - num_subjects : :obj:`int` + n_files : :obj:`int` Number of subjects. - num_items : :obj:`int` + n_elements : :obj:`int` Number of items. item_tile : :obj:`int` Item tile. @@ -79,8 +79,8 @@ def compute_tile_shape_full_subjects( Tile shape. """ tile = storage_utils.compute_full_subject_chunk_shape( - num_subjects=num_subjects, - num_items=num_items, + n_files=n_files, + n_elements=n_elements, item_chunk=item_tile, target_chunk_mb=target_tile_mb, storage_np_dtype=storage_np_dtype, @@ -88,8 +88,8 @@ def compute_tile_shape_full_subjects( logger.debug( 'Computed tile shape: %s (subjects=%d, items=%d, item_tile=%s, target_tile_mb=%.2f)', tile, - num_subjects, - num_items, + n_files, + n_elements, str(item_tile), float(target_tile_mb), ) @@ -148,9 +148,9 @@ def create_scalar_matrix_array( if stacked_values.dtype != storage_np_dtype: stacked_values = stacked_values.astype(storage_np_dtype) - num_subjects, num_items = stacked_values.shape + n_files, n_elements = stacked_values.shape tile_shape = compute_tile_shape_full_subjects( - num_subjects, num_items, tile_voxels, target_tile_mb, storage_np_dtype + n_files, n_elements, tile_voxels, target_tile_mb, storage_np_dtype ) uri = os.path.join(base_uri, dataset_path) @@ -158,10 +158,10 @@ def create_scalar_matrix_array( # Domain and schema dim_subjects = tiledb.Dim( - name='subjects', domain=(0, num_subjects - 1), tile=tile_shape[0], dtype=np.int64 + name='subjects', domain=(0, n_files - 1), tile=tile_shape[0], dtype=np.int64 ) dim_items = tiledb.Dim( - name='items', domain=(0, num_items - 1), tile=tile_shape[1], dtype=np.int64 + name='items', domain=(0, n_elements - 1), tile=tile_shape[1], dtype=np.int64 ) dom = tiledb.Domain(dim_subjects, dim_items) attr_filters = _build_filter_list(compression, compression_level, shuffle) @@ -171,8 +171,8 @@ def create_scalar_matrix_array( logger.info( 'Creating TileDB array %s with shape (%d, %d), dtype=%s, tiles=%s', uri, - num_subjects, - num_items, + n_files, + n_elements, storage_np_dtype, tile_shape, ) @@ -196,8 +196,8 @@ def create_scalar_matrix_array( def create_empty_scalar_matrix_array( base_uri, dataset_path, - num_subjects, - num_items, + n_files, + n_elements, storage_dtype='float32', compression='zstd', compression_level=5, @@ -214,9 +214,9 @@ def create_empty_scalar_matrix_array( Base URI. dataset_path : :obj:`str` Dataset path. - num_subjects : :obj:`int` + n_files : :obj:`int` Number of subjects. - num_items : :obj:`int` + n_elements : :obj:`int` Number of items. storage_dtype : :obj:`str` Storage dtype. @@ -240,17 +240,17 @@ def create_empty_scalar_matrix_array( """ storage_np_dtype = resolve_dtype(storage_dtype) tile_shape = compute_tile_shape_full_subjects( - num_subjects, num_items, tile_voxels, target_tile_mb, storage_np_dtype + n_files, n_elements, tile_voxels, target_tile_mb, storage_np_dtype ) uri = os.path.join(base_uri, dataset_path) _ensure_parent_group(uri) dim_subjects = tiledb.Dim( - name='subjects', domain=(0, num_subjects - 1), tile=tile_shape[0], dtype=np.int64 + name='subjects', domain=(0, n_files - 1), tile=tile_shape[0], dtype=np.int64 ) dim_items = tiledb.Dim( - name='items', domain=(0, num_items - 1), tile=tile_shape[1], dtype=np.int64 + name='items', domain=(0, n_elements - 1), tile=tile_shape[1], dtype=np.int64 ) dom = tiledb.Domain(dim_subjects, dim_items) attr_filters = _build_filter_list(compression, compression_level, shuffle) @@ -260,8 +260,8 @@ def create_empty_scalar_matrix_array( logger.info( 'Creating empty TileDB array %s with shape (%d, %d), dtype=%s, tiles=%s', uri, - num_subjects, - num_items, + n_files, + n_elements, storage_np_dtype, tile_shape, ) @@ -285,18 +285,18 @@ def write_rows_in_column_stripes(uri: str, rows: Sequence[np.ndarray]): Parameters ---------- uri : str - Target array URI with shape (num_subjects, num_elements). + Target array URI with shape (n_files, num_elements). rows : Sequence[np.ndarray] List/sequence of 1D arrays, one per subject, length == num_elements. Each will be cast on write to array attr dtype if needed. """ with tiledb.open(uri, 'r') as Ainfo: dom = Ainfo.schema.domain - num_subjects = dom.dim(0).domain[1] - dom.dim(0).domain[0] + 1 + n_files = dom.dim(0).domain[1] - dom.dim(0).domain[0] + 1 num_elements = dom.dim(1).domain[1] - dom.dim(1).domain[0] + 1 attr_dtype = Ainfo.schema.attr(0).dtype - if len(rows) != num_subjects: + if len(rows) != n_files: raise ValueError('rows length does not match array subjects dimension') # Try to align stripe width to the items tile for best throughput @@ -304,7 +304,7 @@ def write_rows_in_column_stripes(uri: str, rows: Sequence[np.ndarray]): items_tile = Ainfo2.schema.domain.dim(1).tile stripe_width = items_tile if items_tile is not None else max(1, num_elements // 8) - buf = np.empty((num_subjects, stripe_width), dtype=attr_dtype) + buf = np.empty((n_files, stripe_width), dtype=attr_dtype) for start in range(0, num_elements, stripe_width): end = min(start + stripe_width, num_elements) width = end - start diff --git a/src/modelarrayio/storage/utils.py b/src/modelarrayio/storage/utils.py index a011aa2..a2afbed 100644 --- a/src/modelarrayio/storage/utils.py +++ b/src/modelarrayio/storage/utils.py @@ -36,29 +36,29 @@ def resolve_dtype( def compute_full_subject_chunk_shape( - num_subjects: int, - num_items: int, + n_files: int, + n_elements: int, item_chunk: int, target_chunk_mb: float, storage_np_dtype: str | np.dtype | type[np.floating], ) -> tuple[int, int]: """Compute a 2-D chunk/tile shape that spans all subjects.""" - num_subjects = int(num_subjects) - num_items = int(num_items) - if num_subjects <= 0 or num_items <= 0: + n_files = int(n_files) + n_elements = int(n_elements) + if n_files <= 0 or n_elements <= 0: raise ValueError( 'Cannot compute chunk shape with zero-length dimension: ' - f'num_subjects={num_subjects}, num_items={num_items}' + f'n_files={n_files}, n_elements={n_elements}' ) - subjects_per_chunk = num_subjects + subjects_per_chunk = n_files if int(item_chunk) > 0: - items_per_chunk = min(int(item_chunk), num_items) + items_per_chunk = min(int(item_chunk), n_elements) else: bytes_per_value = np.dtype(resolve_dtype(storage_np_dtype)).itemsize target_bytes = float(target_chunk_mb) * 1024.0 * 1024.0 items_per_chunk = max(1, int(target_bytes / (bytes_per_value * subjects_per_chunk))) - items_per_chunk = min(items_per_chunk, num_items) + items_per_chunk = min(items_per_chunk, n_elements) return subjects_per_chunk, items_per_chunk diff --git a/test/test_cifti_cli.py b/test/test_cifti_cli.py index 557b2d3..f2ce2f7 100644 --- a/test/test_cifti_cli.py +++ b/test/test_cifti_cli.py @@ -92,9 +92,9 @@ def test_concifti_cli_creates_expected_hdf5(tmp_path, monkeypatch): # Scalars dataset dset = h5['scalars/THICK/values'] - num_subjects, num_items = dset.shape - assert num_subjects == 2 - assert num_items == n_grayordinates + n_files, n_elements = dset.shape + assert n_files == 2 + assert n_elements == n_grayordinates # Column names exist and match subjects count grp = h5['scalars/THICK'] diff --git a/test/test_h5_storage.py b/test/test_h5_storage.py index 5145ebf..b12f01e 100644 --- a/test/test_h5_storage.py +++ b/test/test_h5_storage.py @@ -46,8 +46,8 @@ def test_resolve_compression_invalid_gzip_level_falls_back() -> None: def test_compute_chunk_shape_full_subjects() -> None: chunk = compute_chunk_shape_full_subjects( - num_subjects=3, - num_items=100, + n_files=3, + n_elements=100, item_chunk=10, target_chunk_mb=2.0, storage_np_dtype=np.float32, @@ -57,8 +57,8 @@ def test_compute_chunk_shape_full_subjects() -> None: def test_compute_chunk_shape_auto_item_chunk() -> None: chunk = compute_chunk_shape_full_subjects( - num_subjects=2, - num_items=50, + n_files=2, + n_elements=50, item_chunk=0, target_chunk_mb=1.0, storage_np_dtype=np.float32, @@ -68,14 +68,14 @@ def test_compute_chunk_shape_auto_item_chunk() -> None: @pytest.mark.parametrize( - ('num_subjects', 'num_items'), + ('n_files', 'n_elements'), [(0, 10), (3, 0), (0, 0)], ) -def test_compute_chunk_shape_rejects_zero_dimension(num_subjects: int, num_items: int) -> None: +def test_compute_chunk_shape_rejects_zero_dimension(n_files: int, n_elements: int) -> None: with pytest.raises(ValueError, match='zero-length'): compute_chunk_shape_full_subjects( - num_subjects, - num_items, + n_files, + n_elements, item_chunk=0, target_chunk_mb=2.0, storage_np_dtype=np.float32, @@ -103,16 +103,12 @@ def test_create_scalar_matrix_dataset_writes_data_and_attrs(tmp_path) -> None: def test_write_rows_in_column_stripes_matches_dense_write(tmp_path) -> None: """Stripe writer should match assigning the full matrix.""" path = tmp_path / 'stripe.h5' - num_subjects, num_elements = 3, 17 - full = np.arange(num_subjects * num_elements, dtype=np.float64).reshape( - num_subjects, num_elements - ) - rows = [full[i].copy() for i in range(num_subjects)] + n_files, num_elements = 3, 17 + full = np.arange(n_files * num_elements, dtype=np.float64).reshape(n_files, num_elements) + rows = [full[i].copy() for i in range(n_files)] with h5py.File(path, 'w') as h5: - dset = h5.create_dataset( - 'm', shape=(num_subjects, num_elements), dtype='f8', chunks=(3, 5) - ) + dset = h5.create_dataset('m', shape=(n_files, num_elements), dtype='f8', chunks=(3, 5)) write_rows_in_column_stripes(dset, rows) with h5py.File(path, 'r') as h5: diff --git a/test/test_storage_utils.py b/test/test_storage_utils.py index 5927dc7..65354b2 100644 --- a/test/test_storage_utils.py +++ b/test/test_storage_utils.py @@ -18,8 +18,8 @@ def test_resolve_dtype_rejects_unknown_values() -> None: def test_compute_full_subject_chunk_shape_auto() -> None: chunk = storage_utils.compute_full_subject_chunk_shape( - num_subjects=4, - num_items=128, + n_files=4, + n_elements=128, item_chunk=0, target_chunk_mb=1.0, storage_np_dtype='float32', diff --git a/test/test_voxels_cli.py b/test/test_voxels_cli.py index 729d82f..403493f 100644 --- a/test/test_voxels_cli.py +++ b/test/test_voxels_cli.py @@ -111,8 +111,8 @@ def test_convoxel_cli_creates_expected_hdf5(tmp_path, monkeypatch): # Scalars dataset dset = h5['scalars/FA/values'] - num_subjects, num_voxels = dset.shape - assert num_subjects == 2 + n_files, num_voxels = dset.shape + assert n_files == 2 assert num_voxels == ijk.shape[1] # Column names exist and match subjects count diff --git a/test/test_voxels_s3.py b/test/test_voxels_s3.py index 16b7520..6c52ffd 100644 --- a/test/test_voxels_s3.py +++ b/test/test_voxels_s3.py @@ -107,13 +107,13 @@ def test_convoxel_s3_parallel(tmp_path, group_mask_path, monkeypatch): with h5py.File(out_h5, 'r') as h5: dset = h5['scalars/alff/values'] - num_subjects, num_voxels = dset.shape + n_files, num_voxels = dset.shape - assert num_subjects == len(OHSU_SUBJECTS) + assert n_files == len(OHSU_SUBJECTS) assert num_voxels > 0 # Each subject should have at least some non-NaN values - for i in range(num_subjects): + for i in range(n_files): assert not np.all(np.isnan(dset[i, :])) # Column names recorded in the file From 5ae093566f428ada8e09c7022057101670b991dd Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 12:52:02 -0400 Subject: [PATCH 18/20] Rename more variables. num_elements --> n_elements num_voxels --> n_voxels --- src/modelarrayio/storage/h5_storage.py | 14 +++++++------- src/modelarrayio/storage/tiledb_storage.py | 12 ++++++------ test/test_h5_storage.py | 6 +++--- test/test_voxels_cli.py | 4 ++-- test/test_voxels_s3.py | 4 ++-- 5 files changed, 20 insertions(+), 20 deletions(-) diff --git a/src/modelarrayio/storage/h5_storage.py b/src/modelarrayio/storage/h5_storage.py index c9bf260..33a0d3c 100644 --- a/src/modelarrayio/storage/h5_storage.py +++ b/src/modelarrayio/storage/h5_storage.py @@ -293,15 +293,15 @@ def write_rows_in_column_stripes(dset, rows): Parameters ---------- dset : h5py.Dataset - Target dataset with shape (n_files, num_elements) and chunking set. + Target dataset with shape (n_files, n_elements) and chunking set. rows : Sequence[np.ndarray] - List/sequence of 1D arrays, one per subject, length == num_elements. + List/sequence of 1D arrays, one per subject, length == n_elements. Each will be cast on write to dset.dtype if needed. """ - n_files, num_elements = dset.shape + n_files, n_elements = dset.shape if len(rows) != n_files: raise ValueError('rows length does not match dataset subjects dimension') - stripe_width = dset.chunks[1] if dset.chunks is not None else max(1, num_elements // 8) + stripe_width = dset.chunks[1] if dset.chunks is not None else max(1, n_elements // 8) logger.info( 'Stripe-writing dataset %s with stripe width=%d (chunks=%s)', dset.name, @@ -312,14 +312,14 @@ def write_rows_in_column_stripes(dset, rows): buf = np.empty((n_files, stripe_width), dtype=dset.dtype) with logging_redirect_tqdm(): for start in tqdm( - range(0, num_elements, stripe_width), + range(0, n_elements, stripe_width), bar_format=( '{percentage:3.0f}% {n_fmt}/{total_fmt} [{elapsed}<{remaining}, {rate_fmt}]' ), ascii=True, - mininterval=max(1, (num_elements / stripe_width) // 200), + mininterval=max(1, (n_elements / stripe_width) // 200), ): - end = min(start + stripe_width, num_elements) + end = min(start + stripe_width, n_elements) width = end - start if width != stripe_width: # resize buffer view on last partial stripe diff --git a/src/modelarrayio/storage/tiledb_storage.py b/src/modelarrayio/storage/tiledb_storage.py index 67ef6a5..c82ae71 100644 --- a/src/modelarrayio/storage/tiledb_storage.py +++ b/src/modelarrayio/storage/tiledb_storage.py @@ -285,15 +285,15 @@ def write_rows_in_column_stripes(uri: str, rows: Sequence[np.ndarray]): Parameters ---------- uri : str - Target array URI with shape (n_files, num_elements). + Target array URI with shape (n_files, n_elements). rows : Sequence[np.ndarray] - List/sequence of 1D arrays, one per subject, length == num_elements. + List/sequence of 1D arrays, one per subject, length == n_elements. Each will be cast on write to array attr dtype if needed. """ with tiledb.open(uri, 'r') as Ainfo: dom = Ainfo.schema.domain n_files = dom.dim(0).domain[1] - dom.dim(0).domain[0] + 1 - num_elements = dom.dim(1).domain[1] - dom.dim(1).domain[0] + 1 + n_elements = dom.dim(1).domain[1] - dom.dim(1).domain[0] + 1 attr_dtype = Ainfo.schema.attr(0).dtype if len(rows) != n_files: @@ -302,11 +302,11 @@ def write_rows_in_column_stripes(uri: str, rows: Sequence[np.ndarray]): # Try to align stripe width to the items tile for best throughput with tiledb.open(uri, 'r') as Ainfo2: items_tile = Ainfo2.schema.domain.dim(1).tile - stripe_width = items_tile if items_tile is not None else max(1, num_elements // 8) + stripe_width = items_tile if items_tile is not None else max(1, n_elements // 8) buf = np.empty((n_files, stripe_width), dtype=attr_dtype) - for start in range(0, num_elements, stripe_width): - end = min(start + stripe_width, num_elements) + for start in range(0, n_elements, stripe_width): + end = min(start + stripe_width, n_elements) width = end - start if width != stripe_width: buf_view = buf[:, :width] diff --git a/test/test_h5_storage.py b/test/test_h5_storage.py index b12f01e..9a03df6 100644 --- a/test/test_h5_storage.py +++ b/test/test_h5_storage.py @@ -103,12 +103,12 @@ def test_create_scalar_matrix_dataset_writes_data_and_attrs(tmp_path) -> None: def test_write_rows_in_column_stripes_matches_dense_write(tmp_path) -> None: """Stripe writer should match assigning the full matrix.""" path = tmp_path / 'stripe.h5' - n_files, num_elements = 3, 17 - full = np.arange(n_files * num_elements, dtype=np.float64).reshape(n_files, num_elements) + n_files, n_elements = 3, 17 + full = np.arange(n_files * n_elements, dtype=np.float64).reshape(n_files, n_elements) rows = [full[i].copy() for i in range(n_files)] with h5py.File(path, 'w') as h5: - dset = h5.create_dataset('m', shape=(n_files, num_elements), dtype='f8', chunks=(3, 5)) + dset = h5.create_dataset('m', shape=(n_files, n_elements), dtype='f8', chunks=(3, 5)) write_rows_in_column_stripes(dset, rows) with h5py.File(path, 'r') as h5: diff --git a/test/test_voxels_cli.py b/test/test_voxels_cli.py index 403493f..3d3bfa9 100644 --- a/test/test_voxels_cli.py +++ b/test/test_voxels_cli.py @@ -111,9 +111,9 @@ def test_convoxel_cli_creates_expected_hdf5(tmp_path, monkeypatch): # Scalars dataset dset = h5['scalars/FA/values'] - n_files, num_voxels = dset.shape + n_files, n_voxels = dset.shape assert n_files == 2 - assert num_voxels == ijk.shape[1] + assert n_voxels == ijk.shape[1] # Column names exist and match subjects count grp = h5['scalars/FA'] diff --git a/test/test_voxels_s3.py b/test/test_voxels_s3.py index 6c52ffd..314dd8e 100644 --- a/test/test_voxels_s3.py +++ b/test/test_voxels_s3.py @@ -107,10 +107,10 @@ def test_convoxel_s3_parallel(tmp_path, group_mask_path, monkeypatch): with h5py.File(out_h5, 'r') as h5: dset = h5['scalars/alff/values'] - n_files, num_voxels = dset.shape + n_files, n_voxels = dset.shape assert n_files == len(OHSU_SUBJECTS) - assert num_voxels > 0 + assert n_voxels > 0 # Each subject should have at least some non-NaN values for i in range(n_files): From 69e23337862637a79ed2af7ba138528ab5de2dd1 Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 12:59:39 -0400 Subject: [PATCH 19/20] Default workers to 1. --- src/modelarrayio/cli/cifti_to_h5.py | 12 +++--------- src/modelarrayio/cli/mif_to_h5.py | 12 +++--------- src/modelarrayio/cli/nifti_to_h5.py | 12 +++--------- src/modelarrayio/cli/parser_utils.py | 5 ++--- 4 files changed, 11 insertions(+), 30 deletions(-) diff --git a/src/modelarrayio/cli/cifti_to_h5.py b/src/modelarrayio/cli/cifti_to_h5.py index 5b8ccd0..2275b75 100644 --- a/src/modelarrayio/cli/cifti_to_h5.py +++ b/src/modelarrayio/cli/cifti_to_h5.py @@ -4,7 +4,6 @@ import argparse import logging -import os from concurrent.futures import ThreadPoolExecutor, as_completed from pathlib import Path @@ -34,7 +33,7 @@ def cifti_to_h5( shuffle=True, chunk_voxels=0, target_chunk_mb=2.0, - workers=None, + workers=1, s3_workers=1, scalar_columns=None, ): @@ -63,7 +62,7 @@ def cifti_to_h5( target_chunk_mb : :obj:`float` Target chunk/tile size in MiB when auto-computing the spatial axis length workers : :obj:`int` - Maximum number of parallel TileDB write workers (``None`` = auto). + Maximum number of parallel TileDB write workers. Default 1. Has no effect when ``backend='hdf5'``. s3_workers : :obj:`int` Number of workers for parallel S3 downloads @@ -173,12 +172,7 @@ def _process_scalar_job(scalar_name, source_files): ) return scalar_name - worker_count = workers if isinstance(workers, int) and workers > 0 else None - if worker_count is None: - cpu_count = os.cpu_count() or 1 - worker_count = min(len(scalar_names), max(1, cpu_count)) - else: - worker_count = min(len(scalar_names), worker_count) + worker_count = min(len(scalar_names), workers) if worker_count <= 1: for scalar_name in scalar_names: diff --git a/src/modelarrayio/cli/mif_to_h5.py b/src/modelarrayio/cli/mif_to_h5.py index b1bff4d..ea753ff 100644 --- a/src/modelarrayio/cli/mif_to_h5.py +++ b/src/modelarrayio/cli/mif_to_h5.py @@ -4,7 +4,6 @@ import argparse import logging -import os from concurrent.futures import ThreadPoolExecutor, as_completed from functools import partial from pathlib import Path @@ -33,7 +32,7 @@ def mif_to_h5( shuffle=True, chunk_voxels=0, target_chunk_mb=2.0, - workers=None, + workers=1, s3_workers=1, scalar_columns=None, ): @@ -66,7 +65,7 @@ def mif_to_h5( target_chunk_mb : :obj:`float` Target chunk/tile size in MiB when auto-computing the spatial axis length workers : :obj:`int` - Maximum number of parallel TileDB write workers. Default 0 (auto). + Maximum number of parallel TileDB write workers. Default 1. Has no effect when ``backend='hdf5'``. s3_workers : :obj:`int` Number of parallel workers for S3 downloads. Default 1. @@ -133,12 +132,7 @@ def mif_to_h5( ) return int(not output.exists()) - worker_count = workers if isinstance(workers, int) and workers > 0 else None - if worker_count is None: - cpu_count = os.cpu_count() or 1 - worker_count = min(len(scalar_names), max(1, cpu_count)) - else: - worker_count = min(len(scalar_names), worker_count) + worker_count = min(len(scalar_names), workers) def _write_scalar_job(scalar_name): scalar_output = ( diff --git a/src/modelarrayio/cli/nifti_to_h5.py b/src/modelarrayio/cli/nifti_to_h5.py index 00c4f4d..12952ad 100644 --- a/src/modelarrayio/cli/nifti_to_h5.py +++ b/src/modelarrayio/cli/nifti_to_h5.py @@ -4,7 +4,6 @@ import argparse import logging -import os from concurrent.futures import ThreadPoolExecutor, as_completed from functools import partial from pathlib import Path @@ -34,7 +33,7 @@ def nifti_to_h5( shuffle=True, chunk_voxels=0, target_chunk_mb=2.0, - workers=None, + workers=1, s3_workers=1, scalar_columns=None, ): @@ -65,7 +64,7 @@ def nifti_to_h5( target_chunk_mb : :obj:`float` Target chunk/tile size in MiB when auto-computing. Default 2.0. workers : :obj:`int` - Maximum number of parallel TileDB write workers. Default 0 (auto). + Maximum number of parallel TileDB write workers. Default 1. Has no effect when ``backend='hdf5'``. s3_workers : :obj:`int` Number of parallel workers for S3 downloads. Default 1. @@ -134,12 +133,7 @@ def nifti_to_h5( ) return int(not output.exists()) - worker_count = workers if isinstance(workers, int) and workers > 0 else None - if worker_count is None: - cpu_count = os.cpu_count() or 1 - worker_count = min(len(scalar_names), max(1, cpu_count)) - else: - worker_count = min(len(scalar_names), worker_count) + worker_count = min(len(scalar_names), workers) def _write_scalar_job(scalar_name): scalar_output = ( diff --git a/src/modelarrayio/cli/parser_utils.py b/src/modelarrayio/cli/parser_utils.py index a083699..5a47836 100644 --- a/src/modelarrayio/cli/parser_utils.py +++ b/src/modelarrayio/cli/parser_utils.py @@ -95,11 +95,10 @@ def add_to_modelarray_args(parser, default_output='output.h5'): type=int, help=( 'Maximum number of parallel TileDB write workers. ' - 'Default 0 (auto, uses CPU count). ' - 'Set to 1 to disable parallel writes. ' + 'Default 1. ' 'Has no effect when --backend=hdf5.' ), - default=0, + default=1, ) s3_group = parser.add_argument_group('S3 arguments') From ae0ddb6fb0a8612f07d130989148ad29f9a4340f Mon Sep 17 00:00:00 2001 From: Taylor Salo Date: Wed, 1 Apr 2026 13:04:55 -0400 Subject: [PATCH 20/20] Update test_parser_utils.py --- test/test_parser_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_parser_utils.py b/test/test_parser_utils.py index 6fa13d4..0e2f75f 100644 --- a/test/test_parser_utils.py +++ b/test/test_parser_utils.py @@ -79,6 +79,6 @@ def test_tiledb_args_group(tmp_path_factory) -> None: args = p.parse_args(['--cohort-file', str(cohort_file), '--backend', 'tiledb']) assert args.output == Path('arrays.tdb') assert args.backend == 'tiledb' - assert args.workers == 0 + assert args.workers == 1 assert args.s3_workers == 1 assert args.log_level == 'INFO'