[Numpy-discussion] Problem looping over numpy array in C

[Sumant S.R. Oemrawsingh]

Hi guys,

I have a problem with looping over numpy arrays in C. I modify the array
in-place (and return None), but after modification, the array doesn't seem
to play nice any more.

Below, I have the C code for a function do_something (a stripped version
of my original function), which has as arguments a 1D numpy array (float
or complex) and either an int or a sequence of ints.

In python, I do the following using only the integer:

>>> a = array([0., 1., 2., 3., 2., 1., 0.])
>>> do_something(a,3)
>>> savetxt('bla',a)
>>>

Which works fine. However, when I do the same, but with a list of any
length larger than 0:

>>> a = array([0., 1., 2., 3., 2., 1., 0.])
>>> do_something(a,[3])
>>> savetxt('bla',a)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python2.5/site-packages/numpy/core/numeric.py", line 767,
in savetxt
    X = asarray(X)
  File "/usr/lib/python2.5/site-packages/numpy/core/numeric.py", line 132,
in asarray
    return array(a, dtype, copy=False, order=order)
TypeError: an integer is required
>>> savetxt('bla',a)
>>>

For some reason, the first time it doesn't work; the asarray fails but
then succeeds on a second try. The resulting file is identical to when I
would have looped over the list in python, and called the do_something
function with each integer separately. For instance:

for i in [3,1,0]:
  do_something(a, i)

works fine as well. So apparently looping in C temporarily leaves the
array in a weird state but manages to automatically be restored after one
exception. I've checked that type(a), a.dtype, a.shape and a.ndim remain
the same before and after calling do_something with a sequence or with an
integer as second argument. That doesn't seem to be the problem.

The reason that I want do the loop in C is that I need some
precalculations, before being able to do the actual loop. But they might
not be time-consuming enough to warrant writing it in C, so I could just
do the loop in python and not have this problem any more. However, if the
loop in C manages to somehow (temporarily) corrupt the array, how can I be
sure that the single-shot call doesn't succeed by accident?

If anyone can help, suggest something to try or spot my mistake, I'd
appreciate it.

Thanks,
Sumant


static PyObject*
do_something(PyObject *self, PyObject *args, PyObject *kwargs)
{
  PyObject *input, *s, *item;
  PyArrayObject *array = NULL;
  npy_intp i,n;
  static char *kwlist[] = {"a", "index", NULL};

  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "Oi", kwlist,
        &input, &i))
    goto sequence;

  if ((array = (PyArrayObject *) PyArray_FromAny(input, NULL, 1, 1,
          NPY_CARRAY | NPY_UPDATEIFCOPY, NULL)) == NULL)
    goto error;

  n = PyArray_DIM(array,0);

  if (PyArray_ISCOMPLEX(array))
    do_something_complex((complex*)(PyArray_DATA(array)), n, i);
  else if (PyArray_ISFLOAT(array))
    do_something_double((double*)(PyArray_DATA(array)), n, i);
  else
    goto error;

  Py_DECREF(array);
  Py_RETURN_NONE;

sequence:
  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OO", kwlist,
        &input, &s))
    goto error;

  if (PySequence_Check(s) == 0)
    goto error;

  if ((array = (PyArrayObject *) PyArray_FromAny(input, NULL, 1, 1,
          NPY_CARRAY | NPY_UPDATEIFCOPY, NULL)) == NULL)
    goto error;

  n = PyArray_DIM(array,0);

  if (PyArray_ISCOMPLEX(array))
  {
    complex* data = (complex*)PyArray_DATA(array);
    for (i = 0; i < PySequence_Size(s); i++)
    {
      item = PySequence_GetItem(s, i);
      do_something_complex(data, n, PyInt_AsSsize_t(item));
      Py_DECREF(item);
    }
  }
  else if (PyArray_ISFLOAT(array))
  {
    double* data = (double*)PyArray_DATA(array);
    for (i = 0; i < PySequence_Size(s); i++)
    {
      item = PySequence_GetItem(s, i);
      do_something_double(data, n, PyInt_AsSsize_t(item));
      Py_DECREF(item);
    }
  }
  else
    goto error;

  Py_DECREF(array);
  Py_RETURN_NONE;

error:
  Py_XDECREF(array);
  return NULL;
}