odeint -- segmentation fault
Dear scipy users and developers, I have the following problem with odeint: Suppose I would call odeint as follows y=odeint(dydt,y,t,....), where dydt is a complicated function that also calls odeint. So dydt is written as def dydt(y,t): <blabla> intermed=odeint(rhs,y0,t,...) <blabla> rhs is written in C and wrapped through SWIG In my case, I get segmentation faults doing this. Unfortunately, I was not able to cook up an artificial example that is simple, and behaves like this. What I can tell is that, 1) if I let dydt return a zero array (or any "direct" computation of the derivative), then the program works, while if dydt uses odeint, it crashes. 2) The return value of odeint does not have to be *used* for the crash to happen: e.g. one can use the complicated routine inside dydt, but still return zeros, and the crash would take place. 3) Using the complicated routine inside a loop that performs e.g. forward Euler, gives no problem at all. The crash really occurs when the value of dydt is passed into odeint. (This was tested by stepping with pdb.) So a segmentation fault occurs when odeint accesses the result of calling dydt if that function also uses odeint in a complicated way. Maybe the fact that inside dydt, odeint is calling a C routine for its time derivative, has something to do with it? Since I do not have access to odeint's code, I cannot tell what is going on, nor how to solve this... How could this problem be solved effectively? Best, Giovanni -- Giovanni Samaey http://www.cs.kuleuven.ac.be/~giovanni/ Katholieke Universiteit Leuven email: giovanni@cs.kuleuven.ac.be Departement Computerwetenschappen phone: +32-16-327081 Celestijnenlaan 200A, B-3001 Heverlee, Belgium fax: +32-16-327996 Office: A04.36
Giovanni Samaey wrote: [snip]
Since I do not have access to odeint's code, I cannot tell what is going on, nor how to solve this...
How could this problem be solved effectively?
First, try an example where rhs is not a SWIGed function. That might be the problem itself. I tried a simple example with a Python rhs (on a Powerbook with CVS SciPy) and I didn't get a segfault. If you want to see odeint's source code to help you debug, just download it: http://www.scipy.org/download/ When you have more information, also please provide your platform and SciPy version so we know where to look. Thanks. -- Robert Kern rkern@ucsd.edu "In the fields of hell where the grass grows high Are the graves of dreams allowed to die." -- Richard Harter
Hi, I have *much* more information now. A simple Python rhs function suffices to get the segmentation fault. A sample code that produces the segfault is attached: The function integrate calls odeint with the "equation being defined through dydt. This dydt initialized another integration, around each point x, and calls a function step, which executes odeint with "equation" rhs. The crucial point to get the segfault is that rhs has a different number of odes than dydt. The sample code added here gives print statements that should help understand execution, and a comment is added that would indicate the segfault. I am working on linux with SciPy version '0.3.0_266.4239', Numeric 23.3 and python 2.3.4 Robert Kern wrote:
Giovanni Samaey wrote:
[snip]
Since I do not have access to odeint's code, I cannot tell what is going on, nor how to solve this...
How could this problem be solved effectively?
First, try an example where rhs is not a SWIGed function. That might be the problem itself. I tried a simple example with a Python rhs (on a Powerbook with CVS SciPy) and I didn't get a segfault.
If you want to see odeint's source code to help you debug, just download it: http://www.scipy.org/download/
When you have more information, also please provide your platform and SciPy version so we know where to look.
Thanks.
-- Giovanni Samaey http://www.cs.kuleuven.ac.be/~giovanni/ Katholieke Universiteit Leuven email: giovanni@cs.kuleuven.ac.be Departement Computerwetenschappen phone: +32-16-327081 Celestijnenlaan 200A, B-3001 Heverlee, Belgium fax: +32-16-327996 Office: A04.36 import scipy from scipy.integrate import odeint,simps def step(y0,t,x): print "In step!" print "y0: ",scipy.shape(y0) y0=scipy.ravel(y0) y=odeint(rhs,y0,scipy.array([t,t+1e-5]),ml=1,mu=1,args=(x,)) print "y: ", scipy.shape(y) ret=scipy.reshape(y[-1,:],(len(y0),)) print "ret: ", scipy.shape(ret) return scipy.reshape(y[-1,:],(len(y0),)) def rhs(y,t,x): print "In rhs!" ny=scipy.size(y) ydot=scipy.zeros(ny,scipy.Float) ydot[1:-1]=y[2:]-2*y[1:-1]+y[:-2] return ydot def dydt(y,t,x): print "In dydt!" y=scipy.reshape(y,(1,len(y))) print "y: ", scipy.shape(y) ynew=scipy.zeros(scipy.shape(y)) for i in scipy.arange(len(x)): print i xx=scipy.arange(x[i]-5e-3,x[i]+5e-3,1e-7) yy=scipy.ones(scipy.shape(xx),scipy.Float) yynew=step(yy,t,xx) print "yynew: ", scipy.shape(yynew) ynew[0,i]=yynew[500] print "ynew: ",scipy.shape(y) print "dydt shape: ", scipy.shape(scipy.ravel((ynew-y)/1e-5)) # the segmentation fault occurs when this return statement is executed! return scipy.ravel((ynew-y)/1e-5) def integrate(y,t,x): print "In integrate!" print "Initial condition shape: ", scipy.shape(scipy.ravel(y)) return odeint(dydt,scipy.ravel(y),t,ml=1,mu=1,args=(x,)) if __name__=="__main__": x=scipy.arange(0,1.05,0.1) y0=1-4*(x-0.5)**2 result=integrate(y0,t=scipy.arange(0,1e-2,5e-3),x=x)
Giovanni Samaey wrote:
Hi,
I have *much* more information now. A simple Python rhs function suffices to get the segmentation fault. A sample code that produces the segfault is attached: The function integrate calls odeint with the "equation being defined through dydt. This dydt initialized another integration, around each point x, and calls a function step, which executes odeint with "equation" rhs.
The crucial point to get the segfault is that rhs has a different number of odes than dydt. The sample code added here gives print statements that should help understand execution, and a comment is added that would indicate the segfault.
I am working on linux with SciPy version '0.3.0_266.4239', Numeric 23.3 and python 2.3.4
Robert Kern wrote:
I think I've fixed this problem. The problem seems to be that lsoda (written in Fortran) made heavy use of common blocks which made it difficult to use in a re-entrant fashion (your dydt function itself calls odeint). ODEPACK provided a function to save and restore the common block which I sprinkled throughout the original fortran code so that the functions called (and jacobians) could use odeint without ruining the needed common block for the rest of the algorithm. Your test script now works on my system and the fixes are in CVS. I hope this helps -Travis
Travis E. Oliphant wrote:
I think I've fixed this problem. The problem seems to be that lsoda (written in Fortran) made heavy use of common blocks which made it difficult to use in a re-entrant fashion (your dydt function itself calls odeint).
This is what I suspected to be the problem (the re-entrant use). I am really thankful for this fix; it helped me a lot. Do you suspect that deeper recursions could cause problems? (Just wondering, I don't need this for now.)
ODEPACK provided a function to save and restore the common block which I sprinkled throughout the original fortran code so that the functions called (and jacobians) could use odeint without ruining the needed common block for the rest of the algorithm.
Your test script now works on my system and the fixes are in CVS.
Thank you. The test script works here too now; and my "real" code works too! On question concerning the CVS version of scipy. If I do regular updates, do I have to rebuild scipy completely, or does the build script take care of checking what changed? And can I break running code by updating regularly? Best, Giovanni
I hope this helps
-Travis
_______________________________________________ SciPy-user mailing list SciPy-user@scipy.net http://www.scipy.net/mailman/listinfo/scipy-user
-- Giovanni Samaey http://www.cs.kuleuven.ac.be/~giovanni/ Katholieke Universiteit Leuven email: giovanni@cs.kuleuven.ac.be Departement Computerwetenschappen phone: +32-16-327081 Celestijnenlaan 200A, B-3001 Heverlee, Belgium fax: +32-16-327996 Office: A04.36
participants (3)
-
Giovanni Samaey -
Robert Kern -
Travis E. Oliphant