Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP

Loading…

ode crashes if rhs returns a tuple instead of a list (Trac #1187) #1714

Closed
scipy-gitbot opened this Issue · 11 comments

1 participant

@scipy-gitbot

Original ticket http://projects.scipy.org/scipy/ticket/1187 on 2010-05-31 by trac user miha, assigned to unknown.

I took the example from

[http://docs.scipy.org/doc/scipy-0.7.x/reference/generated/scipy.integrate.ode.html#scipy.integrate.ode]

and changed f's return type from list to tuple (i removed [ and ] in f's return statement). The example then crashed with the following output:

0-th dimension must be 2 but got 0 (not defined).
rv_cb_arr is NULL
Call-back cb_f_in_zvode__user__routines failed.
Traceback (most recent call last):
  File "proba.py", line 16, in <module>
    r.integrate(r.t+dt)
  File "/usr/lib/python2.6/dist-packages/scipy/integrate/ode.py", line 260, in integrate
    self.f_params,self.jac_params)
  File "/usr/lib/python2.6/dist-packages/scipy/integrate/ode.py", line 554, in run
    y1,t,istate = self.runner(*(args[:5]+tuple(self.call_args)+args[5:]))
SystemError: NULL result without error in PyObject_Call

shell returned 1

I use scipy 0.7.2.

@scipy-gitbot

@WarrenWeckesser wrote on 2011-11-11

This appears to be an f2py bug. f2py generates the file vodemodule.c,
which contains this snippet:

  capi_return = PyObject_CallObject(cb_f_in_dvode__user__routines_capi,(PyObject *)capi_arglist);
#ifdef F2PY_REPORT_ATEXIT
f2py_cb_stop_call_clock();
#endif
  CFUNCSMESSPY("cb:capi_return=",capi_return);
  if (capi_return == NULL) {
    fprintf(stderr,"capi_return is NULL\n");
    goto capi_fail;
  }
  if (capi_return == Py_None) {
    Py_DECREF(capi_return);
    capi_return = Py_BuildValue("()");
  }
  else if (!PyTuple_Check(capi_return)) {
    capi_return = Py_BuildValue("(N)",capi_return);
  }
  capi_j = PyTuple_Size(capi_return);
  capi_i = 0;
/*frompyobj*/
  if (capi_j>capi_i) {
    PyArrayObject *rv_cb_arr = NULL;
    if ((capi_tmp = PyTuple_GetItem(capi_return,capi_i++))==NULL) goto capi_fail;
    rv_cb_arr =  array_from_pyobj(PyArray_DOUBLE,ydot_Dims,1,F2PY_INTENT_IN
|F2PY_INTENT_C
,capi_tmp);
    if (rv_cb_arr == NULL) {
      fprintf(stderr,"rv_cb_arr is NULL\n");
      goto capi_fail;
    }
    MEMCOPY(ydot,rv_cb_arr->data,PyArray_NBYTES(rv_cb_arr));
    if (capi_tmp != (PyObject *)rv_cb_arr) {
      Py_DECREF(rv_cb_arr);
    }
  }
  CFUNCSMESS("cb:cb_f_in_dvode__user__routines:successful\n");
  Py_DECREF(capi_return);

In the first line of the snippet, capi_return is the python object returned
by the call to the the callback function f(). Several checks are made:

  • if capi_return is NONE, something bad happened (not relevant here);

  • if capi_return is Py_None, the user's function returned None. This code
    replaces None with an empty tuple;

  • if capi_return is is not a tuple, it is wrapped in a tuple of length 1:

        capi_return = Py_BuildValue("(N)",capi_return);
    

Note that if the return value is already a tuple, it is unchanged by these
initial checks. Then this line

  capi_j = PyTuple_Size(capi_return);   

results in capi_j being 2, so in the next 'if' statement, capi_j > capi_i, and
the if-block is executed. That code calls PyTuple_GetItem to get the first
element of the tuple, which it then passes to array_from_pyobj. This is the
problem, because if the return value of f() was a tuple, it is only the first
element of the tuple--a scalar--that is being passed to array_from_pyobj in
the next line. The function array_from_pyobj() is defined in fortranobject.c.
It will call check_and_fix_dimensions() (also in fortranobject.c) like this:

        if (check_and_fix_dimensions(arr,rank,dims)) {
            return NULL; /*XXX: set exception */
        }

It is check_and_fix_dimensions() that is printing the error message

0-th dimension must be 2 but got 0 (not defined).

because it got a scalar value, but it expected an array with shape (2,).
Then the code in array_from_pyobj() detects that check_and_fix_dimensions()
has failed, and returns NULL, but--as noted with 'XXX' in the comment--it
does not properly set the error exception. This results in the error message:

SystemError: NULL result without error in PyObject_Call

when the exception propagates up to the user.

I don't know why the code in vodemodule.c wraps the non-tuple result in
a tuple, only to pull out the value later. But this appears to be the source
of the bug. This means, for example, that if f() returns a tuple of length 1
containing the tuple of values, the solver works--but that's just a curiosity,
not a solution or work-around.

@scipy-gitbot

@rgommers wrote on 2011-11-11

note: still present in master.

Are we sure this is an f2py bug and not a result of using f2py the wrong way? If the former, probably good to open a new bug on numpy Trac.

@scipy-gitbot

@WarrenWeckesser wrote on 2011-11-11

Replying to [comment:2 rgommers]:

Are we sure this is an f2py bug and not a result of using f2py the wrong way?

Actually, no (and apologies to Pearu for a possibly premature assessment :). I haven't looked at the input files that lead to the generation of vodemodule.c.

@scipy-gitbot

@pearu wrote on 2011-11-14

f2py supports user-defined functions that have multiple return values.
These return values are expected to be returned as a single tuple object.
So, the issue here is not a result of f2py bug but due to this f2py feature
that is unintentionally triggered by the given example.

As a possible fix, we could have the f2py code to treat the returned tuples
as ordinary sequences when the expected number of returned values is exactly 1.
On the other hand, in the case of f2py misuse (returning 2 values when 1 is expected,
for instance) the values will be silently converted to an array and may result
even a more complicated issues to be analyzed.

I think the best fix is to document this f2py feature in the corresponding
codes, that is, never return tuple as a single return value,
and perhaps improve the error message giving a hint why the crash might have been occurred.

@scipy-gitbot

@rgommers wrote on 2011-11-14

Documenting plus clear error message sounds reasonable.

@scipy-gitbot

trac user RhysU wrote on 2012-12-21

Appears to not just be related to tuple-vs-list. The following

import numpy as np
import numpy.random
import scipy.integrate.ode
import matplotlib.pyplot as plt

# Given state [x, y, z] compute right hand side of the Lorenz equations.
# The default coefficients compute the canonical attractor.
def lorenz(t, xyz, coefficients = [10., 8./3., 28.]):
    x, y, z = xyz
    sigma, beta, rho = coefficients
    return [ sigma*(y-x), x*(rho-z)-y, x*y-beta*z ]

def lorenz_compute(tf = 1000, dt = 0.1, xyz0 = np.random.randn(3)):
    solver = scipy.integrate.ode(lorenz_compute)
    solver.set_integrator('vode')
    solver.set_initial_value(xyz0, 0)
    soln = []
    while solver.successful() and solver.t < tf:
        solver.integrate(solver.t + dt)
        soln.append([solver.t, solver.y[0], solver.y[1], solver.y[2]])
    return np.array(soln)

if __name__ == "__main__":
    a = lorenz_compute()

bombs out with

[4883 rhys@setun lorenz]$ python sample.py
unexpected array size: new_size=3, got array with arr_size=0
rv_cb_arr is NULL
Call-back cb_f_in_dvode__user__routines failed.
Traceback (most recent call last):
  File "sample.py", line 24, in <module>
    a = lorenz_compute()
  File "sample.py", line 19, in lorenz_compute
    solver.integrate(solver.t + dt)
  File "/usr/lib/python2.6/dist-packages/scipy/integrate/ode.py", line 260, in integrate
    self.f_params,self.jac_params)
  File "/usr/lib/python2.6/dist-packages/scipy/integrate/ode.py", line 435, in run
    y1,t,istate = self.runner(*(args[:5]+tuple(self.call_args)+args[5:]))
SystemError: NULL result without error in PyObject_Call

Notice the lack of tuple-vs-list-iness in function lorenz.

@scipy-gitbot

@rgommers wrote on 2013-01-05

The first line of your lorenz_compute function is incorrect. This works as expected:

import numpy as np
from scipy import integrate

# Given state [x, y, z] compute right hand side of the Lorenz equations.
# The default coefficients compute the canonical attractor.
def lorenz(t, xyz, coefficients = [10., 8./3., 28.]):
    x, y, z = xyz
    sigma, beta, rho = coefficients
    return [sigma*(y-x), x*(rho-z)-y, x*y-beta*z]

def lorenz_compute(tf = 1000, dt = 0.1, xyz0 = np.random.randn(3)):
    solver = integrate.ode(lorenz)
    solver.set_integrator('vode')
    solver.set_initial_value(xyz0, 0)
    soln = []
    while solver.successful() and solver.t < tf:
        solver.integrate(solver.t + dt)
        soln.append([solver.t, solver.y[0], solver.y[1], solver.y[2]])

    return np.array(soln)

if __name__ == "__main__":
    a = lorenz_compute()

And then, changing the return of lorenz from list of 3 elements to tuple of elements gives the SystemError again.

@scipy-gitbot

@rgommers wrote on 2013-01-05

Generate an understandable error in #395

That closes this ticket I think; for the f2py change suggested by Pearu a numpy ticket should be opened.

@scipy-gitbot

Milestone changed to 0.12.0 by @rgommers on 2013-01-05

@scipy-gitbot

trac user RhysU wrote on 2013-01-05

Replying to [comment:8 rgommers]:

The first line of your lorenz_compute function is incorrect.

Wow. Indeed. Thanks for pointing out my glaring mistake.

@scipy-gitbot

@rgommers wrote on 2013-02-12

Opened numpy/numpy#2981 for Pearu's suggestion on how to deal with this in f2py. Closing this ticket.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Something went wrong with that request. Please try again.