[cig-commits] r6577 - in cs/cigma/trunk/sandbox: . numpy
luis at geodynamics.org
luis at geodynamics.org
Sun Apr 15 23:16:02 PDT 2007
Author: luis
Date: 2007-04-15 23:16:02 -0700 (Sun, 15 Apr 2007)
New Revision: 6577
Added:
cs/cigma/trunk/sandbox/numpy/
cs/cigma/trunk/sandbox/numpy/Makefile
cs/cigma/trunk/sandbox/numpy/gridloop.c
cs/cigma/trunk/sandbox/numpy/gridloop.f
cs/cigma/trunk/sandbox/numpy/langtangen.py
cs/cigma/trunk/sandbox/numpy/numpy_macros.h
cs/cigma/trunk/sandbox/numpy/tests.py
Log:
Added Python functions which demonstrate how to access a numpy array from C
Added: cs/cigma/trunk/sandbox/numpy/Makefile
===================================================================
--- cs/cigma/trunk/sandbox/numpy/Makefile 2007-04-16 06:09:09 UTC (rev 6576)
+++ cs/cigma/trunk/sandbox/numpy/Makefile 2007-04-16 06:16:02 UTC (rev 6577)
@@ -0,0 +1,19 @@
+CC = gcc
+CFLAGS = -O3 -g
+INCLUDES = -I/usr/include/python2.4 -I/usr/lib/python2.4/site-packages/numpy/core/include
+
+all: ext_gridloop_f.so ext_gridloop_c.so
+
+ext_gridloop_f.so: gridloop.f
+ f2py -m ext_gridloop_f -c gridloop.f
+
+gridloop_c.o: gridloop.c
+ $(CC) $(CFLAGS) $(INCLUDES) -c $< -o $@
+
+ext_gridloop_c.so: gridloop_c.o
+ gcc -shared $^ -o $@
+
+clean:
+ rm -f *.so *.o *.pyc
+
+.PHONY: clean
Added: cs/cigma/trunk/sandbox/numpy/gridloop.c
===================================================================
--- cs/cigma/trunk/sandbox/numpy/gridloop.c 2007-04-16 06:09:09 UTC (rev 6576)
+++ cs/cigma/trunk/sandbox/numpy/gridloop.c 2007-04-16 06:16:02 UTC (rev 6577)
@@ -0,0 +1,213 @@
+#include <Python.h> /* Python as seen from C */
+#include <numpy/arrayobject.h> /* numpy as seen from C */
+
+
+static PyObject *
+gridloop1(PyObject *self, PyObject *args)
+{
+ PyArrayObject *a, *xcoor, *ycoor;
+ PyObject *func1, *arglist, *result;
+ int nx, ny, i, j;
+ double *a_ij, *x_i, *y_j;
+
+ /* arguments: a, xcoor, ycoor, func1 */
+ if (!PyArg_ParseTuple(args, "O!O!O!O:gridloop1",
+ &PyArray_Type, &a,
+ &PyArray_Type, &xcoor,
+ &PyArray_Type, &ycoor,
+ &func1))
+ {
+ return NULL; /* PyArg_ParseTuple has raised an exception */
+ }
+
+ /* alternative parsing without checking the pointer types
+ if (!PyArg_ParseTuple(args, "OOOO", &a, &xcoor, &ycoor, &func1))
+ {
+ return NULL;
+ }
+ // */
+
+ if (a->nd != 2 || a->descr->type_num != PyArray_DOUBLE)
+ {
+ PyErr_Format(PyExc_ValueError,
+ "a array is %d-dimensional or not of type float",
+ a->nd);
+ return NULL;
+ }
+
+ nx = a->dimensions[0];
+ ny = a->dimensions[1];
+
+ if (xcoor->nd != 1 || xcoor->descr->type_num != PyArray_DOUBLE ||
+ xcoor->dimensions[0] != nx)
+ {
+ PyErr_Format(PyExc_ValueError,
+ "xcoor array has wrong dimension (%d), type or length (%d)",
+ xcoor->nd, xcoor->dimensions[0]);
+ return NULL;
+ }
+
+ if (ycoor->nd != 1 || ycoor->descr->type_num != PyArray_DOUBLE ||
+ ycoor->dimensions[0] != ny)
+ {
+ PyErr_Format(PyExc_ValueError,
+ "ycoor array has wrong dimension (%d), type or length (%d)",
+ ycoor->nd, ycoor->dimensions[0]);
+ return NULL;
+ }
+
+ if (!PyCallable_Check(func1))
+ {
+ PyErr_Format(PyExc_TypeError, "func1 is not a callable function");
+ return NULL;
+ }
+
+ for (i = 0; i < nx; i++)
+ {
+ for (j = 0; j < ny; j++)
+ {
+ a_ij = (double *)(a->data + i*a->strides[0] + j*a->strides[1]);
+ x_i = (double *)(xcoor->data + i*xcoor->strides[0]);
+ y_j = (double *)(ycoor->data + j*ycoor->strides[0]);
+ // shorter: xcoor[i], ycoor[j] for one-dim arrays
+
+ /* call Python function pointed to by func1: */
+ arglist = Py_BuildValue("(dd)", *x_i, *y_j);
+ result = PyEval_CallObject(func1, arglist);
+ Py_DECREF(arglist);
+ if (result == NULL)
+ {
+ /* exception in func1 */
+ return NULL;
+ }
+
+ *a_ij = PyFloat_AS_DOUBLE(result);
+
+ Py_DECREF(result);
+
+#ifdef DEBUG
+ printf("a[%d,%d]=func1(%g,%g)=%g\n",i,j,*x_i,*y_j,*a_ij);
+#endif
+ }
+ }
+ return Py_BuildValue(""); /* return None */
+
+ /* alternative (return 0 for success):
+ return Py_BuildValue("i", 0);
+ // */
+}
+
+
+#include "numpy_macros.h"
+
+static PyObject *
+gridloop2(PyObject *self, PyObject *args)
+{
+ PyArrayObject *a, *xcoor, *ycoor;
+ int a_dims[2];
+ PyObject *func1, *arglist, *result;
+ int nx, ny, i, j;
+ double *a_ij, *x_i, *y_j;
+
+ /* arguments: xcoor, ycoor, func1 */
+ if (!PyArg_ParseTuple(args, "O!O!O:gridloop2",
+ &PyArray_Type, &xcoor,
+ &PyArray_Type, &ycoor,
+ &func1))
+ {
+ return NULL; /* PyArg_ParseTuple has raised an exception */
+ }
+
+ nx = xcoor->dimensions[0];
+ ny = ycoor->dimensions[0];
+
+ NDIM_CHECK(xcoor,1); TYPE_CHECK(xcoor, PyArray_DOUBLE);
+ NDIM_CHECK(ycoor,1); TYPE_CHECK(ycoor, PyArray_DOUBLE);
+ CALLABLE_CHECK(func1);
+
+ /* create return array */
+ a_dims[0] = nx;
+ a_dims[1] = ny;
+ a = (PyArrayObject *)PyArray_FromDims(2, a_dims, PyArray_DOUBLE);
+ if (a == NULL)
+ {
+ printf("Creating %dx%d array failed\n", a_dims[0], a_dims[1]);
+ return NULL; /* PyArray_FromDims raised an exception */
+ }
+
+ for (i = 0; i < nx; i++)
+ {
+ for (j = 0; j < ny; j++)
+ {
+ arglist = Py_BuildValue("(dd)", IND1(xcoor,i), IND1(ycoor,j));
+ result = PyEval_CallObject(func1, arglist);
+ Py_DECREF(arglist);
+ if (result == NULL)
+ {
+ return NULL; /* exception in func1 */
+ }
+
+ IND2(a,i,j) = PyFloat_AS_DOUBLE(result);
+
+ Py_DECREF(result);
+
+#ifdef DEBUG
+ printf("a[%d,%d]=func1(%g,%g)=%g\n", i, j,
+ IND1(xcoor,i), IND1(ycoor,j), IND2(a,i,j));
+#endif
+ }
+ }
+
+ return PyArray_Return(a);
+}
+
+
+
+
+/*
+ * doc strings
+ */
+
+static char gridloop1_doc[] = "gridloop1(a, xcoor, ycoor, pyfunc)";
+static char gridloop2_doc[] = "a = gridloop2(xcoor, ycoor, pyfunc)";
+static char module_doc[] = \
+ "module ext_gridloop:\n\
+ gridloop1(a, xcoor, ycoor, pyfunc)\n\
+ a = gridloop2(xcoor, ycoor, pyfunc)";
+
+/*
+ * The method table must always be present - it lists the
+ * functions that should be callable from Python:
+ */
+static PyMethodDef ext_gridloop_methods[] = {
+
+ {"gridloop1", /* name of func when called from Python */
+ gridloop1, /* corresponding C function */
+ METH_VARARGS, /* ordinary (not keyword) arguments */
+ gridloop1_doc}, /* doc string for gridloop1 function */
+
+ {"gridloop2", /* name of func when called from Python */
+ gridloop2, /* corresponding C function */
+ METH_VARARGS, /* ordinary (not keyword) arguments */
+ gridloop2_doc}, /* doc string for gridloop1 function */
+
+ {NULL, NULL} /* required ending of the method table */
+};
+
+PyMODINIT_FUNC
+initext_gridloop_c()
+{
+ /*
+ * Assign the name of the module, the name of
+ * the method table, and optionally a module
+ * doc string:
+ */
+
+ Py_InitModule3("ext_gridloop_c", ext_gridloop_methods, module_doc);
+
+ /* without module doc string:
+ Py_InitModule("ext_gridloop", ext_gridloop_methods);
+ // */
+
+ import_array(); /* required numpy initialization */
+}
Added: cs/cigma/trunk/sandbox/numpy/gridloop.f
===================================================================
--- cs/cigma/trunk/sandbox/numpy/gridloop.f 2007-04-16 06:09:09 UTC (rev 6576)
+++ cs/cigma/trunk/sandbox/numpy/gridloop.f 2007-04-16 06:16:02 UTC (rev 6577)
@@ -0,0 +1,20 @@
+ subroutine gridloop2(a, xcoor, ycoor, nx, ny, func1)
+ integer nx, ny
+ real*8 a(0:nx-1,0:ny-1), xcoor(0:nx-1), ycoor(0:ny-1), func1
+ external func1
+Cf2py intent(out) a
+Cf2py intent(in) xcoor
+Cf2py intent(in) ycoor
+Cf2py depend(nx,ny) a
+
+ integer i,j
+ real*8 x,y
+ do j = 0, ny-1
+ y = ycoor(j)
+ do i = 0, nx-1
+ x = xcoor(i)
+ a(i,j) = func1(x,y)
+ end do
+ end do
+ return
+ end
Added: cs/cigma/trunk/sandbox/numpy/langtangen.py
===================================================================
--- cs/cigma/trunk/sandbox/numpy/langtangen.py 2007-04-16 06:09:09 UTC (rev 6576)
+++ cs/cigma/trunk/sandbox/numpy/langtangen.py 2007-04-16 06:16:02 UTC (rev 6577)
@@ -0,0 +1,99 @@
+import numpy
+import ext_gridloop_c as ext_gridloop
+
+from numpy import sin
+from numpy import zeros, arange
+from numpy import newaxis, size
+
+def sequence(min=0.0, max=None, inc=1.0):
+ """See Langtangen (pp. 150-151)"""
+ return arange(min, max + inc/2.0, inc)
+seq = sequence
+
+class Grid2D(object):
+ """From Langtangen's Python book (pg. 377)"""
+ def __init__(self, xmin=0, xmax=1, dx=0.5,
+ ymin=0, ymax=1, dy=0.5):
+
+ self.xcoor = seq(xmin, xmax, dx)
+ self.ycoor = seq(ymin, ymax, dy)
+
+ self.xcoorv = self.xcoor[:,newaxis]
+ self.ycoorv = self.ycoor[newaxis,:]
+
+ def __call__(self, f):
+ """
+ Treat f either as an expression containing x and y
+ or as a standard Python function f(x,y). Evaluate
+ the formula or function for all grid points and
+ return the corresponding two-dimensional array.
+ """
+ if isinstance(f, basestring):
+ # introduce the names x and y such that
+ # a simple eval(f) will work for the arrays directly
+ x = self.xcoorv
+ y = self.ycoorv
+ a = eval(f)
+ else:
+ a = f(self.xcoorv, self.ycoorv)
+ return a
+
+ def gridloop(self, f):
+ f_is_str = isinstance(f, basestring)
+ if f_is_str:
+ f_compiled = compile(f, '<string>', 'eval')
+ lx,ly = size(self.xcoor), size(self.ycoor)
+ _a = zeros((lx,ly), float)
+ for i in xrange(lx):
+ x = self.xcoor[i]
+ for j in xrange(ly):
+ y = self.ycoor[j]
+ if f_is_str:
+ _a[i,j] = eval(f_compiled)
+ else:
+ _a[i,j] = f(x,y)
+ return _a
+
+class Grid2Deff(Grid2D):
+ """Langtangen (pg. 433)"""
+
+
+ def ext_gridloop1(self, f):
+ """Compute a[i,j] = f(xi,yj) in an external routine."""
+
+ # a is made here, sent to the routine, and then returned
+ lx,ly = size(self.xcoor), size(self.ycoor)
+ a = zeros((lx,ly), float)
+
+ # <Fortran-specific initialization...>
+
+ if isinstance(f, basestring):
+ f_compiled = compile(f, '<string>', 'eval')
+ def f_str(x,y):
+ return eval(f_compiled)
+ func = f_str
+ else:
+ func = f
+
+ ext_gridloop.gridloop1(a, self.xcoor, self.ycoor, func)
+
+ return a
+
+
+ def ext_gridloop2(self, f):
+ """Compute a[i,j] = f(xi,yj) in an external routine."""
+
+ # a is made in the external routine
+ if isinstance(f, basestring):
+ f_compiled = compile(f, '<string>', 'eval')
+ def f_str(x,y):
+ return eval(f_compiled)
+ func = f_str
+ else:
+ func = f
+
+ a = ext_gridloop.gridloop2(self.xcoor, self.ycoor, func)
+
+ return a
+
+
Added: cs/cigma/trunk/sandbox/numpy/numpy_macros.h
===================================================================
--- cs/cigma/trunk/sandbox/numpy/numpy_macros.h 2007-04-16 06:09:09 UTC (rev 6576)
+++ cs/cigma/trunk/sandbox/numpy/numpy_macros.h 2007-04-16 06:16:02 UTC (rev 6577)
@@ -0,0 +1,57 @@
+#ifndef __NUMPY_MACROS_H__
+#define __NUMPY_MACROS_H__
+
+/*
+ * This file defines some macros for programming with
+ * numpy arrays in C extension modules.
+ */
+
+/* define some macros for array dimension/type check and subscripting */
+
+#define QUOTE(s) # s /* turn s into string "s" */
+
+#define NDIM_CHECK(a, expected_ndim) \
+ if (a->nd != expected_ndim) { \
+ PyErr_Format(PyExc_ValueError, \
+ "%s array is %d-dimensional, but expected to be %d-dimensional",\
+ QUOTE(a), a->nd, expected_ndim); \
+ return NULL; \
+ }
+
+#define DIM_CHECK(a, dim, expected_length) \
+ if (dim > a->nd) { \
+ PyErr_Format(PyExc_ValueError, \
+ "%s array has no %d dimension (max dim is %d)", \
+ QUOTE(a), dim, a->nd); \
+ return NULL; \
+ } \
+ if (a->dimensions[dim] != expected_length) { \
+ PyErr_Format(PyExc_ValueError, \
+ "%s array has wrong %d-dimension=%d (expected %d)", \
+ QUOTE(a), dim, a->dimensions[dim], expected_length); \
+ return NULL; \
+ }
+
+#define TYPE_CHECK(a, tp) \
+ if (a->descr->type_num != tp) { \
+ PyErr_Format(PyExc_TypeError, \
+ "%s array is not of correct type (%d)", \
+ QUOTE(a), tp); \
+ return NULL; \
+ }
+
+#define CALLABLE_CHECK(func) \
+ if (!PyCallable_Check(func)) { \
+ PyErr_Format(PyExc_TypeError, \
+ "%s is not a callable function", \
+ QUOTE(func)); \
+ return NULL; \
+ }
+
+
+#define IND1(a,i) *((double *)(a->data + i*a->strides[0]))
+#define IND2(a,i,j) *((double *)(a->data + i*a->strides[0] + j*a->strides[1]))
+#define IND3(a,i,j,k) *((double *)(a->data + i*a->strides[0] + j*a->strides[1] + k*a->strides[2]))
+
+
+#endif
Added: cs/cigma/trunk/sandbox/numpy/tests.py
===================================================================
--- cs/cigma/trunk/sandbox/numpy/tests.py 2007-04-16 06:09:09 UTC (rev 6576)
+++ cs/cigma/trunk/sandbox/numpy/tests.py 2007-04-16 06:16:02 UTC (rev 6577)
@@ -0,0 +1,21 @@
+#!/usr/bin/env python
+
+def f1(x,y):
+ print 'x+2*y =', x + 2*y
+ return x + 2*y
+
+def verify1():
+ """Basic test of the extension module."""
+ from numpy import allclose
+ from langtangen import Grid2Deff
+
+ g = Grid2Deff(dx=0.5, dy=1)
+ f_exact = g(f1) # numpy computation
+
+ f = g.ext_gridloop2(f1)
+ print 'f computed by external gridloop2 function:\n', f
+ if allclose(f, f_exact, atol=1.0e-10, rtol=1.0e-12):
+ print 'f is correct'
+
+if __name__ == '__main__':
+ verify1()
More information about the cig-commits
mailing list