多并行Python进程调用C DLL时返回值异常问题排查

• Python模拟程序通过ctypes调用编译后的C DLL，使用自定义类和ctypes完成数据传递与接收。
• C代码单次执行100,000次随机过程迭代，每个过程需15000次迭代收敛；Python代码循环9次调用该C代码，整体运行在multiprocessing.Pool(processes=12)中，共24次迭代。
• 单实例运行无异常，但同时打开多个控制台并行运行不同模型脚本时，偶尔出现返回值损坏：在9×100,000次迭代（每次返回7个变量）中随机出现错误值（上次检测约2900个），例如本该为a + b的值仅返回b，仿佛a意外变为0（C代码逻辑上不可能发生），疑似内存覆盖或指针问题。
• 目前仅发现某一个值异常，可能存在未被察觉的其他错误，需定位问题根源以确保模拟结果可信。

# -*- coding: utf-8 -*-
import ctypes
import re
import os
import multiprocessing
import numpy as np
import pandas as pd
from scipy.stats import norm

PATH_ = re.match(r"^(.*?fits)", os.path.realpath(__file__)).group(1)
os.chdir(PATH_ + '/GCD_rTru')

myso = ctypes.cdll.LoadLibrary(f'../models/GCDM.dll')
algo = myso.GCDM

class params:
    def __init__(self, x0, v, r, g, Te, Tr, xi, leak, u=0, n_sim_trials=100000):
        self.x0 = x0
        self.v = v
        self.r = r
        self.g = g
        self.Te = Te
        self.Tr = Tr
        self.xi = xi
        self.leak = leak
        self.u = u
        self.drift_guess = 0
        self.sx0 = 0
        self.sv = 0
        self.sTe = 0
        self.sTr = 0
        self.sLeak = 0
        self.s = .1
        self.dt = .05
        self.n = n_sim_trials

        self.resp = np.zeros(self.n)
        self.RT = np.zeros(self.n)
        self.PMT = np.zeros(self.n)
        self.MT = np.zeros(self.n)
        self.firstHit = np.zeros(self.n)
        self.firstHitLoc = np.zeros(self.n)
        self.coactiv = np.zeros(self.n)

        self.maxiter = 15000
        self.randomTable = norm.ppf(np.arange(.0001, 1, .0001))
        self.rangeLow = 0
        self.rangeHigh = len(self.randomTable)

def run_sx(s):
    cls = []
    par = params_df.loc[s]
    for f in [200, 1000, 2000]:
        for mc in [.02, .10, .40]:
            sim = params( x0=0, v=par['k']*mc, r=par[f'r_{f}'], g=par['g'], Te=par['Te'], Tr=par[f'Tr_{f}'], xi=par['xi'], leak=par[f'leak'], u=par[f'u_{f}'] )

            algo(
                ctypes.c_double(sim.x0), ctypes.c_double(sim.v), ctypes.c_double(sim.r), ctypes.c_double(sim.g), ctypes.c_double(sim.Te), ctypes.c_double(sim.Tr),
                ctypes.c_double(sim.xi), ctypes.c_double(sim.leak), ctypes.c_double(sim.u), ctypes.c_double(sim.drift_guess),
                ctypes.c_double(sim.sx0), ctypes.c_double(sim.sv), ctypes.c_double(sim.sTe), ctypes.c_double(sim.sTr), ctypes.c_double(sim.sLeak),
                ctypes.c_double(sim.s), ctypes.c_double(sim.dt),
                ctypes.c_void_p(sim.resp.ctypes.data), ctypes.c_void_p(sim.RT.ctypes.data),
                ctypes.c_void_p(sim.PMT.ctypes.data), ctypes.c_void_p(sim.MT.ctypes.data),
                ctypes.c_void_p(sim.firstHit.ctypes.data), ctypes.c_void_p(sim.firstHitLoc.ctypes.data),
                ctypes.c_int(sim.n), ctypes.c_int(sim.maxiter),
                ctypes.c_int(sim.rangeLow), ctypes.c_int(sim.rangeHigh), ctypes.c_void_p(sim.randomTable.ctypes.data)
            )

            temp = pd.DataFrame(columns=['s', 'condition_sat', 'condition_force', 'condition_mc', 'accuracy', 'RT', 'PMT', 'MT', 'latency_first_pB', 'location_first_pB'])
            temp['response'] = sim.resp
            temp['accuracy'] = [1 if r==1 else 0 for r in sim.resp]
            temp['RT'] = sim.RT
            temp['PMT'] = sim.PMT
            temp['MT'] = sim.MT
            temp['latency_first_pB'] = sim.firstHit
            temp['location_first_pB'] = [1 if loc==1 else 0 for loc in sim.firstHitLoc]
            temp['condition_force'] = np.repeat(f, 100000)
            temp['condition_mc'] = np.repeat(mc, 100000)
            temp['s'] = np.repeat(s, 100000)
            cls.append(temp)

    classifier = pd.concat(cls, ignore_index=True)
    classifier.to_csv(f'processing/data/s{s}.csv')

if __name__=='__main__':
    params_df = pd.read_csv('params.csv', index_col=0)
    n_subjects = len(params_df.index)
    with multiprocessing.Pool(processes=12) as pool:
        res = pool.map(func=run_sx, iterable=params_df.index)
        pool.close()
        pool.join()

#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
#include <math.h>

void GCDM(double x0, double v, double r, double g, double Te, double Tr, double xi, double leak, double u, double drift_guess,
          double sx0, double sv, double sTe, double sTr, double sleak, double s, double dt, 
          double *resp, double *RT, double *PMT, double *MT, double *firstHit, double *firstHitLoc,          
          int n, int maxiter, int rangeLow, int rangeHigh, double *randomTable)
{
    double dv, x, y, urg, zU, zL, randNum;
    int i, iter, randIndex, outOfGate, hit;
    struct timeval t1;
    gettimeofday(&t1, NULL);
    srand(t1.tv_usec * t1.tv_sec);

    for (i=0; i<n; i++) {
        dv = x0;
        x = dv;
        y = dv;
        
        resp[i] = -1.0;
        RT[i] = -1.0;
        PMT[i] = -1.0;
        MT[i] = -1.0;
        firstHit[i] = -1.0;
        firstHitLoc[i] = -1.0;

        iter = 0;
        outOfGate = 0;
        hit = 0;
        do {
            iter = iter+1;
            
            randNum = rand()/(1.0 + RAND_MAX);
            randIndex = (randNum * (rangeHigh - rangeLow + 1)) + rangeLow;
            randNum = randomTable[randIndex];
            dv = dv + (v*dt) + (sqrt(dt)*s*randNum); // decision variable
            
            randNum = rand()/(1.0 + RAND_MAX);
            randIndex = (randNum * (rangeHigh - rangeLow + 1)) + rangeLow;
            randNum = randomTable[randIndex];
            x = dv + (xi/sqrt(dt))*randNum; // if (xi != 0) then noisy transmission else pure transmission
            
            y = y + leak*(x - y)*dt + drift_guess*dt; // motor preparation variable, constant gain
            urg = u*iter*dt; // if (u != 0) then urgency-force signal helping else nothing
            zU = y - g + urg; // neural drives
            zL = -y - g + urg;
            
            if (((zU > 0) || (zL > 0)) && (outOfGate == 0)) {
                if (hit == 0) {
                    hit = 1; // this is the first time the gate is overcome
                    firstHit[i] = ((iter*dt) - (dt/2.0)) + Te; // latency of this first gate overcoming
                    if (zU > 0) { firstHitLoc[i] = 1.0; }
                    if (zL > 0) { firstHitLoc[i] = 2.0; } // location
                }
                outOfGate=1; // the decision variable is now located in the region between an EMG bound and the corresponding response bound
                PMT[i] = ((iter*dt) - (dt/2.0)) + Te; // in case several EMG bound hits occur, PMT will store the latency of the last hit before the response 
            }
            
            if ((outOfGate == 1) && ((zU <= 0) && (zL <= 0))) { outOfGate = 0; }
            
            if (zU >= r) {
                resp[i] = 1.0;
                RT[i] = (iter*dt) - (dt/2.0) + Te + Tr;
                MT[i] = RT[i] - PMT[i];
                break;
            }

            if (zL >= r) {
                resp[i] = 2.0; 
                RT[i] = (iter*dt) - (dt/2.0) + Te + Tr;
                MT[i] = RT[i] - PMT[i];
                break;
            }
        } while (iter<maxiter);
    }
}

内容的提问来源于stack exchange，提问作者Edouard D.