diff --git a/applications/dynaclr/configs/collections/zuben_gut/conform_parquet.py b/applications/dynaclr/configs/collections/zuben_gut/conform_parquet.py new file mode 100644 index 000000000..190421970 --- /dev/null +++ b/applications/dynaclr/configs/collections/zuben_gut/conform_parquet.py @@ -0,0 +1,64 @@ +"""Step 2: conform Zuben's cell-index parquet to the canonical schema (end-to-end). + +Single authoritative step, run after the v3 stores are preprocessed (Step 1): + +1. Repoint ``store_path`` from the v2 originals to the v3 copies (same basename). +2. Fill required/derived columns absent from the source + (``tracks_path``, ``microscope``, ``T_shape``, ``C_shape``). +3. ``preprocess_cell_index`` fills the ``norm_*`` columns from the v3 ``.zattrs``. +4. Seed ``z_focus = z`` LAST — gut has no ``focus_slice`` zattrs, so the bbox-center + ``z`` IS the focus. This must come after step 3 because ``preprocess_cell_index`` + also writes ``z_focus`` (as NaN here, since there is no focus_slice) and would + otherwise clobber it. The datamodule centers the Z window on ``z_focus``. + +Run:: + + uv run --no-sync python applications/dynaclr/configs/collections/zuben_gut/conform_parquet.py +""" + +import os + +import pandas as pd + +from viscy_data.cell_index import preprocess_cell_index, read_cell_index, write_cell_index + +SRC = "/hpc/projects/jacobo_group/zuben/proj/gutCellClassifier/data/dynaclr_cell_index_bbox_center.parquet" +V3_ROOT = "/hpc/projects/organelle_phenotyping/datasets/zuben_gut_development" +OUT = "/hpc/projects/jacobo_group/collab/ed/dynaclr/dynaclr_cell_index_gut_v1.parquet" +N_CHANNELS = 4 + + +def main() -> None: + """Conform the parquet, fill norm stats, and seed z_focus from the bbox-center z.""" + df = pd.read_parquet(SRC) + n_rows_in = len(df) + + # Repoint store_path to the v3 copies (same basename). + df["store_path"] = df["store_path"].astype(str).map(lambda p: f"{V3_ROOT}/{os.path.basename(p)}") + + # Required columns absent from the source parquet. + df["tracks_path"] = "" # ignored by ExperimentRegistry.from_cell_index + df["microscope"] = "" + df["T_shape"] = 1 # static: single timepoint + df["C_shape"] = N_CHANNELS + + write_cell_index(df, OUT) + + # Fill norm_* from the v3 .zattrs (writes z_focus as NaN — no focus_slice). + preprocess_cell_index(OUT, focus_channel="nuclear") + + # Seed z_focus = z LAST so it survives preprocess_cell_index. + out_df = read_cell_index(OUT) + out_df["z_focus"] = out_df["z"].astype("float32") + write_cell_index(out_df, OUT) + + check = read_cell_index(OUT) + print("# Step 2 — conform parquet\n") + print(f"- rows in: **{n_rows_in}**, rows out: **{len(check)}**") + print(f"- output: `{OUT}`") + print(f"- norm_mean NaN: **{int(check['norm_mean'].isna().sum())}**") + print(f"- z_focus NaN: **{int(check['z_focus'].isna().sum())}** (should be 0)") + + +if __name__ == "__main__": + main() diff --git a/applications/dynaclr/configs/collections/zuben_gut/gen_convert_v3.py b/applications/dynaclr/configs/collections/zuben_gut/gen_convert_v3.py new file mode 100644 index 000000000..9a3c34c09 --- /dev/null +++ b/applications/dynaclr/configs/collections/zuben_gut/gen_convert_v3.py @@ -0,0 +1,92 @@ +"""Generate per-store biahub-concatenate configs to copy Zuben's 25 v2 gut stores to v3. + +Each of Zuben's stores is one gut = one experiment = a single position ``A/1/0`` with the +four channels ``[nuclear, septate, brush_border, SuH]``. We convert each store to its own +v3 store under ``OUTPUT_ROOT`` preserving the original basename and the ``A/1/0`` layout, so +the cell-index parquet's ``(store_path, well, fov)`` mapping stays valid — only the +``store_path`` directory changes. + +Emits, per store, ``configs/.yml`` and a single ``submit_all.sh`` that runs +``biahub concatenate`` (conda env ``biautils``) for every store as SLURM jobs. + +Run:: + + uv run --no-sync python applications/dynaclr/configs/collections/zuben_gut/gen_convert_v3.py + bash /convert_v3/submit_all.sh +""" + +import os +from pathlib import Path + +import pandas as pd + +PARQUET = "/hpc/projects/jacobo_group/zuben/proj/gutCellClassifier/data/dynaclr_cell_index_bbox_center.parquet" +OUTPUT_ROOT = "/hpc/projects/organelle_phenotyping/datasets/zuben_gut_development" +CHANNELS = ["nuclear", "septate", "brush_border", "SuH"] +CHUNKS_CZYX = [1, 16, 256, 256] +SHARDS_RATIO = [1, 1, 4, 8, 8] +OME_VERSION = "0.5" +BIAHUB_ENV = "biautils" + +HERE = Path(__file__).resolve().parent +WORKDIR = HERE / "convert_v3" +CONFIGS_DIR = WORKDIR / "configs" + + +def _yaml_config(store_path: str) -> str: + import yaml + + cfg = { + # One input FOV glob per data path: this store's single position A/1/0. + "concat_data_paths": [f"{store_path}/A/1/*"], + "time_indices": "all", + # List-of-lists: one channel list per data path (here, one path). + "channel_names": [list(CHANNELS)], + "X_slice": "all", + "Y_slice": "all", + "Z_slice": "all", + "chunks_czyx": list(CHUNKS_CZYX), + "shards_ratio": list(SHARDS_RATIO), + "output_ome_zarr_version": OME_VERSION, + } + return yaml.safe_dump(cfg, default_flow_style=False, sort_keys=False) + + +def main() -> None: + """Emit one biahub-concatenate config per store plus a SLURM submit driver.""" + df = pd.read_parquet(PARQUET, columns=["store_path"]) + stores = sorted(df["store_path"].astype(str).unique()) + + CONFIGS_DIR.mkdir(parents=True, exist_ok=True) + submit_lines = [ + "#!/bin/bash", + "set -euo pipefail", + "# Auto-generated by gen_convert_v3.py. Submits one biahub concatenate per store.", + f"mkdir -p {OUTPUT_ROOT}", + "", + ] + + for store in stores: + stem = os.path.basename(store) # e.g. AAY6_sox21a_d0_63x_gut1.zarr + cfg_path = CONFIGS_DIR / f"{stem.replace('.zarr', '')}.yml" + cfg_path.write_text(_yaml_config(store)) + out_path = f"{OUTPUT_ROOT}/{stem}" + submit_lines.append(f'echo "=== {stem} ==="') + submit_lines.append( + f"conda run -n {BIAHUB_ENV} biahub concatenate " + f'-c "{cfg_path}" -o "{out_path}" -m -sb "{WORKDIR / "sbatch_overrides.sh"}"' + ) + submit_lines.append("") + + (WORKDIR / "submit_all.sh").write_text("\n".join(submit_lines) + "\n") + (WORKDIR / "sbatch_overrides.sh").write_text( + "#!/bin/bash\n#SBATCH --partition=cpu\n#SBATCH --cpus-per-task=4\n#SBATCH --mem-per-cpu=32G\n" + ) + + print(f"Generated {len(stores)} configs under {CONFIGS_DIR}") + print(f"Submit driver: {WORKDIR / 'submit_all.sh'}") + print(f"Output root: {OUTPUT_ROOT}") + + +if __name__ == "__main__": + main() diff --git a/applications/dynaclr/configs/collections/zuben_gut/gen_preprocess.py b/applications/dynaclr/configs/collections/zuben_gut/gen_preprocess.py new file mode 100644 index 000000000..cb4ee68be --- /dev/null +++ b/applications/dynaclr/configs/collections/zuben_gut/gen_preprocess.py @@ -0,0 +1,56 @@ +"""Generate + submit `viscy preprocess` SLURM jobs for the 25 v3 gut stores. + +Writes per-channel normalization stats (fov/dataset/timepoint) into each store's +``.zattrs``. One SLURM job per store (each store is a single large FOV). + +Run:: + + uv run --no-sync python applications/dynaclr/configs/collections/zuben_gut/gen_preprocess.py + bash /preprocess/submit_all.sh +""" + +import glob +from pathlib import Path + +V3_ROOT = "/hpc/projects/organelle_phenotyping/datasets/zuben_gut_development" +REPO = "/hpc/mydata/eduardo.hirata/repos/viscy" +VENV = f"{REPO}/.venv-dynaclr" +HERE = Path(__file__).resolve().parent +WORKDIR = HERE / "preprocess" + +JOB_TEMPLATE = """#!/bin/bash +#SBATCH --job-name=pp_{stem} +#SBATCH --nodes=1 +#SBATCH --ntasks-per-node=1 +#SBATCH --partition=cpu +#SBATCH --cpus-per-task=32 +#SBATCH --mem-per-cpu=4G +#SBATCH --time=02:00:00 +#SBATCH --output={workdir}/slurm_pp_{stem}_%j.out + +export PYTHONNOUSERSITE=1 +export UV_PROJECT_ENVIRONMENT={venv} + +uv run --project "{repo}" --package dynaclr \\ + viscy preprocess --data_path "{store}" \\ + --channel_names=-1 --num_workers 32 --block_size 32 +""" + + +def main() -> None: + """Emit one `viscy preprocess` SLURM job per v3 store plus a submit driver.""" + WORKDIR.mkdir(parents=True, exist_ok=True) + stores = sorted(glob.glob(f"{V3_ROOT}/*.zarr")) + submit = ["#!/bin/bash", "set -euo pipefail", ""] + for store in stores: + stem = Path(store).name.replace(".zarr", "") + job = WORKDIR / f"pp_{stem}.sh" + job.write_text(JOB_TEMPLATE.format(stem=stem, workdir=WORKDIR, venv=VENV, repo=REPO, store=store)) + submit.append(f"sbatch {job}") + (WORKDIR / "submit_all.sh").write_text("\n".join(submit) + "\n") + print(f"Generated {len(stores)} preprocess jobs under {WORKDIR}") + print(f"Submit: bash {WORKDIR / 'submit_all.sh'}") + + +if __name__ == "__main__": + main() diff --git a/applications/dynaclr/configs/collections/zuben_gut/verify_parquet_zarr.py b/applications/dynaclr/configs/collections/zuben_gut/verify_parquet_zarr.py new file mode 100644 index 000000000..9b16514fe --- /dev/null +++ b/applications/dynaclr/configs/collections/zuben_gut/verify_parquet_zarr.py @@ -0,0 +1,80 @@ +"""Step 0: verify Zuben's cell-index parquet is consistent with its zarr stores. + +Checks, for every store referenced in the parquet: +- the store opens and exposes the expected channel names, +- each (well, fov) referenced resolves to a real position, +- cell centroids fit inside the FOV with room for the requested patch half-width. + +Run:: + + uv run --no-sync python applications/dynaclr/configs/collections/zuben_gut/verify_parquet_zarr.py +""" + +import sys + +import pandas as pd +from iohub import open_ome_zarr + +PARQUET = "/hpc/projects/jacobo_group/zuben/proj/gutCellClassifier/data/dynaclr_cell_index_bbox_center.parquet" +EXPECTED_CHANNELS = ["nuclear", "septate", "brush_border", "SuH"] +YX_PATCH_SIZE = (256, 256) # extraction patch used by the bag-of-channels config + + +def main() -> int: + """Check every store opens, channels match, positions resolve, and report border-OOB cells.""" + df = pd.read_parquet(PARQUET) + y_half = YX_PATCH_SIZE[0] // 2 + x_half = YX_PATCH_SIZE[1] // 2 + + problems: list[str] = [] + n_cells_out_of_bounds = 0 + + for store_path, store_group in df.groupby("store_path", observed=True): + store_path = str(store_path) + try: + with open_ome_zarr(store_path, mode="r") as plate: + channels = list(plate.channel_names) + if channels != EXPECTED_CHANNELS: + problems.append(f"{store_path}: channels {channels} != {EXPECTED_CHANNELS}") + positions = {name for name, _ in plate.positions()} + except Exception as exc: # noqa: BLE001 - surface any open failure + problems.append(f"{store_path}: failed to open ({exc})") + continue + + for (well, fov), fov_group in store_group.groupby(["well", "fov"], observed=True): + pos_key = f"{well}/{fov}" + if pos_key not in positions: + problems.append(f"{store_path}: position {pos_key} not found") + continue + y_shape = int(fov_group["Y_shape"].iloc[0]) + x_shape = int(fov_group["X_shape"].iloc[0]) + y = fov_group["y"].to_numpy() + x = fov_group["x"].to_numpy() + oob = (y < y_half) | (y > y_shape - y_half) | (x < x_half) | (x > x_shape - x_half) + n_cells_out_of_bounds += int(oob.sum()) + + n_stores = df["store_path"].nunique() + n_cells = df["cell_id"].nunique() + + print("# Step 0 — parquet↔zarr consistency\n") + print(f"- stores checked: **{n_stores}**") + print(f"- unique cells: **{n_cells}**") + print(f"- expected channels: `{EXPECTED_CHANNELS}`") + print( + f"- cells whose {YX_PATCH_SIZE} patch would fall out of bounds: " + f"**{n_cells_out_of_bounds}** " + f"({100 * n_cells_out_of_bounds / len(df):.1f}% of rows)" + ) + + if problems: + print(f"\n**{len(problems)} problem(s):**") + for p in problems[:50]: + print(f" - {p}") + return 1 + + print("\n**OK** — all stores open, channels match, positions resolve.") + return 0 + + +if __name__ == "__main__": + sys.exit(main()) diff --git a/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.sh b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.sh new file mode 100644 index 000000000..c3617fb11 --- /dev/null +++ b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.sh @@ -0,0 +1,47 @@ +#!/bin/bash +# DynaCLR-3D-Gut-BagOfChannels (Phase 1) — Zuben gut cells, bag-of-channels SimCLR. +# +# New run: +# sbatch applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.sh +# Resume: +# CKPT_PATH=.../last.ckpt WANDB_RUN_ID= sbatch .../DynaCLR-3D-Gut-BagOfChannels.sh + +#SBATCH --job-name=dynaclr_gut_boc +#SBATCH --nodes=1 +#SBATCH --ntasks-per-node=2 +#SBATCH --gpus-per-node=2 +#SBATCH --partition=gpu +#SBATCH --cpus-per-task=15 +#SBATCH --mem-per-cpu=8G +#SBATCH --time=2-00:00:00 + +# Set WORKSPACE_DIR to YOUR clone of the repo before submitting, e.g. +# WORKSPACE_DIR=/hpc/mydata//repos/VisCy sbatch .sh +# MODEL_ROOT is where checkpoints/configs are written (defaults to your clone's +# models/ dir; override to a shared project path if desired). +WORKSPACE_DIR="${WORKSPACE_DIR:?Set WORKSPACE_DIR to your repo clone path}" + +export PROJECT="DynaCLR-3D-Gut-BagOfChannels" +export RUN_NAME="gut-3d-z24-cellz-64-ntxent-t0p2-self" +export CONFIGS="applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.yml" +export MODEL_ROOT="${MODEL_ROOT:-${WORKSPACE_DIR}/models}" +# Point at the dynaclr-pinned venv in your clone (avoids the shared .venv sync race). +export UV_PROJECT_ENVIRONMENT="${UV_PROJECT_ENVIRONMENT:-${WORKSPACE_DIR}/.venv-dynaclr}" + +# W&B writes sample images to a temp dir under $TMPDIR before upload. On some +# SLURM nodes the default /tmp is per-job and gets swept mid-run, causing +# `FileNotFoundError: .../wandb-media/*.png` at validation image logging. Pin +# TMPDIR + WANDB_DIR to persistent paths we create so they can't disappear. +export TMPDIR="${TMPDIR:-${MODEL_ROOT}/tmp}" +export WANDB_DIR="${WANDB_DIR:-${MODEL_ROOT}/wandb}" +mkdir -p "$TMPDIR" "$WANDB_DIR" + +# The shared trainer recipe logs to the `computational_imaging` W&B entity. Set +# WANDB_ENTITY to your own entity to log there instead (EXTRA_ARGS overrides the +# recipe). Leave unset to keep the default. +if [ -n "${WANDB_ENTITY:-}" ]; then + export EXTRA_ARGS="${EXTRA_ARGS:-} --trainer.logger.init_args.entity=${WANDB_ENTITY}" +fi + +# Absolute path (SLURM spools this script, so $(dirname "$0") would break). +source "${WORKSPACE_DIR}/applications/dynaclr/configs/training/slurm/train.sh" diff --git a/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.yml b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.yml new file mode 100644 index 000000000..b019e8cd7 --- /dev/null +++ b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.yml @@ -0,0 +1,151 @@ +# DynaCLR-3D-Gut-BagOfChannels (Phase 1) +# ======================================= +# 3D bag-of-channels contrastive learning on Zuben's zebrafish gut data. +# One random channel per sample (nuclear | septate | brush_border | SuH), +# 24-slice Z window, 64x64 XY — matches Zuben's 24x64x64 patch convention. +# +# STATIC DATA: single timepoint (t=0), no tracking. SimCLR self-positives +# (same crop → same marker). One marker per batch (batch_group_by=marker) so the +# model can't shortcut on channel identity; batches balanced across the 6 states +# (stratify_by=perturbation) to counter the ~3.5x state imbalance. +# (Phase 3, later: consecutive-state positives to regularize the space temporally.) +# +# Normalization reads per-FOV full-frame stats from the v3 .zattrs written by +# `viscy preprocess`. Single timepoint → timepoint_statistics == fov_statistics. +# +# Z window centered on the parquet `z_focus` column (seeded from bbox-center z +# by conform_parquet.py), 24 slices, no random Z crop. XY extracted at 80 then +# center-cropped to 64. +# Pipeline: extract (24,80,80) → normalize → affine → flip/contrast/noise +# → CenterCrop (24,64,64) [auto-appended]. +# +# Launch: +# sbatch applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-BagOfChannels.sh + +base: + - ../recipes/trainer/fit.yml + - ../recipes/topology/ddp_2gpu.yml + - ../recipes/model/contrastive_encoder_convnext_tiny.yml + +trainer: + precision: bf16-mixed + max_epochs: 150 + logger: + init_args: + project: DynaCLR-3D-Gut-BagOfChannels + name: gut-3d-z24-96to80to64-ntxent-t0p2-self + # Re-list callbacks WITHOUT OnlineEvalCallback (static data → degenerate + # temporal kNN eval). Loss/val + post-hoc PCA/UMAP colored by stage/marker + # is the evaluation signal for Phase 1. + callbacks: + - class_path: lightning.pytorch.callbacks.LearningRateMonitor + init_args: + logging_interval: step + - class_path: lightning.pytorch.callbacks.ModelCheckpoint + init_args: + monitor: loss/val + every_n_epochs: 1 + save_top_k: 5 + save_last: true + +model: + init_args: + encoder: + init_args: + in_stack_depth: 24 + stem_kernel_size: [4, 4, 4] + stem_stride: [4, 4, 4] + projection_dim: 32 + drop_path_rate: 0.1 + loss_function: + init_args: + temperature: 0.2 + lr: 0.00002 + pca_color_keys: "[perturbation,experiment,marker]" + # PCA pairplot cadence (overrides the base recipe's 10). Watch the 6 states + # separate early. + log_embeddings_every_n_epochs: 5 + log_negative_metrics_every_n_epochs: 2 + example_input_array_shape: [1, 1, 24, 64, 64] + +data: + class_path: dynaclr.data.datamodule.MultiExperimentDataModule + init_args: + cell_index_path: /hpc/projects/jacobo_group/collab/ed/dynaclr/dynaclr_cell_index_gut_v1.parquet + focus_channel: null + # The datamodule centers the Z window on the parquet `z_focus` column; + # conform_parquet.py seeds z_focus = bbox-center z (gut cells span the full + # stack within one store, so per-FOV focus won't do). 24-slice window, + # symmetric (offset 0.5), no Z rescale (extraction == window), no random crop. + z_window: 24 + z_extraction_window: 24 + z_focus_offset: 0.5 + # Small XY margin (80 → 64) so the auto-appended CenterCrop trims affine + # rotation zero-fill; no explicit random XY crop. + yx_patch_size: [80, 80] + final_yx_patch_size: [64, 64] + channels_per_sample: 1 + # SimCLR self-positive: anchor == positive (same crop, two augmentations) → + # trivially the SAME marker, and no supervised pull that would fight the + # negatives. State structure emerges unsupervised. + positive_cell_source: self + # One marker per batch: all negatives share the anchor's marker, so the model + # can't take the "tell channels apart" shortcut and must learn within-marker + # (cell-state) structure. + batch_group_by: marker + # Within the single-marker batch, balance the 6 developmental states — the + # dataset is ~3.5x imbalanced (stage_0 1379 vs stage_5 390 per marker). + stratify_by: [perturbation] + # Each gut store is a single FOV → split cells within each gut (not by FOV, + # which would leave 0 val FOVs). All rows of one cell stay on one side. + split_mode: cell + split_ratio: 0.8 + batch_size: 256 + num_workers: 4 + seed: 42 + normalizations: + - class_path: viscy_transforms.NormalizeSampled + init_args: + keys: [channel_0] + level: timepoint_statistics + subtrahend: mean + divisor: std + augmentations: + - class_path: viscy_transforms.BatchedRandAffined + init_args: + keys: [channel_0] + prob: 0.8 + scale_range: [[0.9, 1.1], [0.9, 1.1], [0.9, 1.1]] + rotate_range: [3.14, 0.0, 0.0] + shear_range: [0.05, 0.05, 0.0, 0.05, 0.0, 0.05] + # No random spatial crop — the Z window is already centered on the cell + # plane and XY is tight. The datamodule auto-appends a CenterCrop to + # [24, 64, 64], which also trims affine rotation zero-fill at the edges. + - class_path: viscy_transforms.BatchedRandFlipd + init_args: + keys: [channel_0] + spatial_axes: [1, 2] + prob: 0.5 + - class_path: viscy_transforms.BatchedRandAdjustContrastd + init_args: + keys: [channel_0] + prob: 0.5 + gamma: [0.6, 1.6] + - class_path: viscy_transforms.BatchedRandScaleIntensityd + init_args: + keys: [channel_0] + prob: 0.5 + factors: 0.5 + - class_path: viscy_transforms.BatchedRandGaussianSmoothd + init_args: + keys: [channel_0] + prob: 0.5 + sigma_x: [0.25, 0.50] + sigma_y: [0.25, 0.50] + sigma_z: [0.0, 0.2] + - class_path: viscy_transforms.BatchedRandGaussianNoised + init_args: + keys: [channel_0] + prob: 0.5 + mean: 0.0 + std: 0.1 diff --git a/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.sh b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.sh new file mode 100644 index 000000000..d23dd4b45 --- /dev/null +++ b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.sh @@ -0,0 +1,43 @@ +#!/bin/bash +# DynaCLR-3D-Gut-MultiChannel (Phase 2) — Zuben gut cells, 4-channel input. +# +# New run: +# sbatch applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.sh +# Resume: +# CKPT_PATH=.../last.ckpt WANDB_RUN_ID= sbatch .../DynaCLR-3D-Gut-MultiChannel.sh + +#SBATCH --job-name=dynaclr_gut_4ch +#SBATCH --nodes=1 +#SBATCH --ntasks-per-node=2 +#SBATCH --gpus-per-node=2 +#SBATCH --partition=gpu +#SBATCH --cpus-per-task=15 +#SBATCH --mem-per-cpu=8G +#SBATCH --time=2-00:00:00 + +# Set WORKSPACE_DIR to YOUR clone of the repo before submitting, e.g. +# WORKSPACE_DIR=/hpc/mydata//repos/VisCy sbatch .sh +WORKSPACE_DIR="${WORKSPACE_DIR:?Set WORKSPACE_DIR to your repo clone path}" + +export PROJECT="DynaCLR-3D-Gut-MultiChannel" +export RUN_NAME="gut-3d-4ch-z24-cellz-64-ntxent-t0p2-self" +export CONFIGS="applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.yml" +export MODEL_ROOT="${MODEL_ROOT:-${WORKSPACE_DIR}/models}" +export UV_PROJECT_ENVIRONMENT="${UV_PROJECT_ENVIRONMENT:-${WORKSPACE_DIR}/.venv-dynaclr}" + +# W&B writes sample images to a temp dir under $TMPDIR before upload. On some +# SLURM nodes the default /tmp is per-job and gets swept mid-run, causing +# `FileNotFoundError: .../wandb-media/*.png` at validation image logging. Pin +# TMPDIR + WANDB_DIR to persistent paths we create so they can't disappear. +export TMPDIR="${TMPDIR:-${MODEL_ROOT}/tmp}" +export WANDB_DIR="${WANDB_DIR:-${MODEL_ROOT}/wandb}" +mkdir -p "$TMPDIR" "$WANDB_DIR" + +# The shared trainer recipe logs to the `computational_imaging` W&B entity. Set +# WANDB_ENTITY to your own entity to log there instead. Leave unset for default. +if [ -n "${WANDB_ENTITY:-}" ]; then + export EXTRA_ARGS="${EXTRA_ARGS:-} --trainer.logger.init_args.entity=${WANDB_ENTITY}" +fi + +# Absolute path (SLURM spools this script, so $(dirname "$0") would break). +source "${WORKSPACE_DIR}/applications/dynaclr/configs/training/slurm/train.sh" diff --git a/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.yml b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.yml new file mode 100644 index 000000000..269d8628b --- /dev/null +++ b/applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.yml @@ -0,0 +1,147 @@ +# DynaCLR-3D-Gut-MultiChannel (Phase 2) +# ====================================== +# 4-channel contrastive learning on Zuben's zebrafish gut data. +# All channels stacked as input (nuclear, septate, brush_border, SuH) → +# (B, 4, 24, 64, 64). Per-channel normalization: each channel normalized with +# its own per-FOV full-frame stats from the v3 .zattrs. +# +# STATIC DATA: SimCLR self-positives; batches balanced across the 6 states +# (stratify_by=perturbation) to counter the ~3.5x state imbalance. +# Switch to consecutive-state positives later (Phase 3) to regularize the space +# temporally — needs a stage-adjacency bucket or code change (see plan). +# +# Z window centered on the parquet `z_focus` column (seeded from bbox-center z +# by conform_parquet.py), 24 slices, no random Z crop. XY 80 → CenterCrop 64. +# +# Launch: +# sbatch applications/dynaclr/configs/training/DynaCLR-3D/DynaCLR-3D-Gut-MultiChannel.sh + +base: + - ../recipes/trainer/fit.yml + - ../recipes/topology/ddp_2gpu.yml + - ../recipes/model/contrastive_encoder_convnext_tiny.yml + +trainer: + precision: bf16-mixed + max_epochs: 150 + logger: + init_args: + project: DynaCLR-3D-Gut-MultiChannel + name: gut-3d-4ch-z24-cellz-64-ntxent-t0p2-self + callbacks: + - class_path: lightning.pytorch.callbacks.LearningRateMonitor + init_args: + logging_interval: step + - class_path: lightning.pytorch.callbacks.ModelCheckpoint + init_args: + monitor: loss/val + every_n_epochs: 1 + save_top_k: 5 + save_last: true + +model: + init_args: + encoder: + init_args: + in_channels: 4 + in_stack_depth: 24 + stem_kernel_size: [4, 4, 4] + stem_stride: [4, 4, 4] + projection_dim: 32 + drop_path_rate: 0.1 + loss_function: + init_args: + temperature: 0.2 + lr: 0.00002 + pca_color_keys: "[perturbation,experiment,marker]" + # PCA pairplot cadence (overrides the base recipe's 10). + log_embeddings_every_n_epochs: 5 + log_negative_metrics_every_n_epochs: 2 + example_input_array_shape: [1, 4, 24, 64, 64] + +data: + class_path: dynaclr.data.datamodule.MultiExperimentDataModule + init_args: + cell_index_path: /hpc/projects/jacobo_group/collab/ed/dynaclr/dynaclr_cell_index_gut_v1.parquet + focus_channel: null + # Z window centers on the parquet z_focus column (seeded from bbox-center z). + z_window: 24 + z_extraction_window: 24 + z_focus_offset: 0.5 + yx_patch_size: [80, 80] + final_yx_patch_size: [64, 64] + # All 4 channels stacked per sample → (B, 4, Z, Y, X). No marker shortcut + # here (every sample carries all channels), so no batch_group_by needed. + channels_per_sample: null + # SimCLR self-positive (same crop, two augmentations); state structure emerges + # unsupervised. Batches balanced across the 6 states to counter imbalance. + positive_cell_source: self + stratify_by: [perturbation] + split_mode: cell + split_ratio: 0.8 + batch_size: 64 + num_workers: 4 + seed: 42 + # Per-channel normalization: each channel keyed by its real name, all at + # fov_statistics (== timepoint_statistics for single-timepoint data). + normalizations: + - class_path: viscy_transforms.NormalizeSampled + init_args: + keys: [nuclear] + level: fov_statistics + subtrahend: mean + divisor: std + - class_path: viscy_transforms.NormalizeSampled + init_args: + keys: [septate] + level: fov_statistics + subtrahend: mean + divisor: std + - class_path: viscy_transforms.NormalizeSampled + init_args: + keys: [brush_border] + level: fov_statistics + subtrahend: mean + divisor: std + - class_path: viscy_transforms.NormalizeSampled + init_args: + keys: [SuH] + level: fov_statistics + subtrahend: mean + divisor: std + augmentations: + - class_path: viscy_transforms.BatchedRandAffined + init_args: + keys: [nuclear, septate, brush_border, SuH] + prob: 0.8 + scale_range: [[0.9, 1.1], [0.9, 1.1], [0.9, 1.1]] + rotate_range: [3.14, 0.0, 0.0] + shear_range: [0.05, 0.05, 0.0, 0.05, 0.0, 0.05] + - class_path: viscy_transforms.BatchedRandFlipd + init_args: + keys: [nuclear, septate, brush_border, SuH] + spatial_axes: [1, 2] + prob: 0.5 + - class_path: viscy_transforms.BatchedRandAdjustContrastd + init_args: + keys: [nuclear, septate, brush_border, SuH] + prob: 0.5 + gamma: [0.6, 1.6] + - class_path: viscy_transforms.BatchedRandScaleIntensityd + init_args: + keys: [nuclear, septate, brush_border, SuH] + prob: 0.5 + factors: 0.5 + - class_path: viscy_transforms.BatchedRandGaussianSmoothd + init_args: + keys: [nuclear, septate, brush_border, SuH] + prob: 0.5 + sigma_x: [0.25, 0.50] + sigma_y: [0.25, 0.50] + sigma_z: [0.0, 0.2] + - class_path: viscy_transforms.BatchedRandGaussianNoised + init_args: + keys: [nuclear, septate, brush_border, SuH] + prob: 0.5 + mean: 0.0 + std: 0.1 diff --git a/applications/dynaclr/docs/DAGs/training.md b/applications/dynaclr/docs/DAGs/training.md index ca448dd7c..5a3b7f257 100644 --- a/applications/dynaclr/docs/DAGs/training.md +++ b/applications/dynaclr/docs/DAGs/training.md @@ -30,11 +30,14 @@ dynaclr build-cell-index \ ▼ dynaclr preprocess-cell-index \ /hpc/.../collections/.parquet \ - --focus-channel Phase3D + --focus-channel Phase3D --focus-level fov │ opens each unique FOV once from zarr zattrs: │ norm_mean/std/median/iqr/max/min — per (cell, timepoint, channel) - │ z_focus_mean — per FOV (mean across timepoints) - │ z — per timepoint focus slice index + │ z_focus — focus plane the Z window centers on, + │ from focus_slice at --focus-level + │ (fov = per-FOV mean; per_timepoint = + │ per-timepoint index). `z` (tracked + │ position) is left unchanged. │ drops empty frames (max == 0) ▼ .parquet (ready: self-contained, no zarr reads at training time) @@ -44,8 +47,9 @@ viscy fit --config configs/training/.yml │ OR: sbatch configs/training/.sh (SLURM, recommended) │ MultiExperimentDataModule reads parquet only at init │ tensorstore opens zarr lazily on first batch - │ ExperimentRegistry reads plate.zattrs["focus_slice"] once at startup - │ for z_ranges (z_extraction_window centered on dataset z_focus_mean) + │ Z window centers on the parquet `z_focus` column (per-sample). If z_focus + │ is null, falls back to the per-experiment z_range from zattrs focus_slice, + │ else mid-stack. z_focus_offset sets the below/above split. ▼ checkpoints/ + wandb logs ``` @@ -71,7 +75,7 @@ viscy fit (GPU, hours–days) | Step | Command | Input | Output | | --------------------- | ------------------------------------------------------------------------- | -------------------------------------- | --------------------------------------------------------- | | Build cell index | `dynaclr build-cell-index --num-workers 8` | collection YAML + zarr + tracking CSVs | parquet with TCZYX shape columns | -| Preprocess cell index | `dynaclr preprocess-cell-index --focus-channel Phase3D` | parquet + zarr zattrs | parquet with norm stats, per-timepoint z, empties removed | +| Preprocess cell index | `dynaclr preprocess-cell-index --focus-channel Phase3D --focus-level fov` | parquet + zarr zattrs | parquet with norm stats + z_focus, empties removed | | Train (interactive) | `uv run viscy fit --config configs/training/.yml` | training config + parquet | checkpoints + logs | | Train (SLURM) | `sbatch configs/training/.sh` | training config + parquet | checkpoints + logs | | Resume (SLURM) | `CKPT_PATH=.../last.ckpt sbatch configs/training/.sh` | checkpoint path env var | resumed checkpoints | @@ -85,8 +89,8 @@ viscy fit (GPU, hours–days) | Pixel data (TCZYX arrays) | zarr store on VAST | `prepare run` → concatenate | | Cell tracking (y, x, t, track_id) | tracking.zarr on VAST | `prepare run` → concatenate | | Normalization stats (per FOV/timepoint) | zarr zattrs → parquet `norm_*` columns | `viscy preprocess` → `preprocess-cell-index` | -| Focus slice (per timepoint) | zarr zattrs → parquet `z` column | `viscy preprocess` → `preprocess-cell-index` | -| Focus slice mean (per FOV) | zarr zattrs → parquet `z_focus_mean` | `viscy preprocess` → `preprocess-cell-index` | +| Tracked cell z (per cell) | tracking → parquet `z` column (unchanged by preprocess) | `build-cell-index` | +| Focus plane (Z window center) | zarr zattrs focus_slice → parquet `z_focus` | `viscy preprocess` → `preprocess-cell-index --focus-level` | | TCZYX shape per FOV | parquet columns | `build-cell-index` | | Collection definition | `configs/collections/.yml` in git | manually authored | | Parquet | `/hpc/projects/organelle_phenotyping/models/collections/` | `build-cell-index` | @@ -158,7 +162,8 @@ To reproduce: `build-cell-index` → `preprocess-cell-index` from the same colle ## Notes - `preprocess-cell-index` overwrites the parquet in-place by default. Pass `--output` to write elsewhere. -- `--focus-channel Phase3D` selects which channel's `per_timepoint` focus indices are written to the `z` column. Use the channel that has the sharpest axial contrast (label-free Phase3D for most experiments). -- At training time, `ExperimentRegistry.__post_init__` reads `plate.zattrs["focus_slice"][channel]["dataset_statistics"]["z_focus_mean"]` to compute per-experiment z_ranges for patch extraction. This is the only zarr metadata read at training startup; the parquet is self-contained for all per-cell data. -- The `z` column in the parquet is carried through to embeddings obs during predict — downstream consumers (e.g., visualization) can use it to recover the in-focus plane for each cell at each timepoint. +- `--focus-channel Phase3D` selects which channel's `focus_slice` metadata feeds the `z_focus` column. Use the channel with sharpest axial contrast (label-free Phase3D for most experiments). `--focus-level {fov,per_timepoint}` picks per-FOV mean vs per-timepoint index. +- **Z window centering is parquet-first.** At training time the datamodule centers the Z window on the per-sample `z_focus` column. If `z_focus` is null it falls back to the per-experiment `z_range` (which `ExperimentRegistry` derives from `plate.zattrs["focus_slice"][channel]["dataset_statistics"]["z_focus_mean"]`, else mid-stack). Existing parquets without a `z_focus` column read as all-null → identical to the previous zattrs-only behavior. +- `z_focus` and `z` are distinct: `z_focus` = where the Z window centers; `z` = the tracked cell position (unchanged by preprocess). Datasets without focus_slice (e.g. static bbox-center data) can seed `z_focus = z` in their own prep script. +- The `z` column is carried through to embeddings obs during predict for downstream consumers. - For performance tuning (num_workers, pin_memory, batch_size, augmentation placement), see [profiling.md](profiling.md) — authored after the first validated profiling sweep. diff --git a/applications/dynaclr/src/dynaclr/data/datamodule.py b/applications/dynaclr/src/dynaclr/data/datamodule.py index 1c49aaf82..7462e479b 100644 --- a/applications/dynaclr/src/dynaclr/data/datamodule.py +++ b/applications/dynaclr/src/dynaclr/data/datamodule.py @@ -190,6 +190,7 @@ def __init__( positive_match_columns: list[str] | None = None, positive_channel_source: str = "same", label_columns: dict[str, str] | None = None, + split_mode: str = "fov", max_border_shift: int = -1, shuffle_val: bool = False, pin_memory: bool = True, @@ -203,6 +204,9 @@ def __init__( self.z_window = z_window self.z_extraction_window = z_extraction_window self.z_focus_offset = z_focus_offset + if split_mode not in ("fov", "cell"): + raise ValueError(f"split_mode must be 'fov' or 'cell', got {split_mode!r}") + self.split_mode = split_mode self.yx_patch_size = yx_patch_size self.final_yx_patch_size = final_yx_patch_size self.val_experiments = val_experiments if val_experiments is not None else [] @@ -474,36 +478,44 @@ def _setup_fov_split(self, registry: ExperimentRegistry, cell_index_df: pd.DataF # (81M+ rows for OPS), which hashes a Python tuple per row and # dominates setup-time memory. Per-group isin against a small # Python-set of FOV names is O(group_size) with no object index. - train_fovs_per_exp: dict[str, set[str]] = {} - val_fovs_per_exp: dict[str, set[str]] = {} + # split_key: "fov_name" holds out whole FOVs; "cell_id" holds out cells + # within each experiment (needed when every experiment is a single FOV, + # e.g. one-position-per-store datasets). Splitting on cell_id keeps all + # rows of one cell (its channels) on the same side — no channel leakage. + split_key = "fov_name" if self.split_mode == "fov" else "cell_id" + train_keys_per_exp: dict[str, set[str]] = {} + val_keys_per_exp: dict[str, set[str]] = {} for exp_name, group in full_index.tracks.groupby("experiment"): - fovs = sorted(group["fov_name"].unique()) - n_train = max(1, int(len(fovs) * self.split_ratio)) - rng.shuffle(fovs) - train_fovs_per_exp[exp_name] = set(fovs[:n_train]) - val_fovs_per_exp[exp_name] = set(fovs[n_train:]) - - n_train_fovs = sum(len(s) for s in train_fovs_per_exp.values()) - n_val_fovs = sum(len(s) for s in val_fovs_per_exp.values()) + keys = sorted(group[split_key].unique()) + n_train = max(1, int(len(keys) * self.split_ratio)) + rng.shuffle(keys) + train_keys_per_exp[exp_name] = set(keys[:n_train]) + val_keys_per_exp[exp_name] = set(keys[n_train:]) + + n_train_keys = sum(len(s) for s in train_keys_per_exp.values()) + n_val_keys = sum(len(s) for s in val_keys_per_exp.values()) _logger.info( - "FOV split (ratio=%.2f): %d train FOVs, %d val FOVs", + "%s split (ratio=%.2f): %d train %s, %d val %s", + self.split_mode, self.split_ratio, - n_train_fovs, - n_val_fovs, + n_train_keys, + split_key, + n_val_keys, + split_key, ) def _build_train_mask(df: pd.DataFrame) -> np.ndarray: - """Row-wise boolean mask: True if (experiment, fov_name) is train.""" + """Row-wise boolean mask: True if the split key is in the train set.""" mask = np.zeros(len(df), dtype=bool) # groupby("experiment") returns integer positions in ``df`` via # group.index after reset_index; we rely on the caller passing # reset-indexed frames (which is what MultiExperimentIndex produces). for exp_name, group in df.groupby("experiment", sort=False): - train_fovs = train_fovs_per_exp.get(exp_name, set()) - if not train_fovs: + train_keys = train_keys_per_exp.get(exp_name, set()) + if not train_keys: continue - sub_mask = group["fov_name"].isin(train_fovs).to_numpy() + sub_mask = group[split_key].isin(train_keys).to_numpy() mask[group.index.to_numpy()] = sub_mask return mask diff --git a/applications/dynaclr/src/dynaclr/data/dataset.py b/applications/dynaclr/src/dynaclr/data/dataset.py index a682b744e..aa172a92d 100644 --- a/applications/dynaclr/src/dynaclr/data/dataset.py +++ b/applications/dynaclr/src/dynaclr/data/dataset.py @@ -360,6 +360,7 @@ def _cache_columns(df: pd.DataFrame, columns: list[str]) -> dict: "norm_std", "norm_median", "norm_iqr", + "z_focus", } if self.positive_match_columns: hot_cols.update(self.positive_match_columns) @@ -717,9 +718,13 @@ def _build_norm_meta( ------- NormMeta or None """ - # Parquet path: norm columns present and value is not NA + # Parquet fast-path: one norm_* row = one channel's stats. Only valid in + # bag-of-channels mode (one channel per sample, keyed "channel_0"). In + # all-channels / fixed mode a sample reads multiple channels but the row + # carries only its own channel's stats, so fall through to the zarr + # zattrs path below, which returns every channel's stats. norm_mean_arr = arrays.get("norm_mean") - if norm_mean_arr is not None: + if self._channel_mode == "from_index" and norm_mean_arr is not None: norm_mean = norm_mean_arr[idx] if norm_mean is not None and not (isinstance(norm_mean, float) and np.isnan(norm_mean)): tp_stats = { @@ -728,12 +733,7 @@ def _build_norm_meta( "median": torch.tensor(arrays["norm_median"][idx], dtype=torch.float32), "iqr": torch.tensor(arrays["norm_iqr"][idx], dtype=torch.float32), } - if self._channel_mode == "from_index": - return {"channel_0": {"timepoint_statistics": tp_stats}} - else: - ch_arr = arrays.get("channel_name") - ch_name = ch_arr[idx] if ch_arr is not None else "channel_0" - return {ch_name: {"timepoint_statistics": tp_stats}} + return {"channel_0": {"timepoint_statistics": tp_stats}} # Fallback: read from zarr zattrs (old parquets without norm columns) store_path = arrays["store_path"][idx] @@ -822,13 +822,30 @@ def _slice_patch( channel_names_to_read = exp.channel_names channel_indices = [exp.channel_names.index(name) for name in channel_names_to_read] - # Per-experiment z_range (scale-adjusted window size centered on z_range center) + # Z window sizing comes from the per-experiment z_range; the center is + # the per-sample focus plane. Single source of truth: the parquet + # ``z_focus`` column when populated (written by preprocess-cell-index), + # otherwise fall back to the z_range center (which the registry derived + # from zattrs focus_slice, else mid-stack). z_focus_offset sets the + # fraction of the window placed below the focus plane (0.5 = symmetric). z_start_base, z_end_base = self.index.registry.z_ranges[exp_name] z_window_size = z_end_base - z_start_base z_count = round(z_window_size * scale_z) - z_focus = (z_start_base + z_end_base) // 2 - z_start = z_focus - z_count // 2 + z_focus_arr = arrays.get("z_focus") + z_focus_val = z_focus_arr[idx] if z_focus_arr is not None else None + if z_focus_val is not None and not (isinstance(z_focus_val, float) and np.isnan(z_focus_val)): + z_center = int(round(float(z_focus_val))) + z_below = round(z_count * self.index.registry.z_focus_offset) + z_start = z_center - z_below + else: + z_start = (z_start_base + z_end_base) // 2 - z_count // 2 z_end = z_start + z_count + # Clamp the window inside the image so edge cells still yield a full patch. + z_total = image.shape[2] + if z_start < 0: + z_start, z_end = 0, z_count + elif z_end > z_total: + z_start, z_end = z_total - z_count, z_total patch = image.oindex[ t, [int(c) for c in channel_indices], diff --git a/applications/dynaclr/src/dynaclr/data/preprocess_cell_index.py b/applications/dynaclr/src/dynaclr/data/preprocess_cell_index.py index e49bc1c74..64c87e178 100644 --- a/applications/dynaclr/src/dynaclr/data/preprocess_cell_index.py +++ b/applications/dynaclr/src/dynaclr/data/preprocess_cell_index.py @@ -17,7 +17,14 @@ default=None, help="Channel name for focus_slice lookup (e.g. Phase3D). Default: first channel per FOV.", ) -def main(parquet_path, output, focus_channel): +@click.option( + "--focus-level", + type=click.Choice(["fov", "per_timepoint"]), + default="fov", + show_default=True, + help="focus_slice level written to the z_focus column: per-FOV mean or per-timepoint index.", +) +def main(parquet_path, output, focus_channel, focus_level): """Preprocess a cell index parquet: add normalization stats, focus slice, remove empty frames. Reads precomputed metadata from zarr zattrs and writes them as parquet @@ -27,4 +34,5 @@ def main(parquet_path, output, focus_channel): parquet_path=parquet_path, output_path=output, focus_channel=focus_channel, + focus_level=focus_level, ) diff --git a/packages/viscy-data/src/viscy_data/_typing.py b/packages/viscy-data/src/viscy_data/_typing.py index 0e6baf6cc..0ed9125a0 100644 --- a/packages/viscy-data/src/viscy_data/_typing.py +++ b/packages/viscy-data/src/viscy_data/_typing.py @@ -254,7 +254,7 @@ class TripletSample(TypedDict): "Z_shape", "Y_shape", "X_shape", - "z_focus_mean", + "z_focus", ] CELL_INDEX_NORMALIZATION_COLUMNS = [ diff --git a/packages/viscy-data/src/viscy_data/cell_index.py b/packages/viscy-data/src/viscy_data/cell_index.py index 09a72f64f..ba869067a 100644 --- a/packages/viscy-data/src/viscy_data/cell_index.py +++ b/packages/viscy-data/src/viscy_data/cell_index.py @@ -82,7 +82,7 @@ ("Z_shape", pa.int32()), ("Y_shape", pa.int32()), ("X_shape", pa.int32()), - ("z_focus_mean", pa.float32()), + ("z_focus", pa.float32()), ("norm_mean", pa.float32()), ("norm_std", pa.float32()), ("norm_median", pa.float32()), @@ -238,6 +238,7 @@ def preprocess_cell_index( parquet_path: str | Path, output_path: str | Path | None = None, focus_channel: str | None = None, + focus_level: str = "fov", ) -> None: """Add normalization stats, focus slice, and remove invalid rows. @@ -246,7 +247,8 @@ def preprocess_cell_index( - ``norm_mean``, ``norm_std``, ``norm_median``, ``norm_iqr``, ``norm_max``, ``norm_min`` — per-timepoint, per-channel statistics - - ``z_focus_mean`` — per-FOV focus plane from ``focus_slice`` + - ``z_focus`` — the focus plane the training Z window is centered on, + sourced from ``focus_slice`` at the level chosen by ``focus_level``. Drops rows where timepoint stats are missing or ``norm_max == 0.0`` (empty frames). The processed parquet is written to ``output_path``; @@ -262,12 +264,20 @@ def preprocess_cell_index( focus_channel : str | None Channel name for ``focus_slice`` lookup (e.g. ``"Phase3D"``). When ``None``, uses the first channel_name in each FOV's group. + focus_level : {'fov', 'per_timepoint'} + Which ``focus_slice`` level feeds the ``z_focus`` column: + ``'fov'`` uses the per-FOV ``fov_statistics.z_focus_mean``; + ``'per_timepoint'`` uses the per-timepoint focus index for each + sample's ``t``. The per-cell tracked ``z`` column is left unchanged. Raises ------ ValueError - If a FOV has no normalization metadata (run ``viscy preprocess`` first). + If a FOV has no normalization metadata (run ``viscy preprocess`` first), + or if ``focus_level`` is not one of the accepted values. """ + if focus_level not in ("fov", "per_timepoint"): + raise ValueError(f"focus_level must be 'fov' or 'per_timepoint', got {focus_level!r}") if output_path is None: output_path = parquet_path @@ -315,8 +325,9 @@ def preprocess_cell_index( t_arr = df["t"].astype(int).to_numpy() norm_arrays = {stat: np.full(len(df), float("nan"), dtype=np.float32) for stat in stat_keys} + # z_focus is the plane the training Z window centers on. The per-cell + # tracked ``z`` column is left untouched — it is a distinct quantity. focus_arr = np.full(len(df), float("nan"), dtype=np.float32) - z_arr = df["z"].to_numpy(dtype=np.int16).copy() valid_mask = np.ones(len(df), dtype=bool) for i in range(len(df)): @@ -327,17 +338,18 @@ def preprocess_cell_index( for stat in stat_keys: norm_arrays[stat][i] = float(tp_stats[stat]) fov_key = (store_arr[i], fov_arr[i]) - z_focus = focus_lookup.get(fov_key) - if z_focus is not None: - focus_arr[i] = z_focus - z_t = focus_per_t_lookup.get(fov_key, {}).get(t_arr[i]) - if z_t is not None: - z_arr[i] = z_t + if focus_level == "per_timepoint": + z_t = focus_per_t_lookup.get(fov_key, {}).get(t_arr[i]) + if z_t is not None: + focus_arr[i] = z_t + else: # "fov" + z_focus = focus_lookup.get(fov_key) + if z_focus is not None: + focus_arr[i] = z_focus for stat in stat_keys: df[f"norm_{stat}"] = norm_arrays[stat] - df["z_focus_mean"] = focus_arr - df["z"] = z_arr + df["z_focus"] = focus_arr df = df[valid_mask].reset_index(drop=True) n_dropped = n_before - len(df)