SPT simulations and convergence¶

The repository-only SPT workflow catalogues supplied XDMF/HDF results, reconstructs lobules, serializes quantitative analysis, and generates numerical and geometry figures. It requires access to study data, but does not need FEBio, VTK rendering, or a graphics session. The reusable library remains independent of these study workflows. This analyzes existing solver results; it does not run new simulations.

Input and output directories¶

Keep the network dataset separate from local outputs. The current dataset is available at /data/koeniglab/spt, with this structure:

spt/
├── simulations_sixth/
│   ├── material/
│   │   ├── feb_manipulation_v3.py
│   │   └── lobule_zonation_pattern_bc.feb.template
│   └── xdmf/
│       ├── sim001.xdmf
│       ├── sim001.h5
│       ├── sim001_cells.nc
│       ├── sim001_points.nc
│       └── ...
└── convergence_sixth/
    ├── material/
    │   └── *convergence_*.feb
    └── xdmf/
        ├── cbc_convergence_06250.xdmf
        ├── cbc_convergence_06250.h5
        └── ...

Both names carry the study design. The generation script is versioned, the simulation template names the boundary conditions it was built for, and a convergence run is prefixed with the variant it belongs to. The workflow discovers the script and the template rather than expecting fixed names, and identifies a convergence mesh by the resolution after its last convergence_ marker, so a renamed study is read without changing code.

Use the development environment, or sync just research dependencies with uv sync --locked --group research --no-dev. Commands run from the checkout root:

To run the complete project from catalogue through analysis, figures and manuscript figures, use:

export POROUS_MEDIA_SPT_DATA=/data/koeniglab/spt
export MPLBACKEND=Agg PYVISTA_OFF_SCREEN=true
uv run --no-sync python -m research.spt.pipeline

The complete workflow below describes resuming, inspecting stages, and the results layout.

On PowerShell, set $env:POROUS_MEDIA_SPT_DATA and $env:MPLBACKEND instead. Every command also accepts --input-root /data/koeniglab/spt explicitly. The input default is POROUS_MEDIA_SPT_DATA, falling back to the checkout's data/spt. The output default is POROUS_MEDIA_SPT_RESULTS, falling back to results/ inside the input root, for example /data/koeniglab/spt/results. Results never live inside the repository checkout. Outputs inside a supplied study directory, outputs containing the input root and symlinked output subdirectories are rejected. No command changes the supplied study directories.

Results layout¶

The output root separates deliverables from intermediate data. Every report refresh also writes this overview as README.md:

results/
├── index.html           offline report; open this first
├── README.md            folder guide
├── manuscript/
│   ├── exports/         stable figure-04.svg ... figure-s1-convergence.svg, tissue panels
│   └── <figure>/        versioned generation of each manuscript figure
├── analysis/<study>/    quantitative tables, NetCDF, figures and provenance
├── figures/
│   ├── runs/<study>/<run_id>/   geometry time course of one history
│   └── overviews/<study>/       final-state overview of a study
├── meshes/<study>/<run_id>/     reconstructed lobules and mixed-sector histories
├── catalog/             runs.tsv, optional runs.xlsx, condition figures
├── logs/                pipeline.json and per-stage reports and catalogues
└── cache/
    ├── cells/<study>/<run_id>/      numerical cell, point and summary caches
    ├── analysis/<study>/<run_id>/   per-run analysis results
    ├── plots/<study>/<run_id>/      sampled plot geometry and values
    ├── depth/<study>/<run_id>/      native cap heights for tilted views
    └── previews/                    report thumbnails

Each artifact directory holds a manifest.json naming its current generation-<id>/ folder; previous generations are removed after a successful update. Unused report previews are pruned whenever the report is rebuilt.

Storage profile¶

A complete results tree of 48 simulations and 7 convergence meshes is about 7 GB. The equivalent tree written before these settings was 97 GB. Each setting and its measured effect:

Stacking is what makes the compression effective rather than nominal. Written one step at a time, each array is a chunk of its own and the compressor sees a single frame; a chunk spanning consecutive steps exposes how little changes between them. meshio cannot read HyperSlab references, so these sidecars are read with porous_media.data.xdmf_tools.XDMFTimeSeriesReader, which resolves them and otherwise behaves exactly like meshio's reader. ParaView reads them directly.

The cell histories are stored as float32 by the source solver, so nothing here reduces numerical precision. Everything under cache/ is rebuilt by rerunning the pipeline; meshes/, analysis/, figures/ and manuscript/ are the deliverables. Reconstructing with --all-fields, or requesting SVG or PDF figures, raises these numbers; --all-fields multiplies mesh storage by about five.

Compression and the mesh field selection are recorded in each artifact's manifest.json, and both participate in the stale check: changing either rebuilds the artifacts that depend on it rather than mixing layouts within one tree.

Inspect conditions and available runs¶

information.py reads descriptor and model metadata without loading HDF arrays. After processing, it reuses catalogue metadata when source fingerprints match. Use --refresh to reread descriptors explicitly. It discovers actual convergence IDs and records cell/point counts, available scalar fields, time ranges in seconds, and completion relative to the model's explicit control duration. Models with ambiguous or absent duration metadata have an unknown completion status.

Simulation IDs follow the pattern/flow ordering in the supplied generator. The command parses literal assignments without executing the archived Python script. If no generator is supplied, it labels the repository defaults as unverified. Flow units remain unspecified when the source does not declare them. Convergence filename tokens remain identifiers; no mesh-length units are inferred from them.

uv run --no-sync python -m research.spt.information --dry-run
uv run --no-sync python -m research.spt.information --plots --excel

The first command writes nothing. The second writes catalog/runs.tsv, optional catalog/runs.xlsx, and the condition figures catalog/boundary-flows.png and catalog/condition-colors.png. Figures are saved without opening windows and are closed after use. Metadata failures and missing expected simulations are reported and produce a nonzero exit status.

The dataset inspected on 2026-11-12 contains 56 simulation runs, seven zonation patterns by eight boundary flows, and six convergence runs. cbc_convergence_06250 ends at 53,027.1 s, earlier than the model duration of 86,400 s, and its own model declares no duration, so it is catalogued with an unknown complete. Comparisons must choose a shared time window or explicitly handle different end times.

Pattern 6 is the CYP2E1 zonation pattern. It rises from periportal to perivenous like the linear pattern, but smoothly: measured on sim049 at t = 0, relative protein runs from 0.013 to 0.996 of its range, is flat periportally, steepest at about 0.6 of the position and saturated perivenously. It carries the YlOrBr palette, the row label G CYP2E1 and a sigmoidal profile icon. A pattern index this repository does not name at all is still catalogued and plotted, as pattern <index>, rather than dropped or given an invented label.

A lobule is six rotations of one source sector, so the mixed-pattern comparison view holds six sources and no more. It compares every pattern beside CYP2E1, which is the measured profile rather than one of the synthetic ones. Every pattern still has its own lobule and its own row in the study and manuscript figures. A study that leaves more than six patterns after that exclusion fails rather than silently dropping one, because the sector count fixes how many fit. --no-combinations skips the mixed views entirely.

Build or resume local caches¶

data_processing.py uses existing XDMF/HDF results directly. Its default cache step reads the source sixths. Explicit reconstruction steps prepare full lobules for subsequent analysis, as described below. VTK conversion remains separate.

uv run --no-sync python -m research.spt.data_processing --status
uv run --no-sync python -m research.spt.data_processing --workers 2
uv run --no-sync python -m research.spt.data_processing \
  --study convergence --runs cbc_convergence_06250 --workers 1

--status and --dry-run inspect metadata and cache state without writing. --study accepts both, simulations, or convergence. --runs takes explicit run IDs; unknown selections and missing/empty study directories fail clearly. A complete simulation selection requires every expected run. Use --runs for an intentionally partial study.

The default analysis profile retains the five analysis scalars, cell volumes, position, fluid/solid volumes, pressure, zonation pattern, and cell type. It omits unused point arrays. --profile full retains every scalar cell and point field. Stored field precision is preserved, with type promotion where needed.

First-time adoption of the supplied NetCDF caches validates time coordinates, field dimensions, and first/middle/last HDF frames for every selected field. The manifest labels this as sample validation, not a complete comparison of all values. Choose a stronger policy when needed:

# Compare every selected cached value against HDF before adoption.
uv run --no-sync python -m research.spt.data_processing --legacy-cache full

# Read selected fields directly from HDF, bypassing supplied NetCDF caches.
uv run --no-sync python -m research.spt.data_processing --legacy-cache ignore

# Rebuild a selected local cache from HDF even if its fingerprints match.
uv run --no-sync python -m research.spt.data_processing --runs sim001 --force

Rejected legacy caches fall back to HDF and the report records the reason. Subsequent runs check input sizes and nanosecond modification times, processing options, conditions, and output fingerprints. Unchanged runs do not reopen XDMF or numerical arrays. Changed sources rebuild from HDF, avoiding reuse of stale legacy caches. Changes that deliberately preserve both file size and modification time require --force; metadata fingerprints are not content hashes.

Two spawned workers are the default. --workers 1 runs serially. Each worker owns one run's files and memory. Tune the worker count for the network share rather than the machine's CPU count. An unchanged batch avoids worker startup.

Outputs and failure recovery¶

Each run publishes a complete generation under:

results/cache/cells/<study>/<run_id>/
├── manifest.json
└── generation-<identifier>/
    ├── cells.nc
    ├── points.nc
    ├── summary.nc
    └── limits.json

The manifest identifies the active generation and records provenance and the legacy validation method. Scalar limits are calculated from the selected data without rereading the source. An OS lock prevents concurrent readers/writers from replacing a run during a load. Partial writes never replace the active manifest. A failed run leaves previously completed outputs intact, and the next invocation resumes other completed runs. Source changes during processing cause failure and require a retry after the upload completes.

logs/processing.json records all selected successes and failures, while logs/processing-runs.tsv records the selected successful metadata. These describe the latest invocation and do not overwrite the separate information catalogue. Any failed run gives the processing command a nonzero exit status. Short but readable histories are reported separately from processing failures.

Cache APIs load data into memory and close file handles:

from pathlib import Path

from research.spt.cache import load_processed_run, load_summary
from research.spt.catalog import discover_runs

source = discover_runs(Path("/data/koeniglab/spt"), run_ids=("sim001",))[0]
output = Path("/data/koeniglab/spt/results")
summary = load_summary(source, output)  # Reads only the compact summary.
cells, points, summary = load_processed_run(source, output)

Summaries contain native-time volume-weighted population means and standard deviations, accumulated in float64. Invalid volumes and nonfinite analysis values fail explicitly. Per-cell arrays retain the original time axis for spatial analysis and interpolation. Interpolating a summary is not generally equivalent to interpolating cell values and changing volumes before reduction.

Further analysis and visualization¶

plots_simulations.py and plots_convergence.py are entry points to the shared quantitative analysis below. They use actual run IDs and local homogeneous-lobule caches. Geometry plotting uses its own shared renderer. Numerical cache generation does not invoke either plotting workflow or import PyVista.

Lobule reconstruction is an explicit processing step described below. VTK conversion remains a separate helper for studies without XDMF results. PyVista snapshots and video export require a working rendering setup; the sector comparison figures generated below use Matplotlib and run headlessly. For a portable numerical example, see the library quickstart.

Reconstruct lobules and sector comparisons¶

Create the modified geometries after inspecting the supplied sixths. This step reconstructs existing simulation results by symmetry; it does not run FEBio. The complete workflow then builds numerical caches from the homogeneous lobules:

uv run --no-sync python -m research.spt.data_processing --step all --workers 2

Individual steps can be run or resumed independently:

# Reconstruct every simulation and convergence mesh, preserving native times.
uv run --no-sync python -m research.spt.data_processing --step lobulus

# Compose the mixed convergence view and eight comparisons of zonation patterns.
# This also saves a PNG comparison figure alongside the time series.
uv run --no-sync python -m research.spt.data_processing --step combinations

# Build both sets of geometries and figures, without numerical cache creation.
uv run --no-sync python -m research.spt.data_processing --step geometry

# Analyze reconstructed lobules using the same scalar profiles and cache policy.
uv run --no-sync python -m research.spt.data_processing \
  --step cache --geometry lobulus

# Inspect reconstruction status without writing or starting workers.
uv run --no-sync python -m research.spt.data_processing --step geometry --status

The existing default, --step cache --geometry sixth, remains available for inspecting source-sector results. Subsequent full-lobule analysis should use --geometry lobulus or --step all. Mixed comparisons are excluded from this quantitative analysis path. A selected reconstruction can use --step lobulus --runs sim001; combinations require complete studies and reject --runs. --study simulations and --study convergence work for every step.

The reconstruction writes:

Each lobule places six copies of its source sector counterclockwise around the origin, reflecting alternating sectors to pair matching radial boundaries. Connectivity permutations preserve each source cell's orientation. Scalar fields are copied, Cartesian vectors are transformed with Q, and full stress/strain tensors with Q T Q.T. The Cartesian scalar fields position_x, position_y, and position_z transform together. Unsupported component layouts fail explicitly. Cell volume is preserved per cell; a homogeneous lobule has six times its source's volume and the same volume-weighted means and standard deviations.

Seam points remain separate. This preserves different values on adjoining patterns and permits different convergence resolutions without inventing node correspondence or averaging discontinuities. sector_id identifies sectors 0 through 5, and source_cell_id preserves the original cell identity within each source. These are reconstructed visualization/analysis geometries, not conforming meshes intended as new solver inputs.

The mixed convergence view uses 06250, 09375, 12500, 18750, 25000, and 37250 in that order. 50000 still receives its own homogeneous lobule. The finest input ends at 36,590 s, so the mixed history ends there. Native output times are offset between resolutions. The mixed history uses the first sector's native time grid inside the common interval, including the common endpoints. Continuous scalar/vector/tensor fields interpolate linearly. Integer fields, cell_type, rr_zonation_pattern, rr_necrosis, and gauss_points_TPM retain their previous sampled value. No field is extrapolated. Homogeneous lobules retain every original time sample and do not interpolate any field.

Each dose comparison combines patterns 0 through 5 from the same boundary-flow condition. For example, dose 1 uses sim001, sim009, sim017, sim025, sim033, and sim041. Dose numbers identify the generator's boundary-flow ordering; flow units remain unspecified. Each comparison includes the complete common history. Its figure shows initial protein zonation and final necrosis, with numbered sectors and a pattern legend. The convergence figure shows the final toxicant values of the six numbered resolution sectors. Both figures draw sector outlines without internal mesh edges. Figures render headlessly with Matplotlib, without requiring a VTK display backend.

Geometry generations are stored separately from numerical cache generations:

results/meshes/<study>/<run_id>/
├── manifest.json
└── generation-<identifier>/
    ├── <run_id>.xdmf
    ├── <run_id>.h5
    ├── <run_id>_limits.json
    ├── provenance.json
    ├── model.feb          # Homogeneous lobules, when a source model exists.
    └── <run_id>.png       # Mixed comparisons only.

XDMF references its colocated HDF sidecar by relative filename. Copying the pair is sufficient to open it elsewhere; it does not depend on network HDF files. Each worker stages its input descriptors and HDF sidecars on local storage with sequential copies, avoiding repeated small reads over the network filesystem. Temporary copies are removed after readers close, including on failure. Allow space for one source history per lobule worker or six per mixed-view worker, in addition to generated outputs. Original inputs remain the provenance source and are checked for changes before publication. The HDF writer uses gzip level 4 with the byte-shuffle filter and stacks each attribute into one dataset spanning the series, which is what collapses repeated scalar and time-invariant storage. A sidecar keeps only the fields the pipeline reads back, listed as GEOMETRY_CELL_FIELDS in research/spt/fields.py: the display layers, the analysis fields and the categorical fields. A source run also carries mechanical vectors and 3x3 tensors that nothing downstream reads; measured on one 1.76 GB sidecar, stress, Lagrange_strain and solid_stress_TPM alone were 45% of its storage and the vector fields another 16%. Pass --all-fields to reconstruct the complete mechanical state instead, which multiplies mesh storage by about five. Source fingerprints, conditions, field-transform policy, and temporal alignment are recorded in the manifest and provenance file. Only complete generations become visible through the manifest. Failed reconstruction or figure generation preserves the previous generation, and unchanged runs skip array reads and worker startup. logs/geometry.json describes the latest selected reconstruction invocation. logs/geometry-runs.tsv lists successful outputs, including descriptor and figure paths, so readers do not need to navigate generation directories manually.

Resolve published paths rather than globbing temporary or retired generations:

from pathlib import Path

from research.spt.cache import load_summary
from research.spt.geometry import geometry_path, lobulus_sources

raw = Path("/data/koeniglab/spt")
output = raw / "results"
source = lobulus_sources(raw, output, run_ids=("sim001",))[0]
summary = load_summary(source, output)
lobulus_xdmf = source.xdmf
comparison_xdmf = geometry_path(
    output / "meshes" / "simulations_combinations" / "dose_01"
)

Numerical lobulus caches use the same cache/cells/<study>/<run_id>/ layout as source-sector caches, under their respective *_lobulus study names. Use published lobules for quantitative analysis and sector comparisons for visualization. Mixed-sector histories never enter whole-lobule statistics.

Plot sixths, lobules and combined geometries¶

After reconstruction, create all geometry figures with:

uv run --no-sync python -m research.spt.plot_geometry --workers 2

The input, data and output roots default to the same locations used above. --data-root selects the published reconstructed histories; --output-root selects plot storage. Both default to the results directory. The raw share is read only. Plotting does not change reconstructed histories or numerical caches. These repository commands also work independently:

uv run --no-sync python -m research.spt.plots_simulations_geometry
uv run --no-sync python -m research.spt.plots_convergence_geometry

# Render a selected source sixth, its lobule and the same-dose comparison.
uv run --no-sync python -m research.spt.plot_geometry \
  --study simulations --runs sim048 dose_08 --geometry all

# Inspect pending work without writing files or reading numerical arrays.
uv run --no-sync python -m research.spt.plot_geometry --status

--geometry accepts all, sixth, lobulus, or combined. The complete current dataset produces 119 time-course figures: 55 source sixths, 55 homogeneous lobules, eight dose comparisons and one mixed convergence history. Each of the six study/geometry groups also gets a final-state overview. Simulation overviews follow the generator's pattern and dose ordering. No figure shows the mesh itself; only the triangles of manuscript Figure S1 A display cell edges.

Default figures preserve the manuscript's five snapshots: indices 1, 3, 5, 7 and 9 of ten equally spaced times in the common interval. Explicit --times and other snapshot counts override this selection. The five columns retain extracellular substrate and product, protein, toxic compound and necrosis. The canonical layers in research/spt/fields.py preserve magma for concentrations, hot for protein and binary for necrosis, with the original labels and horizontal scalar bars in every panel. Scalar panels omit internal mesh lines and show only the domain outline, including holes. Homogeneous lobule outlines omit artificial sector seams; independent combined sectors retain their own outlines. Layers can opt into transparent cell edges with show_edges=True in research/spt/fields.py. Boundary-only scalar surfaces use a thin stroke in each cell's own color to close vector-export hairlines, without adding neutral mesh lines. Each panel occupies 4 by 4 inches, matching the original 600 pixel panels at 150 DPI. Scalar bars retain five labels, bold text and one decimal place. Final-state overviews use the same panel sizes and per-panel scalar bars. Small ranges use scientific notation when fixed decimals would repeat labels, so nonzero product concentrations do not appear as five identical zeros.

Times share a common interval within each selected study family. Simulation plots therefore retain 24 h even when convergence ends at 36,590 s. Shared color ranges use validated local numerical-cache limits, then full-history _limits.json metadata beside source descriptors, without scanning network arrays again. Missing or malformed legacy metadata falls back to sampled ranges when no validated cache is available. Explicit layer limits take precedence. Combined views identify each sector. These settings are shared with the historical DataLayer adapter so the two workflows cannot silently choose different colors.

Continuous scalars are linearly interpolated between neighboring native frames. Necrosis, zonation, cell type and other discrete fields retain their previous native value. No extrapolation or averaging between cells is performed. Only requested scalar arrays and neighboring frames are read from HDF; point, vector, tensor and unused scalar arrays are omitted. A compact local sample cache is shared by time courses, overviews and optional videos. Plotting never creates interpolated XDMF/HDF histories.

Choose fields, times and outputs explicitly:

# Final-state figures for only necrosis and toxicant.
uv run --no-sync python -m research.spt.plot_geometry \
  --fields rr_necrosis 'rr_(T)' --snapshots 1 --formats png svg

# Explicit times in seconds must lie within every selected history's interval.
uv run --no-sync python -m research.spt.plot_geometry \
  --study simulations --times 0 21600 43200 86400

# Preserve outliers while showing smaller concentrations on a symlog color scale.
# Use a separate output root so the linear default figures remain unchanged.
uv run --no-sync python -m research.spt.plot_geometry \
  --study convergence --scale symlog --output-root /tmp/spt-convergence-symlog

# Generate only the two study overviews per geometry.
uv run --no-sync python -m research.spt.plot_geometry --only-overviews

# Separate field/time panels are optional, as are videos and GIFs.
uv run --no-sync python -m research.spt.plot_geometry \
  --study simulations --runs dose_08 --geometry combined \
  --panels --video-frames 200 --fps 12 --gif

All geometry figures share an orthographic camera tilted 20 degrees from above, with an azimuth of -60 degrees. Native cap heights and exterior side faces show the mesh thickness without exaggeration. Caps and walls use their owning cell's unchanged scalar value and the canonical colormap; lighting does not alter the colors. The camera applies consistently to sixths, lobules, combined sectors, manuscript figures and optional videos. Quantitative plots and their statistics are unaffected. Colored convergence boundaries and sector numbers follow the same projection, and cropping includes the visible thickness.

SVG and PDF preserve vector cell geometry and text. PNG previews default to 150 DPI; --dpi also controls optional video resolution. FFmpeg is required only when videos are requested. Video artists are reused and frames stream directly to FFmpeg, without writing intermediate PNG panels. GIF conversion is opt-in. --no-overviews omits study overviews. Bounded spawned workers own their files and rendering state; --workers 1 runs serially.

Per-run and study figures default to PNG alone because a dense vector time course dwarfs its raster preview: across one results tree the per-run SVG files averaged 32 MB against 1.0 MB for the same figures as PNG, and SVG accounted for 88% of all figure storage. Request --formats png svg for the few figures that need cell-level vector editing, and --formats png pdf when a vector file has to stay small. Selecting fewer fields or snapshots also reduces rendering time and storage without changing the underlying histories.

The shared Matplotlib renderer runs headlessly without importing PyVista or needing a graphics session. It supports planar XY triangles/quads and vertical wedge/hexahedron extrusions with horizontal caps. Planar elements tilt without inventing thickness. Other geometries fail clearly instead of producing an incorrect projection; use the library's 3D visualization helpers for those. Historical visualize_spt_2d calls remain available, including generator inputs, custom layer labels, unit conversions and fixed color limits. Separate panels and videos now require explicit selection in that helper too.

Outputs and preview caches use the following layout:

results/
├── cache/plots/<study>/<run_id>/manifest.json
├── cache/depth/<study>/<run_id>/manifest.json
├── cache/previews/<fingerprint>.webp
├── figures/runs/<study>/<run_id>/manifest.json
├── figures/overviews/<study>/manifest.json
├── logs/plots-both.tsv
└── logs/plots-both.json

Each manifest points to a complete generation. A failed replacement preserves the previous published figures. Source and output fingerprints enable unchanged runs to skip XML parsing, numerical reads, rendering and worker startup. Changing only a rendering option reuses the existing sampled arrays. New sample caches retain native cap heights. Older planar caches are adopted without resampling scalar histories: one geometry-only read stores the heights separately in cache/depth/. Later rendering changes reuse this small local cache. Mesh-view settings are included in figure provenance and invalidate geometry figure caches. Missing artifacts are regenerated; --force explicitly resamples and rerenders. Fingerprints use size and modification time, so changes preserving both require --force. The catalogue gives the current paths of figures and any requested videos; the final HTML report links them. The plot report records errors, sampled color ranges, elapsed time and the number of selected HDF field arrays read. Files are named by the selected study (both, simulations, or convergence) and describe the latest invocation, including an intentionally partial selection.

--scale symlog is an explicit alternative for concentration ranges containing large outliers. It preserves every value, with a linear region within one millionth of the maximum absolute display limit and logarithmic spacing outside that region. Figures label the transformed scale; necrosis stays linear. The current convergence_50000 has sampled substrate values above 2,000 mM, whereas the six finer histories stay near 1 mM or below. Its product values likewise stretch the shared linear range. Keep those source values visible; use a separate symlog output root when inspecting finer spatial differences.

Time courses, spatial profiles and convergence¶

After numerical cache generation, run both study analyses or select one:

uv run --no-sync python -m research.spt.analysis
uv run --no-sync python -m research.spt.plots_simulations
uv run --no-sync python -m research.spt.plots_convergence

# Change spatial resolution or explicit comparison times, in seconds.
uv run --no-sync python -m research.spt.analysis \
  --bins 30 --spatial-times 0 21600 36000 --time-points 101

# Analyze an explicit subset, retaining its source-derived condition labels.
uv run --no-sync python -m research.spt.analysis \
  --study simulations --runs sim001 sim009 --workers 1

The default --geometry lobulus analyzes every simulation lobule and every convergence lobule. --geometry sixth explicitly analyzes source-sector numerical caches instead. Mixed geometries are excluded. --data-root selects the numerical-cache root; --output-root selects local analysis storage. Neither the raw share nor reconstructed simulation histories are modified.

Manuscript time-course figures use matched-time cell-volume-weighted means. Simulation panels retain the pattern palettes (Greys, Blues, Purples, Reds, Oranges, Greens, YlOrBr), dose shades, black-edged circular markers and population SD error bars. --no-spread explicitly omits the bars; --spread also enables them for convergence. SD describes cell variation, not uncertainty in a mean or a confidence interval. Convergence retains a coarse-to-fine viridis palette, plain lines, the 1.5 h reference line and the legend above the panels. Actual mesh IDs replace obsolete IDs from the archived scripts. Colors remain stable when selecting a subset of doses or meshes.

Panels retain bold 12 point labels, 2.7 by 2.5 inch dimensions, the original axis labels and PP/PV ticks, and 300 DPI PNG exports. The native histories, including every original output and each run's actual endpoint, remain available in native-timecourses.tsv and the per-run native.nc caches.

For comparisons, 51 equally spaced times cover the common interval separately within each selected study family. Continuous cell values and cell volumes are interpolated before weighted reduction. This differs from interpolating summary means when volumes change. Discrete fields retain their previous native value, and no time is extrapolated. Source necrosis fractions are preserved; figures display their volume-weighted percentage while serialized values remain fractions. The cell data are read selectively from local NetCDF caches, without loading complete histories, point arrays, or network HDF data.

Spatial plots default to the first common time, the historical reference time (10 h for simulations and 1.5 h for convergence), and the last common time. If the reference lies outside the common interval, its midpoint is used instead. Figures show individual cell scatters rather than binned curves, including valid fractional necrosis. For quantitative comparisons, twenty equal bins cover the source's normalized rr_position, from periportal 0 to perivenous 1. Each bin stores volume-weighted mean, population SD, minimum, maximum, cell count, and total volume. Position 1 belongs to the last bin. Empty or zero-volume bins retain missing statistics and zero volume; they are never filled by spatial interpolation. Raw cell snapshots are also serialized so later analysis can choose a different reduction.

The initial-properties figure always uses each run's first native frame, even when explicit spatial comparison times start later. Simulation diagnostics show the highest selected dose per pattern; convergence diagnostics show every mesh. Diagnostics retain the manuscript's ten-bin position count histogram and individual cell scatters. Display conversions follow the archived script: mesh cell volume uses m³ to nL (1e12), and fluid/solid volume uses L to nL (1e9). Pressure retains the original Pressure [?] label. Serialized arrays remain in source units; caches do not provide independent unit metadata. The former use of the last frame for initial diagnostics remains corrected.

Simulation figures compare doses by pattern and include matched-time dose responses. Convergence figures use the actual mesh IDs and lobulus cell counts. Resolution panels retain a linear element-count axis, black connecting line and mesh-colored markers at 1.5 h. If that time is absent from explicitly selected snapshots, the middle selected snapshot is used and the report names its time. The mesh with the largest cell count supplies the reference. Exported tables report absolute and relative differences of whole-lobule means at matched times. A zero reference makes relative differences undefined, stored as missing values. Spatial comparisons report RMS differences of bin means, weighted by reference bin volume, with the fraction of reference volume covered by mutually populated bins. These are differences from the finest available mesh, not estimates of error against a continuum solution or an inferred convergence order.

Convergence includes supplemental symlog figures to expose smaller values beside the current coarsest mesh's large substrate/product outliers. --linear-only omits those supplements. All values remain in the serialized results. Default exports are PNG; --formats, --dpi, and --fields select figure output.

Reusable analysis results¶

Each run publishes a compact numerical generation under cache/analysis/<study>/<run_id>/:

Each study publishes an independent generation under analysis/<study>/, containing analysis.nc, native/matched time-course TSV tables, spatial-profile TSV tables, figures, a local HTML gallery, and provenance. Convergence also includes convergence-means.tsv and convergence-profiles.tsv. The study NetCDF combines matched statistics and profiles with run, dose, pattern, cell-count, and native-end-time coordinates; native histories with differing time axes are preserved in the native TSV table and per-run files. No pickle files are used. The study tables and NetCDF files can be moved and read without the raw share.

The figures appear in the analysis section of the offline report. logs/analysis-<study>.tsv supplies current result paths, and logs/analysis-<study>.json records errors, common windows, elapsed time, and the number of native frames read. These describe the latest selection. Unchanged invocations inspect fingerprints without opening numerical arrays or starting workers. Plot styling or format changes reuse numerical reductions. Missing exports rebuild, failed publication retains previous figures, and --force explicitly regenerates the selected analysis. --status writes nothing.

Load a validated per-run result from a checkout:

from pathlib import Path

from research.spt.analysis_data import load_analysis
from research.spt.plot_data import plot_sources

raw = Path("/data/koeniglab/spt")
output = raw / "results"
source = plot_sources(
    raw, output, study="simulations", geometry="lobulus", run_ids=("sim001",)
)[0]
with load_analysis(source, output, "profiles") as profiles:
    print(profiles.sel(field="rr_necrosis"))

Run and verify the complete project¶

uv run --no-sync python -m research.spt.pipeline --status
uv run --no-sync python -m research.spt.pipeline --workers 2

The pipeline runs these stages in order:

1. Catalogue actual inputs and plot simulation conditions.
2. Reconstruct homogeneous lobules and mixed visualization geometries.
3. Build or resume the homogeneous-lobule numerical caches.
4. Serialize quantitative analysis and create its figures.
5. Create or resume sixth, lobulus, and combined geometry figures.
6. Create the manuscript figures and the tissue panels.

It checks every selected result, stops dependent work after a failure, and writes logs/pipeline.json plus the final index.html report at the output root. The JSON names each stage and records its outcome and elapsed time. Repeating the command is a supported resume operation. Existing numerical and geometry outputs are reused only when their provenance and output fingerprints match; they are not rebuilt merely to verify the pipeline. Short source histories are listed separately in truncated_runs; a successful pipeline does not imply that every source solver run reached its model duration.

--study selects a study, and --runs selects an explicit subset. A subset omits mixed comparisons and manuscript figures, which require complete studies. --no-combinations skips the mixed comparison lobules, which place one source per sector and therefore need exactly six zonation patterns or six resolutions; use it for a study whose sources do not fill one lobule. --no-geometry-figures skips the last two stages when only numerical results are needed. --formats selects the per-run and study figure formats and --manuscript-formats those of the manuscript deliverables, while --geometry-snapshots changes the geometry figure panel count. Inspecting a fresh output root reports downstream stages as pending without creating directories or loading HDF arrays.

Portable regressions run every stage from fresh synthetic raw files, including mixed convergence meshes, installed-library imports, serialized results and figures. Run them with uv run --no-sync pytest -q research/tests. Validation of the private current dataset uses the same pipeline against the read-only share.

Shared workflow code¶

SPT modules share field definitions in fields.py, filesystem publication and locking in artifacts.py, and bounded job execution in execution.py. Numerical caches, reconstructed geometries, analysis, and plotting use those same helpers. Serial and parallel batches isolate job failures consistently, and empty batches or a single pending history do not start worker processes. Field ordering, manifest schemas, existing cache generations, and command names remain compatible.

The four plots_* command modules delegate to the shared analysis or geometry workflow. Historical plotting adapters still serve other research projects and delegate to the current renderer. These helpers remain repository-only research code, outside the public library and release artifacts.

Manuscript figures¶

The complete pipeline creates the updated Figures 4-7 of the SPT manuscript, a convergence supplement and the tissue panels in its last stage. Recreate only the manuscript figures after the pipeline has finished with:

MPLBACKEND=Agg PYVISTA_OFF_SCREEN=true uv run --no-sync \
  python -m research.spt.manuscript

This step reads validated local NetCDF summaries, selected cell frames and cached projected meshes. It does not require the raw network drive, rerun reductions or alter existing plots. Unchanged figures use metadata-only cache checks. Missing exports are repaired; --force regenerates selected manuscript figures.

Figures 4 and 5 explicitly sample 0, 6, 12, 18 and 24 h, independently of the historical overview snapshot selection. Labels are derived from the actual times; a nonzero overview frame is never relabeled as zero. Continuous cell fields are interpolated locally and necrosis holds its preceding native frame. Shared scalar limits use full histories of the entire simulation cohort, preserving the mesh palettes. Necrosis is a cell state, so every map shows it on the physical 0-1 scale and float32 round-off just below zero never reaches a colorbar label. Maps of Figures 4, 5 and 7 show lobule outlines without internal mesh edges; outline widths scale with panel size. Only the source meshes in Figure S1 A show their triangles. Figure 6 reuses the existing volume-weighted means, population SD, markers and pattern/dose colors. SD represents spatial heterogeneity, rather than uncertainty across replicate simulations.

Time runs left to right in Figures 4 and 5, with columns for 0, 6, 12, 18 and 24 h. Figure 4 retains all five fields as rows, each with one vertical colorbar on the right. Figure 5 shows only protein, toxic compound and necrosis. Its three bottom colorbars serve every pattern block. Figure 7 uses two swatches beside the complete grid: black means necrosis and white means no necrosis. Its binary legend has no intermediate values. Individual map panels have no colorbars. A simulation whose solver stops before 24 h keeps its place in the pattern-by-dose grid as a dashed panel naming the run and its actual end time; nothing is extrapolated, and the caption and manifest list the affected runs, so a completed rerun invalidates the cached figure. The field palettes and full-cohort limits remain unchanged. Geometry axes use their cropped domain extents, including white living tissue; smaller canvases and closely packed grids reduce whitespace around the meshes.

All manuscript figures use dose labels D1-D8. The representative dose in Figures 4 and 5 is D6. Pattern names and lettering match across Figures 4-7: A Constant, B Linear increase, C Sharp perivenous, D Linear decrease, E Sharp periportal, F Random and G CYP2E1. The row letters and the cohort size follow PATTERN_NAMES in research/spt/manuscript_data.py, so naming a further pattern extends every figure rather than one of them. Figures 6 and 7 share the same row-label format and typography; Figure 5 uses those labels as its pattern-block headings. Every pattern label also carries the same small square black-and-white vector profile panel, including the Figure 4 heading. The horizontal direction runs from periportal (left) to perivenous (right); the vertical direction shows relative initial protein abundance. Flat, rising, falling, sharp-step, sigmoidal and random profiles distinguish the patterns. These schematics explain pattern identity; the random profile is illustrative.

Figures 4, 5, 7 and S1 use the same slightly tilted orthographic geometry view as the study figures, including native thickness and exterior side faces. Figure 6 retains its time-course layout and styling. Captions and provenance record the view, and manuscript scalar ranges, binary necrosis states, pattern icons and dose labels stay unchanged.

The convergence supplement retains every available native time point and actual end time. The 06250 and 09375 histories stop at 10.16 and 21.58 h; no curve is extended to 24 h. The common comparison window ends at 10.16 h. A single row of linear time courses uses the six approved meshes 37250, 25000, 18750, 12500, 09375 and 06250, ordered coarse to fine. Each resolution has a distinct qualitative color, shared by its time courses, source-mesh boundary and mixed-sector boundaries. Source sixths are selected from local lobule meshes, so producing this figure does not reopen the original network meshes. The mixed view numbers those same six resolutions clockwise from upper left, coarse to fine. This display rotation preserves cell values and the archived composition; it applies to the Panel A mesh view and all five readouts. Panel B compares substrate S, product P, protein, toxic compound T and necrosis at 5400 s (1.5 h), the existing convergence reference time. All six sectors use exactly this time, native-frame continuous interpolation and preceding-frame necrosis states. The scalar values come from local homogeneous-lobule cell caches, with composition order and cell counts checked before assembly. These maps show sector outlines without internal mesh edges. Each field uses the canonical palette and full-history limits shared by the six meshes; necrosis uses the fixed range 0-1. Use --comparison-time SECONDS to select another time within every native history; times outside the common interval are rejected. Dashed time-course lines mark the map time; dotted lines mark the common interval end.

Stable manuscript filenames are available under manuscript/exports/, for example figure-04.svg through figure-07.svg and figure-s1-convergence.svg. A format is exported only when it was rendered, so requesting --formats pdf adds the matching .pdf names. These files share storage with validated generations. Regenerating a figure refreshes its stable exports, caption and provenance from the complete generation. Outputs live under manuscript/<figure>/manifest.json in complete generations. Each generation contains the rendered exports, SVG and PNG by default, caption.md and provenance.json, including run conditions, snapshot times and native end times. Dense map panels are rasterized at the chosen export DPI to keep the composite files practical; annotations and time-course curves stay vector in SVG/PDF. Each map is rendered as an independent small tile so PDF export never retains one full-composite raster buffer per panel. All manuscript formats are cropped to displayed content with a one-point safety border. Cropping retains labels, arrows, legends and white living tissue, and preserves the spacing between panels. Hidden map-axis rectangles do not add empty margins; embedded map tiles are cropped geometrically, independently of their scalar colors. SVG/PDF text and time courses remain vector. The default is 300 DPI. Select --dpi 600 when higher-resolution maps are needed. For example, regenerate only Figures 4 and 7 as SVG and PNG:

uv run --no-sync python -m research.spt.manuscript \
  --figures 4 7 --formats svg png --dpi 300

The default representative dose is index 6, matching the original manuscript. --dose changes that selection for Figures 4 and 5 and records the choice in their provenance and captions. Keep the default when reproducing the manuscript. The final inspection report includes these outputs in its Manuscript figures section and uses the same small, cached previews as other figures.

Flat tissue cutouts for linear increase¶

Create additional top-down CYP2E1, toxic-compound and necrosis panels by repeating the linear-increase lobule on a hexagonal lattice:

uv run --no-sync python -m research.spt.tissue_panels

The defaults use D6 (sim014) at 0, 6, 12, 18 and 24 h, with canonical palettes and full-cohort scalar limits. CYP2E1 displays the model's rr_protein field. Continuous fields use native-time interpolation; necrosis retains the preceding native state. The new workflow reads local processed cells and cached XY geometry, without reopening raw HDF histories.

Each panel measures 2 native flat-to-flat lobule widths per side and samples an independent random XY position within a 4 by 4 lobule-width tissue square. Crop centers are sampled uniformly wherever the complete cutout fits. Each field/time panel shows its own tissue location, including in the composite. The default --seed 0 reproduces the crop positions; change --seed to sample another set. Exact centers and the PCG64 generator are recorded in provenance. The hexagonal lattice covers the entire square and its corners, with no exterior margins. The number of intersecting lobules depends on the crop. Border cells are clipped exactly; they keep their original cell values. Native central and portal vessel openings remain visible. Copies repeat the same simulation; this is a tissue visualization rather than a new multi-lobule transport solve.

These panels use a direct XY projection with zero tilt and no visible thickness, independently of the camera used by other geometry figures. New tissue panels omit the internal cell grid by default; --mesh enables the canonical transparent mesh lines. --vessel-color lightgray distinguishes the native lumens from white living tissue without assigning scalar values to the vessels.

Each standalone panel is an exact, margin-free square, 900 by 900 pixels at 300 DPI with the default settings. SVG and PDF versions contain the same bounded map image. An additional annotated composite arranges fields in rows and time from left to right, with shared scalar bars, the pattern icon and a binary necrosis legend. For a final-time selection or a wider square:

uv run --no-sync python -m research.spt.tissue_panels \
  --dose 6 --times-hours 24 --width-lobules 2.5

Widths smaller than 2, negative seeds, invalid dose indices and times outside the source history are rejected. --formats, --dpi, --mesh and --vessel-color control only these new outputs. Existing manuscript figures and numerical results are preserved. The default stable exports include linear-increase-tissue-d6.svg and individual files such as linear-increase-tissue-d6-necrosis-t86400s.png, under manuscript/exports/. Exact fractional seconds use p instead of a decimal point in filenames. Captions and provenance record the source run, time selection, scalar ranges, lattice spacing, vessel coverage and clipping convention.

Outputs publish through one atomic generation under manuscript/linear-increase-tissue-d6/manifest.json. An unchanged invocation checks metadata without opening numerical arrays or rendering. Missing stable exports are restored from that generation. The HTML report lists the composite and each field/time panel separately, with small previews and descriptive titles.

Final inspection report¶

Open index.html in the results directory in a browser after the pipeline finishes. The offline report gathers all current published figures: manuscript figures, numerical analysis, geometry overviews, individual sixths and lobules, mixed sectors, and simulation conditions. Exports of the same figure share one card. Previously generated figures remain included when the latest pipeline invocation selects a subset of the data. README.md beside it explains the folder layout.

Search by run, pattern, dose, or title, and combine section, study family, geometry, pattern, dose, and scale filters. Click a preview for the zoom viewer. Left and right arrow keys move through the filtered selection; Escape closes the viewer. SVG and PDF links open the original exports. Data and provenance links provide the numerical tables, NetCDF results, generation metadata, and execution records. Aggregate figures match every pattern and dose they contain.

Refresh only the report after generating additional figures:

uv run --no-sync python -m research.spt.report

This reads local manifests and output metadata without accessing the raw share, opening numerical arrays, or rerendering figures. It lists missing, modified, or invalid generations as omissions instead of linking incomplete outputs. The pipeline status refers to the latest pipeline invocation; the report does not independently validate raw-source freshness. Condition PNGs are included when present. The gallery uses lazily loaded, lossless WebP previews bounded by 768 by 768 pixels. Opening a figure loads its original resolution; PNG, SVG and PDF downloads remain unchanged. Previews are cached under cache/previews/ using relative source paths and file metadata. An unchanged report refresh reads only metadata and reuses previews without decoding the original images. Missing or changed previews are regenerated locally without rereading simulation data, and previews no longer referenced by the report are removed.

The HTML embeds its styling and interaction code and needs no web server or internet connection. Images and data retain relative links. Move or share the whole results directory, including its generation and cache folders, to keep those links working. The report generator also recognizes relocated generation files from their manifests. Generated reports and private results stay outside the repository.