基于成对二元比较的R代码优化：细菌基因相似度计算提速

Got it, let's cut through the runtime bloat here—1.6 hours for your workload is way too slow, and there are concrete, actionable tweaks you can make to slash that time down to minutes. Here's how to approach the optimization step by step:

1. Ditch Pure Python Loops for Vectorized Operations #

Processing binary strings with Python loops is a performance killer. Instead, convert your binary sequences into numerical arrays (using NumPy) and leverage vectorized operations—these run under the hood in optimized C code, so they're orders of magnitude faster.

• Convert binary strings to numerical arrays: Turn each "1/0" string into a boolean or uint8 NumPy array, or even better, pack the binary data into integers (since bitwise operations are lightning-fast). For example:

  import numpy as np
  
  def str_to_bin_array(bin_str):
      # Convert string to a uint8 array of 0s and 1s
      return np.array(list(bin_str), dtype=np.uint8)
  
  # Or pack into a single integer (great for short sequences)
  def str_to_int(bin_str):
      return int(bin_str, 2)
  

• Calculate similarity/dissimilarity with vectorized functions: For Hamming distance (count of differing bits, perfect for your binary data), use NumPy's bitwise operations instead of manual string comparisons:

  # For two binary arrays a and b
  hamming_distance = np.sum(np.bitwise_xor(a, b))
  # For two integer-packed sequences
  hamming_distance = bin(a ^ b).count('1')
  

2. Precompute & Cache All Model Data #

Stop re-parsing files every time you need to compare models. Load all 250 files into a single in-memory NumPy matrix first—your total dataset size is tiny (250×450=112,500 models; even with 10,000-bit sequences, that's ~130MB of memory) so this is totally feasible.

Example loading script:

import os
import numpy as np

def load_all_models(folder_path):
    all_models = []
    for filename in os.listdir(folder_path):
        with open(os.path.join(folder_path, filename), 'r') as f:
            # Read each line, convert to binary array, and append
            for line in f:
                bin_str = line.strip()
                all_models.append(str_to_bin_array(bin_str))
    # Convert list of arrays to a 2D NumPy matrix
    return np.vstack(all_models)

3. Cut Redundant Calculations in Half #

Your dissimilarity matrix is symmetric (distance between model A and B is the same as B and A) and has zeros on the diagonal. Only compute the upper or lower triangle of the matrix, then mirror it to fill the rest—this cuts your total number of calculations from n² to n(n-1)/2.

Using NumPy to handle this:

n = len(all_models)
diss_matrix = np.zeros((n, n), dtype=np.float32)

# Get indices for the upper triangle (excluding diagonal)
rows, cols = np.triu_indices(n, k=1)

for i, j in zip(rows, cols):
    diss_matrix[i, j] = np.sum(np.bitwise_xor(all_models[i], all_models[j]))
    # Mirror to lower triangle
    diss_matrix[j, i] = diss_matrix[i, j]

4. Use Specialized Libraries for Pairwise Distances #

Don't reinvent the wheel—libraries like SciPy have highly optimized functions for exactly this task. scipy.spatial.distance.pdist computes pairwise distances in C, and it's way faster than any Python loop you can write.

For binary data, the hamming metric gives you the proportion of differing bits (multiply by sequence length to get raw count if needed):

from scipy.spatial.distance import pdist, squareform

# Assume all_models is a 2D NumPy array of 0s and 1s
# Compute pairwise Hamming distances (returns a condensed array)
condensed_dists = pdist(all_models, metric='hamming')
# Convert to full square dissimilarity matrix
diss_matrix = squareform(condensed_dists)

# If you want raw bit count instead of proportion:
seq_length = all_models.shape[1]
diss_matrix = squareform(condensed_dists) * seq_length

This alone can reduce your runtime from hours to minutes.

5. Parallelize the Workload (If Needed) #

Even with vectorized operations, for extremely large datasets, you can split the work across CPU cores using multiprocessing. SciPy's pdist already uses some internal optimizations, but if you need extra speed, you can split the matrix chunks and compute them in parallel with concurrent.futures or multiprocessing.Pool.

Example with ProcessPoolExecutor:

from concurrent.futures import ProcessPoolExecutor

def compute_chunk(chunk_indices):
    chunk_dists = []
    for i, j in chunk_indices:
        chunk_dists.append(np.sum(np.bitwise_xor(all_models[i], all_models[j])))
    return chunk_dists

# Split upper triangle indices into chunks for each core
rows, cols = np.triu_indices(n, k=1)
chunk_size = len(rows) // os.cpu_count()
chunks = [list(zip(rows[i:i+chunk_size], cols[i:i+chunk_size])) for i in range(0, len(rows), chunk_size)]

with ProcessPoolExecutor() as executor:
    results = executor.map(compute_chunk, chunks)

# Flatten results and populate the matrix
flat_results = []
for res in results:
    flat_results.extend(res)
diss_matrix[rows, cols] = flat_results
diss_matrix[cols, rows] = flat_results

6. Optimize File I/O #

If reading files is a bottleneck, use faster I/O methods instead of manual line-by-line reading. For example, use Pandas to read entire files at once:

import pandas as pd

def load_file_fast(filename):
    # Read all lines into a Series, convert each to binary array
    df = pd.read_csv(filename, header=None, dtype=str)
    return df[0].apply(lambda x: np.array(list(x), dtype=np.uint8)).to_numpy()

Putting all these together—loading everything into NumPy, using SciPy's pdist, and leveraging symmetric matrix optimizations—you should see a massive drop in runtime. I've used these exact techniques to cut similar workloads from hours to under 10 minutes.

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