C代码减法运算异常与ti变量计算错误问题排查求助

问题现象 #

在执行pix0与pix1的减法运算、ti变量计算时出现异常：

• pix0和pix1为常量，二者差值应保持恒定，但输出中subst逐渐减小
• ti变量应按固定步长递增，但输出中ti的步长越来越大

输出结果 #

pix0 = 0.000000
pix1 = 10.000000
subst = 10.000000
ti = 0.000000
subst = 10.000000
ti = 0.010000
subst = 9.990000
ti = 0.030000
subst = 9.970000
ti = 0.060000
subst = 9.940000
ti = 0.100000
subst = 9.900000
ti = 0.150000
subst = 9.850000
ti = 0.210000
subst = 9.790000
ti = 0.280000
subst = 9.720000
ti = 0.360000
subst = 9.640000
ti = 0.450000
subst = 9.550000

代码片段 #

C代码 #

#include <stdio.h>
#include <stdlib.h> 
#include <math.h>
#include <gsl/gsl_odeiv2.h>
#include <gsl/gsl_errno.h>

typedef struct {
    double deg_m, deg_p, alpha, alpha_0, beta, n;
} repressilator_params;

int repressilator_func(double t, const double y[], double f[], void *params) {
    (void)(t);
    repressilator_params *p = (repressilator_params *)params;
    double deg_m = p->deg_m;
    double deg_p = p->deg_p;
    double alpha = p->alpha;
    double alpha_0 = p->alpha_0;
    double beta = p->beta;
    double n = p->n;

    f[0] = -deg_m * y[0] + alpha / (1 + pow(y[5], n)) + alpha_0; 
    f[1] = -deg_m * y[1] + alpha / (1 + pow(y[3], n)) + alpha_0; 
    f[2] = -deg_m * y[2] + alpha / (1 + pow(y[4], n)) + alpha_0; 

    f[3] = -deg_p * y[3] + beta * y[0]; 
    f[4] = -deg_p * y[4] + beta * y[1];
    f[5] = -deg_p * y[5] + beta * y[2];

    return GSL_SUCCESS;
}

double* simulate_repressilator(double *params_array, double *initial_conditions, double *time_span, double num_points) {
    gsl_odeiv2_system sys = {repressilator_func, NULL, 6, params_array};

    double pix0 = time_span[0];
    double pix1 = time_span[1];
    printf("pix0 = %f\n", pix0);
    printf("pix1 = %f\n", pix1);
    double h = 1e-6;
    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, h, 1e-8, 0.0);

    double* results = malloc(num_points * 7 * sizeof(double));
    if (!results) {
        fprintf(stderr, "Failed to allocate memory for results\n");
        return NULL;
    }

    int i;
    double ti;
    int numi = (int) num_points;
    for (i = 0; i < numi; i=i+1) {
        double subst = (pix1 - pix0);
        printf("subst = %f\n", subst);
        ti = pix0 + (10 * ((double) i) / (num_points));
        printf("ti = %f\n", ti);
        
        int status = gsl_odeiv2_driver_apply(d, &pix0, ti, initial_conditions);

        if (status != GSL_SUCCESS) {
            printf("Error, return value=%d\n", status);
            free(results);
            return NULL;
        }

        results[i * 7 + 0] = ti;
        results[i * 7 + 1] = initial_conditions[0];
        results[i * 7 + 2] = initial_conditions[1];
        results[i * 7 + 3] = initial_conditions[2];
        results[i * 7 + 4] = initial_conditions[3];
        results[i * 7 + 5] = initial_conditions[4];
        results[i * 7 + 6] = initial_conditions[5];
    }

    gsl_odeiv2_driver_free(d);

    return results;
}

编译命令 #

gcc -shared -o libRepressilator.so -fPIC repressilator.c -lgsl -lgslcblas -lm

Python/ctypes调用代码 #

import ctypes
import numpy as np

c_lib = ctypes.CDLL('./libRepressilator.so')

c_lib.simulate_repressilator.argtypes = [np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags="C_CONTIGUOUS"),  # params_array
                                         np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags="C_CONTIGUOUS"),  # initial_conditions
                                         np.ctypeslib.ndpointer(dtype=np.float64, ndim=1, flags="C_CONTIGUOUS"),  # time_span
                                         ctypes.c_double]  # num_points
c_lib.simulate_repressilator.restype = ctypes.POINTER(ctypes.c_double * 7)


params_array = np.array([1, 0.1, 100, 0.01, 10, 1], dtype=np.float64)  # deg_m, deg_p,alpha, alpha_0, beta, n
initial_conditions = np.array([0.0, 1.0, 0.0, 2.0, 0.0, 3.0], dtype=np.float64)  # m1, m2, m3, p1, p2, p3
time_span = np.array([0, 10], dtype=np.float64)  # Start time, end time
num_points = 1000

result_ptr = c_lib.simulate_repressilator(params_array, initial_conditions, time_span, num_points)


results = np.ctypeslib.as_array(result_ptr.contents, shape=(num_points, 7))

c_lib.free(result_ptr)

问题根源 #

1. pix0被意外修改：gsl_odeiv2_driver_apply函数的第二个参数是当前时间指针，函数执行后会将该指针指向的变量更新为积分完成后的时间。代码中直接传入&pix0，导致每次循环后pix0不再是初始起始时间，而是当前积分结束时间，这就解释了subst = pix1 - pix0的值逐渐减小。

2. ti的计算逻辑错误：代码中用pix0 + (10 * ((double) i) / (num_points))计算目标时间，一是依赖了被修改的pix0，二是硬编码10作为时间跨度，既不符合通用逻辑，也导致ti的步长越来越大。

修复方案 #

修改simulate_repressilator函数，单独维护当前积分时间变量，同时用初始起止时间计算目标时间点：

double* simulate_repressilator(double *params_array, double *initial_conditions, double *time_span, double num_points) {
    gsl_odeiv2_system sys = {repressilator_func, NULL, 6, params_array};

    double t_start = time_span[0];
    double t_end = time_span[1];
    printf("t_start = %f\n", t_start);
    printf("t_end = %f\n", t_end);
    double h = 1e-6;
    gsl_odeiv2_driver *d = gsl_odeiv2_driver_alloc_y_new(&sys, gsl_odeiv2_step_rk8pd, h, 1e-8, 0.0);

    double* results = malloc(num_points * 7 * sizeof(double));
    if (!results) {
        fprintf(stderr, "Failed to allocate memory for results\n");
        return NULL;
    }

    int i;
    double ti;
    int numi = (int) num_points;
    double t_current = t_start; // 单独维护当前积分时间，不修改初始起始时间
    for (i = 0; i < numi; i=i+1) {
        double subst = (t_end - t_start); // 用初始起止时间计算差值，保持恒定
        printf("subst = %f\n", subst);
        ti = t_start + (t_end - t_start) * ((double)i)/num_points; // 均匀生成目标时间点
        printf("ti = %f\n", ti);
        
        int status = gsl_odeiv2_driver_apply(d, &t_current, ti, initial_conditions);

        if (status != GSL_SUCCESS) {
            printf("Error, return value=%d\n", status);
            free(results);
            return NULL;
        }

        results[i * 7 + 0] = ti;
        results[i * 7 + 1] = initial_conditions[0];
        results[i * 7 + 2] = initial_conditions[1];
        results[i * 7 + 3] = initial_conditions[2];
        results[i * 7 + 4] = initial_conditions[3];
        results[i * 7 + 5] = initial_conditions[4];
        results[i * 7 + 6] = initial_conditions[5];
    }

    gsl_odeiv2_driver_free(d);

    return results;
}

修复说明 #

• 用t_start和t_end保存初始起止时间，确保subst始终为恒定值
• 用t_current跟踪当前积分时间点，避免修改初始起始时间
• ti的计算改为基于初始时间跨度的均匀分布，步长固定为(t_end - t_start)/num_points，确保每次递增步长一致

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