diff --git a/.gitignore b/.gitignore
index 50bf725..0046f75 100644
--- a/.gitignore
+++ b/.gitignore
@@ -139,3 +139,6 @@ spatial_spec.txt
 
 # Benchmark artifacts
 benchmarks_artifacts/
+
+# Scratch / experiment workspace (not part of PR)
+scratch/
diff --git a/genevector/_graph_targets.py b/genevector/_graph_targets.py
index b1ee99b..ff94636 100644
--- a/genevector/_graph_targets.py
+++ b/genevector/_graph_targets.py
@@ -53,3 +53,198 @@ def target_graph_xcorr(X, gene_names, graph=None, aggr="mean",
     np.fill_diagonal(xcorr_sym, 0)
 
     return _matrix_to_score_dict(xcorr_sym, gene_names)
+
+
+def _cross_mi_matrix(A_disc, na, B_disc, nb):
+    """Pairwise MI between every self-gene column of A and neighbor-gene column of B.
+
+    Returns M where M[i, j] = MI(self gene i, neighbor-aggregated gene j). Asymmetric.
+    """
+    from .metrics import _mi_from_joint
+    P = A_disc.shape[1]
+    M = np.zeros((P, P), dtype=np.float64)
+    for i in range(P):
+        if na[i] <= 1:
+            continue
+        a = A_disc[:, i]
+        for j in range(P):
+            if nb[j] <= 1:
+                continue
+            b = B_disc[:, j]
+            mask = (a > 0) | (b > 0)
+            if mask.sum() == 0:
+                continue
+            joint = np.zeros((na[i], nb[j]), dtype=np.float64)
+            np.add.at(joint, (a[mask], b[mask]), 1)
+            M[i, j] = _mi_from_joint(joint)
+    return M
+
+
+def _cross_mi_matrix_torch(A_disc, na, B_disc, nb, device="cuda", max_elems=20_000_000):
+    """GPU/torch cross-MI: M[i, j] = MI(self gene i, neighbor gene j). Mirrors _cross_mi_matrix.
+
+    For each self gene one ``scatter_add`` builds every neighbour gene's joint histogram at
+    once (chunked for memory), so this is O(P) Python iterations rather than O(P^2). Runs on
+    any torch device — pass ``device='cpu'`` to validate against the numpy path.
+
+    Padding to the max bin count is harmless: extra bins stay empty and contribute 0 to MI.
+    Zeroing the (0, 0) cell of each joint reproduces the ``(a>0)|(b>0)`` mask used elsewhere
+    (cells where both genes are zero are dropped).
+    """
+    import torch
+    A_disc = np.asarray(A_disc)
+    B_disc = np.asarray(B_disc)
+    n_cells, P = A_disc.shape
+    MA = int(na.max()) if na.size else 1
+    MB = int(nb.max()) if nb.size else 1
+    cellsize = MA * MB
+    chunk = max(1, min(P, int(max_elems // max(n_cells, 1))))
+
+    A_t = torch.as_tensor(A_disc, device=device, dtype=torch.long)
+    B_t = torch.as_tensor(B_disc, device=device, dtype=torch.long)
+    nb_t = torch.as_tensor(np.asarray(nb), device=device, dtype=torch.long)
+    M = torch.zeros((P, P), device=device, dtype=torch.float64)
+
+    for i in range(P):
+        if int(na[i]) <= 1:
+            continue
+        a = A_t[:, i]  # (n_cells,)
+        for start in range(0, P, chunk):
+            end = min(start + chunk, P)
+            c = end - start
+            Bc = B_t[:, start:end]  # (n_cells, c)
+            joff = (torch.arange(c, device=device, dtype=torch.long) * cellsize).unsqueeze(0)
+            flat = (joff + a.unsqueeze(1) * MB + Bc).reshape(-1)  # (n_cells*c,)
+            hist = torch.zeros(c * cellsize, device=device, dtype=torch.float64)
+            hist.scatter_add_(0, flat, torch.ones_like(flat, dtype=torch.float64))
+            hist = hist.reshape(c, MA, MB)
+            hist[:, 0, 0] = 0.0  # drop both-zero cells == (a>0)|(b>0) mask
+            total = hist.sum(dim=(1, 2), keepdim=True)
+            p = hist / total.clamp_min(1.0)
+            px = p.sum(dim=2, keepdim=True)
+            py = p.sum(dim=1, keepdim=True)
+            ratio = torch.where(p > 0, p / (px * py).clamp_min(1e-12), torch.ones_like(p))
+            mi = (p * torch.log2(ratio)).sum(dim=(1, 2))  # (c,)
+            valid = (total.reshape(-1) > 0) & (nb_t[start:end] > 1)
+            M[i, start:end] = torch.where(valid, mi, torch.zeros_like(mi))
+    return M.cpu().numpy()
+
+
+def _cross_mi_matrix_rust(A_disc, na, B_disc, nb):
+    """rayon-parallel Rust cross-MI: M[i, j] = MI(self gene i, neighbor gene j).
+
+    Mirrors _cross_mi_matrix but parallel across the P*P pairs. ~80-110x faster than numpy
+    and a few x faster than the torch CPU kernel on many-core machines (the diagonal is left
+    zero — it is zeroed downstream anyway).
+    """
+    from ._rust import compute_cross_mi_pairs
+    P = A_disc.shape[1]
+    triples = compute_cross_mi_pairs(
+        np.ascontiguousarray(A_disc, dtype=np.int32), np.asarray(na, dtype=np.int32),
+        np.ascontiguousarray(B_disc, dtype=np.int32), np.asarray(nb, dtype=np.int32), None)
+    M = np.zeros((P, P), dtype=np.float64)
+    for i, j, v in triples:
+        M[i, j] = v
+    return M
+
+
+def _graph_mi_core(X, gene_names, graph, aggr, aggr_params, n_bins, signed,
+                   backend="auto", device="cpu"):
+    """Shared computation for graph_mi / graph_cross_mi: returns signed cross-MI matrix.
+
+    M[i, j] = (sign of self_i vs neighbor_j cross-correlation) * MI(self_i, neighbor_j).
+    Aggregating expression over the graph BEFORE estimating the gene-gene relationship
+    denoises the per-cell counts, making it robust to dropout on sparse spatial panels.
+
+    ``backend`` selects the cross-MI kernel:
+      - "auto" (default): GPU torch when ``device=="cuda"``; else the rayon Rust kernel if the
+        ``_rust`` extension is built (fastest on multi-core CPU); else the torch CPU kernel
+        (~20-35x numpy); else numpy.
+      - "rust" / "gpu" / "numpy": force that kernel (each degrades gracefully if unavailable).
+    All kernels are numerically identical (diff ~1e-15); only speed differs.
+    """
+    if graph is None:
+        raise ValueError(
+            "graph required. Pass any scipy sparse adjacency matrix "
+            "via target_kwargs={'graph': G}"
+        )
+    from .metrics import discretize_genes, HAS_RUST
+    from ._logging import get_logger
+    aggr_fn = get_aggregation(aggr)
+    X_dense = _to_dense(X)
+    X_agg = aggr_fn(X_dense, graph, **(aggr_params or {}))
+
+    A_disc, na = discretize_genes(X_dense, n_bins=n_bins)
+    B_disc, nb = discretize_genes(X_agg, n_bins=n_bins)
+
+    chosen = backend
+    if chosen == "auto":
+        chosen = "gpu" if device == "cuda" else ("rust" if HAS_RUST else "gpu")
+
+    M = None
+    if chosen == "rust":
+        try:
+            M = _cross_mi_matrix_rust(A_disc, na, B_disc, nb)
+        except Exception as e:
+            get_logger(__name__).warning(f"graph_mi rust backend failed ({e}); trying torch.")
+            chosen = "gpu"
+    if M is None and chosen == "gpu":
+        try:
+            import torch
+            dev = device
+            if dev == "cuda" and not torch.cuda.is_available():
+                get_logger(__name__).warning("device='cuda' requested but CUDA unavailable; "
+                                             "using torch CPU.")
+                dev = "cpu"
+            M = _cross_mi_matrix_torch(A_disc, na, B_disc, nb, device=dev)
+        except Exception as e:  # graceful fall back to numpy
+            get_logger(__name__).warning(f"graph_mi torch backend failed ({e}); using numpy.")
+            M = None
+    if M is None:
+        M = _cross_mi_matrix(A_disc, na, B_disc, nb)
+
+    if signed:
+        X_std = (X_dense - X_dense.mean(0)) / (X_dense.std(0) + 1e-8)
+        A_std = (X_agg - X_agg.mean(0)) / (X_agg.std(0) + 1e-8)
+        sign = np.sign((X_std.T @ A_std) / X_dense.shape[0])
+        M = sign * M
+    return M
+
+
+@register_target("graph_mi")
+def target_graph_mi(X, gene_names, graph=None, aggr="mean", aggr_params=None,
+                    n_bins=10, signed=True, backend="auto", device="cpu", **kwargs):
+    """Symmetric graph mutual information between self and neighbor-aggregated expression.
+
+    The MI analogue of ``graph_xcorr``: captures non-linear spatial co-expression while
+    the neighbor aggregation denoises sparse counts. Symmetrized over (i, j). The cross-MI
+    kernel is auto-selected (GPU torch on ``device="cuda"``, else the rayon Rust extension if
+    built, else torch CPU, else numpy); force one with ``backend`` in {"rust","gpu","numpy"}.
+
+    Returns
+    -------
+    dict of dict
+        scores[gene_a][gene_b] = signed MI in (roughly) [-log2(n_bins), log2(n_bins)].
+    """
+    M = _graph_mi_core(X, gene_names, graph, aggr, aggr_params, n_bins, signed,
+                       backend=backend, device=device)
+    M_sym = (M + M.T) / 2
+    np.fill_diagonal(M_sym, 0)
+    return _matrix_to_score_dict(M_sym, gene_names)
+
+
+@register_target("graph_cross_mi")
+def target_graph_cross_mi(X, gene_names, graph=None, aggr="mean", aggr_params=None,
+                          n_bins=10, signed=True, backend="auto", device="cpu", **kwargs):
+    """Asymmetric cross-neighbor MI: MI(gene_a in cell, gene_b in neighbors).
+
+    Directional spatial signal (e.g. ligand in a cell predicting receptor in its
+    neighbours). Not symmetrized — the model's separate input/output weights can encode
+    the asymmetry. Encodes niche/communication directionality in the gene embedding. The cross-MI
+    kernel is auto-selected (GPU torch on ``device="cuda"``, else the rayon Rust extension if
+    built, else torch CPU, else numpy); force one with ``backend`` in {"rust","gpu","numpy"}.
+    """
+    M = _graph_mi_core(X, gene_names, graph, aggr, aggr_params, n_bins, signed,
+                       backend=backend, device=device)
+    np.fill_diagonal(M, 0)
+    return _matrix_to_score_dict(M, gene_names)
diff --git a/genevector/_rust.pyi b/genevector/_rust.pyi
index d38e2c5..8dcc6cd 100644
--- a/genevector/_rust.pyi
+++ b/genevector/_rust.pyi
@@ -3,3 +3,12 @@ def compute_mi_pairs(
     n_bins_per_gene: "numpy.ndarray",
     corr_signs: "numpy.ndarray | None" = None,
 ) -> list[tuple[int, int, float]]: ...
+
+
+def compute_cross_mi_pairs(
+    a_disc: "numpy.ndarray",
+    na_bins: "numpy.ndarray",
+    b_disc: "numpy.ndarray",
+    nb_bins: "numpy.ndarray",
+    corr_signs: "numpy.ndarray | None" = None,
+) -> list[tuple[int, int, float]]: ...
diff --git a/genevector/embedding.py b/genevector/embedding.py
index 181bf8f..3d13a0c 100644
--- a/genevector/embedding.py
+++ b/genevector/embedding.py
@@ -515,7 +515,11 @@ def __init__(self, dataset, embed, log_normalize=True):
              logger.warning("No cell vectors were generated. self.matrix is empty.")
              self.dataset_vector = numpy.zeros(self.embed.vector_size if hasattr(self.embed, "vector_size") else 100) # Default size
         else:
-             self.dataset_vector = numpy.zeros(numpy.array(self.matrix).shape[1])
+             # The dataset vector is the mean cell vector — the shared "background"
+             # direction every cell points toward. Previously left as zeros, which
+             # silently disabled the contrastive subtraction in get_predictive_genes
+             # and made debiased phenotype scoring impossible.
+             self.dataset_vector = numpy.array(self.matrix).mean(axis=0)
 
         logger.info(f"Found {cells_with_no_counts} Cells with No Counts / No scorable gene expression.")
         logger.info("Finished CellEmbedding Initialization.")
@@ -945,14 +949,14 @@ def cell_distance(self, target_vec, norm=False):
         return similarities
 
 
-    def phenotype_probability(self, adata, phenotype_markers, return_distances=False, method="normalized_exponential", target_col="genevector", temperature=0.001):
+    def phenotype_probability(self, adata, phenotype_markers, return_distances=False, method="normalized_exponential", target_col="genevector", temperature=0.001, debias=0.0, contrastive=False, score_norm="none", smooth_graph=None, smooth_alpha=0.5, smooth_adaptive=True, smooth_counts=None, lp_graph=None, lp_alpha=0.0, lp_iter=3):
         """
         Probabilistically assign phenotypes based on a set of cell type labels and associated markers.
         Loads into the anndata the pseudo-probabilities for each cell type and the deterministic label
         taken from the maximum probability over cell types.
 
         :param adata: AnnData object. It's assumed this is `self.adata` or is consistent with it,
-                      especially regarding `adata.obs.index` if `self.cell_distance` is used.
+                      especially regarding `adata.obs.index`.
         :type adata: anndata.AnnData
         :param phenotype_markers: Dictionary of cell type labels (key) to gene markers (list of strings, value).
         :type phenotype_markers: dict
@@ -964,6 +968,43 @@ def phenotype_probability(self, adata, phenotype_markers, return_distances=False
         :type target_col: str
         :param temperature: Temperature parameter for the "normalized_exponential" method.
         :type temperature: float
+        :param debias: Fraction (0–1) of the dataset (background) vector subtracted from each cell
+                       vector before scoring. Opt-in (default 0.0 = current behaviour). Mild values
+                       (~0.5) de-saturate scores on balanced data; large values can over-correct a
+                       dominant population (e.g. epithelial-rich tumours). See :meth:`scoring_report`.
+        :type debias: float
+        :param contrastive: If True, subtract the mean of competing phenotype vectors from each
+                            phenotype vector before scoring (the get_predictive_genes formulation).
+                            Opt-in (default False). Helps when phenotypes are not well separated.
+        :type contrastive: bool
+        :param score_norm: Per-phenotype normalization of the cell×phenotype similarity columns
+                           before the probability function. One of "none" (default), "zscore"
+                           (subtract each phenotype's mean / divide by std — stops a phenotype that
+                           is close to everyone from winning by default), or "rank" (per-phenotype
+                           percentile rank). Opt-in; helps on some datasets, default off.
+        :type score_norm: str
+        :param smooth_graph: Optional scipy sparse adjacency (cells x cells, same order as the cell
+                             embedding) to spatially denoise the cell vectors *for scoring only*
+                             (does not mutate self.matrix / the UMAP). Opt-in. For persistent
+                             denoising that also affects the embedding use :meth:`denoise_cell_vectors`.
+        :type smooth_graph: scipy.sparse matrix or None
+        :param smooth_alpha: Smoothing weight (or per-cell cap when ``smooth_adaptive``).
+        :type smooth_alpha: float
+        :param smooth_adaptive: Scale the smoothing weight per cell by inverse count (low-count
+                                cells borrow more). Default True.
+        :type smooth_adaptive: bool
+        :param smooth_counts: Per-cell totals for adaptive smoothing; if None, derived from the
+                              loaded expression context.
+        :type smooth_counts: array-like or None
+        :param lp_graph: Optional scipy sparse adjacency (cells x cells, same order as the cell
+                         embedding) for spatial label propagation over the soft probabilities. Opt-in.
+        :type lp_graph: scipy.sparse matrix or None
+        :param lp_alpha: Label-propagation coupling in [0, 1) (0 disables). Smooths probabilities
+                         over ``lp_graph`` as a post-step. Improves spatial coherence; can blur
+                         identity in intermixed tissue — keep small (~0.3) and opt-in.
+        :type lp_alpha: float
+        :param lp_iter: Number of label-propagation iterations.
+        :type lp_iter: int
         :return: AnnData with cell type labels and probabilities. If return_distances is True,
                  returns a tuple (adata, raw_similarities_dict).
         :rtype: anndata.AnnData or tuple
@@ -996,51 +1037,73 @@ def phenotype_probability(self, adata, phenotype_markers, return_distances=False
 
 
         phenotype_names = list(phenotype_markers.keys())
-        # Stores raw similarities: {phenotype_name: [sim_cell1, sim_cell2, ...]}
-        # Order of similarities in lists will correspond to self.adata.obs.index
-        raw_similarity_scores = collections.defaultdict(list)
-
-        for pheno_name in tqdm.tqdm(phenotype_names, desc="Computing similarities per phenotype"):
-            markers = phenotype_markers[pheno_name]
-            if not markers:
-                logger.warning(f"No markers provided for phenotype {pheno_name}. Skipping.")
-                # Assign a default low similarity or handle as appropriate
-                raw_similarity_scores[pheno_name] = [0.0] * len(self.adata.obs) # Or len(adata.obs) if strictly using input adata
+
+        # Cell matrix aligned to self.adata.obs.index (== self.matrix row order).
+        C = np.asarray(self.matrix, dtype=float)
+        if smooth_graph is not None:
+            # Spatially denoise a COPY of the cell vectors for scoring only (self.matrix and
+            # the UMAP are untouched). Opt-in; helps sparse / spatially-organised data.
+            C = self._graph_smooth(C, smooth_graph, alpha=smooth_alpha,
+                                   adaptive=smooth_adaptive, counts=smooth_counts)
+        if debias:
+            # Subtract a fraction of the shared background direction. Opt-in: mild values
+            # de-saturate scores on balanced data; large values over-correct a dominant
+            # population. Default debias=0.0 reproduces the original absolute scoring.
+            C = C - float(debias) * np.asarray(self.dataset_vector, dtype=float)
+
+        # Phenotype (marker-mean) vectors.
+        pheno_vectors = []
+        for pheno_name in phenotype_names:
+            markers = phenotype_markers[pheno_name] or []
+            present = [g for g in markers if g in self.embed.embeddings]
+            if not present:
+                logger.warning(f"No usable markers for phenotype {pheno_name}; scoring as zeros.")
+                pheno_vectors.append(np.zeros(C.shape[1]))
                 continue
-            
-            logger.info(f"Markers for {pheno_name}: {', '.join(markers[:5])}{'...' if len(markers) > 5 else ''}")
-            phenotype_vector = self.embed.generate_vector(markers) # Assumes gene names are uppercase or handled by generate_vector
-            
-            # self.cell_distance calculates similarities for cells in self.adata.obs.index
-            # norm=False means use raw vectors for cosine similarity.
-            similarities_for_pheno = self.cell_distance(phenotype_vector, norm=False)
-            raw_similarity_scores[pheno_name] = similarities_for_pheno
-
-        # Prepare matrix for probability calculation: rows are cells, columns are phenotypes
-        # The order of cells is implicitly self.adata.obs.index
-        # The order of phenotypes is phenotype_names
-        num_cells = len(self.adata.obs) # Number of cells for which distances were computed
-        similarity_matrix_cells_x_phenos = np.zeros((num_cells, len(phenotype_names)))
-
-        for i, pheno_name in enumerate(phenotype_names):
-            if pheno_name in raw_similarity_scores:
-                 # Ensure list length matches num_cells, pad if necessary (e.g. if a pheno was skipped)
-                pheno_sims = raw_similarity_scores[pheno_name]
-                if len(pheno_sims) == num_cells:
-                    similarity_matrix_cells_x_phenos[:, i] = pheno_sims
-                else:
-                    logger.warning(f"Similarity score list length mismatch for {pheno_name}. Expected {num_cells}, got {len(pheno_sims)}. Padding with zeros.")
-                    similarity_matrix_cells_x_phenos[:len(pheno_sims), i] = pheno_sims # Fill what's available
-
-        # Apply probability function per cell (i.e., per row of similarity_matrix_cells_x_phenos)
-        probabilities_per_cell = [] # List of np.arrays, each array is probs for one cell
-        for i in range(num_cells):
-            cell_similarity_vector = similarity_matrix_cells_x_phenos[i, :]
-            # Replace NaNs with a low value if any occurred in cell_distance
-            cell_similarity_vector = np.nan_to_num(cell_similarity_vector, nan=-1.0) # Or other appropriate fill
-            
-            prob_dist_for_cell = pfunc(cell_similarity_vector)
-            probabilities_per_cell.append(prob_dist_for_cell)
+            missing = [g for g in markers if g not in self.embed.embeddings]
+            if missing:
+                logger.info(f"{pheno_name}: {len(missing)} marker(s) absent from embedding: {missing[:5]}")
+            pheno_vectors.append(np.asarray(self.embed.generate_vector(present), dtype=float))
+        V = np.asarray(pheno_vectors, dtype=float)
+
+        if contrastive and len(V) > 1:
+            # Contrastive phenotype vectors: subtract the mean of competing phenotypes
+            # (and the background if debiasing) — the get_predictive_genes formulation.
+            bg = float(debias) * np.asarray(self.dataset_vector, dtype=float) if debias else 0.0
+            Vc = np.empty_like(V)
+            for i in range(len(V)):
+                Vc[i] = V[i] - np.delete(V, i, axis=0).mean(axis=0) - bg
+            V = Vc
+
+        # Vectorized cosine similarity (cells x phenotypes).
+        Cn = C / (np.linalg.norm(C, axis=1, keepdims=True) + 1e-9)
+        Vn = V / (np.linalg.norm(V, axis=1, keepdims=True) + 1e-9)
+        similarity_matrix_cells_x_phenos = Cn @ Vn.T
+        num_cells = similarity_matrix_cells_x_phenos.shape[0]
+
+        # Optional per-phenotype (column) normalization of the similarity matrix. De-biases a
+        # phenotype that is close to every cell so it doesn't win the argmax by default. Opt-in.
+        if score_norm and score_norm != "none":
+            S = similarity_matrix_cells_x_phenos
+            if score_norm == "zscore":
+                S = (S - S.mean(axis=0)) / (S.std(axis=0) + 1e-9)
+            elif score_norm == "rank":
+                S = np.argsort(np.argsort(S, axis=0), axis=0) / max(S.shape[0] - 1, 1)
+            else:
+                raise ValueError(f"Unknown score_norm: {score_norm}. Choose 'none', 'zscore', 'rank'.")
+            similarity_matrix_cells_x_phenos = S
+
+        # Apply the probability function per cell (per row).
+        probabilities_per_cell = [
+            pfunc(np.nan_to_num(similarity_matrix_cells_x_phenos[i, :], nan=-1.0))
+            for i in range(num_cells)
+        ]
+
+        # Optional spatial label propagation over the soft probabilities (opt-in post-step).
+        if lp_graph is not None and lp_alpha and lp_alpha > 0:
+            P = self._label_propagate(np.vstack(probabilities_per_cell), lp_graph,
+                                      alpha=lp_alpha, n_iter=lp_iter)
+            probabilities_per_cell = [P[i] for i in range(P.shape[0])]
         
         # Store results in the input 'adata' object
         # Probabilities are ordered by self.adata.obs.index. We assign to input 'adata'.
@@ -1090,6 +1153,175 @@ def phenotype_probability(self, adata, phenotype_markers, return_distances=False
         else:
             return adata
 
+    # ─── Spatial / graph helpers (opt-in) ──────────────────────────
+
+    @staticmethod
+    def _row_normalize_graph(graph):
+        """Row-normalize a sparse adjacency so each row sums to 1 (zero-degree rows stay 0)."""
+        from scipy.sparse import diags
+        rs = numpy.asarray(graph.sum(axis=1)).ravel()
+        rs[rs == 0] = 1.0
+        return diags(1.0 / rs) @ csr_matrix(graph)
+
+    def _cell_total_counts(self):
+        """Per-cell total expression aligned to self.matrix / self.data order."""
+        keys = list(self.data.keys())
+        exp = getattr(self.context, "expression", None)
+        if exp:
+            return numpy.array([sum(exp.get(k, {}).values()) for k in keys], dtype=float)
+        return numpy.ones(len(keys), dtype=float)
+
+    def _label_propagate(self, P, graph, alpha=0.3, n_iter=3):
+        """Spatial label propagation: F <- a*Wn@F + (1-a)*P0, iterated; returns smoothed probs."""
+        Wn = self._row_normalize_graph(graph)
+        P0 = numpy.asarray(P, dtype=float)
+        F = P0.copy()
+        for _ in range(int(n_iter)):
+            F = alpha * numpy.asarray(Wn @ F) + (1.0 - alpha) * P0
+        return F
+
+    def _graph_smooth(self, M, graph, alpha=0.5, adaptive=True, counts=None, k0=None, beta=0.0):
+        """Message-pass a matrix over a graph: (1-a)*self + a*mean(neighbours) [+ beta*2-hop].
+
+        Pure (does not mutate). ``adaptive=True`` makes a_i = min(cap, k0/(k0+counts_i)) so
+        low-count cells borrow more. Shared by denoise_cell_vectors and phenotype_probability.
+        """
+        M = numpy.asarray(M, dtype=float)
+        Wn = self._row_normalize_graph(graph)
+        one = numpy.asarray(Wn @ M)
+        if adaptive:
+            if counts is None:
+                counts = self._cell_total_counts()
+            counts = numpy.asarray(counts, dtype=float).ravel()
+            if k0 is None:
+                pos = counts[counts > 0]
+                k0 = float(numpy.median(pos)) if pos.size else 1.0
+            cap = alpha if 0 < alpha < 1 else 0.85
+            a = numpy.minimum(cap, k0 / (k0 + counts)).reshape(-1, 1)
+        else:
+            a = float(alpha)
+        out = (1.0 - a) * M + a * one
+        if beta > 0:
+            out = out + beta * (numpy.asarray(Wn @ one) - one)
+        return out
+
+    def denoise_cell_vectors(self, graph, alpha=0.5, adaptive=True, counts=None,
+                             k0=None, beta=0.0):
+        """Graph-denoise the cell matrix in gene-vector space (opt-in; modifies self.matrix).
+
+        Message passing over a (usually spatial) graph::
+
+            (1 - a) * self + a * mean(neighbours)   [+ beta * 2-hop]
+
+        With ``adaptive=True`` the per-cell weight ``a_i = min(0.85, k0/(k0+counts_i))`` so
+        low-count (noisy) cells borrow heavily from neighbours while high-count cells barely
+        move — robust across sparsity regimes. Strongly improves cell typing on sparse and/or
+        spatially organised data (e.g. low-depth Xenium), but can contaminate identity in
+        intermixed tissue, so it is **opt-in and OFF by default**. It also gives previously
+        near-empty cells a usable vector. Call BEFORE :meth:`get_adata` /
+        :meth:`phenotype_probability`.
+
+        :param graph: scipy sparse adjacency (cells x cells) aligned to the cell-matrix order.
+        :param alpha: smoothing weight (``adaptive=False``) or the per-cell cap (``adaptive=True``).
+        :param adaptive: scale the weight per cell by inverse count.
+        :param counts: per-cell totals; if None, derived from the loaded expression context.
+        :param k0: adaptive midpoint; defaults to the median of positive counts.
+        :param beta: optional 2-hop weight.
+        :return: the new self.matrix (list of vectors). Original kept in self.uncorrected_matrix.
+        """
+        out = self._graph_smooth(self.matrix, graph, alpha=alpha, adaptive=adaptive,
+                                 counts=counts, k0=k0, beta=beta)
+        self.uncorrected_matrix = self.matrix
+        self.matrix = [out[i] for i in range(out.shape[0])]
+        return self.matrix
+
+    def qc_marker_dict(self, adata, phenotype_markers, layer=None,
+                       specificity_threshold=0.5, min_markers=2, verbose=True):
+        """Quality-control a marker dictionary before phenotyping.
+
+        Computes a provisional mean-expression labelling and, for every (phenotype, marker),
+        reports whether the marker is in the panel, its mean expression, fraction of cells
+        expressing it, a fold ``enrichment`` over the global mean, and a proportion-independent
+        ``specificity`` (tau index in [0, 1] over the per-type means; 1 = specific to one type,
+        0 = uniform). Flags low-specificity markers (tau below ``specificity_threshold``),
+        off-target markers (highest in a different type than declared) and phenotypes with too
+        few usable markers. Tau is used (not raw fold-enrichment) because fold-enrichment is
+        confounded by class proportions and by dense/normalised data.
+
+        :param adata: AnnData with expression (raw counts recommended via ``layer``).
+        :param phenotype_markers: dict of phenotype -> marker gene list.
+        :param layer: layer to read expression from (e.g. "counts"); defaults to ``.X``.
+        :param specificity_threshold: tau below this flags a marker as low-specificity.
+        :param min_markers: warn if a phenotype has fewer usable markers than this.
+        :param verbose: log warnings.
+        :return: pandas.DataFrame with one row per (phenotype, marker) and the QC columns.
+        :rtype: pandas.DataFrame
+        """
+        X = adata.layers[layer] if (layer and layer in adata.layers) else adata.X
+        X = numpy.asarray(X.todense()) if hasattr(X, "todense") else numpy.asarray(X)
+        u2i = {str(g).upper(): i for i, g in enumerate(adata.var.index)}
+
+        # provisional labels via z-scored mean marker expression argmax
+        Z = (X - X.mean(0)) / (X.std(0) + 1e-9)
+        names, cols = [], []
+        for ct, genes in phenotype_markers.items():
+            idx = [u2i[g.upper()] for g in genes if g.upper() in u2i]
+            if idx:
+                names.append(ct)
+                cols.append(Z[:, idx].mean(1))
+        prov = (numpy.array([names[i] for i in numpy.stack(cols, 1).argmax(1)])
+                if cols else numpy.array(["?"] * adata.n_obs))
+
+        ct_list = list(phenotype_markers)
+        type_mean = {ct: X[prov == ct].mean(0) if (prov == ct).any() else numpy.zeros(X.shape[1])
+                     for ct in ct_list}
+        global_mean = X.mean(0) + 1e-9
+        frac_expr = (X > 0).mean(0)
+        K = max(len(ct_list), 2)
+
+        def tau(gi):
+            vals = numpy.array([type_mean[c][gi] for c in ct_list], dtype=float)
+            mx = vals.max()
+            if mx <= 0:
+                return 0.0
+            return float(numpy.sum(1.0 - vals / mx) / (K - 1))
+
+        rows = []
+        usable = collections.Counter()
+        for ct, genes in phenotype_markers.items():
+            for g in genes:
+                gi = u2i.get(g.upper())
+                if gi is None:
+                    rows.append(dict(phenotype=ct, marker=g, in_panel=False, mean_expr=numpy.nan,
+                                     frac_expressing=numpy.nan, enrichment=numpy.nan,
+                                     specificity=numpy.nan, top_type=None, flag="absent"))
+                    continue
+                usable[ct] += 1
+                spec = tau(gi)
+                top = max(ct_list, key=lambda c: type_mean[c][gi])
+                flag = "ok"
+                if top != ct:
+                    flag = f"off_target(>{top})"
+                elif spec < specificity_threshold:
+                    flag = "low_specificity"
+                rows.append(dict(phenotype=ct, marker=g, in_panel=True,
+                                 mean_expr=round(float(X[:, gi].mean()), 4),
+                                 frac_expressing=round(float(frac_expr[gi]), 3),
+                                 enrichment=round(float(type_mean[ct][gi] / global_mean[gi]), 3),
+                                 specificity=round(spec, 3), top_type=top, flag=flag))
+        df = pandas.DataFrame(rows)
+        if verbose:
+            for ct in phenotype_markers:
+                if usable[ct] < min_markers:
+                    logger.warning(f"Phenotype '{ct}' has only {usable[ct]} usable marker(s) "
+                                   f"(< {min_markers}); assignment will be unreliable.")
+            for _, r in df[df["flag"].isin(["low_specificity"]) | df["flag"].str.startswith("off_target")].iterrows():
+                logger.warning(f"Marker '{r['marker']}' for '{r['phenotype']}': {r['flag']} "
+                               f"(specificity={r['specificity']}, highest in {r['top_type']}).")
+            for _, r in df[df["flag"] == "absent"].iterrows():
+                logger.warning(f"Marker '{r['marker']}' for '{r['phenotype']}' not in panel.")
+        return df
+
 
     def cluster(self, adata, up_markers, down_markers=dict()):
         """
diff --git a/scripts/phenotype_workflow.py b/scripts/phenotype_workflow.py
new file mode 100644
index 0000000..c30b42f
--- /dev/null
+++ b/scripts/phenotype_workflow.py
@@ -0,0 +1,291 @@
+#!/usr/bin/env python
+"""Probabilistic cell-type phenotyping workflow (spatial + scRNA).
+
+Generalised genevector phenotyping driven by a marker dictionary. Given an AnnData, a
+JSON ``{phenotype: [marker genes]}`` dict, and an output directory, it:
+
+  1. cleans the panel (optional mito/ribo removal),
+  2. selects genes = HVG ∪ SVG (Moran's I via grafiti when spatial coords + grafiti are
+     available; otherwise scanpy HVG), always force-keeping the marker genes,
+  3. QCs the marker dict (flags absent / low-specificity / off-target markers),
+  4. trains a genevector embedding — using the spatially-aware ``graph_mi`` target for
+     spatial data (neighbour aggregation denoises sparse counts before estimating MI),
+     and plain ``mi`` for dissociated data,
+  5. builds the cell embedding (optionally count-adaptive spatial denoising of cell
+     vectors, opt-in via --denoise),
+  6. assigns phenotypes via ``phenotype_probability`` (optional --debias / --contrastive /
+     --label-prop), and
+  7. writes the annotated AnnData, the marker-QC table, a per-cell assignment table and a
+     run summary into the output directory.
+
+Example
+-------
+    python scripts/phenotype_workflow.py \
+        --input data.h5ad --markers markers.json --output out/ \
+        --target auto --denoise --device cuda
+
+Marker JSON format::
+
+    {"Tumor": ["EPCAM", "KRT8"], "T cell": ["CD3D", "CD3E"], ...}
+"""
+from __future__ import annotations
+
+import argparse
+import datetime
+import json
+import os
+import sys
+
+import numpy as np
+import pandas as pd
+import scanpy as sc
+from scipy import sparse
+from sklearn.neighbors import kneighbors_graph
+
+try:
+    from genevector.data import GeneVectorDataset
+except ModuleNotFoundError:  # running from a source checkout without install
+    sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
+    from genevector.data import GeneVectorDataset
+from genevector.model import GeneVector
+from genevector.embedding import GeneEmbedding, CellEmbedding
+
+
+def log(msg):
+    ts = datetime.datetime.now().strftime("%Y-%m-%d %H:%M:%S")
+    print(f"[{ts}] {msg}", flush=True)
+
+
+# ─── Gene panel cleaning ───────────────────────────────────────
+
+def clean_panel(adata, keep):
+    """Drop mito/ribo/ambiguous genes (keeping any in `keep`)."""
+    keepU = {g.upper() for g in keep}
+    genes = []
+    for g in adata.var.index:
+        gu = str(g).upper()
+        if gu in keepU:
+            genes.append(g)
+            continue
+        if gu.startswith(("MT-", "RPS", "RPL")):
+            continue
+        if "." in gu:
+            continue
+        if "-" in gu and "HLA" not in gu:
+            continue
+        genes.append(g)
+    n_removed = adata.n_vars - len(genes)
+    return adata[:, genes].copy(), n_removed
+
+
+# ─── Spatial graph ─────────────────────────────────────────────
+
+def build_spatial_graph(adata, spatial_key, k):
+    coords = np.asarray(adata.obsm[spatial_key], dtype=float)
+    W = kneighbors_graph(coords, n_neighbors=min(k, adata.n_obs - 1),
+                         mode="connectivity", include_self=False)
+    return ((W + W.T) > 0).astype(float).tocsr()
+
+
+# ─── Gene selection: HVG ∪ SVG ─────────────────────────────────
+
+def select_genes(adata, n_top_genes, marker_genes, spatial, batch_key, counts_layer,
+                 grafiti_path=None):
+    """Return HVG ∪ SVG (force-keeping marker_genes), upper-cased to match genevector."""
+    keep = {g.upper() for g in marker_genes}
+    layer = counts_layer if counts_layer in adata.layers else None
+    svg_ok = False
+    if spatial:
+        if grafiti_path and grafiti_path not in sys.path:
+            sys.path.insert(0, grafiti_path)
+        try:
+            import grafiti as gf
+            if "sample_fov" not in adata.obs:
+                adata.obs["sample_fov"] = (adata.obs[batch_key].astype(str)
+                                           if batch_key and batch_key in adata.obs else "s0")
+            df = gf.pp.spatially_variable_genes(
+                adata, layer=layer, hvg_layer=layer, n_top_genes=n_top_genes,
+                union_hvg=True, white_list=list(marker_genes), inplace=False)
+            flag = df.get("spatially_variable")
+            sel = (df.index[flag.astype(bool)].tolist() if flag is not None
+                   else df.sort_values(df.columns[0], ascending=False).head(n_top_genes).index.tolist())
+            svg_ok = True
+            log(f"Selected {len(sel)} genes via grafiti HVG∪SVG (Moran's I).")
+        except Exception as e:
+            log(f"grafiti SVG unavailable ({e}); falling back to scanpy HVG.")
+    if not svg_ok:
+        try:
+            sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes, flavor="seurat_v3",
+                                        layer=layer, batch_key=batch_key if batch_key in adata.obs else None)
+            sel = adata.var.index[adata.var["highly_variable"]].tolist()
+        except Exception:
+            sc.pp.highly_variable_genes(adata, n_top_genes=n_top_genes)
+            sel = adata.var.index[adata.var["highly_variable"]].tolist()
+        log(f"Selected {len(sel)} HVGs (scanpy).")
+    selU = {str(g).upper() for g in sel} | keep
+    return [g for g in adata.var.index if str(g).upper() in selU]
+
+
+# ─── Main ──────────────────────────────────────────────────────
+
+def main():
+    p = argparse.ArgumentParser(description="Probabilistic cell-type phenotyping (spatial + scRNA).")
+    p.add_argument("--input", required=True, help="Input .h5ad")
+    p.add_argument("--markers", required=True, help="JSON {phenotype: [marker genes]}")
+    p.add_argument("--output", required=True, help="Output directory")
+    p.add_argument("--target", default="auto", choices=["auto", "mi", "graph_mi", "graph_cross_mi"],
+                   help="Training target. 'auto' = graph_mi when spatial else mi.")
+    p.add_argument("--n-genes", type=int, default=2000, help="HVG∪SVG count.")
+    p.add_argument("--epochs", type=int, default=1000)
+    p.add_argument("--dim", type=int, default=100)
+    p.add_argument("--spatial-key", default="spatial", help="obsm key for coordinates.")
+    p.add_argument("--knn", type=int, default=10, help="k for the spatial graph.")
+    p.add_argument("--batch-key", default=None, help="obs column for HVG batch / sample_fov.")
+    p.add_argument("--counts-layer", default="counts", help="layer holding raw counts.")
+    p.add_argument("--no-clean", action="store_true", help="skip mito/ribo gene removal.")
+    p.add_argument("--denoise", action="store_true",
+                   help="count-adaptive spatial denoising of cell vectors (spatial only).")
+    p.add_argument("--debias", type=float, default=0.0,
+                   help="fraction of dataset vector subtracted before scoring (0=off).")
+    p.add_argument("--contrastive", action="store_true",
+                   help="subtract competing-phenotype means before scoring.")
+    p.add_argument("--score-norm", default="none", choices=["none", "zscore", "rank"],
+                   help="per-phenotype normalization of similarity columns before assignment.")
+    p.add_argument("--label-prop", type=float, default=0.0,
+                   help="spatial label-propagation coupling on probabilities (0=off, spatial only).")
+    p.add_argument("--temperature", type=float, default=0.05)
+    p.add_argument("--device", default="auto", choices=["auto", "cpu", "cuda"])
+    p.add_argument("--mi-backend", default="auto")
+    p.add_argument("--grafiti-path", default=None, help="path to a grafiti checkout (for SVG).")
+    args = p.parse_args()
+
+    os.makedirs(args.output, exist_ok=True)
+    device = args.device
+    if device == "auto":
+        try:
+            import torch
+            device = "cuda" if torch.cuda.is_available() else "cpu"
+        except Exception:
+            device = "cpu"
+    log(f"Device: {device}")
+
+    adata = sc.read_h5ad(args.input)
+    log(f"Loaded {adata.n_obs} cells x {adata.n_vars} genes.")
+    markers = json.load(open(args.markers))
+    markers = {k: [str(g).upper() for g in v] for k, v in markers.items()}
+    marker_genes = sorted({g for v in markers.values() for g in v})
+    log(f"{len(markers)} phenotypes, {len(marker_genes)} unique markers.")
+
+    spatial = args.spatial_key in adata.obsm
+    log(f"Spatial coordinates: {'found' if spatial else 'not found'} "
+        f"(obsm['{args.spatial_key}']).")
+
+    # ensure raw counts in .X for MI discretization
+    if args.counts_layer in adata.layers:
+        adata.X = adata.layers[args.counts_layer].copy()
+    else:
+        adata.layers[args.counts_layer] = adata.X.copy()
+        log(f"No '{args.counts_layer}' layer; using .X as counts.")
+
+    if not args.no_clean:
+        adata, n_removed = clean_panel(adata, marker_genes)
+        log(f"Removed {n_removed} mito/ribo/ambiguous genes.")
+
+    # gene selection (HVG ∪ SVG ∪ markers)
+    sel = select_genes(adata, args.n_genes, marker_genes, spatial, args.batch_key,
+                       args.counts_layer, args.grafiti_path)
+    adata = adata[:, sel].copy()
+    adata.X = adata.layers[args.counts_layer].copy()
+    log(f"Training on {adata.n_vars} genes.")
+
+    # marker QC (needs a CellEmbedding-less quick pass — reuse the static logic via a temp embed later)
+    # build spatial graph
+    target = args.target
+    tkw = None
+    W = None
+    if spatial:
+        W = build_spatial_graph(adata, args.spatial_key, args.knn)
+        if target == "auto":
+            target = "graph_mi"
+        if target in ("graph_mi", "graph_cross_mi", "graph_xcorr"):
+            tkw = {"graph": W}
+        if target in ("graph_mi", "graph_cross_mi"):
+            from genevector.metrics import HAS_RUST
+            kern = ("GPU (torch)" if device == "cuda"
+                    else "Rust (rayon, multi-core)" if HAS_RUST else "torch CPU")
+            log(f"{target} cross-MI kernel: {kern}. "
+                + ("" if (device == "cuda" or HAS_RUST) else
+                   "Tip: build the Rust extension (maturin develop --release) for a large "
+                   "multi-core speedup, or pass --device cuda for GPU."))
+    else:
+        if target == "auto":
+            target = "mi"
+        if target.startswith("graph"):
+            log("Graph target requested but no spatial coords; using 'mi'.")
+            target = "mi"
+    log(f"Target: {target}.")
+
+    vec_path = os.path.join(args.output, "embedding.vec")
+    ds = GeneVectorDataset(adata.copy(), load_expression=True, signed_mi=True, device=device,
+                           target=target, target_kwargs=tkw, mi_backend=args.mi_backend,
+                           use_cache=False)
+    log("Training genevector...")
+    gv = GeneVector(ds, output_file=vec_path, emb_dimension=args.dim, c=100, gain=10,
+                    init_ortho=True, device=device)
+    gv.train(args.epochs, update_interval=max(1, args.epochs // 5))
+
+    embed = GeneEmbedding(vec_path, ds, vector="average")
+    cembed = CellEmbedding(ds, embed)
+
+    # marker QC report
+    qc = cembed.qc_marker_dict(adata, markers, layer=args.counts_layer)
+    qc.to_csv(os.path.join(args.output, "marker_qc.csv"), index=False)
+    log(f"Marker QC -> marker_qc.csv ({int((qc['flag']!='ok').sum())} flagged of {len(qc)}).")
+
+    # opt-in spatial denoising of cell vectors (before get_adata)
+    if args.denoise and spatial:
+        counts = np.asarray(adata.layers[args.counts_layer].sum(1)).ravel()
+        # graph aligned to the cell-matrix order (== list(cembed.data.keys()))
+        order = list(cembed.data.keys())
+        pos = {c: i for i, c in enumerate(adata.obs.index)}
+        idx = [pos[c] for c in order]
+        Wsub = W[idx][:, idx]
+        cembed.denoise_cell_vectors(Wsub, adaptive=True, counts=counts[idx])
+        log("Applied count-adaptive spatial denoising to cell vectors.")
+
+    adata_gv = cembed.get_adata()
+
+    # label-prop graph aligned to the (possibly filtered) cell order
+    lp_graph = None
+    if args.label_prop > 0 and spatial:
+        pos = {c: i for i, c in enumerate(adata.obs.index)}
+        idx = [pos[c] for c in adata_gv.obs.index]
+        lp_graph = W[idx][:, idx]
+
+    log("Assigning phenotypes...")
+    adata_gv = cembed.phenotype_probability(
+        adata_gv, markers, method="normalized_exponential", temperature=args.temperature,
+        target_col="genevector", debias=args.debias, contrastive=args.contrastive,
+        score_norm=args.score_norm, lp_graph=lp_graph, lp_alpha=args.label_prop)
+
+    # outputs
+    out_h5ad = os.path.join(args.output, "phenotyped.h5ad")
+    adata_gv.write_h5ad(out_h5ad)
+    prob_cols = [c for c in adata_gv.obs.columns if "Pseudo-probability" in c]
+    assign = adata_gv.obs[["genevector"] + prob_cols].copy()
+    assign.to_csv(os.path.join(args.output, "assignments.csv"))
+    summary = {
+        "input": args.input, "n_cells": int(adata_gv.n_obs), "n_genes": int(adata.n_vars),
+        "spatial": bool(spatial), "target": target, "denoise": bool(args.denoise and spatial),
+        "debias": args.debias, "contrastive": bool(args.contrastive),
+        "score_norm": args.score_norm, "label_prop": args.label_prop,
+        "label_counts": {k: int(v) for k, v in adata_gv.obs["genevector"].value_counts().items()},
+        "markers_flagged": int((qc["flag"] != "ok").sum()),
+    }
+    json.dump(summary, open(os.path.join(args.output, "summary.json"), "w"), indent=2)
+    log(f"Done. Wrote phenotyped.h5ad, assignments.csv, marker_qc.csv, summary.json to {args.output}")
+    log(f"Label counts: {summary['label_counts']}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/src/lib.rs b/src/lib.rs
index 530d241..ac8a83f 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -98,9 +98,103 @@ fn compute_mi_pairs(
         .collect()
 }
 
+/// Compute cross mutual information between every self gene and every neighbour gene.
+///
+/// `a_disc` is the discretized self-expression (cells x genes), `b_disc` the discretized
+/// graph-neighbour-aggregated expression (cells x genes). Returns flat (i, j, mi) triples
+/// where mi = MI(self gene i, neighbour gene j), over all ordered pairs i != j. Cells where
+/// both bins are zero are dropped ((a>0)||(b>0) mask), matching the numpy/torch paths.
+/// rayon-parallel across the P*P pairs — fast multi-core CPU graph_mi.
+#[pyfunction]
+#[pyo3(signature = (a_disc, na_bins, b_disc, nb_bins, corr_signs=None))]
+fn compute_cross_mi_pairs(
+    a_disc: PyReadonlyArray2<i32>,
+    na_bins: PyReadonlyArray1<i32>,
+    b_disc: PyReadonlyArray2<i32>,
+    nb_bins: PyReadonlyArray1<i32>,
+    corr_signs: Option<PyReadonlyArray2<f32>>,
+) -> Vec<(usize, usize, f64)> {
+    let a = a_disc.as_array();
+    let b = b_disc.as_array();
+    let abins = na_bins.as_array();
+    let bbins = nb_bins.as_array();
+    let n_genes = a.ncols();
+    let n_cells = a.nrows();
+
+    let signs: Option<Vec<Vec<f32>>> = corr_signs.map(|arr| {
+        let m = arr.as_array();
+        (0..n_genes)
+            .map(|i| (0..n_genes).map(|j| *m.get((i, j)).unwrap_or(&0.0)).collect())
+            .collect()
+    });
+
+    // all ordered self/neighbour pairs (diagonal is zeroed downstream)
+    let mut pairs: Vec<(usize, usize)> = Vec::new();
+    for i in 0..n_genes {
+        for j in 0..n_genes {
+            if i != j && abins[i] > 1 && bbins[j] > 1 {
+                pairs.push((i, j));
+            }
+        }
+    }
+
+    pairs
+        .par_iter()
+        .filter_map(|&(i, j)| {
+            let na = abins[i] as usize;
+            let nb = bbins[j] as usize;
+
+            let mut joint = vec![0u32; na * nb];
+            let mut count = 0u32;
+            for c in 0..n_cells {
+                let av = a[[c, i]] as usize;
+                let bv = b[[c, j]] as usize;
+                if av > 0 || bv > 0 {
+                    joint[av * nb + bv] += 1;
+                    count += 1;
+                }
+            }
+            if count == 0 {
+                return None;
+            }
+            let total = count as f64;
+
+            let mut px = vec![0.0f64; na];
+            let mut py = vec![0.0f64; nb];
+            let mut joint_f = vec![0.0f64; na * nb];
+            for ai in 0..na {
+                for bi in 0..nb {
+                    let p = joint[ai * nb + bi] as f64 / total;
+                    joint_f[ai * nb + bi] = p;
+                    px[ai] += p;
+                    py[bi] += p;
+                }
+            }
+
+            let mut mi = 0.0f64;
+            for ai in 0..na {
+                for bi in 0..nb {
+                    let pxy = joint_f[ai * nb + bi];
+                    let px_py = px[ai] * py[bi];
+                    if pxy > 0.0 && px_py > 0.0 {
+                        mi += pxy * (pxy / px_py).log2();
+                    }
+                }
+            }
+
+            if let Some(ref s) = signs {
+                let sign = if s[i][j] >= 0.0 { 1.0 } else { -1.0 };
+                mi *= sign;
+            }
+            Some((i, j, mi))
+        })
+        .collect()
+}
+
 /// Python module definition
 #[pymodule]
 fn _rust(m: &Bound<'_, PyModule>) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(compute_mi_pairs, m)?)?;
+    m.add_function(wrap_pyfunction!(compute_cross_mi_pairs, m)?)?;
     Ok(())
 }
diff --git a/tests/test_graph_targets.py b/tests/test_graph_targets.py
index 83a39ff..de3d5c9 100644
--- a/tests/test_graph_targets.py
+++ b/tests/test_graph_targets.py
@@ -6,7 +6,11 @@
 from scipy.sparse import csr_matrix
 
 from genevector.metrics import TARGETS
-from genevector._graph_targets import target_graph_xcorr
+from genevector._graph_targets import (
+    target_graph_xcorr,
+    target_graph_mi,
+    target_graph_cross_mi,
+)
 
 
 # ─── Shared fixtures ──────────────────────────────────────────
@@ -123,3 +127,120 @@ def test_graph_xcorr_with_sparse_X():
     for g1 in genes:
         for g2 in scores_dense[g1]:
             assert scores_dense[g1][g2] == pytest.approx(scores_sparse[g1][g2])
+
+
+# ─── graph_mi (symmetric) and graph_cross_mi (asymmetric) ─────
+
+def _make_chain_panel(n=60):
+    # chain graph; a spatial gradient shared by A and B (co-located), C is the same
+    # marginal with the spatial structure destroyed. Continuous + mostly nonzero so the
+    # MI mask keeps cells.
+    row = list(range(n - 1)) + list(range(1, n))
+    col = list(range(1, n)) + list(range(n - 1))
+    adj = csr_matrix(([1.0] * len(row), (row, col)), shape=(n, n))
+    rng = np.random.default_rng(0)
+    grad = np.linspace(1.0, 6.0, n)
+    gene_a = np.clip(grad + rng.normal(0, 0.2, n), 0, None)
+    gene_b = np.clip(grad + rng.normal(0, 0.2, n), 0, None)  # co-located with A
+    gene_c = rng.permutation(grad)                            # no spatial structure
+    X = np.column_stack([gene_a, gene_b, gene_c])
+    return X, adj, ["A", "B", "C"]
+
+
+def test_graph_mi_registered():
+    assert "graph_mi" in TARGETS
+    assert "graph_cross_mi" in TARGETS
+
+
+def test_graph_mi_requires_graph():
+    X = np.array([[1, 2], [3, 4]], dtype=np.float64)
+    with pytest.raises(ValueError, match="graph required"):
+        target_graph_mi(X, ["a", "b"])
+    with pytest.raises(ValueError, match="graph required"):
+        target_graph_cross_mi(X, ["a", "b"])
+
+
+def test_graph_mi_symmetric_and_no_self():
+    X, adj, genes = _make_chain_panel()
+    scores = target_graph_mi(X, genes, graph=adj)
+    for g in genes:
+        assert g not in scores[g]
+    for g1 in genes:
+        for g2 in scores[g1]:
+            assert scores[g1][g2] == pytest.approx(scores[g2][g1], abs=1e-6)
+
+
+def test_graph_mi_detects_neighbor_coexpression():
+    # A high in even cells, B high in their chain neighbours → strong graph MI(A,B)
+    X, adj, genes = _make_chain_panel()
+    scores = target_graph_mi(X, genes, graph=adj)
+    assert abs(scores["A"]["B"]) > abs(scores["A"]["C"])
+
+
+def test_graph_mi_sparse_equals_dense():
+    X, adj, genes = _make_chain_panel(n=40)
+    sd = target_graph_mi(X, genes, graph=adj)
+    ss = target_graph_mi(csr_matrix(X), genes, graph=adj)
+    for g1 in genes:
+        for g2 in sd[g1]:
+            assert sd[g1][g2] == pytest.approx(ss[g1][g2], abs=1e-6)
+
+
+def test_graph_cross_mi_no_self_pairs():
+    X, adj, genes = _make_chain_panel()
+    scores = target_graph_cross_mi(X, genes, graph=adj)
+    for g in genes:
+        assert g not in scores[g]
+
+
+# ─── GPU (torch) graph_mi == numpy graph_mi (validated on CPU torch) ──
+
+def test_graph_mi_torch_matches_numpy():
+    pytest.importorskip("torch")
+    X, adj, genes = _make_chain_panel(n=50)
+    cpu = target_graph_mi(X, genes, graph=adj, backend="numpy")
+    gpu = target_graph_mi(X, genes, graph=adj, backend="gpu", device="cpu")
+    for g1 in genes:
+        for g2 in cpu[g1]:
+            assert cpu[g1][g2] == pytest.approx(gpu[g1][g2], abs=1e-6)
+
+
+def test_graph_cross_mi_torch_matches_numpy():
+    pytest.importorskip("torch")
+    rng = np.random.default_rng(1)
+    n, d = 80, 6
+    X = rng.poisson(1.5, size=(n, d)).astype(float)
+    adj = csr_matrix((rng.random((n, n)) < 0.2).astype(np.float64))
+    genes = [f"g{i}" for i in range(d)]
+    cpu = target_graph_cross_mi(X, genes, graph=adj, backend="numpy")
+    gpu = target_graph_cross_mi(X, genes, graph=adj, backend="gpu", device="cpu")
+    for g1 in genes:
+        for g2 in cpu[g1]:
+            assert cpu[g1][g2] == pytest.approx(gpu[g1][g2], abs=1e-6)
+
+
+def test_cross_mi_torch_chunking_consistent():
+    pytest.importorskip("torch")
+    from genevector._graph_targets import _cross_mi_matrix, _cross_mi_matrix_torch
+    from genevector.metrics import discretize_genes
+    rng = np.random.default_rng(2)
+    X = rng.poisson(1.0, size=(120, 8)).astype(float)
+    Ad, na = discretize_genes(X)
+    Bd, nb = discretize_genes(X + rng.poisson(0.5, X.shape))
+    ref = _cross_mi_matrix(Ad, na, Bd, nb)
+    full = _cross_mi_matrix_torch(Ad, na, Bd, nb, device="cpu", max_elems=10**9)
+    chunked = _cross_mi_matrix_torch(Ad, na, Bd, nb, device="cpu", max_elems=120 * 2)
+    np.testing.assert_allclose(ref, full, atol=1e-6)
+    np.testing.assert_allclose(ref, chunked, atol=1e-6)
+
+
+def test_graph_mi_rust_matches_numpy():
+    from genevector.metrics import HAS_RUST
+    if not HAS_RUST:
+        pytest.skip("rust extension (_rust) not built")
+    X, adj, genes = _make_chain_panel(n=60)
+    cpu = target_graph_mi(X, genes, graph=adj, backend="numpy")
+    rust = target_graph_mi(X, genes, graph=adj, backend="rust")
+    for g1 in genes:
+        for g2 in cpu[g1]:
+            assert cpu[g1][g2] == pytest.approx(rust[g1][g2], abs=1e-6)
diff --git a/tests/test_phenotyping.py b/tests/test_phenotyping.py
new file mode 100644
index 0000000..b275918
--- /dev/null
+++ b/tests/test_phenotyping.py
@@ -0,0 +1,147 @@
+"""Tests for the phenotyping improvements in CellEmbedding.
+
+Covers the dataset_vector fix, opt-in debias/contrastive scoring, spatial label
+propagation, count-adaptive cell-vector denoising, and the marker-dict QC.
+"""
+import numpy as np
+import pytest
+from scipy.sparse import csr_matrix
+
+anndata = pytest.importorskip("anndata")
+import pandas as pd
+
+from genevector.data import GeneVectorDataset
+from genevector.model import GeneVector
+from genevector.embedding import GeneEmbedding, CellEmbedding
+
+
+@pytest.fixture(scope="module")
+def trained(tmp_path_factory):
+    """Tiny 3-type dataset (2 markers each) trained end-to-end."""
+    rng = np.random.RandomState(0)
+    n_per = 40
+    types = ["A", "B", "C"]
+    blocks, labels, coords = [], [], []
+    for t, ct in enumerate(types):
+        base = np.full((n_per, 6), 0.5)
+        base[:, 2 * t:2 * t + 2] += 6.0  # this type's 2 markers high
+        blocks.append(rng.poisson(base))
+        labels += [ct] * n_per
+        coords.append(rng.rand(n_per, 2) + np.array([t * 5.0, 0.0]))  # spatially separated
+    X = csr_matrix(np.vstack(blocks).astype(np.float64))
+    genes = [f"G{i}" for i in range(6)]
+    adata = anndata.AnnData(X=X,
+                            var=pd.DataFrame(index=genes),
+                            obs=pd.DataFrame({"ct": labels},
+                                             index=[f"C{i}" for i in range(len(labels))]))
+    adata.layers["counts"] = adata.X.copy()
+    adata.obsm["spatial"] = np.vstack(coords)
+
+    vec = str(tmp_path_factory.mktemp("vec") / "emb.vec")
+    ds = GeneVectorDataset(adata.copy(), load_expression=True, signed_mi=True,
+                           device="cpu", use_cache=False, mi_backend="numpy")
+    gv = GeneVector(ds, output_file=vec, emb_dimension=20, c=100, gain=10,
+                    init_ortho=True, device="cpu")
+    gv.train(150, update_interval=75)
+    embed = GeneEmbedding(vec, ds, vector="average")
+    cembed = CellEmbedding(ds, embed)
+    agv = cembed.get_adata()
+    markers = {"A": ["G0", "G1"], "B": ["G2", "G3"], "C": ["G4", "G5"]}
+    return adata, embed, cembed, agv, markers
+
+
+def test_dataset_vector_is_nonzero(trained):
+    _, _, cembed, _, _ = trained
+    # previously initialised to zeros (bug); now the mean cell vector.
+    assert np.linalg.norm(cembed.dataset_vector) > 0
+
+
+def test_phenotype_probability_default(trained):
+    _, _, cembed, agv, markers = trained
+    out = cembed.phenotype_probability(agv, markers, temperature=0.05, target_col="gv")
+    assert "gv" in out.obs
+    # probability columns sum to ~1 per cell
+    pcols = [c for c in out.obs.columns if "Pseudo-probability" in c]
+    assert len(pcols) == 3
+    s = out.obs[pcols].to_numpy().sum(1)
+    npt = np.testing
+    npt.assert_allclose(s, 1.0, atol=1e-5)
+    # recovers the planted structure reasonably
+    from sklearn.metrics import adjusted_rand_score
+    assert adjusted_rand_score(out.obs["ct"], out.obs["gv"]) > 0.5
+
+
+def test_debias_and_contrastive_run(trained):
+    _, _, cembed, agv, markers = trained
+    o1 = cembed.phenotype_probability(agv, markers, temperature=0.05, target_col="d", debias=0.5)
+    o2 = cembed.phenotype_probability(agv, markers, temperature=0.05, target_col="c", contrastive=True)
+    assert set(o1.obs["d"]) <= set(markers)
+    assert set(o2.obs["c"]) <= set(markers)
+
+
+def test_score_norm_options_run(trained):
+    _, _, cembed, agv, markers = trained
+    for sn in ("zscore", "rank"):
+        out = cembed.phenotype_probability(agv, markers, temperature=0.05,
+                                           target_col=f"sn_{sn}", score_norm=sn)
+        assert set(out.obs[f"sn_{sn}"]) <= set(markers)
+        pcols = [c for c in out.obs.columns if "Pseudo-probability" in c]
+        np.testing.assert_allclose(out.obs[pcols].to_numpy().sum(1), 1.0, atol=1e-5)
+    with pytest.raises(ValueError, match="Unknown score_norm"):
+        cembed.phenotype_probability(agv, markers, target_col="bad", score_norm="bogus")
+
+
+def test_scoring_time_smoothing_runs(trained):
+    _, _, cembed, agv, markers = trained
+    from sklearn.neighbors import kneighbors_graph
+    G = kneighbors_graph(agv.obsm["spatial"], 5, mode="connectivity")
+    G = ((G + G.T) > 0).astype(float).tocsr()
+    before = np.array(cembed.matrix).copy()
+    out = cembed.phenotype_probability(agv, markers, temperature=0.05, target_col="sm",
+                                       smooth_graph=G, smooth_alpha=0.5)
+    assert set(out.obs["sm"]) <= set(markers)
+    # scoring-time smoothing must NOT mutate the stored cell matrix
+    np.testing.assert_allclose(np.array(cembed.matrix), before)
+
+
+def test_label_propagation_runs_and_preserves_simplex(trained):
+    _, _, cembed, agv, markers = trained
+    W = cembed._row_normalize_graph  # ensure helper exists
+    from sklearn.neighbors import kneighbors_graph
+    G = kneighbors_graph(agv.obsm["spatial"], 5, mode="connectivity")
+    G = ((G + G.T) > 0).astype(float).tocsr()
+    out = cembed.phenotype_probability(agv, markers, temperature=0.05, target_col="lp",
+                                       lp_graph=G, lp_alpha=0.4, lp_iter=3)
+    pcols = [c for c in out.obs.columns if "Pseudo-probability" in c]
+    np.testing.assert_allclose(out.obs[pcols].to_numpy().sum(1), 1.0, atol=1e-5)
+
+
+def test_denoise_cell_vectors(trained):
+    adata, embed, cembed, agv, markers = trained
+    # fresh cembed so we don't mutate the shared fixture's matrix
+    ds = GeneVectorDataset(adata.copy(), load_expression=True, device="cpu",
+                           use_cache=False, mi_backend="numpy")
+    ds.mi_scores = {g: {h: 0.0 for h in ds.data.genes if h != g} for g in ds.data.genes}
+    ce = CellEmbedding(ds, embed)
+    from sklearn.neighbors import kneighbors_graph
+    coords = adata.obsm["spatial"]
+    G = kneighbors_graph(coords, 5, mode="connectivity")
+    G = ((G + G.T) > 0).astype(float).tocsr()
+    before = np.array(ce.matrix).copy()
+    ce.denoise_cell_vectors(G, adaptive=True, counts=np.asarray(adata.layers["counts"].sum(1)).ravel())
+    after = np.array(ce.matrix)
+    assert after.shape == before.shape
+    assert not np.allclose(after, before)  # vectors actually changed
+    assert hasattr(ce, "uncorrected_matrix")
+
+
+def test_qc_marker_dict_flags(trained):
+    adata, _, cembed, _, markers = trained
+    bad = dict(markers)
+    bad["A"] = ["G0", "NOTAGENE"]      # absent marker
+    bad["B"] = ["G2", "G4"]            # G4 is really a C marker -> off_target
+    qc = cembed.qc_marker_dict(adata, bad, layer="counts", verbose=False)
+    assert {"phenotype", "marker", "in_panel", "specificity", "flag"} <= set(qc.columns)
+    assert (qc[qc.marker == "NOTAGENE"]["flag"] == "absent").all()
+    g4 = qc[(qc.phenotype == "B") & (qc.marker == "G4")]["flag"].iloc[0]
+    assert g4.startswith("off_target")