求助：Ising模型蒙特卡洛模拟的Cython优化，Python与Cython性能一致

一、Cython核心优化点 #

你的Cython代码性能提升有限，核心原因是大量关键计算仍在调用Python级别的numpy函数，没有真正利用C的底层效率。以下是针对性优化：

1. 替换Python随机数生成为C级实现

每次循环调用np.random.randint、np.random.random会带来巨大的Python调用开销，改用C标准库的随机函数：

# 在文件顶部添加cimport
from libc.stdlib cimport rand, srand, RAND_MAX
from libc.math cimport log
from libc.time cimport time

# 初始化随机种子（在主函数调用一次即可）
srand(time(NULL))

# 替换随机数生成逻辑：
# 生成0到size-1的随机整数
n = rand() % size
# 生成0-1的随机浮点数
double rand_val = <double>rand() / RAND_MAX

2. 手动计算邻居自旋和，替代np.add.reduce

np.add.reduce(spins[neighbors_list[n]])是Python级别的函数调用，改为C级循环计算：

# 原代码：
delta_E = 2 * spin * np.add.reduce(spins[neighbors_list[n]])

# 优化后：
cdef int neighbor_sum = 0
cdef int idx
for idx in range(4):
    neighbor_sum += spins[neighbors_list[n][idx]]
delta_E = 2 * spin * neighbor_sum

3. 替换np.log为C标准库函数

把np.log(np.random.random(1)[0])换成C的log函数，避免Python调用：

# 原代码：
if delta_E < 0 or -T*np.log(np.random.random(1)[0]) > delta_E:

# 优化后（结合上面的rand_val）：
if delta_E < 0 or (-T * log(rand_val)) > delta_E:

4. 优化自旋数组的定义与访问

将spins定义为memoryview或C数组，进一步降低访问开销：

# 改用memoryview（兼容numpy且访问更快）
cdef int[:] spins = np.random.choice([-1, 1], size=size).astype(np.int32)

# 或者直接用C数组（极致性能，需手动管理内存）
cdef int *spins = <int*>malloc(size * sizeof(int))
if not spins:
    raise MemoryError("Failed to allocate memory")
# 初始化自旋
for i in range(size):
    spins[i] = 1 if (rand() % 2 == 0) else -1
# 函数结束前释放内存
free(spins)

5. 编译选项优化

在编译时开启最高级优化：

# setup.py示例
from setuptools import setup
from Cython.Build import cythonize
import numpy as np

setup(
    ext_modules=cythonize("ising.pyx", compiler_directives={"language_level": "3"}),
    include_dirs=[np.get_include()],
    extra_compile_args=["-O3", "-ffast-math", "-march=native"]
)

二、非Cython的加速方案 #

如果不想折腾Cython，这些方法也能大幅提升性能：

1. Numba即时编译

Numba可以直接对Python函数进行JIT编译，无需手动写Cython语法，对循环密集型代码效果极佳：

import numba
from numba import njit
import numpy as np

@njit(fastmath=True)
def MH_single_flip_numba(neighbors_list, T, iterations, size):
    spins = np.random.choice([-1, 1], size=size)
    for step in range(iterations):
        n = np.random.randint(0, size)
        spin = spins[n]
        neighbor_sum = spins[neighbors_list[n]].sum()
        delta_E = 2 * spin * neighbor_sum
        if delta_E < 0 or -T*np.log(np.random.random()) > delta_E:
            spins[n] = -spin
    return spins

2. 预生成随机数

把所有需要的随机数一次性预生成，避免循环内重复调用随机数函数：

# 预生成所有模拟步骤的随机索引和浮点数
total_steps = len(T_list) * 5000
rand_indices = np.random.randint(0, size, total_steps)
rand_vals = np.random.random(total_steps)

# 在循环中直接取用预生成的随机数
idx = 0
for T in T_list:
    spins = np.random.choice([-1,1], size=size)
    for _ in range(5000):
        n = rand_indices[idx]
        spin = spins[n]
        delta_E = 2 * spin * spins[neighbors_list[n]].sum()
        if delta_E <0 or -T*np.log(rand_vals[idx])>delta_E:
            spins[n] = -spin
        idx +=1

3. 并行化处理温度循环

不同温度的模拟是独立的，可以用并行库同时处理：

from joblib import Parallel, delayed

def run_simulation(T):
    return MH_single_flip_numba(neighbors_list, T, 5000, size)

# 用所有CPU核心并行处理
results = Parallel(n_jobs=-1)(delayed(run_simulation)(T) for T in T_list)

三、优化效果验证 #

优化后，代码执行时间应该能降低一个数量级以上。建议先尝试Numba方案，开发成本最低且性能提升显著；如果追求极致性能，再用Cython的深度优化方案。

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