From 7cd038fb888fc055f88009e3a5af5e3513c1a8c6 Mon Sep 17 00:00:00 2001 From: Eric Cramer <13970720+emcramer@users.noreply.github.com> Date: Tue, 21 Jul 2026 15:33:57 -0700 Subject: [PATCH] fix: PhysiCell reader orientation and label handling Resolve cell matrix orientation from the frame's declared variable count instead of comparing matrix dimensions. The old heuristic assumed every frame holds more cells than the model has variables, silently returning transposed data for sparse initial conditions, early time steps, and simulations approaching extinction. Derive column positions from the MultiCellDS block, which was already parsed into metadata.extra['custom_labels'] but never read. This makes custom and non-indexed columns reachable via a new 'columns' key, and lets models with differing variable counts parse in one session. Also fix to_dataframe() ignoring include_dead_cells, add an orientation guard to parse_microenvironment_mat, return raw_data as loaded from disk, and single-source __version__ from distribution metadata. Reported by the llm-abm-consistency harness (OHSU ChangLab). --- CHANGELOG.md | 72 ++++ spatialtissuepy/__init__.py | 11 +- spatialtissuepy/synthetic/physicell/parser.py | 198 ++++++++++- spatialtissuepy/synthetic/physicell/reader.py | 84 ++++- tests/test_physicell_labels.py | 330 ++++++++++++++++++ 5 files changed, 661 insertions(+), 34 deletions(-) create mode 100644 tests/test_physicell_labels.py diff --git a/CHANGELOG.md b/CHANGELOG.md index 8af0a81..2872ce6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -5,6 +5,78 @@ All notable changes to spatialtissuepy are documented in this file. The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +PhysiCell reader correctness fixes, reported against v0.2.0 by the +llm-abm-consistency harness (OHSU ChangLab) and reproduced against both their +PhysiCell 1.14.2 output and this repository's bundled example simulation. + +### Fixed + +- **Cell matrix orientation is now resolved from the frame's declared variable + count** rather than by comparing the matrix dimensions. The old heuristic + (`if shape[0] > shape[1]: transpose`) assumed every frame holds more cells + than the model has variables, and silently returned transposed garbage for + any frame that does not — sparse initial conditions, early time steps, small + explants, and simulations approaching extinction. A frame that matches + neither orientation now raises `ValueError` instead of returning data. When + no labels are available, the old heuristic still applies but emits a + `UserWarning` naming the ambiguity. +- **`PhysiCellTimeStep.to_dataframe()` now honors `include_dead_cells`.** It + previously ignored the flag while `n_cells`, `positions`, `to_spatial_data()` + and `cell_counts_by_type()` all respected it, so the same object reported + different populations depending on the accessor and row indices did not align + between `to_dataframe()` and `positions`. +- **`parse_microenvironment_mat` validates matrix orientation** against the + declared substrate count when `substrate_names` is supplied, converting a + would-be silent misread into a clear error. +- **`__version__` is now read from the installed distribution metadata**, so it + cannot drift from `pyproject.toml`. It reported `0.1.0` at the `v0.2.0` tag, + which misattributed results in downstream provenance capture. A regression + test asserts the two agree. + +### Added + +- **`parse_cells_mat(..., labels=...)`**: derive column positions from the + MultiCellDS `` block that `parse_physicell_xml` already parses. The + parsed labels were previously stored at `metadata.extra['custom_labels']` and + never read anywhere in the package. Precedence is `index_mapping` > + `labels` > the existing row-count autodetect, so behavior is unchanged when + labels are not supplied. +- **`'columns'` key in the `parse_cells_mat` result**: every labelled column, + including standard PhysiCell fields absent from the hard-coded index tables + (`is_motile`, `migration_speed`) and model-specific custom variables. The + return dict was previously a closed set, so anything else the model wrote was + unreachable. +- **`'orientation'` key** reporting the on-disk layout of `raw_data`, either + `'variables_x_cells'` or `'cells_x_variables'`. +- **`expand_cell_labels()` and `declared_variable_count()`** helpers for + working with `` blocks. Vector labels expand as `{name}_x/_y/_z` for + `size == 3` and `{name}_0 … {name}_{n-1}` otherwise. +- **`to_dataframe(extra_columns=True)`** appends every labelled column. Off by + default, since a full PhysiCell frame carries 150+ variables. +- **`to_dataframe(include_dead_cells=...)`** to override the instance attribute + per call. +- Regression tests (`tests/test_physicell_labels.py`, 25 tests) covering label + expansion, all four orientation cases, cross-model label isolation, the + microenvironment guard, dead-cell filtering, and version consistency. + +### Changed + +- **`raw_data` from `parse_cells_mat` is now returned exactly as loaded from + disk.** It was previously the reoriented matrix, so callers could not use it + to recover the truth when the orientation heuristic misfired. Callers relying + on the transposed `raw_data` should consult the new `'orientation'` key. +- `PhysiCellTimeStep._load_cell_data()` and `to_trajectory_dataframe()` now + resolve columns from each frame's own XML, so models with differing variable + counts parse correctly in the same session. + +### Known gaps + +- `parse_microenvironment_mat` is implemented and exported but still not wired + into `PhysiCellTimeStep` / `PhysiCellSimulation`; substrate fields remain + reachable only through the parser. Planned for v0.3.0. + ## [0.2.0] - 2026-06-26 First beta release. This release adds AI agent access via an MCP server, a diff --git a/spatialtissuepy/__init__.py b/spatialtissuepy/__init__.py index 3e090af..75c116f 100644 --- a/spatialtissuepy/__init__.py +++ b/spatialtissuepy/__init__.py @@ -21,7 +21,16 @@ io : Input/output utilities """ -__version__ = "0.1.0" +from importlib.metadata import PackageNotFoundError +from importlib.metadata import version as _pkg_version + +try: + # Single source of truth: the version declared in pyproject.toml, as + # recorded in the installed distribution metadata. + __version__ = _pkg_version("spatialtissuepy") +except PackageNotFoundError: # pragma: no cover - running from a source tree + __version__ = "0.0.0.dev0" + __author__ = "spatialtissuepy developers" from spatialtissuepy.core.spatial_data import SpatialTissueData diff --git a/spatialtissuepy/synthetic/physicell/parser.py b/spatialtissuepy/synthetic/physicell/parser.py index 874ba67..59fc12f 100644 --- a/spatialtissuepy/synthetic/physicell/parser.py +++ b/spatialtissuepy/synthetic/physicell/parser.py @@ -7,6 +7,7 @@ from __future__ import annotations +import warnings import xml.etree.ElementTree as ET from dataclasses import dataclass from pathlib import Path @@ -287,10 +288,119 @@ def get_cell_type_mapping( } +def expand_cell_labels(labels: Dict[int, Tuple[str, int]]) -> Dict[int, str]: + """ + Expand a MultiCellDS ```` block into a flat column mapping. + + Vector-valued labels occupy consecutive columns in the cell matrix. A label + of ``size == 3`` is expanded to ``{name}_x/_y/_z`` (the convention PhysiCell + uses for spatial vectors); any other ``size > 1`` is expanded positionally + to ``{name}_0 ... {name}_{size-1}``. + + Parameters + ---------- + labels : dict + Mapping of ``{start_index: (name, size)}``, as produced by + :func:`parse_physicell_xml` at ``metadata.extra['custom_labels']``. + + Returns + ------- + dict + Mapping of ``{column_index: column_name}``. + + Examples + -------- + >>> expand_cell_labels({0: ('ID', 1), 1: ('position', 3)}) + {0: 'ID', 1: 'position_x', 2: 'position_y', 3: 'position_z'} + """ + expanded: Dict[int, str] = {} + + for start_index, (name, size) in labels.items(): + if size == 1: + expanded[start_index] = name + elif size == 3: + for offset, axis in enumerate('xyz'): + expanded[start_index + offset] = f'{name}_{axis}' + else: + for offset in range(size): + expanded[start_index + offset] = f'{name}_{offset}' + + return expanded + + +def declared_variable_count(labels: Dict[int, Tuple[str, int]]) -> int: + """ + Number of matrix rows implied by a ```` block. + + This accounts for the width of the final label, so a trailing vector label + is not undercounted the way ``max(labels) + 1`` would be. + + Parameters + ---------- + labels : dict + Mapping of ``{start_index: (name, size)}``. + + Returns + ------- + int + Expected number of variables (rows) per cell. + """ + return max(index + size for index, (_, size) in labels.items()) + + +def _orient_cell_matrix( + cell_matrix: np.ndarray, + n_expected: Optional[int], + mat_path: Path, +) -> Tuple[np.ndarray, str]: + """ + Normalize a cell matrix to ``(n_variables, n_cells)``. + + When the declared variable count is known, orientation is resolved by + matching it against the array shape -- unambiguous even when a frame holds + fewer cells than the model has variables. Without it, fall back to the + legacy magnitude heuristic and warn, since that assumption silently + corrupts small frames. + + Returns + ------- + (np.ndarray, str) + The matrix in ``(n_variables, n_cells)`` orientation, and the + orientation of the input as loaded from disk. + """ + n_rows, n_cols = cell_matrix.shape + + if n_expected is not None: + # Prefer axis 0 as variables, so a square matrix resolves the way + # MultiCellDS actually writes it. + if n_rows == n_expected: + return cell_matrix, 'variables_x_cells' + if n_cols == n_expected: + return cell_matrix.T, 'cells_x_variables' + raise ValueError( + f"cell matrix {cell_matrix.shape} matches neither orientation for " + f"{n_expected} declared variables in {mat_path}" + ) + + if n_rows > n_cols: + warnings.warn( + f"Guessing cell matrix orientation for {mat_path} from its shape " + f"{cell_matrix.shape} because no block was supplied. This " + f"is wrong for any frame holding fewer cells than the model has " + f"variables. Pass labels= to resolve orientation unambiguously.", + UserWarning, + stacklevel=3, + ) + return cell_matrix.T, 'cells_x_variables' + + return cell_matrix, 'variables_x_cells' + + def parse_cells_mat( mat_path: Path, cell_type_mapping: Optional[Dict[int, str]] = None, - index_mapping: Optional[Dict[str, int]] = None + index_mapping: Optional[Dict[str, int]] = None, + labels: Optional[Dict[int, Tuple[str, int]]] = None, ) -> Dict[str, np.ndarray]: """ Parse a PhysiCell cells_physicell.mat file. @@ -303,6 +413,13 @@ def parse_cells_mat( Mapping from cell type IDs to names. index_mapping : dict, optional Custom row index mapping for different PhysiCell versions. + labels : dict, optional + The ```` block from the frame's XML, as + ``{start_index: (name, size)}``. When supplied, column positions and + matrix orientation are resolved from the frame's own declaration rather + than from a hard-coded table, and every labelled column is returned + under ``'columns'``. Precedence for column positions is + ``index_mapping`` > ``labels`` > row-count autodetect. Returns ------- @@ -314,8 +431,18 @@ def parse_cells_mat( - 'volumes': (n_cells,) array of total volumes - 'radii': (n_cells,) array of cell radii - 'phases': (n_cells,) array of cell cycle phases + - 'dead_flags': (n_cells,) array of dead flags - 'ids': (n_cells,) array of cell IDs - - 'raw_data': Full matrix from file + - 'columns': dict of every labelled column, keyed by name. Empty when + no ``labels`` were supplied. + - 'raw_data': Full matrix exactly as loaded from disk + - 'orientation': layout of ``raw_data``, either ``'variables_x_cells'`` + or ``'cells_x_variables'`` + + Notes + ----- + ``raw_data`` is returned unmodified, so it round-trips to what + ``scipy.io.loadmat`` returned; use ``orientation`` to interpret it. """ from scipy.io import loadmat @@ -340,13 +467,10 @@ def parse_cells_mat( if cell_matrix is None: raise ValueError(f"Could not find cell data in {mat_path}") - # Ensure cells are columns - if cell_matrix.shape[0] > cell_matrix.shape[1]: - cell_matrix = cell_matrix.T + raw_data = cell_matrix - n_cells = cell_matrix.shape[1] - - if n_cells == 0: + # An empty frame carries no orientation signal; report it as written. + if 0 in cell_matrix.shape: return { 'positions': np.empty((0, 3)), 'cell_types': np.array([], dtype=str), @@ -354,15 +478,31 @@ def parse_cells_mat( 'volumes': np.array([]), 'radii': np.array([]), 'phases': np.array([], dtype=int), + 'dead_flags': np.array([], dtype=int), 'ids': np.array([], dtype=int), - 'raw_data': cell_matrix, + 'columns': {}, + 'raw_data': raw_data, + 'orientation': 'variables_x_cells', } - # Use provided or default index mapping + label_columns = expand_cell_labels(labels) if labels else {} + n_expected = declared_variable_count(labels) if labels else None + + cell_matrix, orientation = _orient_cell_matrix( + cell_matrix, n_expected, mat_path + ) + + n_cells = cell_matrix.shape[1] + + # Column positions: explicit index_mapping wins, then the frame's own + # labels, then the legacy row-count autodetect. if index_mapping is None: - # Auto-detect based on matrix shape - # PhysiCell 1.10+ typically has 150+ rows per cell - if cell_matrix.shape[0] >= 30: + if label_columns: + index_mapping = { + name: index for index, name in label_columns.items() + } + elif cell_matrix.shape[0] >= 30: + # PhysiCell 1.10+ typically has 150+ rows per cell index_mapping = CELL_DATA_INDICES_V2 else: index_mapping = CELL_DATA_INDICES_LEGACY @@ -423,6 +563,14 @@ def parse_cells_mat( # Fall back to inferring from phase code dead_flags = np.array([1 if p >= 100 else 0 for p in phases]) + # Every labelled column, including ones with no entry in the index tables + # (is_motile, migration_speed) and model-specific custom variables. + columns = { + name: cell_matrix[index, :] + for index, name in sorted(label_columns.items()) + if index < cell_matrix.shape[0] + } + return { 'positions': positions, 'cell_types': cell_types, @@ -432,7 +580,9 @@ def parse_cells_mat( 'phases': phases, 'dead_flags': dead_flags, 'ids': ids, - 'raw_data': cell_matrix, + 'columns': columns, + 'raw_data': raw_data, + 'orientation': orientation, } @@ -479,7 +629,23 @@ def parse_microenvironment_mat( if me_matrix is None: raise ValueError(f"Could not find microenvironment data in {mat_path}") - # Structure: rows 0-2 are x,y,z; row 3 is volume; rows 4+ are substrates + raw_data = me_matrix + + # Structure: rows 0-2 are x,y,z; row 3 is volume; rows 4+ are substrates. + # PhysiCell always writes (4 + n_substrates, n_voxels), but validate it + # rather than assume -- an unchecked transpose here is silent corruption. + if substrate_names is not None: + n_expected = 4 + len(substrate_names) + if me_matrix.shape[0] != n_expected: + if me_matrix.shape[1] == n_expected: + me_matrix = me_matrix.T + else: + raise ValueError( + f"microenvironment matrix {me_matrix.shape} matches neither " + f"orientation for {len(substrate_names)} declared substrates " + f"in {mat_path}" + ) + voxel_positions = me_matrix[:3, :].T # (n_voxels, 3) # Extract substrate concentrations @@ -495,7 +661,7 @@ def parse_microenvironment_mat( return { 'voxel_positions': voxel_positions, 'concentrations': concentrations, - 'raw_data': me_matrix, + 'raw_data': raw_data, } diff --git a/spatialtissuepy/synthetic/physicell/reader.py b/spatialtissuepy/synthetic/physicell/reader.py index 86c52c3..6ca8132 100644 --- a/spatialtissuepy/synthetic/physicell/reader.py +++ b/spatialtissuepy/synthetic/physicell/reader.py @@ -8,6 +8,7 @@ from __future__ import annotations import re +import xml.etree.ElementTree as ET from dataclasses import dataclass, field from pathlib import Path from typing import TYPE_CHECKING, Dict, List, Optional, Tuple, Union @@ -74,9 +75,21 @@ class PhysiCellTimeStep(ABMTimeStep): def _load_cell_data(self) -> Dict[str, np.ndarray]: """Load cell data from MAT file (cached).""" if self._cell_data is None: + # Resolve columns and matrix orientation from this frame's own + # block, so models with differing variable counts each + # parse against their own declaration. + labels = None + try: + labels = self._load_metadata().extra.get('custom_labels') or None + except (OSError, ET.ParseError): + # No readable XML: parse_cells_mat falls back to its row-count + # heuristic and warns. + pass + self._cell_data = parse_cells_mat( self.cells_mat_path, - self.cell_type_mapping + self.cell_type_mapping, + labels=labels, ) return self._cell_data @@ -180,35 +193,72 @@ def to_spatial_data(self) -> SpatialTissueData: metadata=extra_metadata, ) - def to_dataframe(self) -> pd.DataFrame: + def to_dataframe( + self, + include_dead_cells: Optional[bool] = None, + extra_columns: bool = False, + ) -> pd.DataFrame: """ Convert cell data to a pandas DataFrame. + Parameters + ---------- + include_dead_cells : bool, optional + Whether to include dead cells. Defaults to the instance attribute, + so rows stay aligned with ``positions`` and ``to_spatial_data()``. + Pass ``True`` to recover the unfiltered frame. + extra_columns : bool, default False + Append every labelled column from the frame's XML, including + standard fields absent from the index tables (``is_motile``, + ``migration_speed``) and model-specific custom variables. Off by + default because a full PhysiCell frame carries 150+ variables. + Returns ------- pd.DataFrame - DataFrame with cell properties. + DataFrame with cell properties, in the same row order as + ``positions``. """ data = self._load_cell_data() + if include_dead_cells is None: + include_dead_cells = self.include_dead_cells + + if include_dead_cells: + mask = np.ones(len(data['cell_types']), dtype=bool) + else: + mask = data['dead_flags'] == 0 + df = pd.DataFrame({ - 'cell_id': data['ids'], - 'x': data['positions'][:, 0], - 'y': data['positions'][:, 1], - 'z': data['positions'][:, 2], - 'cell_type': data['cell_types'], - 'cell_type_id': data['cell_type_ids'], - 'volume': data['volumes'], - 'radius': data['radii'], - 'phase': data['phases'], - 'is_dead': data['dead_flags'].astype(bool), - 'is_alive': ~data['dead_flags'].astype(bool), + 'cell_id': data['ids'][mask], + 'x': data['positions'][mask, 0], + 'y': data['positions'][mask, 1], + 'z': data['positions'][mask, 2], + 'cell_type': data['cell_types'][mask], + 'cell_type_id': data['cell_type_ids'][mask], + 'volume': data['volumes'][mask], + 'radius': data['radii'][mask], + 'phase': data['phases'][mask], + 'is_dead': data['dead_flags'][mask].astype(bool), + 'is_alive': ~data['dead_flags'][mask].astype(bool), }) - df['time'] = self.time - df['time_index'] = self.time_index + trailing = pd.DataFrame({ + 'time': np.full(len(df), self.time), + 'time_index': np.full(len(df), self.time_index), + }) + + if extra_columns: + # Build in one concat; a full PhysiCell frame adds 150+ columns and + # inserting them one at a time fragments the frame badly. + extra = pd.DataFrame({ + name: values[mask] + for name, values in data['columns'].items() + if name not in df.columns + }) + return pd.concat([df, extra, trailing], axis=1) - return df + return pd.concat([df, trailing], axis=1) def cell_counts_by_type(self) -> Dict[str, int]: """Get cell counts by type.""" diff --git a/tests/test_physicell_labels.py b/tests/test_physicell_labels.py new file mode 100644 index 0000000..a47a58a --- /dev/null +++ b/tests/test_physicell_labels.py @@ -0,0 +1,330 @@ +""" +Regression tests for PhysiCell reader defects fixed in v0.2.1. + +These cover the label-driven column extraction and orientation handling +reported by the llm-abm-consistency harness (OHSU ChangLab, 2026-07-21). + +The bundled example simulation cannot exercise the orientation bug on its own: +it writes 154 variables per frame and its smallest frame holds 906 cells, so +``n_cells > n_variables`` holds everywhere. These tests build small matrices +where that assumption fails. +""" + +import importlib.metadata +import warnings +from pathlib import Path + +import numpy as np +import pandas as pd +import pytest +from scipy.io import loadmat, savemat + +import spatialtissuepy +from spatialtissuepy.synthetic.physicell import ( + discover_physicell_timesteps, + read_physicell_timestep, +) +from spatialtissuepy.synthetic.physicell.parser import ( + declared_variable_count, + expand_cell_labels, + parse_cells_mat, + parse_microenvironment_mat, + parse_physicell_xml, +) + +N_VARS = 87 + + +@pytest.fixture +def example_physicell_dir(): + """Path to example PhysiCell simulation output folder.""" + path = ( + Path(__file__).parent.parent + / 'examples' / 'sample_data' / 'example_physicell_sim' + ) + if not path.exists(): + pytest.skip(f"Example PhysiCell data not found at {path}") + return path + + +@pytest.fixture +def labels(): + """A minimal block in the shape parse_physicell_xml produces.""" + return { + 0: ('ID', 1), + 1: ('position', 3), + 4: ('total_volume', 1), + 5: ('cell_type', 1), + 7: ('current_phase', 1), + 26: ('dead', 1), + 37: ('radius', 1), + 50: ('is_motile', 1), + 52: ('migration_speed', 1), + 86: ('phenotype_state', 1), + } + + +@pytest.fixture +def cell_matrix(): + """A (87, 40) matrix with each cell's values keyed off its column index.""" + matrix = np.zeros((N_VARS, 40), dtype=float) + for cell in range(40): + matrix[:, cell] = np.arange(N_VARS) * 1000.0 + cell + return matrix + + +def _write(tmp_path, name, matrix): + path = tmp_path / name + savemat(str(path), {'cells': matrix}) + return path + + +class TestLabelExpansion: + """Expansion of the MultiCellDS block into flat columns.""" + + def test_scalar_labels_keep_their_name(self): + assert expand_cell_labels({0: ('ID', 1)}) == {0: 'ID'} + + def test_size_three_expands_to_xyz(self): + assert expand_cell_labels({1: ('position', 3)}) == { + 1: 'position_x', 2: 'position_y', 3: 'position_z' + } + + def test_other_vector_sizes_expand_positionally(self): + assert expand_cell_labels({5: ('state', 2)}) == { + 5: 'state_0', 6: 'state_1' + } + + def test_variable_count_accounts_for_trailing_vector_width(self): + # max(labels) + 1 would undercount this by two. + assert declared_variable_count({0: ('ID', 1), 1: ('position', 3)}) == 4 + + +class TestOrientation: + """Issue #2: orientation must come from the label count, not magnitude.""" + + def test_fewer_cells_than_variables(self, tmp_path, labels, cell_matrix): + """An 87x40 matrix with 87 declared labels parses as 40 cells.""" + path = _write(tmp_path, 'small.mat', cell_matrix) + result = parse_cells_mat(path, labels=labels) + + assert result['positions'].shape[0] == 40 + np.testing.assert_array_equal( + result['positions'][:, 0], cell_matrix[1, :] + ) + + def test_genuine_on_disk_transpose(self, tmp_path, labels, cell_matrix): + """A 40x87 matrix also parses as 40 cells.""" + path = _write(tmp_path, 'transposed.mat', cell_matrix.T) + result = parse_cells_mat(path, labels=labels) + + assert result['positions'].shape[0] == 40 + np.testing.assert_array_equal( + result['positions'][:, 0], cell_matrix[1, :] + ) + assert result['orientation'] == 'cells_x_variables' + + def test_square_matrix_prefers_variables_on_axis_zero(self, tmp_path, labels): + square = np.zeros((N_VARS, N_VARS), dtype=float) + for cell in range(N_VARS): + square[:, cell] = np.arange(N_VARS) * 1000.0 + cell + + result = parse_cells_mat(_write(tmp_path, 'sq.mat', square), labels=labels) + + assert result['positions'].shape[0] == N_VARS + np.testing.assert_array_equal(result['positions'][:, 0], square[1, :]) + assert result['orientation'] == 'variables_x_cells' + + def test_neither_orientation_raises(self, tmp_path, labels): + path = _write(tmp_path, 'bad.mat', np.zeros((99, 33))) + + with pytest.raises(ValueError, match='matches neither orientation'): + parse_cells_mat(path, labels=labels) + + def test_raw_data_round_trips(self, tmp_path, labels, cell_matrix): + """raw_data is what loadmat returned, not the reoriented view.""" + path = _write(tmp_path, 'raw.mat', cell_matrix) + result = parse_cells_mat(path, labels=labels) + + np.testing.assert_array_equal(result['raw_data'], loadmat(str(path))['cells']) + assert result['orientation'] == 'variables_x_cells' + + def test_warns_when_orientation_must_be_guessed(self, tmp_path, cell_matrix): + path = _write(tmp_path, 'noguess.mat', cell_matrix) + + with pytest.warns(UserWarning, match='Guessing cell matrix orientation'): + parse_cells_mat(path) + + +class TestLabelledColumns: + """Issue #1: every labelled column must be reachable.""" + + def test_non_indexed_standard_fields(self, tmp_path, labels, cell_matrix): + result = parse_cells_mat( + _write(tmp_path, 'cols.mat', cell_matrix), labels=labels + ) + + np.testing.assert_array_equal( + result['columns']['is_motile'], cell_matrix[50, :] + ) + np.testing.assert_array_equal( + result['columns']['migration_speed'], cell_matrix[52, :] + ) + + def test_model_specific_custom_column(self, tmp_path, labels, cell_matrix): + result = parse_cells_mat( + _write(tmp_path, 'custom.mat', cell_matrix), labels=labels + ) + + np.testing.assert_array_equal( + result['columns']['phenotype_state'], cell_matrix[86, :] + ) + + def test_differing_label_counts_do_not_leak(self, tmp_path, labels, cell_matrix): + """Two models parsed in one session each resolve from their own XML.""" + model_b_labels = dict(labels) + model_b_labels[87] = ('hif', 1) + model_b_labels[88] = ('hypoxia_timer', 1) + model_b_matrix = np.vstack([ + cell_matrix, + np.full((2, cell_matrix.shape[1]), 7.0), + ]) + + result_a = parse_cells_mat( + _write(tmp_path, 'a.mat', cell_matrix), labels=labels + ) + result_b = parse_cells_mat( + _write(tmp_path, 'b.mat', model_b_matrix), labels=model_b_labels + ) + + assert 'hif' not in result_a['columns'] + assert 'hif' in result_b['columns'] + # Index 86 carries a different meaning in each model. + assert result_a['positions'].shape[0] == result_b['positions'].shape[0] + + def test_no_labels_preserves_legacy_output(self, tmp_path): + """With labels=None, output matches v0.2.0 for a default-layout frame.""" + wide = np.zeros((N_VARS, 200), dtype=float) + for cell in range(200): + wide[:, cell] = np.arange(N_VARS) * 1000.0 + cell + path = _write(tmp_path, 'wide.mat', wide) + + result = parse_cells_mat(path) + + assert result['columns'] == {} + np.testing.assert_array_equal(result['positions'][:, 0], wide[1, :]) + np.testing.assert_array_equal(result['dead_flags'], wide[26, :].astype(int)) + + +class TestMicroenvironmentGuard: + """Issue #6: validate the microenvironment matrix orientation.""" + + def test_correct_orientation_parses(self, tmp_path): + me = np.arange(6 * 50, dtype=float).reshape(6, 50) + path = tmp_path / 'me.mat' + savemat(str(path), {'multiscale_microenvironment': me}) + + result = parse_microenvironment_mat(path, ['oxygen', 'glucose']) + + assert result['voxel_positions'].shape == (50, 3) + np.testing.assert_array_equal(result['concentrations']['oxygen'], me[4, :]) + + def test_transposed_input_is_corrected(self, tmp_path): + me = np.arange(6 * 50, dtype=float).reshape(6, 50) + path = tmp_path / 'me_t.mat' + savemat(str(path), {'multiscale_microenvironment': me.T}) + + result = parse_microenvironment_mat(path, ['oxygen', 'glucose']) + + np.testing.assert_array_equal(result['concentrations']['oxygen'], me[4, :]) + + def test_mismatched_substrate_count_raises(self, tmp_path): + path = tmp_path / 'me_bad.mat' + savemat(str(path), {'multiscale_microenvironment': np.zeros((9, 11))}) + + with pytest.raises(ValueError, match='matches neither orientation'): + parse_microenvironment_mat(path, ['oxygen', 'glucose']) + + +class TestToDataFrameDeadCells: + """Issue #4: to_dataframe must honour include_dead_cells.""" + + @pytest.fixture + def timestep_with_dead(self, example_physicell_dir): + for _, xml_path, _ in discover_physicell_timesteps(example_physicell_dir): + timestep = read_physicell_timestep(xml_path) + if timestep.n_dead_cells > 0: + return timestep + pytest.skip("No frame with dead cells in the example simulation") + + def test_excludes_dead_by_default(self, timestep_with_dead): + df = timestep_with_dead.to_dataframe() + + assert len(df) == timestep_with_dead.n_cells + assert not df['is_dead'].any() + + def test_includes_dead_when_requested(self, timestep_with_dead): + df = timestep_with_dead.to_dataframe(include_dead_cells=True) + + assert len(df) == timestep_with_dead.n_cells_total + + def test_row_order_matches_positions(self, timestep_with_dead): + df = timestep_with_dead.to_dataframe() + + np.testing.assert_allclose( + df[['x', 'y', 'z']].to_numpy(), timestep_with_dead.positions + ) + + def test_extra_columns_are_opt_in(self, timestep_with_dead): + base = timestep_with_dead.to_dataframe() + extended = timestep_with_dead.to_dataframe(extra_columns=True) + + assert extended.shape[1] > base.shape[1] + assert 'is_motile' in extended.columns + assert len(extended) == len(base) + + def test_no_fragmentation_warning(self, timestep_with_dead): + with warnings.catch_warnings(record=True) as record: + warnings.simplefilter("always") + timestep_with_dead.to_dataframe(extra_columns=True) + + assert not [ + w for w in record + if issubclass(w.category, pd.errors.PerformanceWarning) + ] + + +class TestReaderUsesFrameLabels: + """The reader resolves columns from each frame's own XML.""" + + def test_custom_column_reachable_through_reader(self, example_physicell_dir): + _, xml_path, _ = discover_physicell_timesteps(example_physicell_dir)[0] + timestep = read_physicell_timestep(xml_path) + + labels = parse_physicell_xml(xml_path).extra['custom_labels'] + expected = expand_cell_labels(labels) + + data = timestep._load_cell_data() + assert data['columns'], "reader should pass labels through to the parser" + assert set(data['columns']) == set(expected.values()) + + def test_reader_emits_no_orientation_warning(self, example_physicell_dir): + _, xml_path, _ = discover_physicell_timesteps(example_physicell_dir)[0] + + with warnings.catch_warnings(record=True) as record: + warnings.simplefilter("always") + read_physicell_timestep(xml_path)._load_cell_data() + + assert not [ + w for w in record + if 'Guessing cell matrix orientation' in str(w.message) + ] + + +class TestVersionSingleSource: + """Issue #5: __version__ must not drift from distribution metadata.""" + + def test_version_matches_package_metadata(self): + assert spatialtissuepy.__version__ == importlib.metadata.version( + "spatialtissuepy" + )