Tableau C vers PyArray

Question

Je suis en train d'écrire une extension C de Python sans utiliser Cython.

Je veux allouer un tableau double en C, l'utiliser dans une fonction interne (qui se trouve être en Fortran) et le retourner. Je fais remarquer que l'interface C-Fortran fonctionne parfaitement en C.

static PyObject *
Py_drecur(PyObject *self, PyObject *args)
{
  // INPUT
  int n;
  int ipoly;
  double al;
  double be;

  if (!PyArg_ParseTuple(args, "iidd", &n, &ipoly, &al, &be))
    return NULL;

  // OUTPUT
  int nd = 1;
  npy_intp dims[] = {n};
  double a[n];
  double b[n];
  int ierr;

  drecur_(n, ipoly, al, be, a, b, ierr);

  // Create PyArray
  PyObject* alpha = PyArray_SimpleNewFromData(nd, dims, NPY_DOUBLE, a);
  PyObject* beta = PyArray_SimpleNewFromData(nd, dims, NPY_DOUBLE, b);

  Py_INCREF(alpha);
  Py_INCREF(beta);

  return Py_BuildValue("OO", alpha, beta);
}

J'ai débogué ce code et j'obtiens un défaut de segmentation lorsque j'essaie de créer un alpha à partir de a. Jusque là, tout fonctionne bien. La fonction drecur_ fonctionne et j'obtiens le même problème si elle est supprimée.

Maintenant, quelle est la manière standard de définir un PyArray autour de données C ? J'ai trouvé de la documentation mais pas de bon exemple. De plus, qu'en est-il des fuites de mémoire ? Est-il correct d'INCREF avant le retour pour que les instances de alpha et beta soient préservées ? Qu'en est-il de la désallocation lorsqu'ils ne sont plus nécessaires ?

EDIT J'ai finalement réussi avec l'approche trouvée dans Livre de recettes NumPy .

static PyObject *
Py_drecur(PyObject *self, PyObject *args)
{
  // INPUT
  int n;
  int ipoly;
  double al;
  double be;
  double *a, *b;
  PyArrayObject *alpha, *beta;

  if (!PyArg_ParseTuple(args, "iidd", &n, &ipoly, &al, &be))
    return NULL;

  // OUTPUT
  int nd = 1;
  int dims[2];
  dims[0] = n;
  alpha = (PyArrayObject*) PyArray_FromDims(nd, dims, NPY_DOUBLE);
  beta = (PyArrayObject*) PyArray_FromDims(nd, dims, NPY_DOUBLE);
  a = pyvector_to_Carrayptrs(alpha);
  b = pyvector_to_Carrayptrs(beta);
  int ierr;

  drecur_(n, ipoly, al, be, a, b, ierr);

  return Py_BuildValue("OO", alpha, beta);
}

double *pyvector_to_Carrayptrs(PyArrayObject *arrayin)  {
  int n=arrayin->dimensions[0];
  return (double *) arrayin->data;  /* pointer to arrayin data as double */
}

N'hésitez pas à faire des commentaires à ce sujet et merci pour les réponses.

Answer 1

Donc la première chose qui semble suspecte, c'est que votre tableau a y b sont dans la portée locale de la fonction. Cela signifie qu'après le retour, vous obtiendrez un accès illégal à la mémoire.

Vous devez donc déclarer les tableaux avec

double *a = malloc(n*sizeof(double));

Ensuite, vous devez vous assurer que la mémoire est ensuite libérée par l'objet que vous avez créé. Voir cette citation de la documentation :

PyObject PyArray_SimpleNewFromData(int nd, npy_intp dims, int typenum, void* data)

Parfois, vous voulez envelopper la mémoire allouée ailleurs dans un objet ndarray pour une utilisation en aval. Cette routine permet de le faire de manière simple. Les trois premiers arguments sont les mêmes que dans PyArray_SimpleNew, le dernier argument est un pointeur vers un bloc de mémoire contiguë que le ndarray doit utiliser comme tampon de données qui sera interprété de façon contiguë comme en C. Une nouvelle référence à un ndarray est retournée, mais le ndarray ne sera pas propriétaire de ses données. Lorsque ce tableau est désalloué, le pointeur ne sera pas libéré.

Vous devez vous assurer que la mémoire fournie n'est pas libérée alors que le tableau retourné existe. La manière la plus simple de gérer cela est de faire en sorte que les données proviennent d'un autre objet Python comptant des références. Le nombre de références de cet objet doit être augmenté après le passage du pointeur, et le membre de base du tableau retourné doit pointer vers l'objet Python qui possède les données. Ensuite, lorsque le ndarray est désalloué, le membre de base sera DECREFé de manière appropriée. Si vous voulez que la mémoire soit libérée dès que le ndarray est désalloué, il suffit de mettre le flag OWNDATA sur le ndarray retourné.

Pour votre deuxième question, le Py_INCREF(alpha); n'est généralement nécessaire que si vous avez l'intention de conserver la référence dans une variable globale ou un membre de classe. Mais comme vous ne faites qu'envelopper une fonction, vous n'avez pas besoin de le faire. Malheureusement, il se peut que la fonction PyArray_SimpleNewFromData ne met pas le compteur de référence à 1, si c'était le cas vous devriez l'augmenter à 1. J'espère que c'était compréhensible ;).

Answer 2

Un problème pourrait être que vos tableaux (a,b) doivent durer au moins aussi longtemps que le numpy-array qui le contient. Vous avez créé vos tableaux dans une portée locale, ils seront donc détruits lorsque vous quitterez la méthode.

Essayez de faire en sorte que python alloue le tableau (par exemple en utilisant PyArray_SimpleNew ), copiez-y votre contenu et passez-lui un pointeur. Vous pouvez également utiliser boost::python pour s'occuper de ces détails, si la construction contre le boost est une option.