Python C module questions

[jeremy.d.brewer@gmail.com]

Hi,

I'm rather new to Python, and I've just written my first Python C
module. I was wondering if some more experience Pythonista would look
over what I've written and given me some pointers (or find some bugs).
I had a few questions while writing this as well.

Also, I know that there are much better ways to compute nearest
neighbors than the brute force n^2 method, but this is plenty fast for
this application (at least the C version is).

1. How do I add docstrings to my module?
2. The Python reference counting and memory managment seemes awfully
repetetive and error prone. Is there a nicer way of doing this? I
know there are some wrappers out there such as Pyrex and SWIG that
might prove useful.
3. This module consisted of turning 4 Python sequences into C double
arrays, performing some calculations, and then turning a C int array
into a Python tuple for return. And it was a lot of work. Are there
any nice wrapper functions out there for turning Python sequences into
C arrays and vice versa?

Thanks,
Jeremy Brewer

The code is below:

#include <Python.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>

/* xmatch: Takes 2 sets of positions (ra1, dec1 of size len1 and ra2,
dec2 of
of size len2) and computes the nearest neighbors. An array
of
integers for the nearest neighbors of ra1, dec1 that
correspond to
the positions in ra2, dec2 is returned, where points further
away
than mindist (in degrees) on the sky are marked as
-(closest index + 1). */

static int *xmatch(double *ra1, double *dec1, int len1, double *ra2,
double *dec2, int len2, double mindist)
{
int i, j, closest, *neighbors;
double mindist2, dist, old_dist;

/* allocate memory for the neighbors */
neighbors = (int *) malloc(len1 * sizeof(int));
if (neighbors == NULL) {
return NULL;
}

/* compare to mindist^2 to save computational time (saves many
calls to
sqrt in inner loop) */
mindist2 = mindist * mindist;

for (i = 0; i < len1; i++) {
/* intial closest match is as far off as possible */
closest = 0;
old_dist = DBL_MAX;

for (j = 0; j < len2; j++) {
/* angular distance on the sky squared */
dist = pow((cos(dec1[i]) * (ra1[i] - ra2[j])), 2) +
pow((dec1[i] - dec2[j]), 2);

/* keep track of only the closest point */
if (dist < old_dist) {
closest = j;
old_dist = dist;
}
}

/* mark points that lie outside the specified distance */
if (old_dist > mindist2) {
neighbors[i] = -(closest + 1);
}
else {
neighbors[i] = closest;
}
}

return neighbors;
}

/* wrap_xmatch: Python wrapper function for xmatch above */
static PyObject *wrap_xmatch(Py Object *self, PyObject *args)
{
PyObject *ra1, *dec1, *ra2, *dec2, *ntup;
int i, len1, len2, len_chk1, len_chk2;
int *neighbors;
double mindist;
double *r1, *d1, *r2, *d2;

/* read 4 sequences (ra1, dec1, ra2, dec2) and a double (mindist)
*/
if (!PyArg_ParseTu ple(args, "OOOOd", &ra1, &dec1, &ra2, &dec2,
&mindist)) {
return NULL;
}

/* check that mindist is > 0 */
if (mindist <= 0.0) {
PyErr_SetString (PyExc_ValueErr or, "mindist must be > 0");
return NULL;
}

/* convert sequences to tuples if necessary */
ra1 = PySequence_Fast (ra1, "ra1 must be iterable");
if (ra1 == NULL) {
return NULL;
}

dec1 = PySequence_Fast (dec1, "dec1 must be iterable");
if (dec1 == NULL) {
return NULL;
}

ra2 = PySequence_Fast (ra2, "ra2 must be iterable");
if (ra2 == NULL) {
return NULL;
}

dec2 = PySequence_Fast (dec2, "dec2 must be iterable");
if (dec2 == NULL) {
return NULL;
}

/* length of input sequences */
len1 = PySequence_Fast _GET_SIZE(ra1);
len_chk1 = PySequence_Fast _GET_SIZE(dec1) ;

if (len1 != len_chk1) {
PyErr_SetString (PyExc_ValueErr or, "ra1 and dec1 must be the
same size");
return NULL;
}

len2 = PySequence_Fast _GET_SIZE(ra2);
len_chk2 = PySequence_Fast _GET_SIZE(dec2) ;

if (len2 != len_chk2) {
PyErr_SetString (PyExc_ValueErr or, "ra2 and dec2 must be the
same size");
return NULL;
}

/* allocate memory for C arrays */
r1 = (double *) malloc(len1 * sizeof(double)) ;
if (r1 == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
return PyErr_NoMemory( );
}

d1 = (double *) malloc(len1 * sizeof(double)) ;
if (d1 == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
return PyErr_NoMemory( );
}

r2 = (double *) malloc(len2 * sizeof(double)) ;
if (r2 == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
return PyErr_NoMemory( );
}

d2 = (double *) malloc(len2 * sizeof(double)) ;
if (d2 == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
return PyErr_NoMemory( );
}

/* copy ra1 and dec1 */
for (i = 0; i < len1; i++) {
PyObject *fitem, *item;

/* copy ra1 */
item = PySequence_Fast _GET_ITEM(ra1, i);
if (item == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
return NULL;
}

fitem = PyNumber_Float( item);
if (fitem == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
PyErr_SetString (PyExc_TypeErro r, "items in ra1 must be
numbers");
return NULL;
}

r1[i] = PyFloat_AS_DOUB LE(fitem);
Py_DECREF(fitem );

/* copy dec1 */
item = PySequence_Fast _GET_ITEM(dec1, i);
if (item == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
return NULL;
}

fitem = PyNumber_Float( item);
if (fitem == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
PyErr_SetString (PyExc_TypeErro r, "items in dec1 must be
numbers");
return NULL;
}

d1[i] = PyFloat_AS_DOUB LE(fitem);
Py_DECREF(fitem );
}

/* copy ra2 and dec2 */
for (i = 0; i < len2; i++) {
PyObject *fitem, *item;

/* copy ra2 */
item = PySequence_Fast _GET_ITEM(ra2, i);
if (item == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
return NULL;
}

fitem = PyNumber_Float( item);
if (fitem == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
PyErr_SetString (PyExc_TypeErro r, "items in ra2 must be
numbers");
return NULL;
}

r2[i] = PyFloat_AS_DOUB LE(fitem);
Py_DECREF(fitem );

/* copy dec1 */
item = PySequence_Fast _GET_ITEM(dec2, i);
if (item == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
return NULL;
}

fitem = PyNumber_Float( item);
if (fitem == NULL) {
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;
free(r1);
free(d1);
free(r2);
free(d2);
PyErr_SetString (PyExc_TypeErro r, "items in dec2 must be
numbers");
return NULL;
}

d2[i] = PyFloat_AS_DOUB LE(fitem);
Py_DECREF(fitem );
}

/* clean up Python objects */
Py_DECREF(ra1);
Py_DECREF(dec1) ;
Py_DECREF(ra2);
Py_DECREF(dec2) ;

/* now compute the cross match */
neighbors = xmatch(r1, d1, len1, r2, d2, len2, mindist);

/* clean up C arrays */
free(r1);
free(d1);
free(r2);
free(d2);

/* build a Python tuple to hold the returned neighbors */
ntup = PyTuple_New(len 1);
if (ntup == NULL) {
return PyErr_NoMemory( );
}

for (i = 0; i < len1; i++) {
PyObject *integer = Py_BuildValue(" i", neighbors[i]);
if (integer == NULL) {
Py_DECREF(ntup) ;
free(neighbors) ;
return PyErr_NoMemory( );
}

PyTuple_SET_ITE M(ntup, i, integer);
}

/* free last C array */
free(neighbors) ;

return ntup;
}

/* functions defined in the module */
static PyMethodDef astro_methods[] = {
{"xmatch", wrap_xmatch, METH_VARARGS, "Cross match 2 sets of
positions."},
{0} /* sentinel */
};

/* initializes the module */
void initastro(void)
{
(void) Py_InitModule(" astro", astro_methods);
}

[Robert Kern]

Re: Python C module questions

jeremy.d.brewer @gmail.com wrote:[color=blue]
> Hi,
>
> I'm rather new to Python, and I've just written my first Python C
> module. I was wondering if some more experience Pythonista would look
> over what I've written and given me some pointers (or find some bugs).
> I had a few questions while writing this as well.
>
> Also, I know that there are much better ways to compute nearest
> neighbors than the brute force n^2 method, but this is plenty fast for
> this application (at least the C version is).
>
> 1. How do I add docstrings to my module?[/color]

You mean a docstring on the module object itself? In initmodule() make a
string object with the contents of the docstring and assign it to
module.__doc__ . There might be a more idiomatic way, but I can't think
of it without trawling through the docs.
[color=blue]
> 2. The Python reference counting and memory managment seemes awfully
> repetetive and error prone. Is there a nicer way of doing this? I
> know there are some wrappers out there such as Pyrex and SWIG that
> might prove useful.[/color]

I'm quite a fan of Pyrex. You may also want to give ctypes a shot. SWIG
doesn't help as much with the reference counting
[color=blue]
> 3. This module consisted of turning 4 Python sequences into C double
> arrays, performing some calculations, and then turning a C int array
> into a Python tuple for return. And it was a lot of work. Are there
> any nice wrapper functions out there for turning Python sequences into
> C arrays and vice versa?[/color]

http://numeric.scipy.org

[color=blue]
> Thanks,
> Jeremy Brewer
>
> The code is below:[/color]

You probably shouldn't post such large pieces of code to the list.

--
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

[Thomas Heller]

Re: Python C module questions

Robert Kern <rkern@ucsd.edu > writes:
[color=blue]
> jeremy.d.brewer @gmail.com wrote:[color=green]
>> Hi,
>>
>> I'm rather new to Python, and I've just written my first Python C
>> module. I was wondering if some more experience Pythonista would look
>> over what I've written and given me some pointers (or find some bugs).
>> I had a few questions while writing this as well.
>>
>> Also, I know that there are much better ways to compute nearest
>> neighbors than the brute force n^2 method, but this is plenty fast for
>> this application (at least the C version is).[/color][/color]

You should give up C with a dumb algorithm running at fast speed.
Implement a better algorithm in Python, maybe you can even outperform
the dumb code.
[color=blue][color=green]
>> 1. How do I add docstrings to my module?[/color]
>
> You mean a docstring on the module object itself? In initmodule() make a
> string object with the contents of the docstring and assign it to
> module.__doc__ . There might be a more idiomatic way, but I can't think
> of it without trawling through the docs.[/color]

Use Py_InitModule3( ) instead.

Thomas

[jbrewer]

Re: Python C module questions

>You probably shouldn't post such large pieces of code to the list.

OK.
[color=blue]
>You mean a docstring on the module object itself?[/color]

Actually, I meant docstrings to the module and the functions, objects,
methods, whatever else in the module. My code was derived from the
Python Cookbook, which left that part out (rather important for
building real C modules).
[color=blue]
>You should give up C with a dumb algorithm running at fast speed.
>Implement a better algorithm in Python, maybe you can even outperform
>the dumb code.[/color]

That would be the next step if I needed even more speed, but the better
algorithm here would be to use KDTrees, which would be overkill (and
would require lots of development time). The C brute force
implementation runs plenty fast for me. It took only took a few hours
to implement and yielded a factor of 7 performance increase. I was
mostly interested in feedback on whether I had done things in a
properly efficient way (there are so many Python list / tuple /
sequence functions). I also find it hard to believe that there's no
standard Python function for converting sequences of one object to
arrays in C (a friend mentioned that Ruby's C API has this).

Another question: how does the distutils package handle version
upgrades? Say for example I find some bugs in my C code and need to
recompile it, will it just overwrite what's present in the
site-packages directory?

Jeremy

[Robert Kern]

Re: Python C module questions

jbrewer wrote:

[I wrote:][color=blue][color=green]
>>You probably shouldn't post such large pieces of code to the list.[/color]
>
> OK.[/color]

BTW, please attribute your quotes.

[Still me:][color=blue][color=green]
>>You mean a docstring on the module object itself?[/color]
>
> Actually, I meant docstrings to the module and the functions, objects,
> methods, whatever else in the module. My code was derived from the
> Python Cookbook, which left that part out (rather important for
> building real C modules).[/color]

In the method definition tables for top-level functions and type methods
there is a place for the docstrings. You had them in your code.

Please read the documentation on writing extensions.

http://docs.python.org/ext/ext.html

Specifically,

http://docs.python.org/ext/methodTable.html

http://docs.python.org/ext/node22.html

[Thomas Heller wrote:][color=blue][color=green]
>>You should give up C with a dumb algorithm running at fast speed.
>>Implement a better algorithm in Python, maybe you can even outperform
>>the dumb code.[/color]
>
> That would be the next step if I needed even more speed, but the better
> algorithm here would be to use KDTrees, which would be overkill (and
> would require lots of development time).[/color]

Not for you, it won't.

google('kdtree python')
[color=blue]
> The C brute force
> implementation runs plenty fast for me. It took only took a few hours
> to implement and yielded a factor of 7 performance increase. I was
> mostly interested in feedback on whether I had done things in a
> properly efficient way (there are so many Python list / tuple /
> sequence functions). I also find it hard to believe that there's no
> standard Python function for converting sequences of one object to
> arrays in C (a friend mentioned that Ruby's C API has this).
>
> Another question: how does the distutils package handle version
> upgrades? Say for example I find some bugs in my C code and need to
> recompile it, will it just overwrite what's present in the
> site-packages directory?[/color]

It will overwrite the files that have been changed.

--
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

[Thomas Heller]

Re: Python C module questions

jeremy.d.brewer @gmail.com writes:

[...][color=blue]
> /* convert sequences to tuples if necessary */
> ra1 = PySequence_Fast (ra1, "ra1 must be iterable");
> if (ra1 == NULL) {
> return NULL;
> }
>
> dec1 = PySequence_Fast (dec1, "dec1 must be iterable");
> if (dec1 == NULL) {
> return NULL;
> }[/color]

You leak a refcount to ra1 here in the case the the second
PySequence_Fast fails.

[...][color=blue]
> /* allocate memory for C arrays */
> r1 = (double *) malloc(len1 * sizeof(double)) ;
> if (r1 == NULL) {
> Py_DECREF(ra1);
> Py_DECREF(dec1) ;
> Py_DECREF(ra2);
> Py_DECREF(dec2) ;
> return PyErr_NoMemory( );
> }
>
> d1 = (double *) malloc(len1 * sizeof(double)) ;
> if (d1 == NULL) {
> Py_DECREF(ra1);
> Py_DECREF(dec1) ;
> Py_DECREF(ra2);
> Py_DECREF(dec2) ;
> return PyErr_NoMemory( );
> }[/color]
and so on, and so on.

You should probably better start your code initializing all local vars
to NULL, like this:

PyObject *ra1 = NULL, *dec1 = NULL, ...
char *r1 = NULL, char *d1 = NULL, ...;

Then in the body of the function

d1 = (double *) malloc(len1 * sizeof(double)) ;
if (d1 == NULL) {
PyErr_NoMemory( );
goto error;
}

and have a cleanup section at the end:

error:
if (d1) free(d1);
....
Py_XDECREF(ra1) ;
Py_XDECREF(dec1 );
return NULL;
}

Reading the sources for a few builtin modules of the Python sources
itself should give you also an idea.
[color=blue]
> I also find it hard to believe that there's no
> standard Python function for converting sequences of one object to
> arrays in C (a friend mentioned that Ruby's C API has this).[/color]

The builtin array module does help with this, you can build the
C-compatible arrays in Python, pass them to your extension, and access
them using the buffer api.

Thomas