multigedi.tools.multigedi

multigedi.tools.multigedi(mdata, modalities, sample_key, *, K=10, seed=None, mode='Bsphere', max_iterations=30, track_interval=1, num_threads=0, verbose=1, use_gpu=False, gpu_device=0, gpu_output_dir=None, batch_size=-1, key_added='multigedi', copy=False)[source]

Run multi-modal GEDI (MultiGEDI) integration.

Trains a hub-and-spoke MultiGEDI model across multiple data modalities, learning shared cell loadings (joint B) and modality-specific metagenes.

Parameters:

  • mdata (MuData) – MuData object. Each key in modalities must be a valid modality name accessible as mdata[name].

  • modalities (dict[str, dict[str, Any]]) –

    Dict mapping modality name → configuration dict. Each config may have: - obs_type : "M" (default), "M_list", or "X" - orthoZ : bool (default True) - fixed_si : scalar fixed cell-effect value, matching R’s

    fixed_si option

    • layer : layer name to use instead of .X (for "M" / "X")

    • layers : 2-tuple of layer names for "M_list" (first may be None to use .X)

    • Z_shrinkage, A_shrinkage, Qi_shrinkage, Rk_shrinkage, oi_shrinkage, o_shrinkage, and si_shrinkage: finite positive scalar regularization strengths (all default to 1, matching R)

  • sample_key (str) – Column in mdata[name].obs containing sample/batch labels. Whole sample blocks may be absent from a modality. For every sample present in more than one modality, those modalities must have the same obs_names in the same within-sample order. The function validates this identity contract before fitting and raises rather than silently coupling different cells.

  • K (int, default: 10) – Number of latent factors.

  • seed (int | None, default: None) – Random seed for reproducible per-modality rSVD initialization.

  • mode (Literal['Bsphere', 'Bl2'], default: 'Bsphere') – B normalization: "Bsphere" (recommended) or "Bl2".

  • max_iterations (int, default: 30) – Number of optimization iterations.

  • track_interval (int, default: 1) – Record sigma² and dZ every this many iterations.

  • num_threads (int, default: 0) – OpenMP thread count (0 = auto).

  • verbose (int, default: 1) – Verbosity level (0 = silent, 1 = normal).

  • use_gpu (bool, default: False) – When True, run optimization on GPU via the CUDA backend. Requires a GPU-enabled build containing the _multigedi_gpu_py extension and a visible CUDA device. Supports 1..N modalities of any combination of obs_type (“M”, “M_list”, “X”), up to the backend’s 16-modality cap.

  • gpu_device (int, default: 0) – CUDA device index (default 0).

  • gpu_output_dir (str | None, default: None) – Reserved (no effect in v1.1; the in-memory pipeline does not write a binary directory). Kept in the signature for backward compatibility with notebooks that still pass it.

  • batch_size (int, default: -1) – GPU solve_Z cell-batch size (Tier-1 BCD batching). Only consulted when use_gpu=True; ignored on the CPU path. -1 (default) selects the single-shot solve_Z below the automatic cell ceiling. A positive value tiles the transient J×N residual scratch into cell-batches. When left at -1 and the total cell count exceeds ~500K, the runner auto-selects a batch size, emits a warning, and records that effective value in the result. Passing 0 or any value < -1 is rejected. See docs/architecture/bcd_batching.md.

  • key_added (str, default: 'multigedi') – Key under which results are stored in mdata.uns and used in each modality’s obsm coordinate name.

  • copy (bool, default: False) – Unused (reserved for API symmetry). Always returns None.

Return type:

None

Returns:

  • None. Results are stored on the MuData container per scverse

  • convention for multi-modal joint analyses

  • - ``mdata.uns[key_added]` — single source` of truth –

    • params — training configuration, including backend/device and requested/effective GPU batch sizes

    • model — joint, plus per-modality Z, D, Bi, o, sigma2, tracking. joint.global_cell_metadata contains row-aligned sample and cell_id lists for joining the shared embedding without relying on row position alone.

    • svd — joint d / u / v from joint SVD

  • - ``mdata[name].obsm[f``”X_{key_added}_pca”:py:class:`] — joint PCA coordinates` – (equal to V * d) sliced to that modality’s local cells, so sc.pp.neighbors(mdata[name], use_rep="X_multigedi_pca") works even when modalities cover different sample subsets.

  • For backwards compatibility a copy of mdata.uns[key_added] is

  • also written to mdata[first_modality].uns[key_added] with a

  • DeprecationWarning. The first-modality location will be removed

  • in a future release; new code should read from ``mdata.uns`.`

Examples

>>> import multigedi as gd
>>> import mudata as md
>>> mdata = md.read_h5mu("data.h5mu")
>>> gd.tl.multigedi(
...     mdata,
...     modalities={
...         "gene":     {"obs_type": "M",      "orthoZ": True},
...         "splicing": {"obs_type": "M_list", "orthoZ": False,
...                      "layers": (None, "M2")},
...     },
...     sample_key="sample",
...     K=20,
...     max_iterations=30,
... )
>>> mdata["gene"].obsm["X_multigedi_pca"]  # joint PCA rows for gene cells

Notes

Reproducibility & relationship to R `multigedi`. The BCD optimizer underlying tl.multigedi implements the R reference model. Given the same exported iter_0 state, local CPU/R comparisons agree to approximately 1e-12. However an independent tl.multigedi run does NOT reproduce R’s exact iter_N factorization — the rSVD initializer draws its random projection matrix from numpy’s RNG, while R draws from R’s Mersenne-Twister. Different initial projections route BCD into different deep local minima of the non-convex objective. Both minima have nearly identical sigma², but the loadings Z, the cell embedding X_multigedi_pca, and the pairwise cell distances differ across backends. Within a single backend, identical mdata and same MultiGEDIModel(seed=...) produce identical output. Pick one backend per dataset; do not mix R-trained and Python-trained outputs downstream. For exact agreement with a specific R reference run to the documented tolerance, load R’s exported iter_0 via MultiGEDIModel.add_modality(..., init_state=ref).

``batch_size`` is part of run identity. On the GPU path, the batch_size argument becomes part of the run’s identity: a run is only reproducible against another run that used the same effective batch_size. batch_size=-1 (the default) uses the single-shot solve_Z below the automatic ceiling. A positive batch_size tiles the solve_Z residual reduction; re-associating that sum in IEEE-754 f64 shifts the trailing ULPs, so two runs with different batch_size values land ≤1e-11 apart — the same way two different seeds would. This drift is expected and joins the existing multigedi numerical-divergence family; it must not be “fixed”. Do not mix outputs produced with different batch_size values downstream. When the GPU runner auto-selects a batch_size (cell count above the ~500K single-GPU solve_Z ceiling), it emits a warning and stores the chosen value as params["effective_batch_size"].