diff --git a/scripts/test.py b/scripts/test.py
index e014e5f..9ee5a90 100644
--- a/scripts/test.py
+++ b/scripts/test.py
@@ -41,7 +41,11 @@ def main() -> None:
         coupling,
         num_nodes,
     )
-    degree = versatility.get_multi_degree(supra, num_layers, num_nodes)
+    degree = versatility.get_multi_degree(
+        supra,
+        nodes=num_nodes,
+        layers=num_layers,
+    )
 
     print(f"nodes={num_nodes}, layers={num_layers}")
     print(f"aggregated degree={degree.tolist()}")
diff --git a/src/MuxVizPy/global_descriptors.py b/src/MuxVizPy/global_descriptors.py
index f73583a..edad792 100644
--- a/src/MuxVizPy/global_descriptors.py
+++ b/src/MuxVizPy/global_descriptors.py
@@ -4,7 +4,7 @@
 import scipy.sparse as sp
 
 
-def compute_average_global_clustering_coefficient(adj: sp.csr_matrix, n:int, l: int, logger: logging.Logger | None = None) -> float:
+def compute_average_global_clustering_coefficient(adj: sp.csr_matrix, *, nodes:int, layers: int, logger: logging.Logger | None = None) -> float:
     """
     Compute the average global clustering coefficient for a multilayer network.
     This implementation matches the R version from the reference, using the formula:
@@ -12,9 +12,9 @@ def compute_average_global_clustering_coefficient(adj: sp.csr_matrix, n:int, l:
 
     Where F is the matrix with 1s everywhere except the diagonal.
     Args:
-        adj: scipy.sparse.csr_matrix - the adjacency matrix of the multilayer network (shape: (n*l, n*l))
-        n: int - number of nodes per layer
-        l: int - number of layers
+        adj: scipy.sparse.csr_matrix - the adjacency matrix of the multilayer network (shape: (nodes*layers, nodes*layers))
+        nodes: int - number of nodes per layer
+        layers: int - number of layers
         logger: Optional[logging.Logger] - logger for debugging (default: None)
 
     Returns:
@@ -61,8 +61,8 @@ def _sparse_pmin(A: sp.csr_matrix, B: sp.csr_matrix) -> sp.csr_matrix:
 
 def compute_average_global_overlap(
     adj: sp.csr_matrix,
-    n: int,
-    l: int,
+    *, nodes: int,
+    layers: int,
     weighted: bool = False,
     logger: logging.Logger | None = None,
 ) -> float:
@@ -79,9 +79,9 @@ def compute_average_global_overlap(
     For undirected networks (symmetric supra-adjacency) the result is halved.
 
     Args:
-        adj: scipy.sparse.csr_matrix - supra-adjacency matrix (n*l, n*l)
-        n: int - number of nodes per layer
-        l: int - number of layers
+        adj: scipy.sparse.csr_matrix - supra-adjacency matrix (nodes*layers, nodes*layers)
+        nodes: int - number of nodes per layer
+        layers: int - number of layers
         weighted: bool - if True use actual edge weights; if False binarize first (default: False)
         logger: Optional[logging.Logger] - logger for debugging (default: None)
 
@@ -92,24 +92,24 @@ def compute_average_global_overlap(
         De Domenico et al. (2015) "Structural reducibility of multilayer networks",
         Nature Communications, 6, 6864.
     """
-    if l < 2:
+    if layers < 2:
         raise ValueError("At least two layers are required.")
 
-    layers = []
-    for alpha in range(l):
-        A = adj[alpha * n:(alpha + 1) * n, alpha * n:(alpha + 1) * n].astype(float)
+    layer_matrices = []
+    for alpha in range(layers):
+        A = adj[alpha * nodes:(alpha + 1) * nodes, alpha * nodes:(alpha + 1) * nodes].astype(float)
         if not weighted:
             A = (A > 0).astype(float)
-        layers.append(A)
+        layer_matrices.append(A)
 
     # Running element-wise minimum: nonzero only where all layers share the edge.
-    overlap = layers[0]
-    for layer in layers[1:]:
+    overlap = layer_matrices[0]
+    for layer in layer_matrices[1:]:
         overlap = _sparse_pmin(overlap, layer)
-    norm_total = sum(float(layer.sum()) for layer in layers)
+    norm_total = sum(float(layer.sum()) for layer in layer_matrices)
 
     sum_overlap = float(overlap.sum())
-    avg = l * sum_overlap / norm_total if norm_total != 0 else 0.0
+    avg = layers * sum_overlap / norm_total if norm_total != 0 else 0.0
 
     # R checks `sum(A - A^T) == 0` before halving.  However sum(A) == sum(A^T)
     # for any matrix (transposing preserves the total), so R's condition is
@@ -117,15 +117,15 @@ def compute_average_global_overlap(
     avg /= 2
 
     if logger and logger.isEnabledFor(logging.DEBUG):
-        logger.debug(f"sum(O): {sum_overlap}, NormTotal: {norm_total}, L: {l}")
+        logger.debug(f"sum(O): {sum_overlap}, NormTotal: {norm_total}, L: {layers}")
 
     return avg
 
 
 def compute_average_global_overlap_matrix(
     adj: sp.csr_matrix,
-    n: int,
-    l: int,
+    *, nodes: int,
+    layers: int,
     weighted: bool = False,
     logger: logging.Logger | None = None,
 ) -> np.ndarray:
@@ -138,35 +138,35 @@ def compute_average_global_overlap_matrix(
 
     This is a Dice-like similarity coefficient in [0, 1].  The diagonal is 1.
     Args:
-        adj: scipy.sparse.csr_matrix - supra-adjacency matrix (n*l, n*l)
-        n: int - number of nodes per layer
-        l: int - number of layers
+        adj: scipy.sparse.csr_matrix - supra-adjacency matrix (nodes*layers, nodes*layers)
+        nodes: int - number of nodes per layer
+        layers: int - number of layers
         weighted: bool - if True use actual edge weights; if False binarize first (default: False)
         logger: Optional[logging.Logger] - logger for debugging (default: None)
 
     Returns:
-        numpy.ndarray of shape (l, l) - symmetric overlap matrix with 1s on diagonal
+        numpy.ndarray of shape (layers, layers) - symmetric overlap matrix with 1s on diagonal
 
     Reference:
         De Domenico et al. (2015) "Structural reducibility of multilayer networks",
         Nature Communications, 6, 6864.
     """
-    if l < 2:
+    if layers < 2:
         raise ValueError("At least two layers are required.")
 
-    layers = []
+    layer_matrices = []
     layer_sums = []
-    for alpha in range(l):
-        A = adj[alpha * n:(alpha + 1) * n, alpha * n:(alpha + 1) * n].astype(float)
+    for alpha in range(layers):
+        A = adj[alpha * nodes:(alpha + 1) * nodes, alpha * nodes:(alpha + 1) * nodes].astype(float)
         if not weighted:
             A = (A > 0).astype(float)
-        layers.append(A)
+        layer_matrices.append(A)
         layer_sums.append(float(A.sum()))
 
-    M = np.eye(l, dtype=float)
-    for l1 in range(l - 1):
-        for l2 in range(l1 + 1, l):
-            O = _sparse_pmin(layers[l1], layers[l2])
+    M = np.eye(layers, dtype=float)
+    for l1 in range(layers - 1):
+        for l2 in range(l1 + 1, layers):
+            O = _sparse_pmin(layer_matrices[l1], layer_matrices[l2])
             denom = layer_sums[l1] + layer_sums[l2]
             val = 2.0 * float(O.sum()) / denom if denom != 0 else 0.0
             M[l1, l2] = val
@@ -180,8 +180,8 @@ def compute_average_global_overlap_matrix(
 
 def compute_average_global_node_overlap_matrix(
     adj: sp.csr_matrix,
-    n: int,
-    l: int,
+    *, nodes: int,
+    layers: int,
     logger: logging.Logger | None = None,
 ) -> np.ndarray:
     """
@@ -195,32 +195,32 @@ def compute_average_global_node_overlap_matrix(
     The diagonal is 1 (all active nodes overlap with themselves).
 
     Args:
-        adj: scipy.sparse.csr_matrix - supra-adjacency matrix (n*l, n*l)
-        n: int - number of nodes per layer
-        l: int - number of layers
+        adj: scipy.sparse.csr_matrix - supra-adjacency matrix (nodes*layers, nodes*layers)
+        nodes: int - number of nodes per layer
+        layers: int - number of layers
         logger: Optional[logging.Logger] - logger for debugging (default: None)
 
     Returns:
-        numpy.ndarray of shape (l, l) - symmetric overlap matrix with 1s on diagonal
+        numpy.ndarray of shape (layers, layers) - symmetric overlap matrix with 1s on diagonal
 
     Reference:
         De Domenico et al. (2015) "Structural reducibility of multilayer networks",
         Nature Communications, 6, 6864.
     """
-    if l < 2:
+    if layers < 2:
         raise ValueError("At least two layers are required.")
 
     active = []
-    for alpha in range(l):
-        A = adj[alpha * n:(alpha + 1) * n, alpha * n:(alpha + 1) * n]
+    for alpha in range(layers):
+        A = adj[alpha * nodes:(alpha + 1) * nodes, alpha * nodes:(alpha + 1) * nodes]
         out_deg = np.asarray(A.sum(axis=1)).ravel()   # row sums
         in_deg  = np.asarray(A.sum(axis=0)).ravel()   # col sums
         active.append(set(np.where((out_deg > 0) | (in_deg > 0))[0]))
 
-    M = np.eye(l, dtype=float)
-    for l1 in range(l - 1):
-        for l2 in range(l1 + 1, l):
-            val = len(active[l1] & active[l2]) / n
+    M = np.eye(layers, dtype=float)
+    for l1 in range(layers - 1):
+        for l2 in range(l1 + 1, layers):
+            val = len(active[l1] & active[l2]) / nodes
             M[l1, l2] = val
             M[l2, l1] = val
 
diff --git a/src/MuxVizPy/mesoscale.py b/src/MuxVizPy/mesoscale.py
index c61b883..5cfff14 100644
--- a/src/MuxVizPy/mesoscale.py
+++ b/src/MuxVizPy/mesoscale.py
@@ -95,8 +95,8 @@ def inter_layer_assortativity(
 
 def compute_local_clustering_coefficient(
     adj: sp.sparse.spmatrix,
-    n: int,
-    l: int,
+    *, nodes: int,
+    layers: int,
     logger: logging.Logger | None = None,
 ) -> np.ndarray:
     """
@@ -110,17 +110,17 @@ def compute_local_clustering_coefficient(
     ----------
     adj : scipy.sparse matrix, shape (N·L, N·L)
         Binary supra-adjacency matrix with supra-node indexing
-        ``row = layer * n + node``.
-    n : int
+        ``row = layer * nodes + node``.
+    nodes : int
         Number of physical nodes.
-    l : int
+    layers : int
         Number of layers.
     logger : optional
         Unused; kept for API consistency with other descriptor functions.
 
     Returns
     -------
-    numpy.ndarray, shape (n,)
+    numpy.ndarray, shape (nodes,)
         Local clustering coefficient in ``[0, 1]`` for each physical node.
         Nodes with no possible triangles (isolated or degree-1 nodes) are
         assigned 0.
@@ -131,13 +131,13 @@ def compute_local_clustering_coefficient(
     networks. *Physical Review X*, 3(4), 041022.
     https://doi.org/10.1103/PhysRevX.3.041022
     """
-    NL = n * l
+    NL = nodes * layers
     A = adj.astype(float)
 
-    # Fold matrix P of shape (NL x n): P[i*n+p, p] = 1 for all layers i.
-    # P.T @ M @ P sums all (l x l) blocks of M into a single (n x n) matrix.
+    # Fold matrix P of shape (NL x nodes): P[i*nodes+p, p] = 1 for all layers i.
+    # P.T @ M @ P sums all (layers x layers) blocks of M into a single (nodes x nodes) matrix.
     idx = np.arange(NL)
-    P = csr_matrix((np.ones(NL), (idx, idx % n)), shape=(NL, n))
+    P = csr_matrix((np.ones(NL), (idx, idx % nodes)), shape=(NL, nodes))
 
     A2 = A @ A
     A3 = A2 @ A
diff --git a/src/MuxVizPy/percolation.py b/src/MuxVizPy/percolation.py
index 7d8ca6a..1a9a89d 100644
--- a/src/MuxVizPy/percolation.py
+++ b/src/MuxVizPy/percolation.py
@@ -10,8 +10,9 @@
 
 def get_percolation(
     g_list: list[gt.Graph],
-    layers: int,
+    *,
     nodes: int,
+    layers: int,
     order: np.ndarray,
 ) -> dict[str, np.ndarray | float]:
     """
@@ -21,10 +22,10 @@ def get_percolation(
     ----------
     g_list : list of graph_tool.Graph
         List of graphs, one per layer of the multilayer network.
-    layers : int
-        Number of layers in the network.
     nodes : int
         Number of physical nodes (used for normalization).
+    layers : int
+        Number of layers in the network.
     order : np.ndarray
         Array specifying the order in which nodes are removed during percolation.
 
diff --git a/src/MuxVizPy/utils/parsing.py b/src/MuxVizPy/utils/parsing.py
index 0bfe110..d0714fe 100644
--- a/src/MuxVizPy/utils/parsing.py
+++ b/src/MuxVizPy/utils/parsing.py
@@ -820,10 +820,9 @@ def build_tensor_from_supra_adjacency_matrix(a, *, nodes, layers):
 
 def build_transition_matrix_from_adjacency_matrix(
         adj: sp.csr_matrix,
-        n: int,
-        l: int,
+        *, nodes: int,
+        layers: int,
         kind: str,
-        *,
         alpha: float | None = None,
         logger: logging.Logger | None = None,
 ) -> sp.csr_matrix:
@@ -838,11 +837,11 @@ def build_transition_matrix_from_adjacency_matrix(
     Notes:
         - Teleportation term makes the matrix dense; returned as CSR for consistency.
     """
-    NL = n * l
+    NL = nodes * layers
     if not sp.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     if adj.shape != (NL, NL):
-        raise ValueError(f"Incompatible shape {adj.shape} for n={n}, l={l} (expected {(NL, NL)})")
+        raise ValueError(f"Incompatible shape {adj.shape} for n={nodes}, l={layers} (expected {(NL, NL)})")
 
     kind = kind.lower().strip()
 
@@ -982,7 +981,7 @@ def get_aggregate_network(
 
 
 def supra_adjacency_to_block_tensor(
-    supra: sp.spmatrix, layers: int, nodes: int
+    supra: sp.spmatrix, *, nodes: int, layers: int
 ) -> list[list[sp.spmatrix]]:
     """
     Convert supra-adjacency matrix to block tensor format.
@@ -991,10 +990,10 @@ def supra_adjacency_to_block_tensor(
     ----------
     supra : scipy.sparse matrix
         Supra-adjacency matrix.
-    layers : int
-        Number of layers.
     nodes : int
         Number of nodes.
+    layers : int
+        Number of layers.
 
     Returns
     -------
@@ -1009,7 +1008,7 @@ def supra_adjacency_to_block_tensor(
 
 
 def node_tensor_to_network_list(
-    tensor: list[sp.spmatrix], layers: int, nodes: int
+    tensor: list[sp.spmatrix], *, nodes: int, layers: int
 ) -> list[gt.Graph]:
     """
     Convert node-layer adjacency matrices into graph_tool graphs.
@@ -1018,10 +1017,10 @@ def node_tensor_to_network_list(
     ----------
     tensor : list of scipy.sparse matrices
         List of per-layer adjacency matrices.
-    layers : int
-        Number of layers.
     nodes : int
         Number of nodes.
+    layers : int
+        Number of layers.
 
     Returns
     -------
@@ -1078,6 +1077,7 @@ def supra_idx(node, layer):
 def create_supra_transition_matrix_virus(
     supra: sp.spmatrix,
     node_tensor: list[sp.spmatrix],
+    *,
     nodes: int,
     layers: int,
     p_intra: float = 1,
diff --git a/src/MuxVizPy/versatility.py b/src/MuxVizPy/versatility.py
index bdcf8a5..f9c4100 100644
--- a/src/MuxVizPy/versatility.py
+++ b/src/MuxVizPy/versatility.py
@@ -30,13 +30,13 @@ def get_largest_eigenvalue(
 # Block accumulation and aggregation helpers (integrated from hornet/node_based)
 # ---------------------------------------------------------------------------
 
-def is_in_diagonal_block(i: int, j: int, n: int, l: int) -> bool:
+def is_in_diagonal_block(i: int, j: int, *, nodes: int, layers: int) -> bool:
     """Check if nodes i and j are in the same diagonal block (same layer)."""
-    return (i // n) == (j // n)
+    return (i // nodes) == (j // nodes)
 
 
 def _accumulate_on_diagonal_blocks(
-    adj: sps.csr_matrix, n: int, l: int,
+    adj: sps.csr_matrix, *, nodes: int, layers: int,
     is_out_of_diagonal: bool, weights: np.ndarray | None = None,
 ) -> np.ndarray:
     """
@@ -45,10 +45,10 @@ def _accumulate_on_diagonal_blocks(
     Parameters
     ----------
     adj : scipy.sparse.csr_matrix
-        Supra-adjacency matrix of shape (n*l, n*l).
-    n : int
+        Supra-adjacency matrix of shape (nodes*layers, nodes*layers).
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     is_out_of_diagonal : bool
         If True, accumulate off-diagonal (inter-layer) blocks.
@@ -58,22 +58,22 @@ def _accumulate_on_diagonal_blocks(
     Returns
     -------
     np.ndarray
-        Shape (n, l) array of accumulated values.
+        Shape (nodes, layers) array of accumulated values.
     """
-    NL = n * l
+    NL = nodes * layers
     if adj.shape != (NL, NL):
         raise ValueError(f"Adjacency matrix shape {adj.shape} does not match expected shape {(NL, NL)}")
 
     coo = adj.tocoo(copy=False)
-    row_layers = coo.row // n
-    col_layers = coo.col // n
+    row_layers = coo.row // nodes
+    col_layers = coo.col // nodes
     diagonal_mask = row_layers == col_layers
     if is_out_of_diagonal:
         diagonal_mask = ~diagonal_mask
 
-    tgt_nodes = coo.col[diagonal_mask] % n
+    tgt_nodes = coo.col[diagonal_mask] % nodes
     tgt_layers = col_layers[diagonal_mask]
-    flat_idx = tgt_layers * n + tgt_nodes
+    flat_idx = tgt_layers * nodes + tgt_nodes
 
     if weights is None:
         w = coo.data[diagonal_mask]
@@ -81,7 +81,7 @@ def _accumulate_on_diagonal_blocks(
         w = weights[diagonal_mask]
 
     accum = np.bincount(flat_idx, weights=w, minlength=NL)
-    return accum.reshape(l, n).T  # shape (n, l)
+    return accum.reshape(layers, nodes).T  # shape (nodes, layers)
 
 
 def aggregate_metrics_over_layers(metrics: np.ndarray, method: str = "sum") -> np.ndarray:
@@ -116,188 +116,188 @@ def aggregate_metrics_over_layers(metrics: np.ndarray, method: str = "sum") -> n
 # Per-layer degree / strength metrics (integrated from hornet/node_based)
 # ---------------------------------------------------------------------------
 
-def compute_indegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute indegree of each node stratified by layer. Returns shape (n, l)."""
+def compute_indegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute indegree of each node stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = np.ones_like(adj.data, dtype=np.float64)
-    indegree = _accumulate_on_diagonal_blocks(adj, n, l, is_out_of_diagonal=False, weights=weights)
+    indegree = _accumulate_on_diagonal_blocks(adj, nodes=nodes, layers=layers, is_out_of_diagonal=False, weights=weights)
     if logger and logger.isEnabledFor(logging.DEBUG):
         logger.debug(f"Indegree shape: {indegree.shape}, mean: {indegree.mean()}, max: {indegree.max()}, min: {indegree.min()}")
     return indegree
 
 
-def compute_aggregated_indegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated indegree over layers. Returns shape (n,)."""
-    indegree = compute_indegree(adj, n, l, logger=logger)
+def compute_aggregated_indegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated indegree over layers. Returns shape (nodes,)."""
+    indegree = compute_indegree(adj, nodes=nodes, layers=layers, logger=logger)
     return aggregate_metrics_over_layers(indegree, method=method)
 
 
-def compute_instrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute instrength of each node stratified by layer. Returns shape (n, l)."""
+def compute_instrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute instrength of each node stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = adj.data.astype(np.float64)
-    instrength = _accumulate_on_diagonal_blocks(adj, n, l, is_out_of_diagonal=False, weights=weights)
+    instrength = _accumulate_on_diagonal_blocks(adj, nodes=nodes, layers=layers, is_out_of_diagonal=False, weights=weights)
     if logger and logger.isEnabledFor(logging.DEBUG):
         logger.debug(f"Instrength shape: {instrength.shape}, mean: {instrength.mean()}, max: {instrength.max()}, min: {instrength.min()}")
     return instrength
 
 
-def compute_aggregated_instrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", *, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated instrength over layers. Returns shape (n,)."""
-    instrength = compute_instrength(adj, n, l, logger=logger)
+def compute_aggregated_instrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated instrength over layers. Returns shape (nodes,)."""
+    instrength = compute_instrength(adj, nodes=nodes, layers=layers, logger=logger)
     return aggregate_metrics_over_layers(instrength, method=method)
 
 
-def compute_outdegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute outdegree of each node stratified by layer. Returns shape (n, l)."""
+def compute_outdegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute outdegree of each node stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = np.ones_like(adj.data, dtype=np.float64)
-    outdegree = _accumulate_on_diagonal_blocks(adj.T, n, l, is_out_of_diagonal=False, weights=weights)
+    outdegree = _accumulate_on_diagonal_blocks(adj.T, nodes=nodes, layers=layers, is_out_of_diagonal=False, weights=weights)
     if logger and logger.isEnabledFor(logging.DEBUG):
         logger.debug(f"Outdegree shape: {outdegree.shape}, mean: {outdegree.mean()}, max: {outdegree.max()}, min: {outdegree.min()}")
     return outdegree
 
 
-def compute_aggregated_outdegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", *, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated outdegree over layers. Returns shape (n,)."""
-    outdegree = compute_outdegree(adj, n, l, logger=logger)
+def compute_aggregated_outdegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated outdegree over layers. Returns shape (nodes,)."""
+    outdegree = compute_outdegree(adj, nodes=nodes, layers=layers, logger=logger)
     return aggregate_metrics_over_layers(outdegree, method=method)
 
 
-def compute_outstrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute outstrength of each node stratified by layer. Returns shape (n, l)."""
+def compute_outstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute outstrength of each node stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = adj.data.astype(np.float64)
-    outstrength = _accumulate_on_diagonal_blocks(adj.T, n, l, is_out_of_diagonal=False, weights=weights)
+    outstrength = _accumulate_on_diagonal_blocks(adj.T, nodes=nodes, layers=layers, is_out_of_diagonal=False, weights=weights)
     if logger and logger.isEnabledFor(logging.DEBUG):
         logger.debug(f"Outstrength shape: {outstrength.shape}, mean: {outstrength.mean()}, max: {outstrength.max()}, min: {outstrength.min()}")
     return outstrength
 
 
-def compute_aggregated_outstrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated outstrength over layers. Returns shape (n,)."""
-    outstrength = compute_outstrength(adj, n, l, logger=logger)
+def compute_aggregated_outstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated outstrength over layers. Returns shape (nodes,)."""
+    outstrength = compute_outstrength(adj, nodes=nodes, layers=layers, logger=logger)
     return aggregate_metrics_over_layers(outstrength, method=method)
 
 
-def compute_multiindegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute multi-indegree (inter-layer blocks) stratified by layer. Returns shape (n, l)."""
+def compute_multiindegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute multi-indegree (inter-layer blocks) stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = np.ones_like(adj.data, dtype=np.float64)
-    return _accumulate_on_diagonal_blocks(adj, n, l, is_out_of_diagonal=True, weights=weights)
+    return _accumulate_on_diagonal_blocks(adj, nodes=nodes, layers=layers, is_out_of_diagonal=True, weights=weights)
 
 
-def compute_aggregated_multiindegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated multi-indegree over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_multiindegree(adj, n, l, logger=logger), method=method)
+def compute_aggregated_multiindegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated multi-indegree over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_multiindegree(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_multiinstrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute multi-instrength (inter-layer blocks) stratified by layer. Returns shape (n, l)."""
+def compute_multiinstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute multi-instrength (inter-layer blocks) stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = adj.data.astype(np.float64)
-    return _accumulate_on_diagonal_blocks(adj, n, l, is_out_of_diagonal=True, weights=weights)
+    return _accumulate_on_diagonal_blocks(adj, nodes=nodes, layers=layers, is_out_of_diagonal=True, weights=weights)
 
 
-def compute_aggregated_multiinstrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", *, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated multi-instrength over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_multiinstrength(adj, n, l, logger=logger), method=method)
+def compute_aggregated_multiinstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated multi-instrength over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_multiinstrength(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_multioutdegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute multi-outdegree (inter-layer blocks) stratified by layer. Returns shape (n, l)."""
+def compute_multioutdegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute multi-outdegree (inter-layer blocks) stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = np.ones_like(adj.data, dtype=np.float64)
-    return _accumulate_on_diagonal_blocks(adj.T, n, l, is_out_of_diagonal=True, weights=weights)
+    return _accumulate_on_diagonal_blocks(adj.T, nodes=nodes, layers=layers, is_out_of_diagonal=True, weights=weights)
 
 
-def compute_aggregated_multioutdegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated multi-outdegree over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_multioutdegree(adj, n, l, logger=logger), method=method)
+def compute_aggregated_multioutdegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated multi-outdegree over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_multioutdegree(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_multioutstrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute multi-outstrength (inter-layer blocks) stratified by layer. Returns shape (n, l)."""
+def compute_multioutstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute multi-outstrength (inter-layer blocks) stratified by layer. Returns shape (nodes, layers)."""
     if not sps.isspmatrix_csr(adj):
         adj = adj.tocsr(copy=False)
     weights = adj.data.astype(np.float64)
-    return _accumulate_on_diagonal_blocks(adj.T, n, l, is_out_of_diagonal=True, weights=weights)
+    return _accumulate_on_diagonal_blocks(adj.T, nodes=nodes, layers=layers, is_out_of_diagonal=True, weights=weights)
 
 
-def compute_aggregated_multioutstrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated multi-outstrength over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_multioutstrength(adj, n, l, logger=logger), method=method)
+def compute_aggregated_multioutstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated multi-outstrength over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_multioutstrength(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_total_indegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute total indegree (intra + inter layer) stratified by layer. Returns shape (n, l)."""
-    return compute_indegree(adj, n, l, logger=logger) + compute_multiindegree(adj, n, l, logger=logger)
+def compute_total_indegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute total indegree (intra + inter layer) stratified by layer. Returns shape (nodes, layers)."""
+    return compute_indegree(adj, nodes=nodes, layers=layers, logger=logger) + compute_multiindegree(adj, nodes=nodes, layers=layers, logger=logger)
 
 
-def compute_aggregated_total_indegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated total indegree over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_total_indegree(adj, n, l, logger=logger), method=method)
+def compute_aggregated_total_indegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated total indegree over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_total_indegree(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_total_instrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute total instrength (intra + inter layer) stratified by layer. Returns shape (n, l)."""
-    return compute_instrength(adj, n, l, logger=logger) + compute_multiinstrength(adj, n, l, logger=logger)
+def compute_total_instrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute total instrength (intra + inter layer) stratified by layer. Returns shape (nodes, layers)."""
+    return compute_instrength(adj, nodes=nodes, layers=layers, logger=logger) + compute_multiinstrength(adj, nodes=nodes, layers=layers, logger=logger)
 
 
-def compute_aggregated_total_instrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated total instrength over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_total_instrength(adj, n, l, logger=logger), method=method)
+def compute_aggregated_total_instrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated total instrength over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_total_instrength(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_total_outdegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute total outdegree (intra + inter layer) stratified by layer. Returns shape (n, l)."""
-    return compute_outdegree(adj, n, l, logger=logger) + compute_multioutdegree(adj, n, l, logger=logger)
+def compute_total_outdegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute total outdegree (intra + inter layer) stratified by layer. Returns shape (nodes, layers)."""
+    return compute_outdegree(adj, nodes=nodes, layers=layers, logger=logger) + compute_multioutdegree(adj, nodes=nodes, layers=layers, logger=logger)
 
 
-def compute_aggregated_total_outdegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated total outdegree over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_total_outdegree(adj, n, l, logger=logger), method=method)
+def compute_aggregated_total_outdegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated total outdegree over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_total_outdegree(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_total_outstrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute total outstrength (intra + inter layer) stratified by layer. Returns shape (n, l)."""
-    return compute_outstrength(adj, n, l, logger=logger) + compute_multioutstrength(adj, n, l, logger=logger)
+def compute_total_outstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute total outstrength (intra + inter layer) stratified by layer. Returns shape (nodes, layers)."""
+    return compute_outstrength(adj, nodes=nodes, layers=layers, logger=logger) + compute_multioutstrength(adj, nodes=nodes, layers=layers, logger=logger)
 
 
-def compute_aggregated_total_outstrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated total outstrength over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_total_outstrength(adj, n, l, logger=logger), method=method)
+def compute_aggregated_total_outstrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated total outstrength over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_total_outstrength(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_multidegree(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute multidegree (total indegree + total outdegree) stratified by layer. Returns shape (n, l)."""
-    return compute_total_indegree(adj, n, l, logger=logger) + compute_total_outdegree(adj, n, l, logger=logger)
+def compute_multidegree(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute multidegree (total indegree + total outdegree) stratified by layer. Returns shape (nodes, layers)."""
+    return compute_total_indegree(adj, nodes=nodes, layers=layers, logger=logger) + compute_total_outdegree(adj, nodes=nodes, layers=layers, logger=logger)
 
 
-def compute_aggregated_multidegree(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated multidegree over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_multidegree(adj, n, l, logger=logger), method=method)
+def compute_aggregated_multidegree(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated multidegree over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_multidegree(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
-def compute_multistrength(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute multistrength (total instrength + total outstrength) stratified by layer. Returns shape (n, l)."""
-    return compute_total_instrength(adj, n, l, logger=logger) + compute_total_outstrength(adj, n, l, logger=logger)
+def compute_multistrength(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute multistrength (total instrength + total outstrength) stratified by layer. Returns shape (nodes, layers)."""
+    return compute_total_instrength(adj, nodes=nodes, layers=layers, logger=logger) + compute_total_outstrength(adj, nodes=nodes, layers=layers, logger=logger)
 
 
-def compute_aggregated_multistrength(adj: sps.csr_matrix, n: int, l: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
-    """Compute aggregated multistrength over layers. Returns shape (n,)."""
-    return aggregate_metrics_over_layers(compute_multistrength(adj, n, l, logger=logger), method=method)
+def compute_aggregated_multistrength(adj: sps.csr_matrix, *, nodes: int, layers: int, method: str = "sum", logger: logging.Logger | None = None) -> np.ndarray:
+    """Compute aggregated multistrength over layers. Returns shape (nodes,)."""
+    return aggregate_metrics_over_layers(compute_multistrength(adj, nodes=nodes, layers=layers, logger=logger), method=method)
 
 
 def compute_multi_degree(
-    adj: sps.csr_matrix, n: int, l: int,
+    adj: sps.csr_matrix, *, nodes: int, layers: int,
     is_directed: bool = True,
     logger: logging.Logger | None = None,
 ) -> np.ndarray:
@@ -311,10 +311,10 @@ def compute_multi_degree(
     Parameters
     ----------
     adj : scipy.sparse matrix
-        Supra-adjacency matrix of shape (n*l, n*l).
-    n : int
+        Supra-adjacency matrix of shape (nodes*layers, nodes*layers).
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     is_directed : bool
         If True, returns in + out. If False, returns (in + out) / 2.
@@ -323,10 +323,10 @@ def compute_multi_degree(
     Returns
     -------
     np.ndarray
-        Multi-degree per physical node, shape (n,).
+        Multi-degree per physical node, shape (nodes,).
     """
-    indeg = compute_aggregated_indegree(adj, n, l, method="sum", logger=logger)
-    outdeg = compute_aggregated_outdegree(adj, n, l, method="sum", logger=logger)
+    indeg = compute_aggregated_indegree(adj, nodes=nodes, layers=layers, method="sum", logger=logger)
+    outdeg = compute_aggregated_outdegree(adj, nodes=nodes, layers=layers, method="sum", logger=logger)
     if is_directed:
         return indeg + outdeg
     else:
@@ -337,26 +337,26 @@ def compute_multi_degree(
 # Centrality implementations (integrated from hornet/node_based)
 # ---------------------------------------------------------------------------
 
-def compute_eigenvector_centrality(adj: sps.csr_matrix, n: int, l: int, logger: logging.Logger | None = None) -> np.ndarray:
+def compute_eigenvector_centrality(adj: sps.csr_matrix, *, nodes: int, layers: int, logger: logging.Logger | None = None) -> np.ndarray:
     """
     Compute multi-layer eigenvector centrality using the dominant eigenvector of A^T.
 
     Parameters
     ----------
     adj : scipy.sparse matrix
-        Supra-adjacency matrix of shape (n*l, n*l).
-    n : int
+        Supra-adjacency matrix of shape (nodes*layers, nodes*layers).
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     logger : logging.Logger, optional
 
     Returns
     -------
     np.ndarray
-        Max-normalized eigenvector centrality per physical node, shape (n,).
+        Max-normalized eigenvector centrality per physical node, shape (nodes,).
     """
-    NL = n * l
+    NL = nodes * layers
     if not sps.isspmatrix(adj):
         raise TypeError("adj must be a SciPy sparse matrix")
     if adj.shape != (NL, NL):
@@ -365,7 +365,7 @@ def compute_eigenvector_centrality(adj: sps.csr_matrix, n: int, l: int, logger:
     AT = adj.transpose().tocsc()
     lam, lvec = get_perron_eigenpair(AT, logger=logger)
 
-    X = np.reshape(lvec, (n, l), order="F")
+    X = np.reshape(lvec, (nodes, layers), order="F")
     eig_centrality = X.sum(axis=1)
 
     maxv = eig_centrality.max() if eig_centrality.size else 0.0
@@ -377,8 +377,7 @@ def compute_eigenvector_centrality(adj: sps.csr_matrix, n: int, l: int, logger:
 
 
 def compute_katz_centrality(
-    adj: sps.csr_matrix, n: int, l: int,
-    *,
+    adj: sps.csr_matrix, *, nodes: int, layers: int,
     alpha: float | None = None,
     solver: str = "direct",
     maxiter: int = 10000,
@@ -396,7 +395,7 @@ def compute_katz_centrality(
     By default, ``alpha = (1 - EPS) / rho(A)``, where ``rho(A)`` is the
     spectral radius of the supra-adjacency matrix ``A``. Callers can pass an
     explicit ``alpha`` for networks whose spectral radius is zero. The vector
-    x is then reshaped to ``(n, l)`` in column-major order, summed across
+    x is then reshaped to ``(nodes, layers)`` in column-major order, summed across
     layers to aggregate replicas, and max-normalized so that the largest score
     equals 1.
 
@@ -420,10 +419,10 @@ def compute_katz_centrality(
     Parameters
     ----------
     adj : scipy.sparse matrix
-        Supra-adjacency matrix of shape ``(n*l, n*l)``.
-    n : int
+        Supra-adjacency matrix of shape ``(nodes*layers, nodes*layers)``.
+    nodes : int
         Number of physical nodes.
-    l : int
+    layers : int
         Number of layers.
     alpha : float, optional
         Katz attenuation factor. It must be nonnegative and, when the
@@ -448,7 +447,7 @@ def compute_katz_centrality(
 
     Returns
     -------
-    katz : np.ndarray, shape (n,)
+    katz : np.ndarray, shape (nodes,)
         Max-normalized Katz centrality per physical node, ``float32``.
     eigenvalue : float, optional
         Rayleigh quotient estimate. Returned only when
@@ -459,14 +458,14 @@ def compute_katz_centrality(
     TypeError
         If ``adj`` is not a SciPy sparse matrix.
     ValueError
-        If ``adj.shape != (n*l, n*l)`` or if ``solver`` is not one of the
+        If ``adj.shape != (nodes*layers, nodes*layers)`` or if ``solver`` is not one of the
         supported values. Also raised when the automatic attenuation factor
         is undefined or an explicit factor is outside the stable range.
 
     The ``"neumann"`` branch uses a deterministic starting vector.
     """
     EPS = 1e-5
-    NL = n * l
+    NL = nodes * layers
     if not sps.isspmatrix(adj):
         raise TypeError("adj must be a SciPy sparse matrix")
     if adj.shape != (NL, NL):
@@ -536,7 +535,7 @@ def compute_katz_centrality(
 
     eigenvalue = float((x.T @ adj @ x) / (x.T @ x))
 
-    X = np.reshape(x, (n, l), order="F")
+    X = np.reshape(x, (nodes, layers), order="F")
     katz_centrality = X.sum(axis=1)
     maxv = katz_centrality.max() if katz_centrality.size else 0.0
     katz = np.zeros_like(katz_centrality, dtype=np.float32) if maxv == 0 else (katz_centrality / maxv).astype(np.float32)
@@ -552,8 +551,7 @@ def compute_katz_centrality(
 
 
 def compute_multi_rw_centrality(
-    adj: sps.csr_matrix, n: int, l: int, kind: str,
-    *,
+    adj: sps.csr_matrix, *, nodes: int, layers: int, kind: str,
     alpha: float = 0.85, tol: float = 1e-12, max_iter: int = 10000,
     logger: logging.Logger | None = None,
 ) -> np.ndarray:
@@ -569,10 +567,10 @@ def compute_multi_rw_centrality(
     Parameters
     ----------
     adj : scipy.sparse matrix
-        Supra-adjacency matrix of shape (n*l, n*l).
-    n : int
+        Supra-adjacency matrix of shape (nodes*layers, nodes*layers).
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     kind : str
         ``"classical"`` or ``"pagerank"``.
@@ -587,14 +585,14 @@ def compute_multi_rw_centrality(
     Returns
     -------
     np.ndarray
-        Max-normalized RW centrality per physical node, shape (n,).
+        Max-normalized RW centrality per physical node, shape (nodes,).
     """
-    NL = n * l
+    NL = nodes * layers
     kind = kind.lower().strip()
 
     if kind == "classical":
         tran_matrix = parsing_utils.build_transition_matrix_from_adjacency_matrix(
-            adj, n, l, kind="classical", logger=logger,
+            adj, nodes=nodes, layers=layers, kind="classical", logger=logger,
         )
         if NL == 0:
             return np.empty(0, dtype=np.float32)
@@ -603,9 +601,9 @@ def compute_multi_rw_centrality(
 
         total = vec.sum()
         if total == 0:
-            return np.zeros(n, dtype=np.float32)
+            return np.zeros(nodes, dtype=np.float32)
         x = vec / total
-        x = np.reshape(x, (n, l), order="F")
+        x = np.reshape(x, (nodes, layers), order="F")
         x = x.sum(axis=1)
 
         maxv = x.max() if x.size else 0.0
@@ -628,7 +626,7 @@ def compute_multi_rw_centrality(
             raise ValueError("max_iter must be a positive integer")
 
         P = parsing_utils.build_transition_matrix_from_adjacency_matrix(
-            adj, n, l, kind="classical", logger=logger,
+            adj, nodes=nodes, layers=layers, kind="classical", logger=logger,
         ).tocsr()
         if NL == 0:
             return np.empty(0, dtype=np.float32)
@@ -665,7 +663,7 @@ def compute_multi_rw_centrality(
         if logger and logger.isEnabledFor(logging.DEBUG):
             logger.debug("PageRank converged in %d iterations with error %.3e", it + 1, err)
 
-        X = np.reshape(x, (n, l), order="F")
+        X = np.reshape(x, (nodes, layers), order="F")
         mpc = X.sum(axis=1)
         maxv = mpc.max() if mpc.size else 0.0
         mpc_norm = np.zeros_like(mpc, dtype=np.float32) if maxv == 0 else (mpc / maxv).astype(np.float32)
@@ -680,7 +678,7 @@ def compute_multi_rw_centrality(
 
 
 def compute_multipagerank_centrality(
-    adj: sps.csr_matrix, n: int, l: int,
+    adj: sps.csr_matrix, *, nodes: int, layers: int,
     alpha: float = 0.85, tol: float = 1e-12, max_iter: int = 10000,
     logger: logging.Logger | None = None,
 ) -> np.ndarray:
@@ -693,10 +691,10 @@ def compute_multipagerank_centrality(
     Parameters
     ----------
     adj : scipy.sparse matrix
-        Supra-adjacency matrix of shape (n*l, n*l).
-    n : int
+        Supra-adjacency matrix of shape (nodes*layers, nodes*layers).
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     alpha : float
         Damping factor. Default 0.85.
@@ -709,16 +707,16 @@ def compute_multipagerank_centrality(
     Returns
     -------
     np.ndarray
-        Max-normalized PageRank per physical node, shape (n,).
+        Max-normalized PageRank per physical node, shape (nodes,).
     """
     return compute_multi_rw_centrality(
-        adj, n, l, kind="pagerank",
+        adj, nodes=nodes, layers=layers, kind="pagerank",
         alpha=alpha, tol=tol, max_iter=max_iter, logger=logger,
     )
 
 
 def compute_multi_hub_centrality(
-    adj: sps.csr_matrix, n: int, l: int,
+    adj: sps.csr_matrix, *, nodes: int, layers: int,
     eps: float = 1e-16, max_attempts: int = 10,
     approx: bool = False, approx_args: dict | None = None,
     logger: logging.Logger | None = None,
@@ -734,9 +732,9 @@ def compute_multi_hub_centrality(
     ----------
     adj : scipy.sparse matrix
         Supra-adjacency matrix.
-    n : int
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     eps : float
         Positive perturbation used by the approximate solver. The exact
@@ -753,9 +751,9 @@ def compute_multi_hub_centrality(
     Returns
     -------
     np.ndarray
-        Max-normalized hub centrality per physical node, shape (n,).
+        Max-normalized hub centrality per physical node, shape (nodes,).
     """
-    NL = n * l
+    NL = nodes * layers
     if not sps.isspmatrix(adj):
         raise TypeError("adj must be a SciPy sparse matrix")
     if adj.shape != (NL, NL):
@@ -778,7 +776,7 @@ def compute_multi_hub_centrality(
         eigenval, eigenvec = get_largest_real_eigenvalue(AA, logger=logger)
 
     eigenvec = np.abs(np.asarray(np.real_if_close(eigenvec), dtype=np.float64))
-    hc = eigenvec.reshape((l, n)).sum(axis=0)
+    hc = eigenvec.reshape((layers, nodes)).sum(axis=0)
     maxv = hc.max() if hc.size else 0.0
     hc = (
         np.zeros_like(hc, dtype=np.float32)
@@ -795,7 +793,7 @@ def compute_multi_hub_centrality(
 
 
 def compute_multi_authority_centrality(
-    adj: sps.csr_matrix, n: int, l: int,
+    adj: sps.csr_matrix, *, nodes: int, layers: int,
     eps: float = 1e-16, max_attempts: int = 10,
     approx: bool = False, approx_args: dict | None = None,
     logger: logging.Logger | None = None,
@@ -811,9 +809,9 @@ def compute_multi_authority_centrality(
     ----------
     adj : scipy.sparse matrix
         Supra-adjacency matrix.
-    n : int
+    nodes : int
         Number of nodes.
-    l : int
+    layers : int
         Number of layers.
     eps : float
         Positive perturbation used by the approximate solver. The exact
@@ -830,9 +828,9 @@ def compute_multi_authority_centrality(
     Returns
     -------
     np.ndarray
-        Max-normalized authority centrality per physical node, shape (n,).
+        Max-normalized authority centrality per physical node, shape (nodes,).
     """
-    NL = n * l
+    NL = nodes * layers
     if not sps.isspmatrix(adj):
         raise TypeError("adj must be a SciPy sparse matrix")
     if adj.shape != (NL, NL):
@@ -855,7 +853,7 @@ def compute_multi_authority_centrality(
         eigenval, eigenvec = get_largest_real_eigenvalue(AAT, logger=logger)
 
     eigenvec = np.abs(np.asarray(np.real_if_close(eigenvec), dtype=np.float64))
-    X = np.reshape(eigenvec, (n, l), order="F")
+    X = np.reshape(eigenvec, (nodes, layers), order="F")
     authority_centrality = X.sum(axis=1)
     maxv = authority_centrality.max() if authority_centrality.size else 0.0
     ac = (
@@ -921,8 +919,8 @@ def get_multi_degree(
         directed = True if is_directed is None else is_directed
         return compute_multi_degree(
             supra,
-            nodes,
-            layers,
+            nodes=nodes,
+            layers=layers,
             is_directed=directed,
             logger=logger,
         )
diff --git a/tests/test_dimension_arguments.py b/tests/test_dimension_arguments.py
index 30dc758..131de11 100644
--- a/tests/test_dimension_arguments.py
+++ b/tests/test_dimension_arguments.py
@@ -1,30 +1,59 @@
-"""Node and layer counts must be passed by name.
+"""Node and layer counts must be passed by name, everywhere in the package.
 
 Both are plain integers, so a transposed positional call cannot be detected at
-runtime. These functions previously took them in the order (layers, nodes), which
-is the opposite of the rest of the package.
+runtime. Rather than list the functions or the modules, these tests state the rule
+and check it against whatever the package currently defines, so a new function that
+takes both counts positionally fails here rather than shipping.
 """
 
+import importlib
+import inspect
+
 import pytest
-import scipy.sparse as sp
 
-from MuxVizPy import topology, versatility
-from MuxVizPy.utils import parsing
+from MuxVizPy import _LAZY_MODULES
+
+UTILS_MODULES = ["parsing", "io", "decomposition_utils", "katz_utils"]
+
+
+def analysis_modules():
+    names = [f"MuxVizPy.{n}" for n in _LAZY_MODULES if n != "utils"]
+    names += [f"MuxVizPy.utils.{n}" for n in UTILS_MODULES]
+    return [importlib.import_module(name) for name in names]
+
+
+MODULES = analysis_modules()
+
+
+def dimension_functions():
+    found = []
+    for module in MODULES:
+        for name, obj in vars(module).items():
+            if not inspect.isfunction(obj) or obj.__module__ != module.__name__:
+                continue
+            params = inspect.signature(obj).parameters
+            if "nodes" in params and "layers" in params:
+                found.append(pytest.param(obj, id=f"{module.__name__}.{name}"))
+    return found
+
+
+FUNCTIONS = dimension_functions()
+
+
+def test_the_package_defines_such_functions():
+    assert len(FUNCTIONS) > 40, "the discovery above found almost nothing"
 
-SUPRA = sp.eye(4, format="csr")
 
-FUNCTIONS = [
-    topology.get_connected_components,
-    topology.get_multi_path_statistics,
-    topology.get_SP_similarity_matrix,
-    versatility.get_multi_degree,
-    versatility.get_multi_Kcore_centrality,
-    parsing.supra_adjacency_to_network_list,
-    parsing.build_tensor_from_supra_adjacency_matrix,
-]
+@pytest.mark.parametrize("function", FUNCTIONS)
+def test_both_counts_are_keyword_only(function):
+    params = inspect.signature(function).parameters
+    for name in ("nodes", "layers"):
+        assert params[name].kind is inspect.Parameter.KEYWORD_ONLY, (
+            f"{function.__name__} takes `{name}` positionally"
+        )
 
 
-@pytest.mark.parametrize("function", FUNCTIONS, ids=lambda f: f.__name__)
-def test_dimensions_cannot_be_passed_positionally(function):
-    with pytest.raises(TypeError):
-        function(SUPRA, 2, 2)
+@pytest.mark.parametrize("function", FUNCTIONS)
+def test_no_short_dimension_aliases_remain(function):
+    params = inspect.signature(function).parameters
+    assert "n" not in params and "l" not in params
diff --git a/tests/test_edge_case_regressions.py b/tests/test_edge_case_regressions.py
index bebde97..38682da 100644
--- a/tests/test_edge_case_regressions.py
+++ b/tests/test_edge_case_regressions.py
@@ -41,7 +41,7 @@ def test_two_layer_global_overlap_uses_both_layers_in_denominator():
     supra = sp.block_diag([layer0, layer1], format="csr")
 
     overlap = global_descriptors.compute_average_global_overlap(
-        supra, n=3, l=2
+        supra, nodes=3, layers=2
     )
 
     assert overlap == pytest.approx(1.0 / 3.0)
@@ -54,7 +54,7 @@ def test_eigenvector_centrality_supports_tiny_nonnegative_networks():
     with warnings.catch_warnings():
         warnings.simplefilter("ignore", RuntimeWarning)
         centrality = versatility.compute_eigenvector_centrality(
-            adjacency, n=2, l=1
+            adjacency, nodes=2, layers=1
         )
 
     assert centrality.shape == (2,)
@@ -201,7 +201,7 @@ def test_empty_adjacency_has_zero_average_global_clustering():
     """An edgeless graph has a defined global clustering coefficient of zero."""
     result = (
         global_descriptors.compute_average_global_clustering_coefficient(
-            sp.csr_matrix((3, 3)), n=3, l=1
+            sp.csr_matrix((3, 3)), nodes=3, layers=1
         )
     )
 
@@ -432,8 +432,8 @@ def test_pagerank_rejects_an_invalid_damping_factor():
     with pytest.raises(ValueError, match="alpha"):
         versatility.compute_multi_rw_centrality(
             adjacency,
-            n=3,
-            l=1,
+            nodes=3,
+            layers=1,
             kind="pagerank",
             alpha=1.5,
         )
@@ -447,7 +447,7 @@ def test_katz_rejects_undefined_automatic_alpha_for_zero_spectral_radius():
 
     with pytest.raises(ValueError, match="spectral radius|alpha"):
         versatility.compute_katz_centrality(
-            adjacency, n=10, l=1, solver="direct"
+            adjacency, nodes=10, layers=1, solver="direct"
         )
 
 
@@ -467,10 +467,10 @@ def signed_eigenpair(_matrix, logger=None):
     )
 
     hubs = versatility.compute_multi_hub_centrality(
-        adjacency, n=8, l=1, max_attempts=1
+        adjacency, nodes=8, layers=1, max_attempts=1
     )
     authorities = versatility.compute_multi_authority_centrality(
-        adjacency, n=8, l=1, max_attempts=1
+        adjacency, nodes=8, layers=1, max_attempts=1
     )
 
     assert np.all(hubs >= 0)
diff --git a/tests/test_global_descriptors.py b/tests/test_global_descriptors.py
index 53331d6..bd273b4 100644
--- a/tests/test_global_descriptors.py
+++ b/tests/test_global_descriptors.py
@@ -26,13 +26,13 @@ class TestGlobalDescriptorsCorrectness:
 
     def test_agcc_returns_float(self, net_adjacency, net_n, net_l):
         result = global_descriptors.compute_average_global_clustering_coefficient(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert isinstance(result, float)
 
     def test_agcc_non_negative(self, net_adjacency, net_n, net_l):
         result = global_descriptors.compute_average_global_clustering_coefficient(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert result >= 0.0
 
@@ -40,27 +40,27 @@ def test_agcc_non_negative(self, net_adjacency, net_n, net_l):
 
     def test_agov_returns_float(self, net_adjacency, net_n, net_l):
         result = global_descriptors.compute_average_global_overlap(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert isinstance(result, float)
 
     def test_agov_non_negative(self, net_adjacency, net_n, net_l):
         result = global_descriptors.compute_average_global_overlap(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert result >= 0.0
 
     def test_agov_raises_single_layer(self, net_adjacency, net_n):
         with pytest.raises(ValueError):
-            global_descriptors.compute_average_global_overlap(net_adjacency, net_n, 1)
+            global_descriptors.compute_average_global_overlap(net_adjacency, nodes=net_n, layers=1)
 
     def test_agov_binary_equals_weighted_on_binary_adj(self, net_adjacency, net_n, net_l):
         """For a binary adjacency matrix, weighted=True and weighted=False must agree."""
         r_bin = global_descriptors.compute_average_global_overlap(
-            net_adjacency, net_n, net_l, weighted=False
+            net_adjacency, nodes=net_n, layers=net_l, weighted=False
         )
         r_wt = global_descriptors.compute_average_global_overlap(
-            net_adjacency, net_n, net_l, weighted=True
+            net_adjacency, nodes=net_n, layers=net_l, weighted=True
         )
         assert r_bin == pytest.approx(r_wt, abs=1e-10)
 
@@ -68,25 +68,25 @@ def test_agov_binary_equals_weighted_on_binary_adj(self, net_adjacency, net_n, n
 
     def test_agov_mat_shape(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert M.shape == (net_l, net_l)
 
     def test_agov_mat_diagonal_is_one(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         np.testing.assert_array_equal(np.diag(M), np.ones(net_l))
 
     def test_agov_mat_symmetric(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         np.testing.assert_array_almost_equal(M, M.T)
 
     def test_agov_mat_off_diagonal_in_range(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         off_diag = M[~np.eye(net_l, dtype=bool)]
         assert np.all(off_diag >= 0.0)
@@ -94,31 +94,31 @@ def test_agov_mat_off_diagonal_in_range(self, net_adjacency, net_n, net_l):
 
     def test_agov_mat_raises_single_layer(self, net_adjacency, net_n):
         with pytest.raises(ValueError):
-            global_descriptors.compute_average_global_overlap_matrix(net_adjacency, net_n, 1)
+            global_descriptors.compute_average_global_overlap_matrix(net_adjacency, nodes=net_n, layers=1)
 
     # --- average global node overlap matrix (L×L) -------------------------
 
     def test_agnov_mat_shape(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_node_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert M.shape == (net_l, net_l)
 
     def test_agnov_mat_diagonal_is_one(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_node_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         np.testing.assert_array_equal(np.diag(M), np.ones(net_l))
 
     def test_agnov_mat_symmetric(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_node_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         np.testing.assert_array_almost_equal(M, M.T)
 
     def test_agnov_mat_in_range(self, net_adjacency, net_n, net_l):
         M = global_descriptors.compute_average_global_node_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert np.all(M >= 0.0)
         assert np.all(M <= 1.0)
@@ -126,7 +126,7 @@ def test_agnov_mat_in_range(self, net_adjacency, net_n, net_l):
     def test_agnov_mat_raises_single_layer(self, net_adjacency, net_n):
         with pytest.raises(ValueError):
             global_descriptors.compute_average_global_node_overlap_matrix(
-                net_adjacency, net_n, 1
+                net_adjacency, nodes=net_n, layers=1
             )
 
     # --- cross-function invariants ----------------------------------------
@@ -136,10 +136,10 @@ def test_agov_mat_consistent_with_scalar(self, net_adjacency, net_n, net_l):
         directed networks with the R-compatible normalisation, but agov_mat
         (Dice) is always in [0, 1] — verify they are both non-negative."""
         scalar = global_descriptors.compute_average_global_overlap(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         mat = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         assert scalar >= 0.0
         assert np.all(mat >= 0.0)
@@ -148,10 +148,10 @@ def test_agnov_mat_le_agov_mat_off_diagonal(self, net_adjacency, net_n, net_l):
         """Node overlap >= edge overlap off-diagonal: if two layers share an edge
         they must share both endpoint nodes, so edge overlap <= node overlap."""
         edge_mat = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         node_mat = global_descriptors.compute_average_global_node_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         mask = ~np.eye(net_l, dtype=bool)
         # node overlap is per-node (fraction of N); edge overlap is Dice (fraction of edges)
@@ -176,7 +176,7 @@ def test_agcc_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
         if "agcc" not in net_muxviz_results:
             pytest.skip("agcc not in reference results")
         computed = global_descriptors.compute_average_global_clustering_coefficient(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         # R as.vector() on a scalar gives a length-1 list in JSON
         expected = float(np.array(net_muxviz_results["agcc"]).ravel()[0])
@@ -187,7 +187,7 @@ def test_agov_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
         if "agov" not in net_muxviz_results:
             pytest.skip("agov not in reference results")
         computed = global_descriptors.compute_average_global_overlap(
-            net_adjacency, net_n, net_l, weighted=False
+            net_adjacency, nodes=net_n, layers=net_l, weighted=False
         )
         expected = float(np.array(net_muxviz_results["agov"]).ravel()[0])
         np.testing.assert_allclose(computed, expected, rtol=1e-4, atol=1e-4,
@@ -197,7 +197,7 @@ def test_agov_mat_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_result
         if "agov_mat" not in net_muxviz_results:
             pytest.skip("agov_mat not in reference results")
         computed = global_descriptors.compute_average_global_overlap_matrix(
-            net_adjacency, net_n, net_l, weighted=False
+            net_adjacency, nodes=net_n, layers=net_l, weighted=False
         )
         # R as.vector() on a matrix is column-major (Fortran order)
         r_flat = np.array(net_muxviz_results["agov_mat"], dtype=float)
@@ -211,7 +211,7 @@ def test_agnov_mat_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_resul
         if "agnov_mat" not in net_muxviz_results:
             pytest.skip("agnov_mat not in reference results")
         computed = global_descriptors.compute_average_global_node_overlap_matrix(
-            net_adjacency, net_n, net_l
+            net_adjacency, nodes=net_n, layers=net_l
         )
         # R as.vector() on a matrix is column-major (Fortran order)
         r_flat = np.array(net_muxviz_results["agnov_mat"], dtype=float)
diff --git a/tests/test_mesoscale.py b/tests/test_mesoscale.py
index a756baa..4cd6c1a 100644
--- a/tests/test_mesoscale.py
+++ b/tests/test_mesoscale.py
@@ -53,30 +53,30 @@ class TestLocalClusteringCorrectness:
     """compute_local_clustering_coefficient returns sane types/shapes/values."""
 
     def test_returns_ndarray(self, net_adjacency, net_n, net_l):
-        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, net_n, net_l)
+        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, nodes=net_n, layers=net_l)
         assert isinstance(result, np.ndarray)
 
     def test_shape(self, net_adjacency, net_n, net_l):
-        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, net_n, net_l)
+        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n,)
 
     def test_range_non_negative(self, net_adjacency, net_n, net_l):
-        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, net_n, net_l)
+        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, nodes=net_n, layers=net_l)
         assert np.all(result >= 0.0), f"Negative values found: {result[result < 0]}"
 
     def test_range_at_most_one(self, net_adjacency, net_n, net_l):
-        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, net_n, net_l)
+        result = mesoscale.compute_local_clustering_coefficient(net_adjacency, nodes=net_n, layers=net_l)
         assert np.all(result <= 1.0), f"Values > 1 found: {result[result > 1]}"
 
     def test_zero_adjacency_gives_zero_clustering(self):
         n, l = 5, 3
         zero_adj = sp.csr_matrix((n * l, n * l))
-        result = mesoscale.compute_local_clustering_coefficient(zero_adj, n, l)
+        result = mesoscale.compute_local_clustering_coefficient(zero_adj, nodes=n, layers=l)
         assert np.allclose(result, 0.0)
 
     def test_sample_network(self, sample_adjacency, n_nodes, n_layers):
         result = mesoscale.compute_local_clustering_coefficient(
-            sample_adjacency, n_nodes, n_layers
+            sample_adjacency, nodes=n_nodes, layers=n_layers
         )
         assert result.shape == (n_nodes,)
         assert np.all(result >= 0.0)
@@ -98,7 +98,7 @@ def test_local_clus_vs_muxviz(
         if "local_clus" not in net_muxviz_results:
             pytest.skip(f"'local_clus' not in reference results for '{network_config}'")
 
-        computed = mesoscale.compute_local_clustering_coefficient(net_adjacency, net_n, net_l)
+        computed = mesoscale.compute_local_clustering_coefficient(net_adjacency, nodes=net_n, layers=net_l)
         expected = np.asarray(net_muxviz_results["local_clus"], dtype=np.float64).ravel()
 
         compare_metrics(
diff --git a/tests/test_regressions.py b/tests/test_regressions.py
index a1259ec..d641bbc 100644
--- a/tests/test_regressions.py
+++ b/tests/test_regressions.py
@@ -118,7 +118,7 @@ def test_pagerank_transition_matrix_is_stochastic():
     )
 
     transition = parsing.build_transition_matrix_from_adjacency_matrix(
-        adjacency, n=4, l=1, kind="pagerank", alpha=0.85
+        adjacency, nodes=4, layers=1, kind="pagerank", alpha=0.85
     )
 
     np.testing.assert_allclose(
diff --git a/tests/test_utils_parsing.py b/tests/test_utils_parsing.py
index ef43016..bdc9227 100644
--- a/tests/test_utils_parsing.py
+++ b/tests/test_utils_parsing.py
@@ -206,7 +206,7 @@ def test_virus_transition_rejects_invalid_intralayer_weight(
 
     def test_classical_row_stochastic(self, sample_adjacency, n_nodes, n_layers):
         T = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="classical"
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="classical"
         )
         row_sums = np.asarray(T.sum(axis=1)).ravel()
         # Rows with outgoing edges should sum to 1; zero-degree rows sum to 0
@@ -217,25 +217,25 @@ def test_classical_row_stochastic(self, sample_adjacency, n_nodes, n_layers):
 
     def test_classical_nonnegative(self, sample_adjacency, n_nodes, n_layers):
         T = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="classical"
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="classical"
         )
         assert T.min() >= 0.0
 
     def test_classical_shape(self, sample_adjacency, n_nodes, n_layers):
         T = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="classical"
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="classical"
         )
         assert T.shape == (NL, NL)
 
     def test_pagerank_shape(self, sample_adjacency, n_nodes, n_layers):
         T = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="pagerank", alpha=0.85
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="pagerank", alpha=0.85
         )
         assert T.shape == (NL, NL)
 
     def test_pagerank_nonnegative(self, sample_adjacency, n_nodes, n_layers):
         T = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="pagerank", alpha=0.85
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="pagerank", alpha=0.85
         )
         assert T.min() >= 0.0
 
@@ -244,10 +244,10 @@ def test_pagerank_includes_teleportation_and_dangling_mass(
     ):
         alpha = 0.85
         T_class = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="classical"
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="classical"
         )
         T_pr = parsing.build_transition_matrix_from_adjacency_matrix(
-            sample_adjacency, n_nodes, n_layers, kind="pagerank", alpha=alpha
+            sample_adjacency, nodes=n_nodes, layers=n_layers, kind="pagerank", alpha=alpha
         )
 
         expected = alpha * T_class.toarray()
@@ -264,24 +264,24 @@ def test_pagerank_includes_teleportation_and_dangling_mass(
     def test_pagerank_invalid_alpha_raises(self, sample_adjacency, n_nodes, n_layers):
         with pytest.raises(ValueError):
             parsing.build_transition_matrix_from_adjacency_matrix(
-                sample_adjacency, n_nodes, n_layers, kind="pagerank", alpha=0.0
+                sample_adjacency, nodes=n_nodes, layers=n_layers, kind="pagerank", alpha=0.0
             )
         with pytest.raises(ValueError):
             parsing.build_transition_matrix_from_adjacency_matrix(
-                sample_adjacency, n_nodes, n_layers, kind="pagerank", alpha=1.5
+                sample_adjacency, nodes=n_nodes, layers=n_layers, kind="pagerank", alpha=1.5
             )
 
     def test_unknown_kind_raises(self, sample_adjacency, n_nodes, n_layers):
         with pytest.raises(NotImplementedError):
             parsing.build_transition_matrix_from_adjacency_matrix(
-                sample_adjacency, n_nodes, n_layers, kind="bogus"
+                sample_adjacency, nodes=n_nodes, layers=n_layers, kind="bogus"
             )
 
     @pytest.mark.parametrize("kind", ["diffusive", "maxent", "physical", "relaxed-physical"])
     def test_unimplemented_kinds_raise(self, sample_adjacency, n_nodes, n_layers, kind):
         with pytest.raises(NotImplementedError):
             parsing.build_transition_matrix_from_adjacency_matrix(
-                sample_adjacency, n_nodes, n_layers, kind=kind
+                sample_adjacency, nodes=n_nodes, layers=n_layers, kind=kind
             )
 
 
diff --git a/tests/test_versatility.py b/tests/test_versatility.py
index 283e647..ee844a9 100644
--- a/tests/test_versatility.py
+++ b/tests/test_versatility.py
@@ -53,26 +53,26 @@ def test_approximate_largest_eigenvalue_returns_tuple(self, net_adjacency, net_n
     # --- block accumulation helpers ---------------------------------------
 
     def test_is_in_diagonal_block(self):
-        assert versatility.is_in_diagonal_block(0, 1, n=10, l=3) is True
-        assert versatility.is_in_diagonal_block(0, 10, n=10, l=3) is False
-        assert versatility.is_in_diagonal_block(15, 19, n=10, l=3) is True
+        assert versatility.is_in_diagonal_block(0, 1, nodes=10, layers=3) is True
+        assert versatility.is_in_diagonal_block(0, 10, nodes=10, layers=3) is False
+        assert versatility.is_in_diagonal_block(15, 19, nodes=10, layers=3) is True
 
     def test_accumulate_diagonal_blocks_shape(self, net_adjacency, net_n, net_l):
         result = versatility._accumulate_on_diagonal_blocks(
-            net_adjacency, net_n, net_l, is_out_of_diagonal=False,
+            net_adjacency, nodes=net_n, layers=net_l, is_out_of_diagonal=False,
         )
         assert result.shape == (net_n, net_l)
 
     def test_accumulate_off_diagonal_blocks_shape(self, net_adjacency, net_n, net_l):
         result = versatility._accumulate_on_diagonal_blocks(
-            net_adjacency, net_n, net_l, is_out_of_diagonal=True,
+            net_adjacency, nodes=net_n, layers=net_l, is_out_of_diagonal=True,
         )
         assert result.shape == (net_n, net_l)
 
     def test_accumulate_wrong_shape_raises(self, net_n, net_l):
         wrong = sp.eye(5, format="csr")
         with pytest.raises(ValueError, match="does not match"):
-            versatility._accumulate_on_diagonal_blocks(wrong, net_n, net_l, is_out_of_diagonal=False)
+            versatility._accumulate_on_diagonal_blocks(wrong, nodes=net_n, layers=net_l, is_out_of_diagonal=False)
 
     # --- aggregate_metrics_over_layers ------------------------------------
 
@@ -100,57 +100,57 @@ def test_aggregate_unknown_raises(self):
     # --- per-layer degree / strength shapes -------------------------------
 
     def test_compute_indegree_shape(self, net_adjacency, net_n, net_l):
-        result = versatility.compute_indegree(net_adjacency, net_n, net_l)
+        result = versatility.compute_indegree(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_aggregated_indegree_shape(self, net_adjacency, net_n, net_l):
-        result = versatility.compute_aggregated_indegree(net_adjacency, net_n, net_l)
+        result = versatility.compute_aggregated_indegree(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n,)
 
     def test_compute_outdegree_shape(self, net_adjacency, net_n, net_l):
-        result = versatility.compute_outdegree(net_adjacency, net_n, net_l)
+        result = versatility.compute_outdegree(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_aggregated_outdegree_shape(self, net_adjacency, net_n, net_l):
-        result = versatility.compute_aggregated_outdegree(net_adjacency, net_n, net_l)
+        result = versatility.compute_aggregated_outdegree(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n,)
 
     def test_compute_instrength_shape(self, net_interaction, net_n, net_l):
-        result = versatility.compute_instrength(net_interaction, net_n, net_l)
+        result = versatility.compute_instrength(net_interaction, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_outstrength_shape(self, net_interaction, net_n, net_l):
-        result = versatility.compute_outstrength(net_interaction, net_n, net_l)
+        result = versatility.compute_outstrength(net_interaction, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_multiindegree_shape(self, net_adjacency, net_n, net_l):
-        result = versatility.compute_multiindegree(net_adjacency, net_n, net_l)
+        result = versatility.compute_multiindegree(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_multioutdegree_shape(self, net_adjacency, net_n, net_l):
-        result = versatility.compute_multioutdegree(net_adjacency, net_n, net_l)
+        result = versatility.compute_multioutdegree(net_adjacency, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_multiinstrength_shape(self, net_interaction, net_n, net_l):
-        result = versatility.compute_multiinstrength(net_interaction, net_n, net_l)
+        result = versatility.compute_multiinstrength(net_interaction, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     def test_compute_multioutstrength_shape(self, net_interaction, net_n, net_l):
-        result = versatility.compute_multioutstrength(net_interaction, net_n, net_l)
+        result = versatility.compute_multioutstrength(net_interaction, nodes=net_n, layers=net_l)
         assert result.shape == (net_n, net_l)
 
     # --- degree consistency: aggregated indegreesum = indegree + multiindegree
 
     def test_indegree_plus_multiindegree_equals_total(self, net_adjacency, net_n, net_l):
-        intra = versatility.compute_aggregated_indegree(net_adjacency, net_n, net_l)
-        inter = versatility.compute_aggregated_multiindegree(net_adjacency, net_n, net_l)
+        intra = versatility.compute_aggregated_indegree(net_adjacency, nodes=net_n, layers=net_l)
+        inter = versatility.compute_aggregated_multiindegree(net_adjacency, nodes=net_n, layers=net_l)
         total = np.asarray(net_adjacency.sum(axis=0)).ravel()
         total_per_node = total.reshape(net_l, net_n).sum(axis=0)
         np.testing.assert_allclose(intra + inter, total_per_node, atol=1e-10)
 
     def test_outdegree_plus_multioutdegree_equals_total(self, net_adjacency, net_n, net_l):
-        intra = versatility.compute_aggregated_outdegree(net_adjacency, net_n, net_l)
-        inter = versatility.compute_aggregated_multioutdegree(net_adjacency, net_n, net_l)
+        intra = versatility.compute_aggregated_outdegree(net_adjacency, nodes=net_n, layers=net_l)
+        inter = versatility.compute_aggregated_multioutdegree(net_adjacency, nodes=net_n, layers=net_l)
         total = np.asarray(net_adjacency.sum(axis=1)).ravel()
         total_per_node = total.reshape(net_l, net_n).sum(axis=0)
         np.testing.assert_allclose(intra + inter, total_per_node, atol=1e-10)
@@ -158,13 +158,13 @@ def test_outdegree_plus_multioutdegree_equals_total(self, net_adjacency, net_n,
     # --- centrality shapes and ranges -------------------------------------
 
     def test_compute_eigenvector_centrality_shape(self, net_adjacency, net_n, net_l):
-        ec = versatility.compute_eigenvector_centrality(net_adjacency, net_n, net_l)
+        ec = versatility.compute_eigenvector_centrality(net_adjacency, nodes=net_n, layers=net_l)
         assert ec.shape == (net_n,)
         assert ec.max() == pytest.approx(1.0, abs=1e-6)
 
     def test_compute_katz_centrality_exact_shape(self, net_interaction, net_n, net_l):
         katz, eigenvalue = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="direct", return_eigenvalue=True,
+            net_interaction, nodes=net_n, layers=net_l, solver="direct", return_eigenvalue=True,
         )
         assert katz.shape == (net_n,)
         assert katz.max() == pytest.approx(1.0, abs=1e-6)
@@ -173,7 +173,7 @@ def test_compute_katz_centrality_exact_shape(self, net_interaction, net_n, net_l
     def test_compute_katz_centrality_neumann_shape(self, net_interaction, net_n, net_l):
         np.random.seed(42)
         katz, eigenvalue = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="neumann",
+            net_interaction, nodes=net_n, layers=net_l, solver="neumann",
             maxiter=100000, tol=1e-4,
             return_eigenvalue=True,
         )
@@ -182,36 +182,36 @@ def test_compute_katz_centrality_neumann_shape(self, net_interaction, net_n, net
         assert isinstance(eigenvalue, float)
 
     def test_compute_multi_rw_centrality_classical_shape(self, net_adjacency, net_n, net_l):
-        rc = versatility.compute_multi_rw_centrality(net_adjacency, net_n, net_l, kind="classical")
+        rc = versatility.compute_multi_rw_centrality(net_adjacency, nodes=net_n, layers=net_l, kind="classical")
         assert rc.shape == (net_n,)
         assert rc.max() == pytest.approx(1.0, abs=1e-6)
 
     def test_compute_multi_rw_centrality_pagerank_shape(self, net_adjacency, net_n, net_l):
-        pr = versatility.compute_multi_rw_centrality(net_adjacency, net_n, net_l, kind="pagerank")
+        pr = versatility.compute_multi_rw_centrality(net_adjacency, nodes=net_n, layers=net_l, kind="pagerank")
         assert pr.shape == (net_n,)
         assert pr.max() == pytest.approx(1.0, abs=1e-6)
 
     def test_compute_multipagerank_centrality_shape(self, net_adjacency, net_n, net_l):
-        pr = versatility.compute_multipagerank_centrality(net_adjacency, net_n, net_l)
+        pr = versatility.compute_multipagerank_centrality(net_adjacency, nodes=net_n, layers=net_l)
         assert pr.shape == (net_n,)
         assert pr.max() == pytest.approx(1.0, abs=1e-6)
 
     def test_multipagerank_wrapper_matches_direct(self, net_adjacency, net_n, net_l):
         """compute_multipagerank_centrality is a wrapper and must match compute_multi_rw_centrality(kind='pagerank')."""
-        direct = versatility.compute_multi_rw_centrality(net_adjacency, net_n, net_l, kind="pagerank")
-        wrapper = versatility.compute_multipagerank_centrality(net_adjacency, net_n, net_l)
+        direct = versatility.compute_multi_rw_centrality(net_adjacency, nodes=net_n, layers=net_l, kind="pagerank")
+        wrapper = versatility.compute_multipagerank_centrality(net_adjacency, nodes=net_n, layers=net_l)
         np.testing.assert_array_equal(direct, wrapper)
 
     def test_compute_multi_rw_centrality_unknown_kind_raises(self, net_adjacency, net_n, net_l):
         with pytest.raises(ValueError, match="Unknown RW kind"):
-            versatility.compute_multi_rw_centrality(net_adjacency, net_n, net_l, kind="bogus")
+            versatility.compute_multi_rw_centrality(net_adjacency, nodes=net_n, layers=net_l, kind="bogus")
 
     def test_compute_multi_hub_centrality_shape(self, net_adjacency, net_n, net_l):
-        hc = versatility.compute_multi_hub_centrality(net_adjacency, net_n, net_l)
+        hc = versatility.compute_multi_hub_centrality(net_adjacency, nodes=net_n, layers=net_l)
         assert hc.shape == (net_n,)
 
     def test_compute_multi_authority_centrality_shape(self, net_interaction, net_n, net_l):
-        ac = versatility.compute_multi_authority_centrality(net_interaction, net_n, net_l)
+        ac = versatility.compute_multi_authority_centrality(net_interaction, nodes=net_n, layers=net_l)
         assert ac.shape == (net_n,)
 
     # --- public API shapes ------------------------------------------------
@@ -221,7 +221,7 @@ def test_get_multi_degree_shape(self, net_adjacency, net_n, net_l):
         assert deg.shape == (net_n,)
 
     def test_compute_eigenvector_centrality_public_shape(self, net_adjacency, net_n, net_l):
-        ec = versatility.compute_eigenvector_centrality(net_adjacency, net_n, net_l)
+        ec = versatility.compute_eigenvector_centrality(net_adjacency, nodes=net_n, layers=net_l)
         assert ec.shape == (net_n,)
 
 
@@ -235,32 +235,32 @@ class TestKatzCentralityAgreement:
     def test_katz_invalid_solver_raises(self, net_interaction, net_n, net_l):
         with pytest.raises(ValueError, match="Unknown solver"):
             versatility.compute_katz_centrality(
-                net_interaction, net_n, net_l, solver="invalid_solver",
+                net_interaction, nodes=net_n, layers=net_l, solver="invalid_solver",
             )
 
     def test_katz_neumann_warns_when_not_converged(self, net_interaction, net_n, net_l):
         with pytest.warns(UserWarning, match="did not converge"):
             versatility.compute_katz_centrality(
-                net_interaction, net_n, net_l, solver="neumann",
+                net_interaction, nodes=net_n, layers=net_l, solver="neumann",
                 maxiter=1, tol=1e-10,
             )
 
     def test_katz_gmres_matches_direct(self, net_interaction, net_n, net_l):
         exact = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="direct",
+            net_interaction, nodes=net_n, layers=net_l, solver="direct",
         )
         gmres_result = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="gmres",
+            net_interaction, nodes=net_n, layers=net_l, solver="gmres",
             maxiter=1000, tol=1e-8,
         )
         np.testing.assert_allclose(gmres_result, exact, atol=1e-3)
 
     def test_katz_bicgstab_matches_direct(self, net_interaction, net_n, net_l):
         exact = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="direct",
+            net_interaction, nodes=net_n, layers=net_l, solver="direct",
         )
         bicgstab_result = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="bicgstab",
+            net_interaction, nodes=net_n, layers=net_l, solver="bicgstab",
             maxiter=1000, tol=1e-8,
         )
         np.testing.assert_allclose(bicgstab_result, exact, atol=5e-4)
@@ -302,7 +302,7 @@ def test_pagerank_rejects_invalid_alpha(self, alpha):
 
         with pytest.raises(ValueError, match="alpha"):
             versatility.compute_multi_rw_centrality(
-                adjacency, n=3, l=1, kind="pagerank", alpha=alpha
+                adjacency, nodes=3, layers=1, kind="pagerank", alpha=alpha
             )
 
     @pytest.mark.parametrize(
@@ -317,7 +317,7 @@ def test_pagerank_rejects_invalid_iteration_parameters(
 
         with pytest.raises(ValueError, match=parameter):
             versatility.compute_multi_rw_centrality(
-                adjacency, n=3, l=1, kind="pagerank", **kwargs
+                adjacency, nodes=3, layers=1, kind="pagerank", **kwargs
             )
 
     def test_pagerank_raises_when_iteration_does_not_converge(self):
@@ -328,8 +328,8 @@ def test_pagerank_raises_when_iteration_does_not_converge(self):
         with pytest.raises(RuntimeError, match="converge"):
             versatility.compute_multi_rw_centrality(
                 adjacency,
-                n=3,
-                l=1,
+                nodes=3,
+                layers=1,
                 kind="pagerank",
                 tol=1e-16,
                 max_iter=1,
@@ -337,7 +337,7 @@ def test_pagerank_raises_when_iteration_does_not_converge(self):
 
     def test_pagerank_empty_network_returns_empty_result(self):
         result = versatility.compute_multi_rw_centrality(
-            sp.csr_matrix((0, 0)), n=0, l=1, kind="pagerank"
+            sp.csr_matrix((0, 0)), nodes=0, layers=1, kind="pagerank"
         )
 
         assert result.shape == (0,)
@@ -351,7 +351,7 @@ def test_eigenvector_is_deterministic_on_bipartite_network(self):
         )
 
         results = [
-            versatility.compute_eigenvector_centrality(adjacency, n=4, l=1)
+            versatility.compute_eigenvector_centrality(adjacency, nodes=4, layers=1)
             for _ in range(3)
         ]
 
@@ -373,7 +373,7 @@ def test_classical_random_walk_is_nonnegative_and_deterministic(self):
 
         results = [
             versatility.compute_multi_rw_centrality(
-                adjacency, n=4, l=1, kind="classical"
+                adjacency, nodes=4, layers=1, kind="classical"
             )
             for _ in range(3)
         ]
@@ -386,14 +386,14 @@ def test_empty_networks_return_empty_centralities(self):
         adjacency = sp.csr_matrix((0, 0))
 
         results = [
-            versatility.compute_eigenvector_centrality(adjacency, n=0, l=1),
-            versatility.compute_katz_centrality(adjacency, n=0, l=1),
+            versatility.compute_eigenvector_centrality(adjacency, nodes=0, layers=1),
+            versatility.compute_katz_centrality(adjacency, nodes=0, layers=1),
             versatility.compute_multi_rw_centrality(
-                adjacency, n=0, l=1, kind="classical"
+                adjacency, nodes=0, layers=1, kind="classical"
             ),
-            versatility.compute_multi_hub_centrality(adjacency, n=0, l=1),
+            versatility.compute_multi_hub_centrality(adjacency, nodes=0, layers=1),
             versatility.compute_multi_authority_centrality(
-                adjacency, n=0, l=1
+                adjacency, nodes=0, layers=1
             ),
         ]
 
@@ -405,10 +405,10 @@ def test_explicit_katz_alpha_supports_zero_radius_network(self):
         )
 
         first = versatility.compute_katz_centrality(
-            adjacency, n=10, l=1, alpha=0.5, solver="neumann"
+            adjacency, nodes=10, layers=1, alpha=0.5, solver="neumann"
         )
         second = versatility.compute_katz_centrality(
-            adjacency, n=10, l=1, alpha=0.5, solver="neumann"
+            adjacency, nodes=10, layers=1, alpha=0.5, solver="neumann"
         )
 
         assert np.all(first >= 0)
@@ -420,7 +420,7 @@ def test_katz_rejects_invalid_explicit_alpha(self, alpha):
 
         with pytest.raises(ValueError, match="alpha"):
             versatility.compute_katz_centrality(
-                adjacency, n=2, l=1, alpha=alpha
+                adjacency, nodes=2, layers=1, alpha=alpha
             )
 
     def test_exact_hits_balances_tied_components_deterministically(self):
@@ -433,12 +433,12 @@ def test_exact_hits_balances_tied_components_deterministically(self):
         )
 
         hub_results = [
-            versatility.compute_multi_hub_centrality(adjacency, n=8, l=1)
+            versatility.compute_multi_hub_centrality(adjacency, nodes=8, layers=1)
             for _ in range(3)
         ]
         authority_results = [
             versatility.compute_multi_authority_centrality(
-                adjacency, n=8, l=1
+                adjacency, nodes=8, layers=1
             )
             for _ in range(3)
         ]
@@ -543,7 +543,7 @@ class TestVersatilityReference:
 
     def test_katz_exact_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
         computed = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="direct", return_eigenvalue=False,
+            net_interaction, nodes=net_n, layers=net_l, solver="direct", return_eigenvalue=False,
         )
         # The Python and R eigensolvers differ slightly in their numerical results.
         compare_metrics(computed, net_muxviz_results["katz"], "Katz exact (vs muxViz R)",
@@ -552,7 +552,7 @@ def test_katz_exact_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_re
     def test_katz_neumann_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
         np.random.seed(42)
         computed = versatility.compute_katz_centrality(
-            net_interaction, net_n, net_l, solver="neumann",
+            net_interaction, nodes=net_n, layers=net_l, solver="neumann",
             maxiter=100000, tol=1e-4,
             return_eigenvalue=False,
         )
@@ -560,56 +560,56 @@ def test_katz_neumann_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_
                          rtol=0.05, atol=0.05)
 
     def test_pagerank_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_multipagerank_centrality(net_interaction, net_n, net_l)
+        computed = versatility.compute_multipagerank_centrality(net_interaction, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["pagerank"], "PageRank (vs muxViz R)")
 
     def test_hub_exact_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_multi_hub_centrality(net_adjacency, net_n, net_l)
+        computed = versatility.compute_multi_hub_centrality(net_adjacency, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["hub"], "Hub (vs muxViz R)")
 
     def test_hub_approx_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
         np.random.seed(42)
         computed = versatility.compute_multi_hub_centrality(
-            net_adjacency, net_n, net_l, approx=True,
+            net_adjacency, nodes=net_n, layers=net_l, approx=True,
             approx_args={"maxiter": 10000, "tol": 1e-10},
         )
         compare_metrics(computed, net_muxviz_results["hub"], "Hub approx (vs muxViz R)",
                         rtol=0.05, atol=0.05)
 
     def test_auth_exact_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_multi_authority_centrality(net_interaction, net_n, net_l)
+        computed = versatility.compute_multi_authority_centrality(net_interaction, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["auth"], "Authority (vs muxViz R)")
 
     def test_auth_approx_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
         np.random.seed(42)
         computed = versatility.compute_multi_authority_centrality(
-            net_interaction, net_n, net_l, approx=True,
+            net_interaction, nodes=net_n, layers=net_l, approx=True,
             approx_args={"maxiter": 10000, "tol": 1e-10},
         )
         compare_metrics(computed, net_muxviz_results["auth"], "Authority approx (vs muxViz R)",
                         rtol=0.05, atol=0.05)
 
     def test_eigenvector_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_eigenvector_centrality(net_interaction, net_n, net_l)
+        computed = versatility.compute_eigenvector_centrality(net_interaction, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["eigenvector"], "Eigenvector (vs muxViz R)",
                         rtol=5e-4, atol=5e-4)
 
     # --- degree / strength reference tests --------------------------------
 
     def test_indegree_sum_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_aggregated_indegree(net_adjacency, net_n, net_l)
+        computed = versatility.compute_aggregated_indegree(net_adjacency, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["indegree"], "Indegree (vs muxViz R)")
 
     def test_outdegree_sum_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_aggregated_outdegree(net_adjacency, net_n, net_l)
+        computed = versatility.compute_aggregated_outdegree(net_adjacency, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["outdegree"], "Outdegree (vs muxViz R)")
 
     def test_instrength_sum_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_aggregated_instrength(net_interaction, net_n, net_l)
+        computed = versatility.compute_aggregated_instrength(net_interaction, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["instrength"], "Instrength (vs muxViz R)")
 
     def test_outstrength_sum_vs_muxviz(self, net_interaction, net_n, net_l, net_muxviz_results):
-        computed = versatility.compute_aggregated_outstrength(net_interaction, net_n, net_l)
+        computed = versatility.compute_aggregated_outstrength(net_interaction, nodes=net_n, layers=net_l)
         compare_metrics(computed, net_muxviz_results["outstrength"], "Outstrength (vs muxViz R)")
 
     # --- multi-degree derived tests (indegreesum - indegree = multiindegree)
@@ -617,28 +617,28 @@ def test_outstrength_sum_vs_muxviz(self, net_interaction, net_n, net_l, net_muxv
     def test_multiindegree_vs_muxviz_derived(self, net_adjacency, net_n, net_l, net_muxviz_results):
         if "indegreesum" not in net_muxviz_results:
             pytest.skip("indegreesum not in reference results")
-        computed = versatility.compute_aggregated_multiindegree(net_adjacency, net_n, net_l)
+        computed = versatility.compute_aggregated_multiindegree(net_adjacency, nodes=net_n, layers=net_l)
         expected = np.array(net_muxviz_results["indegreesum"]) - np.array(net_muxviz_results["indegree"])
         compare_metrics(computed, expected, "MultiIndegree (vs muxViz derived)")
 
     def test_multioutdegree_vs_muxviz_derived(self, net_adjacency, net_n, net_l, net_muxviz_results):
         if "outdegreesum" not in net_muxviz_results:
             pytest.skip("outdegreesum not in reference results")
-        computed = versatility.compute_aggregated_multioutdegree(net_adjacency, net_n, net_l)
+        computed = versatility.compute_aggregated_multioutdegree(net_adjacency, nodes=net_n, layers=net_l)
         expected = np.array(net_muxviz_results["outdegreesum"]) - np.array(net_muxviz_results["outdegree"])
         compare_metrics(computed, expected, "MultiOutdegree (vs muxViz derived)")
 
     def test_multiinstrength_vs_muxviz_derived(self, net_interaction, net_n, net_l, net_muxviz_results):
         if "instrengthsum" not in net_muxviz_results:
             pytest.skip("instrengthsum not in reference results")
-        computed = versatility.compute_aggregated_multiinstrength(net_interaction, net_n, net_l)
+        computed = versatility.compute_aggregated_multiinstrength(net_interaction, nodes=net_n, layers=net_l)
         expected = np.array(net_muxviz_results["instrengthsum"]) - np.array(net_muxviz_results["instrength"])
         compare_metrics(computed, expected, "MultiInstrength (vs muxViz derived)")
 
     def test_multioutstrength_vs_muxviz_derived(self, net_interaction, net_n, net_l, net_muxviz_results):
         if "outstrengthsum" not in net_muxviz_results:
             pytest.skip("outstrengthsum not in reference results")
-        computed = versatility.compute_aggregated_multioutstrength(net_interaction, net_n, net_l)
+        computed = versatility.compute_aggregated_multioutstrength(net_interaction, nodes=net_n, layers=net_l)
         expected = np.array(net_muxviz_results["outstrengthsum"]) - np.array(net_muxviz_results["outstrength"])
         compare_metrics(computed, expected, "MultiOutstrength (vs muxViz derived)")
 
@@ -647,7 +647,7 @@ def test_multioutstrength_vs_muxviz_derived(self, net_interaction, net_n, net_l,
     def test_multi_degree_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):
         """R GetMultiDegree = GetMultiInDegree + GetMultiOutDegree (directed)."""
         expected = np.array(net_muxviz_results["indegree"]) + np.array(net_muxviz_results["outdegree"])
-        computed = versatility.compute_multi_degree(net_adjacency, net_n, net_l, is_directed=True)
+        computed = versatility.compute_multi_degree(net_adjacency, nodes=net_n, layers=net_l, is_directed=True)
         compare_metrics(computed, expected, "MultiDegree (vs muxViz R derived)")
 
     def test_multi_degree_hornet_vs_muxviz(self, net_adjacency, net_n, net_l, net_muxviz_results):