Numba/NumbaLSODA中如何使用eval()或传入代码生成的ydot数组

Numba的cfunc要求代码必须能在nopython模式下编译为静态机器码，eval是运行时动态解析字符串的内置函数，Numba无法对其做类型推断和编译，因此完全不可能在@cfunc修饰的导数函数中使用eval，无需尝试适配。

以下两种方案都可以保留你动态生成计算逻辑的需求，同时完全兼容Numba的编译要求，性能和硬编码的导数函数一致。

方案1：提前生成完整导数函数代码再编译 #

你已经有生成每个ydot[i]表达式字符串的逻辑，可以直接把这些表达式拼接为完整的导数函数代码字符串，再通过exec执行代码完成Numba编译，完全不需要在函数内部调用eval。
示例代码如下：

from numba import cfunc
from numbaLSODA import lsoda_sig

# 你已有的逻辑生成的每个ydot项的表达式列表
expr_list = [
    "coeff[0]*y[0]*y[1] + coeff[1]*y[2] - coeff[2]*y[3]",
    "coeff[3]*y[0] - coeff[4]*y[1]**2",
    # 剩下的所有ydot表达式依次填入
]

# 拼接生成完整的导数函数代码字符串
func_code_str = """
@cfunc(lsoda_sig)
def deriv(t, y, ydot, coeff):
"""
for idx, expr in enumerate(expr_list):
    func_code_str += f"    ydot[{idx}] = {expr}\n"

# 执行代码字符串，完成编译得到可用的deriv函数
local_namespace = {}
exec(func_code_str, globals(), local_namespace)
deriv = local_namespace["deriv"]

该方案完全复用你现有的表达式生成逻辑，不需要额外改造成本。

方案2：结构化存储计算规则，循环计算 #

如果你的表达式有统一的结构（比如常见的反应动力学模型，每个项都是系数乘部分状态量的乘积），可以把计算规则预存在数组中，在导数函数内用循环完成计算。Numba对循环的优化能力极强，性能不会有损失。
示例代码如下：

import numpy as np
from numba import cfunc
from numbaLSODA import lsoda_sig

# 预定义计算规则数组，每行对应一个计算项，格式为：
# [目标ydot下标, coeff数组下标, 第一个y的下标, 第二个y的下标]
# 第二个y的下标填-1表示该为单状态量项，正负号直接存在coeff数组中即可
term_rules = np.array([
    [0, 0, 0, 1],  # 对应 ydot[0] += coeff[0] * y[0] * y[1]
    [0, 1, 2, -1], # 对应 ydot[0] += coeff[1] * y[2]
    [0, 2, 3, -1], # 对应 ydot[0] -= coeff[2] * y[3]（负号存在coeff[2]中）
    # 所有计算规则依次填入
], dtype=np.int64)
N_TERMS = term_rules.shape[0]
N_STATES = 5 # 替换为你的状态量总数

@cfunc(lsoda_sig)
def deriv(t, y, ydot, coeff):
    # 初始化ydot为0
    for i in range(N_STATES):
        ydot[i] = 0.0
    # 遍历所有规则计算累加
    for k in range(N_TERMS):
        target_idx = term_rules[k, 0]
        c_idx = term_rules[k, 1]
        y1_idx = term_rules[k, 2]
        y2_idx = term_rules[k, 3]
        if y2_idx < 0:
            term_val = coeff[c_idx] * y[y1_idx]
        else:
            term_val = coeff[c_idx] * y[y1_idx] * y[y2_idx]
        ydot[target_idx] += term_val

该方案不需要动态生成代码，仅需修改term_rules即可调整计算逻辑，适合规则统一的场景。

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