Mailman 3 February 2007 - NumPy-Discussion

f2py and Fortran90 gfortran_filename error
by Tyler Hayes 28 Feb '07

28 Feb '07

Hello All: I was suggested to post this question here with the f2py experts from comp.lang.python. I have been able to use the example F77 files suggested in the f2py User Manual, but the problems happened when I tried my own F90 file (attached below. I am able to create a <module_name>.so file using the following method: (1) Created a Fortran90 program matsolve.f90 Note: The program compiles fine and prints the proper output for the simple matrix specified. (2) f2py matsolve.f90 -m matsolve2 -h matsolve2.pyf This created the matsolve2.pyf fine (3) f2py -c matsolve2.pyf --f90exec=/usr/bin/gfortran matsolve.f90 Note: I had to specify the f90exec path as f2py did not automatically find it. This was the only way I could generate a *.so file. The rub: When I import it into Python, I receive the following message: >>> import matsolve2 Traceback (most recent call last): File "<stdin>", line 1, in ? ImportError: ./matsolve2.so: undefined symbol: _gfortran_filename Any suggestions are greatly appreciated. (I almost started smoking again last night! ARGH!) Cheers, t. PS - Below I have attached the F90 program for LU decomposition with a simple test case. ! MATSOLVE.f90 ! ! Start main program PROGRAM MATSOLVE IMPLICIT NONE INTEGER,PARAMETER :: n=3 INTEGER :: i,j REAL,DIMENSION(n) :: x,b REAL,DIMENSION(n,n) :: A,L,U ! Initialize the vectors and matrices with a test case from text ! Using the one given in Appendix A from Thompson. ! Known vector "b" b(1) = 12. b(2) = 11. b(3) = 2. ! Known coefficient matrix "A", and initialize L and U DO i=1,n DO j=1,n L(i,j) = 0. U(i,j) = 0. END DO END DO ! Create matrix A A(1,1) = 3. A(1,2) = -1. A(1,3) = 2. A(2,1) = 1. A(2,2) = 2. A(2,3) = 3. A(3,1) = 2. A(3,2) = -2. A(3,3) = -1. ! Call subroutine to create L and U matrices from A CALL lumake(L,U,A,n) ! Print results PRINT *, '-----------------------' DO i=1,n DO j=1,n PRINT *, i, j, A(i,j), L(i,j), U(i,j) END DO END DO PRINT *, '-----------------------' ! Call subroutine to solve for "x" using L and U CALL lusolve(x,L,U,b,n) ! Print results PRINT *, '-----------------------' DO i=1,n PRINT *, i, x(i) END DO PRINT *, '-----------------------' END PROGRAM MATSOLVE ! Create subroutine to make L and U matrices SUBROUTINE lumake(LL,UU,AA,n1) IMPLICIT NONE INTEGER,PARAMETER :: n=3 INTEGER :: i,j,k REAL :: LUSUM INTEGER,INTENT(IN) :: n1 REAL,DIMENSION(n,n),INTENT(IN) :: AA REAL,DIMENSION(n,n),INTENT(OUT) :: LL,UU ! We first note that the diagonal in our UPPER matrix is ! going to be UU(j,j) = 1.0, this allows us to initialize ! the first set of expressions UU(1,1) = 1. ! Find first column of LL DO i = 1,n1 LL(i,1) = AA(i,1)/UU(1,1) END DO ! Now find first row of UU DO j = 2,n1 UU(1,j) = AA(1,j)/LL(1,1) END DO ! Now find middle LL elements DO j = 2,n1 DO i = j,n1 LUSUM = 0. DO k = 1,j-1 LUSUM = LUSUM + LL(i,k)*UU(k,j) END DO LL(i,j) = AA(i,j) - LUSUM END DO ! Set Diagonal UU UU(j,j) = 1. ! Now find middle UU elements DO i = j+1,n1 LUSUM = 0. DO k = 1,j-1 LUSUM = LUSUM + LL(j,k)*UU(k,i) END DO UU(j,i) = (AA(j,i) - LUSUM)/LL(j,j) END DO END DO END SUBROUTINE lumake ! Make subroutine to solve for x SUBROUTINE lusolve(xx,L2,U2,bb,n2) IMPLICIT NONE INTEGER,PARAMETER :: n=3 INTEGER :: i,j,k REAL :: LYSUM,UXSUM REAL,DIMENSION(n):: y INTEGER,INTENT(IN) :: n2 REAL,DIMENSION(n),INTENT(IN) :: bb REAL,DIMENSION(n,n),INTENT(IN) :: L2,U2 REAL,DIMENSION(n),INTENT(OUT) :: xx ! Initialize DO i=1,n2 y(i) = 0. xx(i) = 0. END DO ! Solve L.y = b y(1) = bb(1)/L2(1,1) DO i = 2,n2 LYSUM = 0. DO k = 1,i-1 LYSUM = LYSUM + L2(i,k)*y(k) END DO y(i) = (bb(i) - LYSUM)/L2(i,i) END DO ! Now do back subsitution for U.x = y xx(n2) = y(n2)/U2(n2,n2) DO j = n2-1,1,-1 UXSUM = 0. DO k = j+1,n2 UXSUM = UXSUM + U2(j,k)*xx(k) END DO xx(j) = y(j) - UXSUM END DO END SUBROUTINE lusolve -- Tyler Joseph Hayes 600 Talbot St. -- Apt. 812 London, Ontario N6A 5L9 Tel : 519.435.0967 Fax : 519.661.3198 Cell : 416.655.7897 email: thayes(a)uwo.ca GPG Key ID# 0x3AE05130

3 9

[ANN] mlabrap-1.0b: a high level python to matlab
by Alexander Schmolck 28 Feb '07

28 Feb '07

URL --- <http://mlabwrap.sourceforge.net> Description ----------- Mlabwrap-1.0 is a high-level python to matlab(tm) bridge that makes calling matlab functions from python almost as convenient as using a normal python library. It is available under a very liberal license (BSD/MIT) and should work on all major platforms and (non-ancient) python and matlab versions and either numpy or Numeric (Numeric support will be dropped in the future). News ---- version 1.0b brings python 2.5 compatibility and various small fixes (improved error handling for 7.3, improvements to setup.py etc.). Provided I don't get any bug reports within the next two weeks or so the only difference between this version and 1.0 final will be cheeseshop support. Since mlabwrap 1.0 will be the last version that offers Numeric support anyone who wants to put off the switch to numpy a bit longer and is interested in using mlabwrap is strongly encouraged to download, test and possibly submit a bug report now. Examples -------- Creating a simple line plot: >>> from mlabwrap import mlab; mlab.plot([1,2,3],'-o') Creating a surface plot: >>> from mlabwrap import mlab; from numpy import * >>> xx = arange(-2*pi, 2*pi, 0.2) >>> mlab.surf(subtract.outer(sin(xx),cos(xx))) Creating a neural network and training it on the xor problem (requires netlab) >>> net = mlab.mlp(2,3,1,'logistic') >>> net = mlab.mlptrain(net, [[1,1], [0,0], [1,0], [0,1]], [0,0,1,1], 1000) Thanks go to Taylor Berg and many others who tested and provided feedback for the upcoming 1.0 release; more acknowledgements can be found on the website. cheers, 'as

2 1

Re: [Numpy-discussion] How to tell numpy to use gfortran as a, compiler
by Jon Wright 28 Feb '07

28 Feb '07

>Date: Sun, 25 Feb 2007 15:43:03 +0100 (CET) >From: "Sturla Molden" <sturla(a)molden.no> >Subject: Re: [Numpy-discussion] How to tell numpy to use gfortran as a > compiler ? >Content-Type: text/plain;charset=iso-8859-1 > > [...] >There is a third Fortran compiler based on the GCC backend called "g95". >It is not a part of GCC and uses copyrighted source code illegally >(numerous GPL violations). The head developer is rumored to have serious >bad karma. Apart from that, g95 is an ok Fortran 95 compiler. Always >ensure that ambiguous switches like "gnu95" means gfortran and not g95. I >don't know what NumPy does. > > These comments are out of order - even for the internet. I remember gfortran was a fork of g95 which occured when the developers wanted to go in different directions. It does not violate the gpl, source is available for download. Jon

1 0

Draft PEP for the new buffer interface to be in Python 3000
by Travis Oliphant 27 Feb '07

27 Feb '07

PEP: <unassigned> Title: Revising the buffer protocol Version: $Revision: $ Last-Modified: $Date: $ Author: Travis Oliphant <oliphant(a)ee.byu.edu> Status: Draft Type: Standards Track Created: 28-Aug-2006 Python-Version: 3000 Abstract This PEP proposes re-designing the buffer API (PyBufferProcs function pointers) to improve the way Python allows memory sharing in Python 3.0 In particular, it is proposed that the multiple-segment and character buffer portions of the buffer API are eliminated and additional function pointers are provided to allow sharing any multi-dimensional nature of the memory and what data-format the memory contains. Rationale The buffer protocol allows different Python types to exchange a pointer to a sequence of internal buffers. This functionality is '''extremely''' useful for sharing large segments of memory between different high-level objects, but it's too limited and has issues. 1. There is the little (never?) used "sequence-of-segments" option (bf_getsegcount) 2. There is the apparently redundant character-buffer option (bf_getcharbuffer) 3. There is no way for a consumer to tell the buffer-API-exporting object it is "finished" with its view of the memory and therefore no way for the exporting object to be sure that it is safe to reallocate the pointer to the memory that it owns (the array object reallocating its memory after sharing it with the buffer object which held the original pointer led to the infamous buffer-object problem). 4. Memory is just a pointer with a length. There is no way to describe what's "in" the memory (float, int, C-structure, etc.) 5. There is no shape information provided for the memory. But, several array-like Python types could make use of a standard way to describe the shape-interpretation of the memory (!wxPython, GTK, pyQT, CVXOPT, !PyVox, Audio and Video Libraries, ctypes, !NumPy, data-base interfaces, etc.) There are two widely used libraries that use the concept of discontiguous memory: PIL and NumPy. Their view of discontiguous arrays is a bit different, though. NumPy uses the notion of constant striding in each dimension as it's basic concept of an array. In this way a simple sub-region of a larger array can be described without copying the data. Strided memory is a common way to describe data to many computing libraries (such as the BLAS and LAPACK). The PIL uses a more opaque memory representation. Sometimes an image is contained in a contiguous segment of memory, but sometimes it is contained in an array of pointers to the contiguous segments (usually lines) of the image. This allows the image to not be loaded entirely into memory. The PIL is where the idea of multiple buffer segments in the original buffer interface came from, I believe. The buffer interface should allow discontiguous memory areas to share standard striding information. However, consumers that do not want to deal with strided memory should also be able to request a contiguous segment easily. Proposal Overview * Eliminate the char-buffer and multiple-segment sections of the buffer-protocol. * Unify the read/write versions of getting the buffer. * Add a new function to the protocol that should be called when the consumer object is "done" with the view. * Add a new function to allow the protocol to describe what is in memory (unifying what is currently done now in struct and array) * Add a new function to allow the protocol to share shape information * Fix all objects in core and standard library to conform to the new interface * Extend the struct module to handle more format specifiers Specification Change the PyBufferProcs structure to typedef struct { getbufferproc bf_getbuffer releasebufferproc bf_releasebuffer formatbufferproc bf_getbufferformat shapebufferproc bf_getbuffershape } typedef PyObject *(*getbufferproc)(PyObject *obj, void **buf, Py_ssize_t *len, int requires) Return a pointer to memory in buf and the length of that memory buffer in buf. Requirements for the memory are provided in requires (PYBUFFER_WRITE, PYBUFFER_ONESEGMENT). NULL is returned and an error raised if the object cannot return a view with those requirements. Otherwise, an object-specific "view" object is returned (which can just be a borrowed reference to obj). This view object should be used in the other API calls and does not need to be decref'd. It should be "released" if the interface exporter provides the bf_releasebuffer function. typedef int (*releasebufferproc)(PyObject *view) This function is called when a view of memory previously acquired from the object is no longer needed. It is up to the exporter of the API to make sure all views have been released before eliminating a reference to a previously returned pointer. It is up to consumers of the API to call this function on the object whose view is obtained when it is no longer needed. A -1 is returned on error and 0 on success. typedef char *(*formatbufferproc)(PyObject *view, int *itemsize) Get the format-string of the memory using the struct-module string syntax (see below for proposed additions to that syntax). Also, there is never an alignment assumption in this string---the full byte-layout is always required. If the implied size of this string is smaller than the length of the buffer then it is assumed that the string is repeated. If itemsize is not NULL, then return the size implied by the format string. This could be the entire length of the buffer or just the length of each element. It is equivalent to *itemsize = PyObject_SizeFromFormat(ret) if ret is the returned string. However, very often objects already know the itemsize without having to compute it separately. typedef PyObject *(*shapebufferproc)(PyObject *view) Return a 2-tuple of lists containing shape information: (shape, strides). The strides object can be None if the memory is C-style contiguous) otherwise it provides the striding in each dimension. All of these routines are optional for a type object (but the last three make no sense unless the first one is implemented). New C-API calls are proposed int PyObject_CheckBuffer(PyObject *obj) return 1 if the getbuffer function is available otherwise 0 PyObject * PyObject_GetBuffer(PyObject *obj, void **buf, Py_ssize_t *len, int requires) return a borrowed reference to a "view" object of memory for the object. Requirements for the memory should be given in requires (PYBUFFER_WRITE, PYBUFFER_ONESEGMENT). The memory pointer is in *buf and its length in *len. Note, the memory is not considered a single segment of memory unless PYBUFFER_ONESEGMENT is used in requires. Get possible striding using PyObject_GetBufferShape on the view object. int PyObject_ReleaseBuffer(PyObject *view) call this function to tell obj that you are done with your "view" This is a no-op if the object doesn't implement a release function. Only call this after a previous PyObject_GetBuffer has succeeded. Return -1 on error. char * PyObject_GetBufferFormat(PyObject *view, int *itemsize) Return a NULL-terminated string indicating the data-format of the memory buffer. The string is in struct-module syntax with the exception that there is never an alignment assumption (all bytes must be accounted for). If the length of the buffer indicated by this string is smaller than the total length of the buffer, then a repeat of the string is implied to fill the length of the buffer. If itemsize is not NULL, then return the implied size of each item (this could be calculated from the format string but it is often known by the view object anyway). PyObject * PyObject_GetBufferShape(PyObject *view) Return a 2-tuple of lists (shape, stride) providing the multi-dimensional shape of the memory area. The stride shows how many bytes to skip in each dimension to move in that dimension from the start of the array. Memory that is not a single contiguous-buffer can be represented with the pointer returned from GetBuffer and the shape and strides returned from GetBufferShape. int PyObject_SizeFromFormat(char *) Return the implied size of the data-format area from a struct-style description. Additions to the struct string-syntax The struct string-syntax is missing some characters to fully implement data-format descriptions already available elsewhere (in ctypes and NumPy for example). Here are the proposed additions: Character Description ================================== '1' bit (number before states how many bits) '?' platform _Bool type 'g' long double 'F' complex float 'D' complex double 'G' complex long double 'c' ucs-1 (latin-1) encoding 'u' ucs-2 'w' ucs-4 'O' pointer to Python Object 'T{}' structure (detailed layout inside {}) '(k1,k2,...,kn)' multi-dimensional array of whatever follows ':name:' optional name of the preceeding element '&' specific pointer (prefix before another charater) 'X{}' pointer to a function (optional function signature inside {}) The struct module will be changed to understand these as well and return appropriate Python objects on unpacking. Un-packing a long-double will return a c-types long_double. Unpacking 'u' or 'w' will return Python unicode. Unpacking a multi-dimensional array will return a list of lists. Un-packing a pointer will return a ctypes pointer object. Un-packing a bit will return a Python Bool. Endian-specification ('=','>','<') is also allowed inside the string so that it can change if needed. The previously-specified endian string is enforce at all times. The default endian is '='. According to the struct-module, a number can preceed a character code to specify how many of that type there are. The (k1,k2,...,kn) extension also allows specifying if the data is supposed to be viewed as a (C-style contiguous, last-dimension varies the fastest) multi-dimensional array of a particular format. Functions should be added to ctypes to create a ctypes object from a struct description, and add long-double, and ucs-2 to ctypes. Code to be affected All objects and modules in Python that export or consume the old buffer interface will be modified. Here is a partial list. * buffer object * bytes object * string object * array module * struct module * mmap module * ctypes module anything else using the buffer API Issues and Details The proposed locking mechanism relies entirely on the objects implementing the buffer interface to do their own thing. Ideally an object that implements the buffer interface should keep at least a number indicating how many releases are extant. The handling of discontiguous memory is new and can be seen as a modification of the multiple-segment interface. It is motivated by NumPy (used to be Numeric). NumPy objects should be able to share their strided memory with code that understands how to manage strided memory. Code should also be able to request contiguous memory if needed and objects exporting the buffer interface should be able to handle that either by raising an error (or constructing a read-only contiguous object and returning that as the view). Currently the struct module does not allow specification of nested structures. It seems like specifying a nested structure should be specified as several ways of viewing memory areas (ctypes and NumPy) already allow this. Copyright This PEP is placed in the public domain

3 7

Buffer Interface for Python 3.0
by Travis Oliphant 27 Feb '07

27 Feb '07

An update for those of you who did not get the chance to come to PyCon. PyCon was very well attended this year and there were some excellent discussions and presentations. From PyCon I learned that Python 3000 is closer than I had previously thought. What this means for me, is that I am now focussing toward getting a re-vamped buffer interface into Python 3.0 This will help us interact with the Python developers more effectively. Once we have the buffer interface hammered out for Python 3.0 we can back-port the result to Python 2.6 Thus my buffer PEP is being re-vamped. Anybody who would like to comment or contribute to the design of the new buffer interface is welcome to voice their opinion. Basically, what we are going to do now is 1) Return the data-format specification in an extended struct-style string 2) Return the shape information in a tuple of lists: (shape, strides) There are two questions I'm grappling with right now: 1) Do we propose the inclusion of offsets in the shape information? NumPy does not use offsets internally but simply has a pointer to the start of the array. 2) The buffer interface needs to understand the idea of discontiguous arrays. If the shape/stride information is separate from the pointer-to-data call, then the user needs to know if that pointer-to-data is a "contiguous chunk" or just the beginning of a strided memory area (and so should not be treated as a single-segment). 3) If we support strided memory areas, then we should probably also allow some way for PIL-like objects to report their buffer sequence (I'm sure this was the origin of the multi-segment buffer protocol to begin with). Or we could just ignore that possibility. The PIL would have to copy memory in order to share it's images. Anybody with ideas is welcome to participate. What I have so far is at http://wiki.python.org/moin/ArrayInterface Thanks, -Travis

2 1

building an array using smaller ones
by Rudolf Sykora 27 Feb '07

27 Feb '07

Hello everybody, I wonder how I could most easily accomplish the following: Say I have sth like: a = array( [1, 2] ) and I want to use this array to build another array in the following sence: b = array( [[1, 2, 3, a], [5, a, 6, 7], [0, 2-a, 3, 4]]) # this doesn't work I would like to obtain b = array( [[1, 2, 3, 1, 2], [5 ,1 ,2 ,6 ,7], [0, 1, 0, 3, 4]] ) I know a rather complicated way but believe there must be an easy one. Thank you very much. Ruda

1 0

inconsistent mgrid results
by Andrew Corrigan 27 Feb '07

27 Feb '07

I'm confused about the following: >>> print mgrid[2.45:2.6:0.05, 0:5:1] [[[ 2.45 2.45 2.45 2.45 2.45] [ 2.5 2.5 2.5 2.5 2.5 ]] [[ 0. 1. 2. 3. 4. ] [ 0. 1. 2. 3. 4. ]]] >>> print mgrid[2.45:2.6:0.05] [ 2.45 2.5 2.55] In the first case in the first dimension I get 2.45, 2.5. In the second case in the first dimension I get 2.45, 2.5, 2.55 In both cases I'm using 2.45:2.6:0.05 to specify the grid in the first dimension. Why is there a difference? Am I using mgrid incorrectly? I'm using the latest numpy release from sourceforge. Thanks, Andrew

1 0

nan, warning and error modes
by David Cournapeau 26 Feb '07

26 Feb '07

Hi, I am developing some numpy code, which sometimes fail because of nan. This is likely to be due to some bad coding on my side, and as such any NaN is a bug for this particular piece of code. Is there a way to get a warning when the first Nan is detected in the code (or even a faulty assertion) ? It looks like there are some variables/functions related to that in numpy, but I didn't find any useful document on the matter, either in the numpy book nor on the scipy website. cheers, David

2 3

How to tell numpy to use gfortran as a compiler ?
by David Cournapeau 26 Feb '07

26 Feb '07

Hi, I try to compile numpy using gfortran, using: python setup.py config --fcompiler=gnu But this does not work. Whatever option I try, numpy build system uses g77, and as a result, I have problems with my ATLAS library compiled with gfortran. What should I do to compiler numpy with blas/lapack compiler with gfortran ? cheers, David

3 6

Advice on masked array implementation
by Fred Clare 26 Feb '07

26 Feb '07

We would like some advice on how to proceed with implementing masked array capabilities in a large package of climate-related analysis functions. We are in the initial stages of trying to duplicate functionality from an existing package written in a locally-developed scripting language. The existing functionality depends heavily on masked arrays (actually on variables with attributes, with "fill_value" being one such). We have polled our user base and, while all responders plan to convert to NumPy, many have not started the conversion or are in transition. It is our experience that converting users to new ways is usually a multi-year undertaking, so it seems that we may need to support both Numeric and NumPy installations for some time to come. Even if people have converted to the NumPy version of our package, they may still be importing packages that have not been converted. Our initial design attempts to deal with the possibility of users potentially mixing (or using exclusively) Numeric arrays, Numeric masked arrays, NumPy arrays, and NumPy masked arrays. For example, suppose you have a function: result = func(arg0, arg1) where the two arguments and return variable can be any one of the four types of arrays mentioned. Currently we are testing to see if either argument is a NumPy or Numeric masked array. If just Numeric masked arrays, then we return a Numeric masked array. If just NumPy masked arrays, then we return a NumPy masked array. If one is a Numeric masked array and the other is a NumPy masked array, then we return a NumPy masked array. Similar checking is done for using just Numeric and/or NumPy non-masked arrays. Does this seem like a reasonable approach? It is tempting just to go with NumPy, but then we will have a large class of users who cannot access the new functionality. We have followed the discussion on the development of the new maskedarray module, but have not used it. I went to http://projects.scipy.org/scipy/numpy/attachment/wiki/MaskedArray/ maskedarray.py as referenced in a posting from Pierre GM, but I got "Internal Error." Has there been any decision as to whether maskedarray will be in NumPy version 1.1? Any estimate as to when 1.1 would be out? If we commit to numpy.core.ma now, how much trouble will it be to convert to the new maskedarray? Is there any user documentation on maskedarray and details on the differences between it and numpy.core.ma? Thanks, Fred Clare

2 1