Implement mdl metric

examples/io_and_data/example_compat_conversion.py

@@ -10,7 +10,7 @@
 from py3plex.compat import convert, to_ir
 from py3plex.compat.schema import infer_schema
 from py3plex.io.schema import MultiLayerGraph, Node, Layer, Edge
-
+from py3plex.compat.exceptions import CompatibilityError
 
 def create_sample_graph():
     """Create a sample multilayer network."""
@@ -82,7 +82,6 @@ def main():
     print("\n3. Error Handling (Strict Mode)")
     print("-" * 60)
     try:
-        from py3plex.compat.exceptions import CompatibilityError
         matrix = convert(graph, "scipy_sparse", strict=True)
         print(f"    Conversion succeeded")
     except CompatibilityError as e:

py3plex/algorithms/community_detection/autocommunity_executor.py

@@ -455,6 +455,7 @@ def _evaluate_algorithms(
     from py3plex.algorithms.community_detection.multilayer_quality_metrics import (
         replica_consistency,
         layer_entropy,
+        mdl_score,
     )
 
     import inspect
@@ -519,9 +520,7 @@ def _evaluate_algorithms(
                         value = 0.0
 
                 elif metric_name == "mdl":
-                    # Description length (if available)
-                    # Placeholder for now
-                    value = 0.0
+                    value = mdl_score(partition, network)
 
                 elif metric_name == "replica_consistency":
                     # Multilayer coherence metric

py3plex/algorithms/community_detection/multilayer_quality_metrics.py

@@ -9,6 +9,8 @@
   layers are assigned to the same community (multilayer coherence)
 - layer_entropy: Measures the balance of community sizes within each layer,
   averaged across layers (degeneracy guardrail)
+- mdl_score: Minimum Description Length based on SBM two-part description
+  length (lower is better)
 """
 
 from __future__ import annotations
@@ -318,8 +320,237 @@ def layer_entropy(
         return 0.0
 
     mean_entropy = float(np.mean(layer_entropies))
-    
+
     # Apply clipping
     clipped_entropy = np.clip(mean_entropy, clip[0], clip[1])
-    
+
     return float(clipped_entropy)
+
+
+def _mdl_single_layer(
+    layer_partition: Dict[Any, Any],
+    layer_edges: List[tuple],
+    directed: bool = False,
+    include_model_cost: bool = True,
+) -> float:
+    """Compute SBM two-part description length for a flat edge set.
+
+    Despite the name, this is used for two things in `mdl_score`: the
+    per-layer intra-layer term (with `include_model_cost=True`, using a
+    single layer's local node -> community assignment), and the global
+    inter-layer term (with `include_model_cost=False`, using the full
+    across-all-layers node-layer -> community assignment, restricted to
+    edges that cross layers). See `mdl_score` for why model cost is only
+    charged once, per layer.
+
+    Args:
+        layer_partition: Dict mapping node -> community_id.
+            Community ids need only be hashable (e.g. singleton sentinel
+            objects for unassigned nodes are supported).
+        layer_edges: List of (u, v) edge pairs to account for. Both
+            endpoints must be keys of `layer_partition`.
+        directed: Whether the underlying graph is directed. Controls the
+            maximum-edge normalization for block pairs. Defaults to False
+            (undirected), preserving backward-compatible behavior for callers
+            that do not pass the flag.
+        include_model_cost: Whether to add the n * log2(k) model-cost term.
+            Set to False when `layer_partition`/`layer_edges` describe a
+            block structure whose assignment cost was already charged
+            elsewhere (e.g. the inter-layer term reuses assignments already
+            paid for by the per-layer model cost).
+
+    Returns:
+        Description length in bits.
+    """
+    communities: Dict[Any, List[Any]] = defaultdict(list)
+    for node, comm in layer_partition.items():
+        communities[comm].append(node)
+
+    k = len(communities)
+    n = len(layer_partition)
+
+    if k == 0 or n == 0:
+        return 0.0
+
+    # Count edges between/within each block pair. Community ids are only
+    # required to be hashable (not orderable) -- singleton sentinel ids used
+    # for unassigned nodes are plain objects, so pairs are keyed with a
+    # frozenset rather than a min/max tuple.
+    # Self-loops are skipped: they do not fit the n_r*(n_r-1) pair model
+    edge_counts: Dict[frozenset, int] = defaultdict(int)
+    for u, v in layer_edges:
+        if u == v:
+            continue
+        r = layer_partition.get(u)
+        s = layer_partition.get(v)
+        if r is None or s is None:
+            continue
+        edge_counts[frozenset((r, s))] += 1
+
+    # Part 1: model cost - n * log2(k) bits to encode node assignments
+    model_cost = n * np.log2(k) if (include_model_cost and k > 1) else 0.0
+
+    # Part 2: data cost - binary entropy per block pair. Block pairs with no
+    # observed edges contribute exactly 0 (H(p=0) == 0), so it suffices to
+    # iterate over the pairs that actually have edges rather than every pair
+    # of communities. This keeps the cost O(|E|) instead of O(k^2), which
+    # matters once unassigned nodes are represented as singleton communities
+    # (k can then be as large as n).
+    data_cost = 0.0
+    for pair, e_rs in edge_counts.items():
+        if len(pair) == 1:
+            (r,) = pair
+            s = r
+        else:
+            r, s = tuple(pair)
+        n_r = len(communities[r])
+        n_s = len(communities[s])
+
+        # Maximum possible edges between blocks r and s.
+        # Directed graphs allow both (u->v) and (v->u): the off-diagonal
+        # capacity doubles and the on-diagonal uses n_r*(n_r-1).
+        if r == s:
+            m_rs = n_r * (n_r - 1) if directed else n_r * (n_r - 1) / 2
+        else:
+            m_rs = 2 * n_r * n_s if directed else n_r * n_s
+        if m_rs == 0:
+            continue
+
+        # Clamp to [0, 1]. Parallel/directional edges can push p above 1,
+        # which makes log2(1 - p) return nan and silently corrupt the score.
+        p = min(e_rs / m_rs, 1.0)
+
+        # Binary entropy H(p); 0 when p == 0 or p == 1
+        if 0.0 < p < 1.0:
+            data_cost += m_rs * (-p * np.log2(p) - (1.0 - p) * np.log2(1.0 - p))
+
+    return model_cost + data_cost
+
+
+def mdl_score(
+    partition: Dict[Any, int],
+    network: Any,
+) -> float:
+    """Compute MDL (Minimum Description Length) score for a multilayer partition.
+
+    Computes the SBM two-part description length independently per layer and
+    sums the results, then adds one more term for inter-layer (cross-layer)
+    edges. Operating per layer for the intra-layer term avoids the flattening
+    problem where pooling all layers into core_network inflates block sizes
+    with cross-layer node-layer pairs that share no edges, corrupting the
+    edge density estimates that drive the data cost.
+
+    Per-layer computation (intra-layer edges only):
+      - Part 1 (model cost): n_ℓ * log2(k_ℓ) bits to encode the community
+        assignment of each node in layer ℓ, where n_ℓ is the number of
+        node-layer pairs in that layer and k_ℓ is the number of distinct
+        communities present in that layer.
+      - Part 2 (data cost): for every block pair (r, s) within layer ℓ,
+        the binary entropy of the observed intra-layer edge density times the
+        maximum possible intra-layer edges between those blocks.
+
+    Inter-layer computation (edges whose endpoints are in different layers,
+    e.g. multiplex coupling edges or general cross-layer edges):
+      - No additional model cost is charged -- every node-layer pair's
+        community assignment was already paid for by the per-layer model
+        cost above.
+      - Data cost only: block sizes N_r/N_s are the *global* community
+        sizes (summed across all layers), and the block pair (r, s) is
+        scored using the observed inter-layer edge density between them,
+        the same binary-entropy formula as the intra-layer case. This
+        slightly overstates the space of possible inter-layer edges (it
+        does not subtract same-layer node pairs from N_r * N_s, since
+        inter-layer edges in py3plex are not restricted to same-node
+        couplings and can connect arbitrary node-layer pairs), which is a
+        conservative approximation rather than an exact block model.
+
+    Total MDL = Σ_ℓ (model_cost_ℓ + data_cost_ℓ) + data_cost_inter_layer
+
+    Lower score is better. Ignoring inter-layer edges entirely would let two
+    partitions with identical intra-layer structure but very different (e.g.
+    incoherent vs. replica-consistent) cross-layer coupling score identically,
+    which is misleading for a *multilayer* MDL metric.
+
+    Partial partitions: every node present in the network is accounted for,
+    not just the ones covered by `partition`. A node the partition omits is
+    scored as its own singleton community, so both n_ℓ (model cost) and its
+    incident edges (data cost) still count. Without this, an algorithm could
+    lower its MDL score simply by leaving hard nodes and their edges out of
+    the partition.
+
+    Args:
+        partition: Dict mapping (node_id, layer) -> community_id.
+        network: Multilayer network with a core_network NetworkX graph.
+
+    Returns:
+        Total description length in bits (float).  Returns 0.0 for empty
+        inputs.
+
+    Examples:
+        >>> # Perfect 2-community split on a clique pair
+        >>> partition = {('A', 'L'): 0, ('B', 'L'): 0, ('C', 'L'): 1, ('D', 'L'): 1}
+        >>> score = mdl_score(partition, net)
+        >>> assert score >= 0.0
+    """
+    G = getattr(network, "core_network", None)
+    node_set = set(G.nodes()) if G is not None else set()
+    all_nodes = node_set | set(partition.keys())
+    if not all_nodes:
+        return 0.0
+
+    directed = G.is_directed() if G is not None else False
+
+    # Collect layers and per-layer node assignments. Every node in the
+    # network is included -- not just the ones covered by `partition` -- so
+    # that an algorithm returning a partial partition cannot shrink n_ℓ or
+    # skip edges just by omitting nodes. Nodes missing from `partition` get a
+    # unique sentinel object as their "community", i.e. each is treated as
+    # its own singleton community.
+    layer_partitions: Dict[Any, Dict[Any, Any]] = defaultdict(dict)
+    n_unassigned = 0
+    for node in all_nodes:
+        layer = node[1] if isinstance(node, tuple) and len(node) >= 2 else None
+        if node in partition:
+            layer_partitions[layer][node] = partition[node]
+        else:
+            layer_partitions[layer][node] = object()
+            n_unassigned += 1
+
+    if n_unassigned:
+        warnings.warn(
+            f"Partition covers {len(partition)} of {len(all_nodes)} nodes; "
+            f"{n_unassigned} unassigned node(s) are scored as singleton "
+            "communities so they cannot be omitted to lower the MDL score.",
+            stacklevel=2
+        )
+
+    # Split edges into intra-layer (same layer on both ends) and inter-layer
+    # (endpoints in different layers, e.g. multiplex coupling edges).
+    layer_edges: Dict[Any, List[tuple]] = defaultdict(list)
+    inter_layer_edges: List[tuple] = []
+    for u, v in (G.edges() if G is not None else []):
+        u_layer = u[1] if isinstance(u, tuple) and len(u) >= 2 else None
+        v_layer = v[1] if isinstance(v, tuple) and len(v) >= 2 else None
+        if u_layer == v_layer:
+            layer_edges[u_layer].append((u, v))
+        else:
+            inter_layer_edges.append((u, v))
+
+    # Sum description length across layers (intra-layer model + data cost)
+    total = 0.0
+    global_partition: Dict[Any, Any] = {}
+    for layer, lpartition in layer_partitions.items():
+        total += _mdl_single_layer(lpartition, layer_edges[layer], directed=directed)
+        global_partition.update(lpartition)
+
+    # Add the inter-layer data cost: community assignments were already paid
+    # for above, so only the observed cross-layer edge density is charged.
+    if inter_layer_edges:
+        total += _mdl_single_layer(
+            global_partition,
+            inter_layer_edges,
+            directed=directed,
+            include_model_cost=False,
+        )
+
+    return total