Menu
▾
▴
SF.net SVN: matplotlib: [4655] trunk/py4science/examples
|
From: <jd...@us...> - 2007-12-06 20:01:55
|
Revision: 4655
http://matplotlib.svn.sourceforge.net/matplotlib/?rev=4655&view=rev
Author: jdh2358
Date: 2007-12-06 12:01:48 -0800 (Thu, 06 Dec 2007)
Log Message:
-----------
added pyrex simple sum demo
Modified Paths:
--------------
trunk/py4science/examples/lotka_volterra.py
Added Paths:
-----------
trunk/py4science/examples/pyrex/
trunk/py4science/examples/pyrex/Makefile
trunk/py4science/examples/pyrex/c_numpy.pxd
trunk/py4science/examples/pyrex/c_python.pxd
trunk/py4science/examples/pyrex/numpyx.pyx
trunk/py4science/examples/pyrex/run_test.py
trunk/py4science/examples/pyrex/setup.py
trunk/py4science/examples/pyrex/sums.pyx
trunk/py4science/examples/pyrex/sums_test.py
Modified: trunk/py4science/examples/lotka_volterra.py
===================================================================
--- trunk/py4science/examples/lotka_volterra.py 2007-12-06 19:18:49 UTC (rev 4654)
+++ trunk/py4science/examples/lotka_volterra.py 2007-12-06 20:01:48 UTC (rev 4655)
@@ -34,24 +34,22 @@
f = y[:,1] # extract the foxes vector
p.figure()
-p.plot(t, r, label='rabbits')
-p.plot(t, f, label='foxes')
+p.subplots_adjust(hspace=0.3)
+p.subplot(211)
+p.plot(t, r, color='blue', label='rabbits', lw=2)
+p.plot(t, f, color='green', label='foxes', lw=2)
p.xlabel('time (years)')
p.ylabel('population')
-p.title('population trajectories')
+p.title('population trajectories and phase plane')
p.grid()
p.legend()
-p.savefig('lotka_volterra.png', dpi=150)
-p.savefig('lotka_volterra.eps')
-p.figure()
-p.plot(r, f)
+p.subplot(212, aspect='equal')
+p.plot(r, f, 'k', lw=2)
p.xlabel('rabbits')
p.ylabel('foxes')
-p.title('phase plane')
-
# make a direction field plot with quiver
rmax = 1.1 * r.max()
fmax = 1.1 * f.max()
@@ -65,13 +63,13 @@
dR = dr(R, F)
dF = df(R, F)
-p.contour(R, F, dR, levels=[0], linewidths=3, colors='black')
-p.contour(R, F, dF, levels=[0], linewidths=3, colors='black')
+CSR = p.contour(R, F, dR, levels=[0], linewidths=3, colors='blue')
+CSF = p.contour(R, F, dF, levels=[0], linewidths=3, colors='green')
+
p.ylabel('foxes')
-p.title('trajectory, direction field and null clines')
-p.savefig('lotka_volterra_pplane.png', dpi=150)
-p.savefig('lotka_volterra_pplane.eps')
+p.savefig('lotka_volterra.png', dpi=150)
+p.savefig('lotka_volterra.eps')
p.show()
Added: trunk/py4science/examples/pyrex/Makefile
===================================================================
--- trunk/py4science/examples/pyrex/Makefile (rev 0)
+++ trunk/py4science/examples/pyrex/Makefile 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,9 @@
+all:
+ python setup.py build_ext --inplace
+
+test: all
+ python sums_test.py
+
+.PHONY: clean
+clean:
+ rm -rf *~ *.so *.c *.o build
Added: trunk/py4science/examples/pyrex/c_numpy.pxd
===================================================================
--- trunk/py4science/examples/pyrex/c_numpy.pxd (rev 0)
+++ trunk/py4science/examples/pyrex/c_numpy.pxd 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,125 @@
+# :Author: Travis Oliphant
+
+cdef extern from "numpy/arrayobject.h":
+
+ cdef enum NPY_TYPES:
+ NPY_BOOL
+ NPY_BYTE
+ NPY_UBYTE
+ NPY_SHORT
+ NPY_USHORT
+ NPY_INT
+ NPY_UINT
+ NPY_LONG
+ NPY_ULONG
+ NPY_LONGLONG
+ NPY_ULONGLONG
+ NPY_FLOAT
+ NPY_DOUBLE
+ NPY_LONGDOUBLE
+ NPY_CFLOAT
+ NPY_CDOUBLE
+ NPY_CLONGDOUBLE
+ NPY_OBJECT
+ NPY_STRING
+ NPY_UNICODE
+ NPY_VOID
+ NPY_NTYPES
+ NPY_NOTYPE
+
+ cdef enum requirements:
+ NPY_CONTIGUOUS
+ NPY_FORTRAN
+ NPY_OWNDATA
+ NPY_FORCECAST
+ NPY_ENSURECOPY
+ NPY_ENSUREARRAY
+ NPY_ELEMENTSTRIDES
+ NPY_ALIGNED
+ NPY_NOTSWAPPED
+ NPY_WRITEABLE
+ NPY_UPDATEIFCOPY
+ NPY_ARR_HAS_DESCR
+
+ NPY_BEHAVED
+ NPY_BEHAVED_NS
+ NPY_CARRAY
+ NPY_CARRAY_RO
+ NPY_FARRAY
+ NPY_FARRAY_RO
+ NPY_DEFAULT
+
+ NPY_IN_ARRAY
+ NPY_OUT_ARRAY
+ NPY_INOUT_ARRAY
+ NPY_IN_FARRAY
+ NPY_OUT_FARRAY
+ NPY_INOUT_FARRAY
+
+ NPY_UPDATE_ALL
+
+ cdef enum defines:
+ # Note: as of Pyrex 0.9.5, enums are type-checked more strictly, so this
+ # can't be used as an integer.
+ NPY_MAXDIMS
+
+ ctypedef struct npy_cdouble:
+ double real
+ double imag
+
+ ctypedef struct npy_cfloat:
+ double real
+ double imag
+
+ ctypedef int npy_intp
+
+ ctypedef extern class numpy.dtype [object PyArray_Descr]:
+ cdef int type_num, elsize, alignment
+ cdef char type, kind, byteorder, hasobject
+ cdef object fields, typeobj
+
+ ctypedef extern class numpy.ndarray [object PyArrayObject]:
+ cdef char *data
+ cdef int nd
+ cdef npy_intp *dimensions
+ cdef npy_intp *strides
+ cdef object base
+ cdef dtype descr
+ cdef int flags
+
+ ctypedef extern class numpy.flatiter [object PyArrayIterObject]:
+ cdef int nd_m1
+ cdef npy_intp index, size
+ cdef ndarray ao
+ cdef char *dataptr
+
+ ctypedef extern class numpy.broadcast [object PyArrayMultiIterObject]:
+ cdef int numiter
+ cdef npy_intp size, index
+ cdef int nd
+ # These next two should be arrays of [NPY_MAXITER], but that is
+ # difficult to cleanly specify in Pyrex. Fortunately, it doesn't matter.
+ cdef npy_intp *dimensions
+ cdef void **iters
+
+ object PyArray_ZEROS(int ndims, npy_intp* dims, NPY_TYPES type_num, int fortran)
+ object PyArray_EMPTY(int ndims, npy_intp* dims, NPY_TYPES type_num, int fortran)
+ dtype PyArray_DescrFromTypeNum(NPY_TYPES type_num)
+ object PyArray_SimpleNew(int ndims, npy_intp* dims, NPY_TYPES type_num)
+ int PyArray_Check(object obj)
+ object PyArray_ContiguousFromAny(object obj, NPY_TYPES type,
+ int mindim, int maxdim)
+ npy_intp PyArray_SIZE(ndarray arr)
+ npy_intp PyArray_NBYTES(ndarray arr)
+ void *PyArray_DATA(ndarray arr)
+ object PyArray_FromAny(object obj, dtype newtype, int mindim, int maxdim,
+ int requirements, object context)
+ object PyArray_FROMANY(object obj, NPY_TYPES type_num, int min,
+ int max, int requirements)
+ object PyArray_NewFromDescr(object subtype, dtype newtype, int nd,
+ npy_intp* dims, npy_intp* strides, void* data,
+ int flags, object parent)
+
+ void PyArray_ITER_NEXT(flatiter it)
+
+ void import_array()
Added: trunk/py4science/examples/pyrex/c_python.pxd
===================================================================
--- trunk/py4science/examples/pyrex/c_python.pxd (rev 0)
+++ trunk/py4science/examples/pyrex/c_python.pxd 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,20 @@
+# -*- Mode: Python -*- Not really, but close enough
+
+# Expose as much of the Python C API as we need here
+
+cdef extern from "stdlib.h":
+ ctypedef int size_t
+
+cdef extern from "Python.h":
+ ctypedef int Py_intptr_t
+ void* PyMem_Malloc(size_t)
+ void* PyMem_Realloc(void *p, size_t n)
+ void PyMem_Free(void *p)
+ char* PyString_AsString(object string)
+ object PyString_FromString(char *v)
+ object PyString_InternFromString(char *v)
+ int PyErr_CheckSignals()
+ object PyFloat_FromDouble(double v)
+ void Py_XINCREF(object o)
+ void Py_XDECREF(object o)
+ void Py_CLEAR(object o) # use instead of decref
Added: trunk/py4science/examples/pyrex/numpyx.pyx
===================================================================
--- trunk/py4science/examples/pyrex/numpyx.pyx (rev 0)
+++ trunk/py4science/examples/pyrex/numpyx.pyx 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,128 @@
+# -*- Mode: Python -*- Not really, but close enough
+
+cimport c_python
+cimport c_numpy
+import numpy
+
+# Numpy must be initialized
+c_numpy.import_array()
+
+def print_array_info(c_numpy.ndarray arr):
+ cdef int i
+
+ print '-='*10
+ print 'printing array info for ndarray at 0x%0lx'%(<c_python.Py_intptr_t>arr,)
+ print 'print number of dimensions:',arr.nd
+ print 'address of strides: 0x%0lx'%(<c_python.Py_intptr_t>arr.strides,)
+ print 'strides:'
+ for i from 0<=i<arr.nd:
+ # print each stride
+ print ' stride %d:'%i,<c_python.Py_intptr_t>arr.strides[i]
+ print 'memory dump:'
+ print_elements( arr.data, arr.strides, arr.dimensions,
+ arr.nd, sizeof(double), arr.dtype )
+ print '-='*10
+ print
+
+
+
+def sum_elements(c_numpy.ndarray arr):
+ cdef int i
+ cdef double x, val
+
+ x = 0.
+ val = 0.
+ for i from 0<=i<arr.dimensions[0]:
+ val = (<double*>(arr.data + i*arr.strides[0]))[0]
+ x = x + val
+
+ return x
+
+
+def scale_elements(int N):
+ cdef int i
+ cdef double x, val
+
+ x = 0.
+ val = 0.
+ for i from 0<=i<N:
+ val = 2.5 * i
+ x = x + val
+ return x
+
+
+cdef print_elements(char *data,
+ c_python.Py_intptr_t* strides,
+ c_python.Py_intptr_t* dimensions,
+ int nd,
+ int elsize,
+ object dtype):
+ cdef c_python.Py_intptr_t i,j
+ cdef void* elptr
+
+ if dtype not in [numpy.dtype(numpy.object_),
+ numpy.dtype(numpy.float64)]:
+ print ' print_elements() not (yet) implemented for dtype %s'%dtype.name
+ return
+
+ if nd ==0:
+ if dtype==numpy.dtype(numpy.object_):
+ elptr = (<void**>data)[0] #[0] dereferences pointer in Pyrex
+ print ' ',<object>elptr
+ elif dtype==numpy.dtype(numpy.float64):
+ print ' ',(<double*>data)[0]
+ elif nd == 1:
+ for i from 0<=i<dimensions[0]:
+ if dtype==numpy.dtype(numpy.object_):
+ elptr = (<void**>data)[0]
+ print ' ',<object>elptr
+ elif dtype==numpy.dtype(numpy.float64):
+ print ' ',(<double*>data)[0]
+ data = data + strides[0]
+ else:
+ for i from 0<=i<dimensions[0]:
+ print_elements(data, strides+1, dimensions+1, nd-1, elsize, dtype)
+ data = data + strides[0]
+
+
+def test_methods(c_numpy.ndarray arr):
+ """Test a few attribute accesses for an array.
+
+ This illustrates how the pyrex-visible object is in practice a strange
+ hybrid of the C PyArrayObject struct and the python object. Some
+ properties (like .nd) are visible here but not in python, while others
+ like flags behave very differently: in python flags appears as a separate,
+ object while here we see the raw int holding the bit pattern.
+
+ This makes sense when we think of how pyrex resolves arr.foo: if foo is
+ listed as a field in the c_numpy.ndarray struct description, it will be
+ directly accessed as a C variable without going through Python at all.
+ This is why for arr.flags, we see the actual int which holds all the flags
+ as bit fields. However, for any other attribute not listed in the struct,
+ it simply forwards the attribute lookup to python at runtime, just like
+ python would (which means that AttributeError can be raised for
+ non-existent attributes, for example)."""
+
+ print 'arr.any() :',arr.any()
+ print 'arr.nd :',arr.nd
+ print 'arr.flags :',arr.flags
+
+def test():
+ """this function is pure Python"""
+ arr1 = numpy.array(-1e-30,dtype=numpy.float64)
+ arr2 = numpy.array([1.0,2.0,3.0],dtype=numpy.float64)
+
+ arr3 = numpy.arange(9,dtype=numpy.float64)
+ arr3.shape = 3,3
+
+ four = 4
+ arr4 = numpy.array(['one','two',3,four],dtype=numpy.object_)
+
+ arr5 = numpy.array([1,2,3]) # int types not (yet) supported by print_elements
+
+ for arr in [arr1,arr2,arr3,arr4,arr5]:
+ print_array_info(arr)
+
+
+ arr = numpy.arange(10.0)
+ print 'sum el', sum_elements(arr)
Added: trunk/py4science/examples/pyrex/run_test.py
===================================================================
--- trunk/py4science/examples/pyrex/run_test.py (rev 0)
+++ trunk/py4science/examples/pyrex/run_test.py 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,3 @@
+#!/usr/bin/env python
+from numpyx import test
+test()
Property changes on: trunk/py4science/examples/pyrex/run_test.py
___________________________________________________________________
Name: svn:executable
+ *
Added: trunk/py4science/examples/pyrex/setup.py
===================================================================
--- trunk/py4science/examples/pyrex/setup.py (rev 0)
+++ trunk/py4science/examples/pyrex/setup.py 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,42 @@
+#!/usr/bin/env python
+"""Install file for example on how to use Pyrex with Numpy.
+
+For more details, see:
+http://www.scipy.org/Cookbook/Pyrex_and_NumPy
+http://www.scipy.org/Cookbook/ArrayStruct_and_Pyrex
+"""
+
+from distutils.core import setup
+from distutils.extension import Extension
+
+# Make this usable by people who don't have pyrex installed (I've committed
+# the generated C sources to SVN).
+try:
+ from Pyrex.Distutils import build_ext
+ has_pyrex = True
+except ImportError:
+ has_pyrex = False
+
+import numpy
+
+# Define a pyrex-based extension module, using the generated sources if pyrex
+# is not available.
+if has_pyrex:
+ pyx_sources = ['sums.pyx']
+ cmdclass = {'build_ext': build_ext}
+else:
+ pyx_sources = ['numpyx.c']
+ cmdclass = {}
+
+
+pyx_ext = Extension('sums',
+ pyx_sources,
+ include_dirs = [numpy.get_include()])
+
+# Call the routine which does the real work
+setup(name = 'sums',
+ description = 'Small example on using Pyrex to write a Numpy extension',
+ url = 'http://www.scipy.org/Cookbook/Pyrex_and_NumPy',
+ ext_modules = [pyx_ext],
+ cmdclass = cmdclass,
+ )
Added: trunk/py4science/examples/pyrex/sums.pyx
===================================================================
--- trunk/py4science/examples/pyrex/sums.pyx (rev 0)
+++ trunk/py4science/examples/pyrex/sums.pyx 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,39 @@
+# -*- Mode: Python -*- Not really, but close enough
+
+cimport c_python
+cimport c_numpy
+import numpy
+
+# Numpy must be initialized
+c_numpy.import_array()
+
+def sum_elements(c_numpy.ndarray arr):
+ cdef int i
+ cdef double x, val
+
+ x = 0.
+ val = 0.
+ for i from 0<=i<arr.dimensions[0]:
+ val = (<double*>(arr.data + i*arr.strides[0]))[0]
+ x = x + val
+
+ return x
+
+
+
+def sum_elements2(c_numpy.ndarray arr):
+ cdef int i
+ cdef double x, val
+
+ arr = numpy.asarray(arr, numpy.float_)
+
+ if arr.nd!=1:
+ raise RuntimeError('only 1D arrays supported; found shape=%s'%str(arr.shape))
+ assert(arr.nd==1)
+ x = 0.
+ val = 0.
+ for i from 0<=i<arr.dimensions[0]:
+ val = (<double*>(arr.data + i*arr.strides[0]))[0]
+ x = x + val
+
+ return x
Added: trunk/py4science/examples/pyrex/sums_test.py
===================================================================
--- trunk/py4science/examples/pyrex/sums_test.py (rev 0)
+++ trunk/py4science/examples/pyrex/sums_test.py 2007-12-06 20:01:48 UTC (rev 4655)
@@ -0,0 +1,29 @@
+import numpy
+import sums
+
+def sumpy(arr):
+
+ total = 0.
+ for val in x:
+ total += val
+ return total
+
+x = numpy.arange(10)
+y = numpy.random.rand(10,10)
+
+print 'sum(x)', sums.sum_elements(x)
+print 'sum2(x)', sums.sum_elements2(x)
+print 'sum(y)', sums.sum_elements(y)
+#print 'sum2(y)', sums.sum_elements2(y)
+
+x = numpy.arange(1e6)
+import time
+start = time.time()
+s1 = sums.sum_elements2(x)
+now = time.time()
+print 'pyrex time', now - start
+
+start = time.time()
+s1 = sumpy(x)
+now = time.time()
+print 'python time', now - start
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
