hensing / pydde Goto Github PK

This project deserves to be brought to current python standards. There are still references to Numeric!
I'll try to contribute a pull request when I have a bit of free time.

Hi Henning,

when using PyDDE with more than 10 state variables I get a segmentation fault possibly due to memory problems. Before digging into the problem myself, I was wondering if you are aware of this issue and can offer some quick help.

Best,
Jannis

Whenever a grad function raises an exception instead of displaying the python exception this program generates a segfaut (this is an aritificial example I cooked up from a real dde that was segfaulting for the same reason after iteration 100):

from __future__ import division
import numpy as np
import scipy as sp
import time
from scipy.special import binom
from PyDDE import pydde

tau = 40
d0 = 1 / 100
G = 10 ** (-2)
a = 0.8
A = 100
# NB: B * x_0 < ln(2) for boundaries in U_x to be defined
B = 0.06
x_0 = 3
t_delay = 50
ddecons = np.array([tau, d0, G, a, A, B, x_0, t_delay])

# Duration of the simulation:
stoptime = 1000.0
numpoints = 1000
timestep = 1.0 # this is the timesteps on which the output is given
dt = 0.1 # this is the integration timestep of the solver


# Initial Conditions
N0 = 0.01
No0 = 10
N_infty0 = 0
ddeist = np.array([N0, No0, N_infty0])
# Set state-scaling array for use in error control when values are very close to 0
ddestsc = np.array([0,0,0])

#  Create a dde object
dde_eg = pydde.dde()

def ddegrad(s, c, t):
    if int(s[0]) != 1 :
        return np.array([1,1,1])
    raise BaseException()


ts = np.arange(0.0, stoptime, timestep)

dde_eg.dde(y=ddeist, times=ts,
           func=ddegrad, parms=ddecons,
           tol=0.000005, dt=dt, hbsize=3000, nlag=0, ssc=ddestsc)

looking at the source (inside wrapper.c) the missing of the exception seems to take place here (pyarraycheck is not working it seems):

    result = PyEval_CallObject(grad_func,arglist);
    //printf("Got grad_func result.\n");
    assert(result);
    assert(PyArray_Check(result));
    //printf("result = %p\n", result);
    currentdata = dblArray_from_PyArray(result);

Whenever using PyEval_CallObject it is good to check that it does not return NULL. If it returns NULL this means that an exception occurred within grad_func in python and we should EXPLICITLY call the exception (i.e. https://docs.python.org/3.0/extending/extending.html ) otherwise this will most probably lead to a segfault which it does dblArray_from_PyArray is expecting a pyobject sequence and when passing it to numpy to cast to a contiguous array it raises a segfault since it is receiving a NULL object and not a sequence.

This makes it very hard for users to debug what is going on as it my not be obvious why their gradient function is raising an exception since it jumps to a segfault and not the exception message.