Re: [Swig-user] How to output 2D C arrays to Python as NumPy arrays?

[Alan W. I. <ir...@be...>]

On 2012-01-10 11:05-0800 Alan W. Irwin wrote:

> I already have a solution for outputting 2D C arrays to Python as
> NumPy arrays.  That method is to treat the NumPy arrays as input
> arrays from the SWIG perspective, but simply allow the C code to
> overwrite the contents of those preexisting NumPy arrays that appear
> in Python argument lists. Obviously that is a massive hack that
> violates the sanctity of function arguments that Python users expect.
> Therefore, I would like to replace what I have got with something more
> appropriate to Python, but I need some help to do that.

To finish off this topic, I got that help from Nickolas Fotopoulos who
referred me to the extremely useful numpy.i documentation at
http://docs.scipy.org/doc/numpy/reference/swig.html, and I found the
numpy.i file itself (which for the record is available at
https://github.com/numpy/numpy/blob/55472ca381c472d86ebb899cb2d786e8e9be8ba5/doc/swig/numpy.i)
was also extremely useful.  However, the 2D ARGOUT numpy.i typemaps
were somewhat limited for my needs.  For example, I wanted to create a
1-D NumPy array of strings and fill that with the output char** matrix
results rather than a 2-D NumPY array of single characters. Also, for
this particular use case (unlike the use case I have discussed in
another current thread where the code block order is so important) the
created NumPy arrays had known maximum dimensions of reasonable size.
So in this case, the approach I took was to create the arrays with
those known maximum sizes, and then after the call to the C function,
resize those arrays to the correct shape (known only at that stage)
before returning them to the Python environment.

Here are the resulting typemaps for the two cases I discussed before
which were previously handled by overwriting the contents
of pre-existing NumPy arrays that appeared in argument lists.

// Output 1D NumPy array of strings corresponding to contents of output
// 2D char ** array from C routine.
%typemap( in, numinputs = 0 ) char **OutMatrixCk7( PyObject * array = NULL )
{
     // Must take an initial guess at dimension because this code must
     // be executed before the actual dimension is known as a result of
     // the c_ephcom_read_constants call.
     // N.B. If EPHCOM_MAXCON is ever exceeded by the number of
     // constants/values in an ephemeris (not likely in the immediate
     // future since the latest JPL ephemeris DE423 only has 222
     // constants at the moment) the core ephcom library would not be
     // able to even read or write a binary ephemeris until this
     // ephcom.h header macro was increased.
     Alen = EPHCOM_MAXCON;
     npy_intp     dims[] = { Alen };
     int          i;
     PyArray_Descr* dtype_desc = PyArray_DescrFromType( NPY_STRING );
     // Force equivalent of zeros(EPHCOM_MAXCON, dtype='SX'), where
     // X is the numerical value of EPHCOM_MAX_CNAM_LENGTH;
     dtype_desc->elsize = EPHCOM_MAX_CNAM_LENGTH;
     array = PyArray_Zeros( 1, dims, dtype_desc, 0 );
     $1    = (char **) malloc( sizeof ( char* ) * Alen );
     for ( i = 0; i < Alen; i++ )
     {
         $1[i] = ( ( (PyArrayObject *) array )->data + i * ( (PyArrayObject *) array )->strides[0] );
     }
}
%typemap( argout ) char **OutMatrixCk7
{
     npy_intp     dims[] = { result };
     PyArray_Dims newshape_location;
     PyArray_Dims * newshape = &newshape_location;
     newshape->len = 1;
     newshape->ptr = dims;
     PyObject * dummy;

     // resize array$argnum to actual size needed.
     // Follow directions at http://www.mail-archive.com/num...@sc.../msg13013.html
     dummy = PyArray_Resize( (PyArrayObject *) array$argnum, newshape, 0, NPY_CORDER );
     if ( !dummy )
         SWIG_fail;
     Py_DECREF( dummy );
     $result = SWIG_Python_AppendOutput( $result, array$argnum );
}
%typemap( freearg ) char **OutMatrixCk7 { free( $1 );}

// Output 2D array whose first dimension is determined by previous
// vector dimension and whose second dimension is EPHCOM_MAXPVCOORD
%typemap( in, numinputs = 0 ) double **OutMatrixCkCk6( PyObject * array = NULL )
{
     // Must take an initial guess at dimension because swig deploys
     // this code before Alen is defined.
     // N.B. If EPHCOM_MAXOBJECTS is ever exceeded by the number of
     // solar system objects that have interpolation data in an
     // ephemeris (say when detailed asteroid motions are included) the
     // core ephcom library would not be able to even read or write a
     // binary ephemeris until this ephcom.h header macro was
     // increased.
     int      i, size, guessed_dimension = EPHCOM_MAXOBJECTS;
     npy_intp dims[] = { guessed_dimension, EPHCOM_MAXPVCOORD };
     array = PyArray_SimpleNew( 2, dims, NPY_DOUBLE );
     if ( !array )
         return NULL;

     size = (int) ( sizeof ( double ) * EPHCOM_MAXPVCOORD );
     $1   = (double **) malloc( sizeof ( double * ) * guessed_dimension );
     for ( i = 0; i < guessed_dimension; i++ )
         $1[i] = (double *) ( ( (PyArrayObject *) array )->data + i * size );
}
%typemap( freearg ) double **OutMatrixCkCk6 { free( $1 ); }
%typemap( argout ) double **OutMatrixCkCk6
{
     npy_intp     dims[] = { Alen, EPHCOM_MAXPVCOORD };
     PyArray_Dims newshape_location;
     PyArray_Dims * newshape = &newshape_location;
     newshape->len = 2;
     newshape->ptr = dims;
     PyObject * dummy;

     // resize array$argnum to actual size needed.
     // Follow directions at http://www.mail-archive.com/num...@sc.../msg13013.html
     dummy = PyArray_Resize( (PyArrayObject *) array$argnum, newshape, 0, NPY_CORDER );
     if ( !dummy )
         SWIG_fail;
     Py_DECREF( dummy );
     $result = SWIG_Python_AppendOutput( $result, array$argnum );
}

To see how these typemaps fit together with the ephcomc API and other
typemaps for that API, see

http://timeephem.svn.sourceforge.net/viewvc/timeephem/trunk/jpl_ephemerides/ephcom2/bindings/swig_support/ephcomcapi.i?view=log

and

http://timeephem.svn.sourceforge.net/viewvc/timeephem/trunk/jpl_ephemerides/ephcom2/bindings/python/ephcomcmodule.i?view=log

As a result of these recent changes I can call functions in
this API from python as follows without relying on any overwriting of
the contents of NumPy arguments to these functions:

actual_number, cnames, values, sss = ephcomc.ephcom_read_constants("../../jpleph_bindata/de405/JPLEPH")

where cnames (known as char **OutMatrixCk7 above), values, and sss are
all output arrays that are filled with constant data from a particular
JPL ephemeris of the planets.  Furthermore, the shapes of these arrays
are either known in advance (sss always has just 3 elements) or are
determined by the value (e.g., actual_number) returned by the call to
the C library routine, ephcom_read_constants.

et2 = array([2305447.5, 0.])
nlist = 15
list = ones(nlist, dtype=int32)
return_code, pv = ephcomc.ephcom_interpolate_list(
"../../jpleph_bindata/de405/JPLEPH", et2, 0, list )

where in this case pv (known as double **OutMatrixCkCk6 above)
contains data interpolated at et2 from the JPL ephemeris for the solar
system objects whose corresponding list values are unity (i.e., all of
the first nlist objects).  pv has a 2D shape with the first dimension
determined by the size of list and a known second dimension of
EPHCOM_MAXPVCOORD = 6.

Following up on the previous topic of not being able to control the
order of the code blocks associated with typemaps, there is no way to
pass dimension information from the code block associated with the
list input array typemap to the codeblock associated with the pv
(a.k.a. double **OutMatrixCkCk6) typemap because the pv codeblock
appears before the list codeblock in the SWIG-generated wrapper code
because of the numinputs=0 qualifier of the former.  As discussed in
the other thread, one workaround in this case would have been to have
a combined typemap for the C equivalent of list and pv.  Instead, for
this particular case I guessed at the size of pv (since that was
fairly convenient to do with known reasonable maximum sizes for the
dimensions) and refined it afterwards once the Alen information was
available in the code concerning the length of list.

I am going to regret all these workarounds if it turns out there
is a way to control the order of the codeblocks in the C wrapper
code generated by SWIG so if somebody could give me the definitive
answer on that question, I would appreciate it.

Alan
__________________________
Alan W. Irwin

Astronomical research affiliation with Department of Physics and Astronomy,
University of Victoria (astrowww.phys.uvic.ca).

Programming affiliations with the FreeEOS equation-of-state
implementation for stellar interiors (freeeos.sf.net); the Time
Ephemerides project (timeephem.sf.net); PLplot scientific plotting
software package (plplot.sf.net); the libLASi project
(unifont.org/lasi); the Loads of Linux Links project (loll.sf.net);
and the Linux Brochure Project (lbproject.sf.net).
__________________________

Linux-powered Science
__________________________