求Python中可批量采样q态Potts模型的MCMC代码（用于蛋白生成）

针对你需要从q=21态Potts模型（对应20种氨基酸+空位）生成MCMC样本、支持批量模型处理及频繁调用的需求，以下是推荐的工具与实现方案：

1. pyDCA：原生DCA/Potts模型采样工具 #

pyDCA是专门用于直接耦合分析（DCA）的Python库，内置针对q态Potts模型的MCMC采样模块，可直接接收fields和couplings参数，适配频繁调用场景。

示例代码 #

import pydca
from pydca.mcmc_sampler import mcmc_sampler

def get_samples(fields, couplings, N_samples):
    # 初始化采样器，指定q=21对应蛋白序列的状态数
    sampler = mcmc_sampler.MCMCSampler(
        fields=fields,
        couplings=couplings,
        q=21,
        seq_len=fields.shape[0]
    )
    # 执行烧录（消除初始状态影响）
    sampler.burn_in(n_burn=1000)
    # 生成指定数量的样本
    samples = sampler.sample(n_samples=N_samples)
    return samples

核心优势 #

• 原生支持DCA训练输出的fields和couplings格式，无需额外转换
• 采样逻辑经过优化，效率接近bmDCA且支持单次调用生成多样本
• 可通过循环批量处理多组模型参数，适配训练流程中的频繁调用

2. 轻量自定义MCMC采样器 #

若现有库无法满足性能或灵活性需求，可自行实现Metropolis-Hastings采样逻辑，无第三方依赖，适合频繁调用场景（可避免库的额外开销）。

示例代码（优化版） #

import numpy as np

def get_samples(fields, couplings, N_samples):
    N = fields.shape[0]
    q = 21
    samples = np.zeros((N_samples, N), dtype=int)
    current_seq = np.random.randint(0, q, size=N)
    
    # 预计算耦合矩阵的对称项，优化能量差计算效率
    sym_couplings = couplings + couplings.transpose(1,0,3,2)
    
    for step in range(1000 + N_samples):
        # 随机选择突变位置
        pos = np.random.randint(N)
        old_state = current_seq[pos]
        new_state = np.random.randint(0, q)
        if new_state == old_state:
            continue
        
        # 仅计算突变位置相关的能量差（避免全序列计算，提升效率）
        delta_E = fields[pos, new_state] - fields[pos, old_state]
        delta_E += np.sum(sym_couplings[pos, :, new_state, current_seq])
        
        # Metropolis接受准则
        if delta_E >= 0 or np.random.rand() < np.exp(delta_E):
            current_seq[pos] = new_state
        
        # 烧录完成后保存样本
        if step >= 1000:
            samples[step - 1000] = current_seq.copy()
    
    return samples

核心优势 #

• 完全自定义，可针对蛋白序列的q=21态场景做针对性优化
• 无第三方依赖，调用开销极低，适合训练流程中的高频调用
• 可结合multiprocessing实现多组模型的并行批量处理

3. GPU加速采样（PyTorch/TensorFlow） #

若需处理大量模型或生成海量样本，可使用深度学习框架实现GPU加速的并行采样，大幅提升批量处理效率。

PyTorch示例代码 #

import torch

def get_samples(fields, couplings, N_samples, device='cuda'):
    # 将参数转移到GPU
    fields = torch.tensor(fields, dtype=torch.float32, device=device)
    couplings = torch.tensor(couplings, dtype=torch.float32, device=device)
    N = fields.shape[0]
    q = 21
    
    samples = torch.zeros((N_samples, N), dtype=torch.int32, device=device)
    current_seq = torch.randint(0, q, (N,), device=device)
    # 预计算对称耦合矩阵
    sym_couplings = couplings + couplings.permute(1,0,3,2)
    
    for step in range(1000 + N_samples):
        pos = torch.randint(0, N, (1,), device=device).item()
        old_state = current_seq[pos]
        new_state = torch.randint(0, q, (1,), device=device).item()
        if new_state == old_state:
            continue
        
        delta_E = fields[pos, new_state] - fields[pos, old_state]
        delta_E += torch.sum(sym_couplings[pos, :, new_state, current_seq])
        
        accept_prob = torch.min(torch.tensor(1.0, device=device), torch.exp(delta_E))
        if torch.rand(1, device=device) < accept_prob:
            current_seq[pos] = new_state
        
        if step >= 1000:
            samples[step - 1000] = current_seq.clone()
    
    return samples.cpu().numpy()

核心优势 #

• GPU加速可同时处理多组模型参数或生成海量样本
• 适合大规模批量处理场景，效率远高于CPU实现
• 可结合DataLoader实现模型参数的批量输入与采样

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