On 06/17/11 00:36, Ondrej Certik wrote:
On Thu, Jun 16, 2011 at 3:02 AM, Robert Cimrmancimr...@ntc.zcu.cz wrote:
On 06/15/11 22:21, Ondrej Certik wrote:
On Wed, Jun 15, 2011 at 5:36 AM, Robert Cimrmancimr...@ntc.zcu.cz wrote:
Hi
any Cython expert around (Ondrej?).
Our build system already supports cython extension modules and I succeeded in using it for simple stuff, like the sfepy/fem/physics/extmods/cdft.pyx, see [1].
Now I would like to proceed further, but am only a beginner, so an advice would come handy:
I need to map somehow numpy arrays to the FMField struct defined in sfepy/fem/extmods/fmfield.h to be able to export C functions (e.g. term evaluation functions) to Python and pass numpy arrays in. What is the best way to do this? Also being able to call Python from C would be nice (e.g. some numpy array functions).
Here is an example how to get to the pointer of the C array in numpy:
def atom_lda2(int Z, ndarray[double, mode="c"] R not None, ndarray[double, mode="c"] Rp not None, ndarray[int, mode="c"] no not None, ndarray[int, mode="c"] lo not None, ndarray[double, mode="c"] fo not None, ndarray[double, mode="c"] V_tot not None, ndarray[double, mode="c"] ks_energies not None, double reigen_eps, int reigen_max_iter, double mixing_eps, double mixing_alpha, int mixing_max_iter): cdef int n = len(R) cdef int n_orb = len(no) assert len(Rp) == n assert len(V_tot) == n assert len(lo) == n_orb assert len(fo) == n_orb assert len(ks_energies) == n_orb cdef double E_tot cdef ndarray[double, mode="c"] density = empty(n) cdef ndarray[double, ndim=2, mode="fortran"] orbitals = empty((n, n_orb), order="F")
rdirac.c_atom_lda2(&n,&n_orb,&Z,&R[0],&Rp[0], &no[0],&lo[0],&fo[0],&V_tot[0],&ks_energies[0],&E_tot, &density[0],&orbitals[0, 0], &reigen_eps,&reigen_max_iter,&mixing_eps,&mixing_alpha, &mixing_max_iter) return E_tot, density, orbitalsAnd C the signature of the rdirac.c_atom_lda2 function is:
cdef extern void c_atom_lda2(int *n, int *n_orb, int *Z, double *R, double *Rp, int *no, int *lo, double *fo, double *V_tot, double *ks_energies, double *E_tot, double *density, double *orbitals, double *reigen_eps, int *reigen_max_iter, double *mixing_eps, double *mixing_alpha, int *mixing_max_iter)
This function is implemented directly in Fortran, but in your case, you can implement it in C. In C, you might want to pass some variables directly (not by reference), so instead of doing "&n" above, you would just do "n". Makes sense?
ok, so&array[0], when used on something declared as ndarray, gives the underlying pointer to the first item. nice!
That's right, it gives the pointer with the right type. Other option is to use the internal data structures of NumPy, but I like the above approach, because Cython handles the rest.
Well, I am considering to get rid of FMField completely later, when/if the numpy refactoring effort by the Enthought guys succeds. For the moment, I will stick with it and try as you suggest.
So here is the field:
typedef struct FMField { int32 nCell; int32 nLev; int32 nRow; int32 nCol; float64 *val0; float64 *val; int32 nAlloc; int32 cellSize;
int32 offset; int32 nColFull;
int32 stride; } FMField;
So tell me more how you want to use it. Who owns the memory? The easiest (see above) is when Python allocates all the memory, and then you just pass a reference to it to C/Fortran, and it fills it in.
Yes, that's what we do now with swig too - allocate the memory in Python (using zeros, empty, whatever...), and then make a FMField that points to that memory (required C-contiguous for simplicity) and has properly initialized the shape attributes etc. FMField is essentially a 4D array, (n_elements, n_quad_points, n_row, n_col) that enables setting a particular element data to be active (by setting the *val pointer) and performing basic linear algebra on it - for each quadrature point it's a matrix.
Right. So can you just pass in there the pointer to the NumPy array, using the above approach?
I could, I think, write a driving routine for each C function that converts arrays to FMFields, but prefer not to :) So I will probably try to expose FMField to Python, so that I could call, for example, a C function
c_fun(FMField *a, FMField *b);
from Python directly as
c_fun(FMField(array_a), FMField(array_b))
where array_{a, b} are numpy arrays. Or is it even possible to avoid the explicit conversion of array arguments, like we have currently with swig typemaps?
The other thing I need is to be able to pass e.g. a 2D array to a function with signature like function(int *data, int n_row, int n_col); so that the pointer goes to *data argument, and the shape to n_row and n_col.
Just do&some_numpy_array[0,0]. See the "orbitals" array above. You need to make sure the ordering is correct. I use Fortran ordering, you will use C ordering. Also in all these things you have to make sure things are contiguous.
The best way to do it is to declare all the numpy arrays as C (or Fortran) contiguous, and Cython raises an exception if the user passes in something else. That way errors will never propagate silently.
Yes, that's cool, that Cython can check contiguousness and "not None".
Thanks for you response!
Code the above up in some branch, and ping me if you need further help resolving some technical stuff.
Sure, I will make a branch on github.
Thanks, r.