无需大型二维数组求平均尺寸分布及随机矩阵聚类标记技术咨询

Hey there! Let's break down your two problems with practical, efficient solutions—these are common tasks in random lattice analysis, so I've got you covered.

1. Generating Random 0-1 Matrices and Labeling Connected Clusters #

We have two reliable approaches here: using scientific computing libraries for quick implementation, or manually implementing the Union-Find (Disjoint Set Union, DSU) algorithm for full control over connectivity rules.

Using SciPy & NumPy (Quick, Minimal Code) #

NumPy lets us generate the random matrix in one line, and SciPy's ndimage.label handles connected component labeling out of the box. You can choose between 4-neighbor (up/down/left/right) or 8-neighbor (including diagonals) connectivity.

import numpy as np
from scipy.ndimage import label

def generate_and_label_cluster(L, p, connectivity=8):
    # Generate 0-1 matrix where 1s have probability p
    matrix = np.random.binomial(n=1, p=p, size=(L, L))
    # Define connectivity structure: 3x3 for 8-neighbor, cross for 4-neighbor
    struct = np.ones((3,3)) if connectivity ==8 else np.array([[0,1,0],[1,1,1],[0,1,0]])
    # Label clusters and count total clusters
    labeled_matrix, num_clusters = label(matrix, structure=struct)
    return matrix, labeled_matrix, num_clusters

# Example usage
L = 20
p = 0.3
matrix, labeled_clusters, cluster_count = generate_and_label_cluster(L, p)
print(f"Generated {L}x{L} matrix with {cluster_count} unique clusters")

Manual Union-Find Implementation (Customizable) #

If you want to avoid external dependencies or tweak the connectivity logic, the DSU algorithm is perfect. It uses path compression and union-by-size to ensure fast cluster merging and lookup.

import numpy as np

def generate_and_label_dsu(L, p, connectivity=8):
    # Generate random 0-1 matrix
    matrix = np.random.binomial(1, p, (L, L))
    # DSU structures: parent array maps index to root, size array tracks cluster size
    parent = np.arange(L*L)
    size = np.ones(L*L, dtype=int)

    def find(u):
        # Path compression to speed up future lookups
        while parent[u] != u:
            parent[u] = parent[parent[u]]
            u = parent[u]
        return u

    def union(u, v):
        # Union by size: attach smaller cluster to larger one's root
        u_root = find(u)
        v_root = find(v)
        if u_root == v_root:
            return
        if size[u_root] < size[v_root]:
            u_root, v_root = v_root, u_root
        parent[v_root] = u_root
        size[u_root] += size[v_root]

    # Flatten (i,j) matrix index to 1D: i*L + j
    for i in range(L):
        for j in range(L):
            if matrix[i,j] == 1:
                current_idx = i*L + j
                # Check left neighbor (current row)
                if j > 0 and matrix[i,j-1] == 1:
                    union(current_idx, i*L + (j-1))
                # Check upper neighbor (previous row)
                if i > 0 and matrix[i-1,j] == 1:
                    union(current_idx, (i-1)*L + j)
                # Add diagonal checks for 8-neighbor connectivity
                if connectivity ==8:
                    if i>0 and j>0 and matrix[i-1,j-1]==1:
                        union(current_idx, (i-1)*L + (j-1))
                    if i>0 and j<L-1 and matrix[i-1,j+1]==1:
                        union(current_idx, (i-1)*L + (j+1))

    # Create labeled matrix from DSU roots
    labeled_matrix = np.zeros((L, L), dtype=int)
    cluster_id = 1
    root_to_id = {}
    for i in range(L):
        for j in range(L):
            if matrix[i,j] ==1:
                root = find(i*L +j)
                if root not in root_to_id:
                    root_to_id[root] = cluster_id
                    cluster_id +=1
                labeled_matrix[i,j] = root_to_id[root]
    return matrix, labeled_matrix, cluster_id-1, size

# Example usage
L = 20
p = 0.3
matrix, labeled_clusters, cluster_count, cluster_sizes = generate_and_label_dsu(L, p)
print(f"Generated {L}x{L} matrix with {cluster_count} unique clusters")

2. Efficient Average Cluster Size Distribution Without Large 2D Arrays #

When working with large L, storing a full L×L matrix eats up too much memory. Instead, we use memory-efficient methods that process data incrementally or focus only on non-zero elements.

Approach 1: Row-by-Row Incremental DSU #

This method only stores the current and previous rows' cluster data, keeping memory usage at O(L) instead of O(L²). We update cluster information as we process each row, and track size distributions on the fly.

import numpy as np
from collections import defaultdict

def compute_avg_cluster_distribution(L, p, num_trials=100, connectivity=4):
    # Track total cluster counts per size across all trials
    size_distribution = defaultdict(int)

    for _ in range(num_trials):
        prev_row_clusters = {}  # Maps column index to cluster root ID for the previous row
        current_dsu = {}        # DSU for current row clusters
        current_size = {}       # Tracks size of each cluster in current DSU

        def find(u):
            while current_dsu[u] != u:
                current_dsu[u] = current_dsu[current_dsu[u]]
                u = current_dsu[u]
            return u

        def union(u, v):
            u_root = find(u)
            v_root = find(v)
            if u_root == v_root:
                return
            if current_size[u_root] < current_size[v_root]:
                u_root, v_root = v_root, u_root
            current_dsu[v_root] = u_root
            current_size[u_root] += current_size[v_root]

        # Define neighbor offsets based on connectivity
        offsets = [(-1,0), (0,-1)]  # 4-neighbor: check upper and left
        if connectivity ==8:
            offsets.extend([(-1,-1), (-1,1)])  # Add diagonals for 8-neighbor

        for i in range(L):
            current_row = np.random.binomial(1, p, L)
            current_row_clusters = {}

            for j in range(L):
                if current_row[j] ==1:
                    # Create new cluster for current element
                    cluster_id = f"row_{i}_col_{j}"
                    current_dsu[cluster_id] = cluster_id
                    current_size[cluster_id] =1

                    # Merge with left neighbor in current row
                    if j>0 and current_row[j-1]==1:
                        left_id = current_row_clusters[j-1]
                        union(cluster_id, left_id)

                    # Merge with relevant neighbors from previous row
                    if i>0:
                        for di, dj in offsets:
                            ni, nj = i+di, j+dj
                            if 0<=ni<L and 0<=nj<L and nj in prev_row_clusters:
                                neighbor_id = prev_row_clusters[nj]
                                union(cluster_id, neighbor_id)

                    # Map current column to its cluster root
                    current_row_clusters[j] = find(cluster_id)

            # Add finished clusters (not connected to current row) to distribution
            for root in prev_row_clusters.values():
                if root not in current_dsu:
                    size_distribution[current_size[root]] +=1

            # Track unique clusters in current row
            seen_roots = set()
            for root in current_row_clusters.values():
                if root not in seen_roots:
                    seen_roots.add(root)
                    size_distribution[current_size[root]] +=1

            # Prepare for next row: update previous row cluster data
            prev_row_clusters = {j: find(root) for j, root in current_row_clusters.items()}
            # Clean up to save memory
            current_dsu.clear()
            current_size.clear()

        # Add remaining clusters from the last row
        seen_roots = set()
        for root in prev_row_clusters.values():
            if root not in seen_roots:
                seen_roots.add(root)
                size_distribution[current_size.get(root,1)] +=1

    # Compute average distribution by dividing by number of trials
    avg_dist = {size: count/num_trials for size, count in size_distribution.items()}
    return dict(sorted(avg_dist.items()))

# Example usage
L = 100
p = 0.3
avg_dist = compute_avg_cluster_distribution(L, p, num_trials=10)
print("Average Cluster Size Distribution:")
for size, avg_count in avg_dist.items():
    print(f"Size {size}: {avg_count:.2f} clusters per trial")

Approach 2: Sparse Event-Driven Monte Carlo #

For low values of p (most elements are 0), we only generate and process coordinates of 1s. This reduces memory usage to O(N) where N ≈ p*L²—far smaller than O(L²) for large L.

import numpy as np
from collections import defaultdict

def sparse_cluster_distribution(L, p, num_trials=100, connectivity=4):
    size_distribution = defaultdict(int)

    for _ in range(num_trials):
        # Generate all coordinates where element is 1
        total_elements = L*L
        num_ones = np.random.binomial(total_elements, p)
        # Sample unique indices and convert to (i,j) coordinates
        indices = np.random.choice(total_elements, num_ones, replace=False)
        ones_coords = [(idx//L, idx%L) for idx in indices]
        ones_coords.sort()  # Process in row-major order

        # DSU setup for sparse coordinates
        parent = {}
        size = {}

        def find(u):
            while parent[u] != u:
                parent[u] = parent[parent[u]]
                u = parent[u]
            return u

        def union(u, v):
            u_root = find(u)
            v_root = find(v)
            if u_root == v_root:
                return
            if size[u_root] < size[v_root]:
                u_root, v_root = v_root, u_root
            parent[v_root] = u_root
            size[u_root] += size[v_root]

        # Define neighbor offsets
        offsets = [(-1,0), (0,-1)]
        if connectivity ==8:
            offsets.extend([(-1,-1), (-1,1)])

        for (i,j) in ones_coords:
            parent[(i,j)] = (i,j)
            size[(i,j)] =1
            # Check all relevant neighbors that have already been processed
            for di, dj in offsets:
                ni, nj = i+di, j+dj
                if 0<=ni<L and 0<=nj<L and (ni, nj) in parent:
                    union((i,j), (ni, nj))

        # Count unique cluster sizes
        seen_roots = set()
        for coord in ones_coords:
            root = find(coord)
            if root not in seen_roots:
                seen_roots.add(root)
                size_distribution[size[root]] +=1

    # Compute average distribution
    avg_dist = {size: count/num_trials for size, count in size_distribution.items()}
    return dict(sorted(avg_dist.items()))

# Example usage
L = 200
p = 0.1
avg_dist = sparse_cluster_distribution(L, p, num_trials=50)
print("Average Sparse Cluster Size Distribution:")
for size, avg_count in avg_dist.items():
    print(f"Size {size}: {avg_count:.2f} clusters per trial")

内容的提问来源于stack exchange，提问作者user9713570