diff --git a/src/modelarrayio/cli/export_results.py b/src/modelarrayio/cli/export_results.py index 3e16c88..4cc13b8 100644 --- a/src/modelarrayio/cli/export_results.py +++ b/src/modelarrayio/cli/export_results.py @@ -14,6 +14,7 @@ from modelarrayio.cli.h5_to_cifti import h5_to_cifti from modelarrayio.cli.h5_to_mif import h5_to_mif from modelarrayio.cli.h5_to_nifti import h5_to_nifti +from modelarrayio.cli.h5_to_odx import h5_to_odx from modelarrayio.cli.parser_utils import _is_file, add_log_level_arg from modelarrayio.utils.misc import detect_modality_from_path @@ -87,7 +88,7 @@ def export_results( f'The template file appears to be MIF/fixel ({template_path!r}). ' 'For MIF/fixel results, supply --index-file and --directions-file.' ) - modality = 'cifti' + modality = 'odx' if detected == 'odx' else 'cifti' else: raise ValueError( 'Cannot determine modality. Provide --mask (NIfTI), ' @@ -128,6 +129,21 @@ def export_results( ) return 0 + if modality == 'odx': + if cohort_file is None and example_file is None: + raise ValueError('One of --cohort-file or --example-file is required for ODX results.') + if example_file is None: + logger.warning('No example ODX provided; using first source_file from cohort.') + example_file = pd.read_csv(cohort_file)['source_file'].iloc[0] + h5_to_odx( + example_odx=example_file, + in_file=in_file, + analysis_name=analysis_name, + compress=compress, + output_dir=output_path, + ) + return 0 + # cifti if cohort_file is None and example_file is None: raise ValueError('One of --cohort-file or --example-file is required for CIFTI results.') diff --git a/src/modelarrayio/cli/h5_to_odx.py b/src/modelarrayio/cli/h5_to_odx.py new file mode 100644 index 0000000..006d87a --- /dev/null +++ b/src/modelarrayio/cli/h5_to_odx.py @@ -0,0 +1,70 @@ +"""Convert ModelArray HDF5 results back to an ODX file. + +The ODX analogue of :mod:`modelarrayio.cli.h5_to_mif`. Reads the +``results//results_matrix`` written by ModelArray and paints each +result metric (e.g. ``.estimate``, ``.p.value.fdr``) onto a template +ODX's group-fixel geometry as a per-fixel (DPF) array, producing a single ODX +that can be visualized (e.g. in trxviz). For each ``p.value`` metric a +``1m.p.value`` (1 - p) companion is also written, matching ``h5_to_mif``. +""" + +from __future__ import annotations + +import logging +from pathlib import Path + +import h5py +import numpy as np + +from modelarrayio.cli import utils as cli_utils +from modelarrayio.utils.odx import write_odx_results + +logger = logging.getLogger(__name__) + + +def h5_to_odx(example_odx, in_file, analysis_name, compress, output_dir): + """Write ModelArray results from an HDF5 file onto a template ODX. + + Parameters + ---------- + example_odx : path-like + Template ODX whose group-fixel geometry the results are painted onto + (any of the per-subject ODX used to build the HDF5, or the group ODX). + in_file : path-like + HDF5 file containing ``results//results_matrix``. + analysis_name : str + Name of the ModelArray results group inside the HDF5 file. + compress : bool + Accepted for signature parity with the other ``h5_to_*`` writers; the + ``.odx`` archive is already compressed. + output_dir : path-like + Directory where ``.odx`` is written. + + Returns + ------- + status : :obj:`int` + 0 if successful. + """ + del compress # ODX archives are self-compressed; kept for API parity. + output_path = Path(output_dir) + output_path.mkdir(parents=True, exist_ok=True) + + with h5py.File(in_file, 'r') as h5_data: + results_matrix = h5_data[f'results/{analysis_name}/results_matrix'] + results_names = cli_utils.read_result_names( + h5_data, analysis_name, results_matrix, logger=logger + ) + results: dict[str, np.ndarray] = {} + for result_col, result_name in enumerate(results_names): + safe = cli_utils.sanitize_result_name(result_name) + column = np.asarray(results_matrix[result_col, :], dtype=np.float32) + results[f'{analysis_name}_{safe}'] = column + if 'p.value' in safe: + results[f'{analysis_name}_{safe.replace("p.value", "1m.p.value")}'] = ( + 1.0 - column + ).astype(np.float32) + + out_odx = output_path / f'{analysis_name}.odx' + write_odx_results(example_odx, results, out_odx) + logger.info('Wrote %d result DPF arrays to %s', len(results), out_odx) + return 0 diff --git a/src/modelarrayio/cli/odx_to_h5.py b/src/modelarrayio/cli/odx_to_h5.py new file mode 100644 index 0000000..eabde09 --- /dev/null +++ b/src/modelarrayio/cli/odx_to_h5.py @@ -0,0 +1,117 @@ +"""Convert per-subject ODX fixel data to an HDF5 (or TileDB) ModelArray file. + +The ODX analogue of :mod:`modelarrayio.cli.mif_to_h5`. Each cohort row points +to one per-subject ODX (single-column DPF) produced by +``odx combine --per-subject-odx DIR``; the shared group-fixel geometry is read +directly from the first ODX, so no separate index/directions files are needed. +""" + +from __future__ import annotations + +import logging +from pathlib import Path + +import h5py + +from modelarrayio.cli import utils as cli_utils +from modelarrayio.utils.odx import gather_fixels_from_odx, load_cohort_odx + +logger = logging.getLogger(__name__) + + +def odx_to_h5( + cohort_long, + backend='hdf5', + output=Path('fixelarray.h5'), + storage_dtype='float32', + compression='gzip', + compression_level=4, + shuffle=True, + chunk_voxels=0, + target_chunk_mb=2.0, + workers=1, + s3_workers=1, + split_outputs=False, +): + """Load per-subject ODX fixel data and write a ModelArray store. + + Parameters mirror :func:`modelarrayio.cli.mif_to_h5.mif_to_h5` minus the + ``index_file``/``directions_file`` arguments — ODX carries its own geometry. + + Returns + ------- + status : :obj:`int` + 0 if successful, non-zero if an expected output was not written. + """ + if cohort_long.empty: + raise ValueError('Cohort file does not contain any ODX scalar entries.') + + # Every template-space ODX shares the group-fixel geometry; read it once. + first_source = str(cohort_long.iloc[0]['source_file']) + fixel_table, voxel_table = gather_fixels_from_odx(first_source) + + logger.info('Extracting ODX data...') + scalars, sources_lists = load_cohort_odx(cohort_long, s3_workers) + if not sources_lists: + raise ValueError('Unable to derive scalar sources from cohort file.') + + scalar_names = list(sources_lists.keys()) + + if backend == 'hdf5': + if split_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) + cli_utils.write_table_dataset(h5_file, 'voxels', voxel_table) + cli_utils.write_hdf5_scalar_matrices( + h5_file, + 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()) + + # tiledb backend + for scalar_name in scalar_names: + scalar_output = ( + cli_utils.prefixed_output_path(output, scalar_name) if split_outputs else output + ) + cli_utils.write_tiledb_scalar_matrices( + scalar_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, + ) + return 0 diff --git a/src/modelarrayio/cli/to_modelarray.py b/src/modelarrayio/cli/to_modelarray.py index 6b39e83..76d1cf8 100644 --- a/src/modelarrayio/cli/to_modelarray.py +++ b/src/modelarrayio/cli/to_modelarray.py @@ -11,6 +11,7 @@ from modelarrayio.cli.cifti_to_h5 import cifti_to_h5 from modelarrayio.cli.mif_to_h5 import mif_to_h5 from modelarrayio.cli.nifti_to_h5 import nifti_to_h5 +from modelarrayio.cli.odx_to_h5 import odx_to_h5 from modelarrayio.cli.parser_utils import _is_file, add_log_level_arg from modelarrayio.utils.misc import load_and_normalize_cohort @@ -37,7 +38,7 @@ def to_modelarray( ): """Load neuroimaging data and write to an HDF5 or TileDB modelarray file. - The modality (NIfTI, CIFTI, or MIF/fixel) is autodetected from the source + The modality (NIfTI, CIFTI, MIF/fixel, or ODX) is autodetected from the source file extensions listed in the cohort file. Parameters @@ -91,6 +92,10 @@ def to_modelarray( ) return mif_to_h5(index_file=index_file, directions_file=directions_file, **common_kwargs) + if modality == 'odx': + # ODX carries its own group-fixel geometry, so no index/directions files. + return odx_to_h5(**common_kwargs) + # cifti return cifti_to_h5(**common_kwargs) @@ -106,7 +111,7 @@ def _parse_to_modelarray(): parser = argparse.ArgumentParser( description=( 'Convert neuroimaging data to a modelarray HDF5 file. ' - 'The modality (NIfTI, CIFTI, or MIF/fixel) is autodetected from ' + 'The modality (NIfTI, CIFTI, MIF/fixel, or ODX) is autodetected from ' 'the source file extensions in the cohort file.' ), formatter_class=argparse.ArgumentDefaultsHelpFormatter, diff --git a/src/modelarrayio/utils/misc.py b/src/modelarrayio/utils/misc.py index 7530cf6..e88092a 100644 --- a/src/modelarrayio/utils/misc.py +++ b/src/modelarrayio/utils/misc.py @@ -14,7 +14,7 @@ def detect_modality_from_path(path: str) -> str: - """Return ``'cifti'``, ``'mif'``, or ``'nifti'`` based on file extension. + """Return ``'cifti'``, ``'mif'``, ``'odx'``, or ``'nifti'`` based on file extension. Parameters ---------- @@ -36,11 +36,13 @@ def detect_modality_from_path(path: str) -> str: return 'cifti' if path.endswith(('.mif.gz', '.mif')): return 'mif' + if path.rstrip('/').endswith('.odx'): + return 'odx' if path.endswith(('.nii.gz', '.nii')): return 'nifti' raise ValueError( f'Cannot detect modality from file extension: {path!r}. ' - 'Expected .mif, .nii, .nii.gz, or a CIFTI compound extension ' + 'Expected .mif, .odx, .nii, .nii.gz, or a CIFTI compound extension ' '(e.g. .dscalar.nii, .pscalar.nii).' ) diff --git a/src/modelarrayio/utils/odx.py b/src/modelarrayio/utils/odx.py new file mode 100644 index 0000000..99efacf --- /dev/null +++ b/src/modelarrayio/utils/odx.py @@ -0,0 +1,177 @@ +"""Utility functions for ODX fixel data (the PennLINC odx-rs Python bindings). + +ODX files are the native container produced by ``odx combine``. Each +template-space ODX carries its own group-fixel geometry (offsets + directions + +``compact_to_ijk``) plus per-fixel scalar arrays (DPF), so — unlike the MIF +path — no separate ``index``/``directions`` files are needed. + +The cohort consumes **one per-subject ODX per row** (a single-column DPF), the +direct analogue of one ``.mif`` per subject in :mod:`modelarrayio.utils.mif`. +Produce them with ``odx combine --per-subject-odx DIR``. +""" + +from __future__ import annotations + +from collections import defaultdict +from pathlib import Path + +import numpy as np +import pandas as pd +from tqdm import tqdm + +__all__ = ['gather_fixels_from_odx', 'load_cohort_odx', 'write_odx_results'] + + +def _import_odx(): + """Import the optional ``odx`` package with a helpful error if missing.""" + try: + import odx + except ImportError as exc: # pragma: no cover - exercised only without odx + raise ImportError( + "The 'odx' package is required to read ODX fixel data. Install the " + 'odx-rs Python bindings (e.g. `pip install odx`) to use the odx modality.' + ) from exc + return odx + + +def gather_fixels_from_odx(odx_path): + """Build ``(fixel_table, voxel_table)`` from an ODX file's geometry. + + Mirrors :func:`modelarrayio.utils.mif.gather_fixels`, but sources the fixel + directions and per-voxel fixel layout straight from an ODX (``offsets`` + + ``directions`` + ``compact_to_ijk``). ODX lays fixels out in compact-voxel + order with monotonic ``offsets``, so the voxel/fixel ids need no resorting. + + Parameters + ---------- + odx_path : path-like + Path to an ODX file (archive or directory). + + Returns + ------- + fixel_table : :obj:`pandas.DataFrame` + Columns ``fixel_id, voxel_id, x, y, z``. + voxel_table : :obj:`pandas.DataFrame` + Columns ``voxel_id, i, j, k``. + """ + odx = _import_odx() + obj = odx.load(str(odx_path)) + offsets = np.asarray(obj.offsets, dtype=np.int64) # (n_vox + 1,) + directions = np.asarray(obj.directions, dtype=np.float32) # (n_fixels, 3) + ijk = np.asarray(obj.compact_to_ijk, dtype=np.int64) # (n_vox, 3) + n_vox = int(ijk.shape[0]) + counts = np.diff(offsets) if offsets.size else np.zeros(0, dtype=np.int64) + n_fixels = int(offsets[-1]) if offsets.size else 0 + + voxel_table = pd.DataFrame( + { + 'voxel_id': np.arange(n_vox, dtype=np.int32), + 'i': ijk[:, 0], + 'j': ijk[:, 1], + 'k': ijk[:, 2], + } + ) + fixel_voxel_ids = np.repeat(np.arange(n_vox, dtype=np.int32), counts) + fixel_table = pd.DataFrame( + { + 'fixel_id': np.arange(n_fixels, dtype=np.int32), + 'voxel_id': fixel_voxel_ids, + 'x': directions[:, 0], + 'y': directions[:, 1], + 'z': directions[:, 2], + } + ) + return fixel_table, voxel_table + + +def _read_odx_column(obj, scalar_name, source): + """Read one per-subject scalar column (a single-column DPF) from an ODX.""" + names = set(obj.dpf_names()) + if scalar_name not in names: + raise ValueError( + f"ODX '{source}' has no per-fixel scalar '{scalar_name}'. " + f'Available DPF arrays: {sorted(names)}' + ) + arr = np.asarray(obj.dpf(scalar_name), dtype=np.float32) # (n_fixels, ncols) + if arr.ndim != 2 or arr.shape[1] != 1: + raise ValueError( + f"ODX '{source}' scalar '{scalar_name}' has shape {arr.shape}; the odx " + 'modality expects one per-subject ODX per cohort row (a single-column ' + 'DPF). Produce per-subject ODX files with `odx combine --per-subject-odx DIR`.' + ) + return np.ascontiguousarray(arr[:, 0]) + + +def load_cohort_odx(cohort_long, s3_workers=1): + """Load all ODX scalar rows from the cohort. + + Each cohort row points to one per-subject ODX with a single-column DPF, + exactly as the MIF path uses one ``.mif`` per subject. Returns the same + structures as :func:`modelarrayio.utils.mif.load_cohort_mif`. + + Parameters + ---------- + cohort_long : :obj:`pandas.DataFrame` + Long-format cohort dataframe with columns ``scalar_name`` and + ``source_file``. + s3_workers : :obj:`int` + Accepted for signature parity with the MIF loader; ODX reads are + memory-mapped and run serially. + + 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 (HDF5 column metadata). + """ + odx = _import_odx() + scalars: dict[str, list[np.ndarray]] = defaultdict(list) + sources_lists: dict[str, list[str]] = defaultdict(list) + for row in tqdm(list(cohort_long.itertuples(index=False)), desc='Loading ODX data'): + scalar_name = str(row.scalar_name) + source = str(row.source_file) + obj = odx.load(source) + scalars[scalar_name].append(_read_odx_column(obj, scalar_name, source)) + sources_lists[scalar_name].append(source) + return scalars, sources_lists + + +def write_odx_results(template_odx, results, out_path): + """Write per-fixel result arrays onto a template ODX's geometry. + + Rebuilds an ODX from ``template_odx``'s group-fixel geometry (mask, affine, + offsets, directions) and attaches each result as a single-column DPF, so the + statistics map back onto the same fixels (e.g. for visualization in trxviz). + + Parameters + ---------- + template_odx : path-like + An ODX whose geometry defines the output fixels (typically one of the + per-subject ODX files, or the group ODX, used to build the HDF5). + results : Mapping[str, numpy.ndarray] + Maps a DPF name to a 1-D array of length ``n_fixels``. + out_path : path-like + Output ODX path (``.odx`` archive or directory). + + Returns + ------- + out_path : :obj:`pathlib.Path` + """ + odx = _import_odx() + tmpl = odx.load(str(template_odx)) + dims = tuple(int(d) for d in tmpl.dimensions) + affine = np.asarray(tmpl.affine, dtype=np.float64) + mask = np.ascontiguousarray(np.asarray(tmpl.mask, dtype=np.uint8).reshape(-1)) + offsets = np.asarray(tmpl.offsets, dtype=np.int64) + directions = np.ascontiguousarray(np.asarray(tmpl.directions, dtype=np.float32)) + n_fixels = int(offsets[-1]) if offsets.size else 0 + + builder = odx.OdxBuilder(affine, dims, mask) + for v in range(offsets.size - 1): + builder.push_voxel_peaks(directions[offsets[v]:offsets[v + 1]]) + for name, arr in results.items(): + col = np.ascontiguousarray(np.asarray(arr, dtype=np.float32).reshape(n_fixels, 1)) + builder.set_dpf(name, col) + builder.finalize().save(str(out_path)) + return Path(out_path) diff --git a/test/test_h5_to_odx_unit.py b/test/test_h5_to_odx_unit.py new file mode 100644 index 0000000..cda9064 --- /dev/null +++ b/test/test_h5_to_odx_unit.py @@ -0,0 +1,77 @@ +"""Unit tests for the ODX results write-back (h5_to_odx + export-results routing). + +These monkeypatch the ODX writer so they run without the optional ``odx`` +package, mirroring ``test_mif_to_h5_unit.py``. +""" + +from __future__ import annotations + +from pathlib import Path + +import h5py +import numpy as np + +from modelarrayio.cli import export_results as export_results_mod +from modelarrayio.cli import h5_to_odx as h5_to_odx_mod + + +def _results_h5(path: Path) -> None: + with h5py.File(path, 'w') as f: + grp = f.create_group('results/myana') + # 2 metrics x 3 fixels + rm = grp.create_dataset( + 'results_matrix', + data=np.array([[1.0, 2.0, 3.0], [0.01, 0.5, 0.9]], dtype=np.float32), + ) + rm.attrs['colnames'] = ['ses2.estimate', 'ses2.p.value'] + + +def test_h5_to_odx_builds_result_dpf_dict(monkeypatch, tmp_path: Path) -> None: + h5 = tmp_path / 'results.h5' + _results_h5(h5) + captured = {} + + def fake_write(template_odx, results, out_path): + captured['template'] = template_odx + captured['results'] = results + captured['out'] = out_path + return Path(out_path) + + monkeypatch.setattr(h5_to_odx_mod, 'write_odx_results', fake_write) + status = h5_to_odx_mod.h5_to_odx( + example_odx='tmpl.odx', + in_file=h5, + analysis_name='myana', + compress=False, + output_dir=tmp_path, + ) + assert status == 0 + assert captured['template'] == 'tmpl.odx' + assert captured['out'] == tmp_path / 'myana.odx' + res = captured['results'] + # estimate metric carried through; p.value gets a 1m.p.value companion + assert np.allclose(res['myana_ses2.estimate'], [1.0, 2.0, 3.0]) + assert np.allclose(res['myana_ses2.p.value'], [0.01, 0.5, 0.9]) + assert np.allclose(res['myana_ses2.1m.p.value'], [0.99, 0.5, 0.1]) + + +def test_export_results_routes_odx(monkeypatch, tmp_path: Path) -> None: + h5 = tmp_path / 'results.h5' + _results_h5(h5) + calls = {} + + def fake_h5_to_odx(**kwargs): + calls.update(kwargs) + return 0 + + monkeypatch.setattr(export_results_mod, 'h5_to_odx', fake_h5_to_odx) + status = export_results_mod.export_results( + in_file=h5, + analysis_name='myana', + output_dir=tmp_path / 'out', + example_file=str(tmp_path / 'group.odx'), # .odx → odx modality + ) + assert status == 0 + # routed to the ODX writer (not mis-routed to CIFTI), with the .odx template + assert calls['example_odx'] == str(tmp_path / 'group.odx') + assert calls['analysis_name'] == 'myana' diff --git a/test/test_odx_to_h5_unit.py b/test/test_odx_to_h5_unit.py new file mode 100644 index 0000000..71d740c --- /dev/null +++ b/test/test_odx_to_h5_unit.py @@ -0,0 +1,94 @@ +"""Unit tests for the ODX modality (detection + odx_to_h5 branching). + +These tests monkeypatch the ODX readers so they run without the optional +``odx`` package installed, mirroring ``test_mif_to_h5_unit.py``. +""" + +from __future__ import annotations + +from pathlib import Path + +import h5py +import numpy as np +import pandas as pd +import pytest + +from modelarrayio.cli import odx_to_h5 +from modelarrayio.utils.misc import detect_modality_from_path + + +def _cohort() -> pd.DataFrame: + return pd.DataFrame( + { + 'scalar_name': ['angle_deg', 'angle_deg'], + 'source_file': ['sub-01.odx', 'sub-02.odx'], + 'subject': ['sub-01', 'sub-02'], + } + ) + + +def _fixtures(): + fixel_table = pd.DataFrame( + {'fixel_id': [0, 1], 'voxel_id': [0, 0], 'x': [1.0, 0.0], 'y': [0.0, 1.0], 'z': [0.0, 0.0]} + ) + voxel_table = pd.DataFrame({'voxel_id': [0], 'i': [0], 'j': [0], 'k': [0]}) + scalars = {'angle_deg': [np.array([0.0, 5.0], np.float32), np.array([3.0, 7.0], np.float32)]} + sources = {'angle_deg': ['sub-01.odx', 'sub-02.odx']} + return fixel_table, voxel_table, scalars, sources + + +def test_detect_modality_odx() -> None: + assert detect_modality_from_path('group/sub-01.odx') == 'odx' + assert detect_modality_from_path('group/sub-01.odx/') == 'odx' # directory layout + + +def test_odx_to_h5_raises_when_sources_missing(monkeypatch, tmp_path: Path) -> None: + monkeypatch.setattr( + odx_to_h5, 'gather_fixels_from_odx', lambda _src: (pd.DataFrame(), pd.DataFrame()) + ) + monkeypatch.setattr(odx_to_h5, 'load_cohort_odx', lambda cohort_long, s3_workers: ({}, {})) + with pytest.raises(ValueError, match='Unable to derive scalar sources'): + odx_to_h5.odx_to_h5(_cohort(), output=tmp_path / 'out.h5') + + +def test_odx_to_h5_writes_modelarray_schema(monkeypatch, tmp_path: Path) -> None: + fixel_table, voxel_table, scalars, sources = _fixtures() + monkeypatch.setattr( + odx_to_h5, 'gather_fixels_from_odx', lambda *_a, **_k: (fixel_table, voxel_table) + ) + monkeypatch.setattr(odx_to_h5, 'load_cohort_odx', lambda *_a, **_k: (scalars, sources)) + + out = tmp_path / 'angle.h5' + status = odx_to_h5.odx_to_h5(_cohort(), backend='hdf5', output=out) + assert status == 0 + + with h5py.File(out, 'r') as f: + values = f['scalars/angle_deg/values'] + # ModelArray's per-scalar values matrix: (n_subjects, n_fixels) + assert values.shape == (2, 2) + assert np.allclose(values[0], [0.0, 5.0]) + assert np.allclose(values[1], [3.0, 7.0]) + cols = [c.decode() if isinstance(c, bytes) else str(c) for c in f['scalars/angle_deg/column_names'][()]] + assert cols == ['sub-01.odx', 'sub-02.odx'] + # fixel geometry is carried for mapping results back + assert 'fixels' in f + assert 'voxels' in f + assert list(f['fixels'].attrs['column_names']) == ['fixel_id', 'voxel_id', 'x', 'y', 'z'] + + +def test_odx_to_h5_split_outputs(monkeypatch, tmp_path: Path) -> None: + fixel_table, voxel_table, _scalars, _sources = _fixtures() + scalars = { + 'angle_deg': [np.array([0.0, 5.0], np.float32), np.array([3.0, 7.0], np.float32)], + 'afd': [np.array([0.1, 0.2], np.float32), np.array([0.3, 0.4], np.float32)], + } + sources = {'angle_deg': ['sub-01.odx', 'sub-02.odx'], 'afd': ['sub-01.odx', 'sub-02.odx']} + monkeypatch.setattr( + odx_to_h5, 'gather_fixels_from_odx', lambda *_a, **_k: (fixel_table, voxel_table) + ) + monkeypatch.setattr(odx_to_h5, 'load_cohort_odx', lambda *_a, **_k: (scalars, sources)) + + status = odx_to_h5.odx_to_h5(_cohort(), backend='hdf5', output=tmp_path / 'fixels.h5', split_outputs=True) + assert status == 0 + assert (tmp_path / 'angle_deg_fixels.h5').exists() + assert (tmp_path / 'afd_fixels.h5').exists()