Mailman 3 April 2005 - NumPy-Discussion

% and fmod
by Arnd Baecker April 14, 2005

April 14, 2005

Dear all, I encountered the following puzzling behaviour of the modulo operator %: In [1]: import Numeric In [2]: print Numeric.__version__ 23.8 In [3]: x=Numeric.arange(10.0) In [4]: print x%4 [ 0. 1. 2. 3. 0. 1. 2. 3. 0. 1.] In [5]: print 3.0%4 3.0 In [6]: print (-x)%4 [-0. -1. -2. -3. -0. -1. -2. -3. -0. -1.] # <====== In [7]: print (-3.0)%4 # vs. 1.0 # <====== (OK) In [8]: print Numeric.fmod(x,4) [ 0. 1. 2. 3. 0. 1. 2. 3. 0. 1.] In [9]: print Numeric.fmod(-x,4) [-0. -1. -2. -3. -0. -1. -2. -3. -0. -1.] So it seems that for arrays % behaves like fmod! This seems in contrast to what one finds in the python 2.3 documentation: "5.6. Binary arithmetic operations" """The % (modulo) operator yields the remainder from the division of the first argument by the second. [...] The arguments may be floating point numbers, e.g., 3.14%0.7 equals 0.34 (since 3.14 equals 4*0.7 + 0.34.) The modulo operator always yields a result with the same sign as its second operand (or zero); the absolute value of the result is strictly smaller than the absolute value of the second operand.""" I am presently teaching a course on computational physics with python and the students have huge difficulties with % behaving differently for arrays and scalars. I am aware that (according to Kernighan/Ritchie) the C standard does not define the result of % when any of the operands is negative. So can someone help me: is the different behaviour of % for scalars and arrays a bug, a feature, or what should I tell my students ? ;-). Many thanks, Arnd P.S.: BTW: the documentation for fmod and remainder is pretty short on this: In [3]:fmod? Type: ufunc String Form: <ufunc 'fmod'> Namespace: Interactive Docstring: fmod(x,y) is remainder(x,y) In [4]:remainder? Type: ufunc String Form: <ufunc 'remainder'> Namespace: Interactive Docstring: returns remainder of division elementwise Are contributions of more detailed doc-strings welcome ? P.P.S.: for numarray one gets even less information: In [1]: import numarray In [2]: numarray.fmod? Type: _BinaryUFunc Base Class: <class 'numarray.ufunc._BinaryUFunc'> String Form: <UFunc: 'remainder'> Namespace: Interactive Docstring: Class for ufuncs with 2 input and 1 output arguments In [3]: numarray.remainder? Type: _BinaryUFunc Base Class: <class 'numarray.ufunc._BinaryUFunc'> String Form: <UFunc: 'remainder'> Namespace: Interactive Docstring: Class for ufuncs with 2 input and 1 output arguments In [4]: print numarray.__version__ 1.1.1 P^3.S: scipy's mod seems to be an alternative: In [1]: import scipy In [2]: scipy.mod? Type: function Base Class: <type 'function'> String Form: <function mod at 0x40383994> Namespace: Interactive File: /usr/lib/python2.3/site-packages/scipy_base/function_base.py Definition: scipy.mod(x, y) Docstring: x - y*floor(x/y) For numeric arrays, x % y has the same sign as x while mod(x,y) has the same sign as y. In [3]: x=-scipy.arange(10) In [4]: x%4 Out[4]: array([ 0, -1, -2, -3, 0, -1, -2, -3, 0, -1]) In [5]: scipy.mod(x,4) Out[5]: array([ 0., 3., 2., 1., 0., 3., 2., 1., 0., 3.]) In [6]: scipy.mod?? Type: function Base Class: <type 'function'> String Form: <function mod at 0x40383994> Namespace: Interactive File: /usr/lib/python2.3/site-packages/scipy_base/function_base.py Definition: scipy.mod(x, y) Source: def mod(x,y): """ x - y*floor(x/y) For numeric arrays, x % y has the same sign as x while mod(x,y) has the same sign as y. """ return x - y*Numeric.floor(x*1.0/y)

1 0

ANN: matplotlib-0.80
by John Hunter April 14, 2005

April 14, 2005

A lot of development has gone into matplotlib since the last major release, which I'll summarize here. For details, see the notes for the incremental releases at http://matplotlib.sf.net/whats_new.html. Improvements since 0.70 -- contouring: Lots of new contour funcitonality with line and polygon contours provided by contour and contourf. Automatic inline contour labeling with clabel. See http://matplotlib.sourceforge.net/screenshots.html#pcolor_demo -- QT backend Sigve Tjoraand, Ted Drain and colleagues at the JPL collaborated on a QTAgg backend -- Unicode strings are rendered in the agg and postscript backends. Currently, all the symbols in the unicode string have to be in the active font file. In later releases we'll try and support symbols from multiple ttf files in one string. See examples/unicode_demo.py -- map and projections A new release of the basemap toolkit - See http://matplotlib.sourceforge.net/screenshots.html#plotmap -- Auto-legends The automatic placement of legends is now supported with loc='best'; see examples/legend_auto.py. We did this at the matplotlib sprint at pycon -- Thanks John Gill and Phil! Note that your legend will move if you interact with your data and you force data under the legend line. If this is not what you want, use a designated location code. -- Quiver (direction fields) Ludovic Aubry contributed a patch for the matlab compatible quiver method. This makes a direction field with arrows. See examples/quiver_demo.py -- Performance optimizations Substantial optimizations in line marker drawing in agg -- Robust log plots Lots of work making log plots "just work". You can toggle log y Axes with the 'l' command -- nonpositive data are simply ignored and no longer raise exceptions. log plots should be a lot faster and more robust -- Many more plotting functions, bugfixes, and features, detailed in the 0.71, 0.72, 0.73 and 0.74 point release notes at http://matplotlib.sourceforge.net/whats_new.html http://matplotlib.sourceforge.net JDH

1 0

ANN: numarray-1.3.0
by Todd Miller April 14, 2005

April 14, 2005

Release Notes for numarray-1.3.0 Numarray is an array processing package designed to efficiently manipulate large multi-dimensional arrays. Numarray is modelled after Numeric and features c-code generated from python template scripts, the capacity to operate directly on arrays in files, arrays of heterogeneous records, string arrays, and in-place operation on memory mapped files. I. ENHANCEMENTS 1. Migration of NumArray.__del__ to C (tp_dealloc). Overall performance. 2. Removal of dictionary update from array view creation improves performance of view/slice/subarray creation. This should e.g. improve the performance of wxPython sequence protocol access to Nx2 arrays. Subclasses now need to do a.flags |= numarray.generic._UPDATEDICT to ensure that dictionary based attributes are inherited by views. NumArrays no longer do this by default. 2. Modifications to support scipy.special. 3. Removal of an unnecessary getattr() from ufunc calling sequence. Ufunc performance. II. BUGS FIXED / CLOSED 1179355 average() broken in numarray 1.2.3 1167184 Floating point exception in numarray's dot() 1151892 Bug in matrixmultiply with zero size arrays 1160184 RecArray reversal 1156172 Incorect error message for shape incompatability 1155538 Incorrect error message when multiplying arrays See http://sourceforge.net/tracker/?atid=450446&group_id=1369&func=browse for more details. III. CAUTIONS This release should be backward binary compatible with numarray 1.1.1 and 1.2.3. WHERE ----------- Numarray-1.3.0 windows executable installers, source code, and manual is here: http://sourceforge.net/project/showfiles.php?group_id=1369 Numarray is hosted by Source Forge in the same project which hosts Numeric: http://sourceforge.net/projects/numpy/ The web page for Numarray information is at: http://stsdas.stsci.edu/numarray/index.html Trackers for Numarray Bugs, Feature Requests, Support, and Patches are at the Source Forge project for NumPy at: http://sourceforge.net/tracker/?group_id=1369 REQUIREMENTS ------------------------------ numarray-1.3.0 requires Python 2.2.2 or greater. Python-2.3.4 or Python-2.4.1 is recommended. AUTHORS, LICENSE ------------------------------ Numarray was written by Perry Greenfield, Rick White, Todd Miller, JC Hsu, Paul Barrett, Phil Hodge at the Space Telescope Science Institute. We'd like to acknowledge the assitance of Francesc Alted, Paul Dubois, Sebastian Haase, Chuck Harris, Tim Hochberg, Nadav Horesh, Edward C. Jones, Eric Jones, Jochen Kuepper, Travis Oliphant, Pearu Peterson, Peter Verveer, Colin Williams, Rory Yorke, and everyone else who has contributed with comments and feedback. Numarray is made available under a BSD-style License. See LICENSE.txt in the source distribution for details. -- Todd Miller jmiller(a)stsci.edu

1 0

numarray, nd_image transforms, and multi-channel images
by Alexandre Guimond April 14, 2005

April 14, 2005

Hi all. I've been looking at numarray to do some image processing. A lot of the work I do deal with transforming images, either with affine transformations, or vector field. Numarray seems somewhat well equiped to address these issues, but I am concerned about one aspect. It seems that the transformation code (affine_transforrm and geometric_transform) computes input coordonates for every output coordinate in the resulting array. If I have an RGB image for which the transformation is the same for all 3 RGB channels, I would assume that this will triple the workload unncessarily. It might have a dramatic effect for the geometric transformation which will most often be slower then affine. Is there any way around this, e.g. is it possible to specify numarray to use the same interpolation coefficients for the last "n" dimention of the array, or to tell numarray to only compute interpolation coefficients and apply those seperatly for each channel? thx for any help / info. alex.

2 1

scipy.base
by Pearu Peterson April 13, 2005

April 13, 2005

Hi Travis, I have committed scipy.{distutils,base} to Numeric3 CVS repository. scipy.distutils is a reviewed version of scipy_distutils and as one of its new features there is Configuration class that allows one to write much simpler setup.py files for subpackages. See setup.py files under Numeric3/scipy directory for examples. scipy.base is a very minimal copy of scipy_base plus ndarray modules. When using setup_scipy.py for building, the ndarray package is installed as scipy.base and from scipy.base import * should work equivalently to from ndarray import * for instance. I have used information from Numeric3/setup.py to implement Numeric3/scipy/base/setup.py and it should be updated whenever Numeric3/setup.py is changed. However, I would recommend start using scipy.base instead of ndarray as using both may cause unexpected behaviour when installed ndarray is older than scipy.base installation (see [*]). In Numeric3 CVS repository that would mean replacing setup.py with setup_scipy.py and any modification to ndarray setup scripts should be done in scipy/base/setup.py. We can apply this step whenever you feel confident with new setup.py files. Let me know if you have any troubles with them. To clean up Numeric3 CVS repository completely then Include, Src, Lib, CodeGenerators directories should be moved under the scipy/base directory. However, this step can be omitted if you would prefer working with files at the top directory of Numeric3. Current setup.py scripts fully support this approach as well. There are also few open issues and questions. First, how to name Numeric3 project when it installs scipy.base, scipy.distutils, Numeric packages, etc? This name will be used when creating source distributions and also as part of the path where header files will be installed. At the moment setup_scipy.py uses the name 'ndarray'. And so `setup_scipy.py sdist`, for example, produces ndarray-30.0.tar.gz file; `setup_scipy.py install` installs header files under <prefix>/include/ndarray/ directory. Though this is fine with me, I am not sure that this is an ideal situation. I think we should choose the name now and stick to it forever, especially since 3rd party extension modules need to know where to look for ndarray header files. This name cannot be 'numarray', obviously, but there are options like 'ndarray', 'numpy', and may be others. In fact, 'Numeric' (with version 3x.x) would be also an option but that would be certainly cause some problems when one wants both Numeric 2x.x and Numeric 3x.x to be installed in the system, the header files would end up in the same directory, for instance. As a workaround, we could force installing Numeric3 header files to <prefix>/include/Numeric/3/ or something. I acctually like this idea but I wonder what other think about this. Second, is it already possible to use ndarray C/API as a replacement of Numeric C/API, i.e. would simple replacement of #include "Numeric/arrayobject.h" with #include "ndarray/arrayobject.h" work? And if not, will it ever be? This would be interesting to know as an extension writer. [*] Due to keeping changes to Numeric3 sources minimal, scipy.base multiarray and umath modules first try to import ndarray and then scipy.base whenever ndarray is missing. One should remove ndarray installation from the system before using scipy.base. Regards, Pearu

3 2

Subclassing and metadata
by Travis Oliphant April 12, 2005

April 12, 2005

I think I've found a possible solution for subclasses that want to handle metadata. Essentially, any subclass that defines the method _update_meta(self, other) will get that method called when an array is sliced, or subscripted. Anytime an array is created where a subtype is the caller, this method will be called if it is available. Here is a simple example: import ndarray class subclass(ndarray.ndarray): def __new__(self, shape, *args, **kwds): self = ndarray.ndarray.__new__(subclass, shape, 'V4') return self def __init__(self, shape, *args, **kwds): self.dict = kwds return def _update_meta(self, obj): self.dict = obj.dict Comments? -Travis

1 0

Silent overflow of Int32 array
by Tim Churches April 11, 2005

April 11, 2005

I just got caught by code equivalent to this (with NumPy 23.8 on 32 bit Linux): >>> import Numeric as N >>> a = N.array((2000000000,1000000000),typecode=N.Int32) >>> N.add.reduce(a) -1294967296 OK, it is an elementary mistake, but the silent overflow caught me unawares. casting the array to Float64 before summing it avoids the error, but in my instance the actual data is a rank-1 array of 21 million integers with a mean value of about 140 (which adds up more than sys.maxint), and casting to Float64 will use quite a lot of memory (as well as taking some time). Any advice for catching or avoiding such overflow without always incurring a performance and memory hit by always casting to Float64? Shouldn't add.reduce() be checking for overflow and raising an error? Then it would be easy to upcast only when overflow (or underflow) occurs, rather than always. Tim C

2 2

Re: [Numpy-discussion] Questions about the array interface.
by Scott Gilbert April 11, 2005

April 11, 2005

--- Scott Gilbert <xscottg(a)yahoo.com> wrote: > > --- Andrew Straw <strawman(a)astraw.com> wrote: > > > > I feared it was something like that. (No platform independent way to > > represent special values like nan, inf, and so on.) So I think if > > we're going to require an encoding character such as '<' or '>' we > > should also include one that means native which CAN handle these > > special values... And document why it's needed and why it may get one > > into trouble. > > > > Let's say we used ']' for big endian native, and '[' for little endian > native? Does that just indicate the possible presence of NaNs for Infs > in the data? > > Adding those codes doesn't have any affect on whether or not libraries > can deal with them. I guess I'm not understanding something. > I think I'm understanding my problem in understanding :-). There IS a platform independant way to represent NaNs and Infs. It's pretty clearly spelled out in IEEE-754: http://stevehollasch.com/cgindex/coding/ieeefloat.html I think something we've been assuming is that the array data is basically IEEE-754 compliant (maybe it needs to be byteswapped). If that's not true, then we're going to need some new typecodes. We're not supporting the ability to pass VAX floating point around (Are we????). The problem is that you can't make any safe assumptions about whether your current platform will deal with IEEE-754 data in any predictable way if it contains NaNs or Infs. So additional typecodes won't really solve anything. Tim Peter's explanation is a good representation of Python's official position regarding floating point issues, but a much simpler explanation is possible... The struct module in "standard mode" decodes the data one character at a time and builds a float from them. You can see this in the _PyFloat_Unpack8 function in the floatobject.c file. In other words, this routine probably works on a VAX too (taking a IEEE-754 double and building a VAX floating point as it goes). You can also see the comment in there that says it doesn't handle NaNs or Infs. I don't think we need another indicator for '>' big-endian or '<' for little-endian. Cheers, -Scott

5 5

Re: [Numpy-discussion] Questions about the array interface.
by Scott Gilbert April 9, 2005

April 9, 2005

--- konrad.hinsen(a)laposte.net wrote: > > This discussion has been coming up regularly for a few years. Until now > the concensus has always been that Python should make no assumptions > that go beyond what a C compiler can promise. Which means no > assumptions about floating-point representation. > > Of course the computing world is changing, and IEEE format may well be > ubiquitous by now. Vaxes must be in the museum by now. But how about > mainframes? IBM mainframes didn't use IEEE when I used them (last time > 15 years ago), and they are still around, possibly compatible to their > ancestors. > I've been following this mailing list for a few years now, but I skip a lot of threads. I almost certainly skipped this topic in the past since it wasn't relevant to me. I'm only interested in it now since it's relevant to this data interchange business, so I'm sorry if this is a rehash... Trying to stay portable is a good goal, and I can understand why Python proper would try to adhere to the restrictions it does. Despite the claim, Python makes plenty of assumptions that a standards conformant C compiler could break. If numpy doesn't make some assumptions about floating point representation, it's going to kill the possibity of passing data across machines, and that's pretty unacceptable. I'm not comfortable saying "ubiquitous" since I don't know what the mainframe or super computing community is making use of, and I don't know what sort of little machines Python is running on. The closest thing to a mainframe that I've ever used was a Convex, and I never knew what it's floating point representation was. However, I know that x86, PPC, AMD-64, IA64, Alpha, Sparc, and whatever HPUX and SGIs are running on all use IEEE-754 format. That's probably 99.999% of all machines capable of running Python, and at least that percentage of users. It would be a shame to gum up this typecode thing for situations that don't occur in practice. If it has to be done, then I recommend we use the '@' code in place of the '<' or '>' for platforms that are out of the ordinary. It's important to specify that '@' is only to be used on floating point data that is not IEEE-754. In this case it doesn't mean "native" like it does in the struct module, it means "weird" :-). > > Another detail to consider is that although most machines use the IEEE > representation, hardly any respects the IEEE rules for floating point > operations in all detail. In particular, trusting that Inf and NaN will > be treated as IEEE postulates is a risky business. > See that's the thing. Why burden how you label the data with the restrictions of the current machine? You can take the data off the machine. Whether or not I can rely on what NaN*Inf will give me, I know that I can take NaN and Inf to another machine and get the same interpretation of the data. This whole thread started because Andrew Straw showed that struct.pack('<d',nan) causes an exception, but that's just a limitation in the struct module. He was definitely running it on a machine that was capable of representing an 8 byte little-endian NaN. He doesn't need a new typecode until he tries to transport data from some esoteric mainframe. Cheers, -Scott

1 0

Alternate C-only array protocol for speed?
by David M. Cooke April 8, 2005

April 8, 2005

It seems that people are worried about speed of the attribute-based array interface when using small arrays in C. Here's an alternative: Define some attribute (for now, call it __array_c__), which returns a CObject whose value (which you get with PyCObject_GetVoidPtr) would be a pointer to a struct describing the array. It would look something like typedef struct { int version; int nd; Py_LONG_LONG *shape; char typecode; Py_LONG_LONG *strides; Py_LONG_LONG offset; void *data; } SimpleCArray; (The order here follows that of the array interface spec; if somebody's got any comments on what mixing int's, Py_LONG_LONG, and char's in a struct does to the packing and potential alignment problems I'd like to know.) version is there as a sanity check: I'd say for this version it's something like 0xDECAF ('cause it's lightweight, see ;-). It's primarily a check that you've got the right thing (sinc CObjects are intrinsically opaque types). Then: - the array object guarantees that the data, etc. remains alive, probably by passing itself as the desc parameter to the CObject. The array data would have to stay at the same location and the same size while the reference is held. - typecode follows that of the __array_typestr__ attribute - shape and strides are pointers to arrays of at least nd elements. - this doesn't handle byteswapped as-is. Maybe a flags, or endian, attribute could be added. - you can still have the full attribute-base array interface (__array_strides__, etc.) to fall back on. If the typecode is 'V', you'll have to look at __array_descr__. Creating one from a Numeric PyArrayObject would go like this: PyObject *create_SimpleCArray(PyArrayObject *a) { SimpleCArray *ca = PyMem_New(SimpleCArray, 1); ca->version = 0xDECAF; ca->nd = a->nd; ca->shape = PyMem_New(Py_LONG_LONG, ca->nd); for (i = 0; i < ca->nd; i++) { ca->shape[i] = a->dimensions[i]; } ca->strides = PyMem_New(Py_LONG_LONG, ca->nd); for (i = 0; i < ca->nd; i++) { ca->strides[i] = a->strides[i]; } ca->offset = 0; ca->data = &my_data; Py_INCREF(a); PyObject *co = PyCObject_FromVoidPtrAndDesc(ca, a, free_numeric_simplecarray); return co; } where void free_numeric_simplecarray(SimpleCArray *ca, PyArrayObject *a) { PyMem_Free(ca->shape); PyMem_Free(ca->strides); PyMem_Free(ca); Py_DECREF(a); } Some points: - you have to keep the CObject around: destroying it will potentially destroy the array you're looking at. - I was thinking that maybe adding a PyObject *owner could make it easier to keep track of the owner; I'm not sure, as the descr argument in CObjects can easily play that role. - The creator of the SimpleCArray is free to add elements to the end (as long as they don't affect the padding/alignment of the previous ones: haven't thought about this). You could put the real owner of the array data there, for example (say, if it was wrapping a Blitz++ array). Or have a small _strides[30] array at the end, and strides would point to that (saving you a memory allocation). This simple C interface would, I think, alleviate much worries about speed for small arrays, and even for large arrays. -- |>|\/|< /--------------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm(a)physics.mcmaster.ca

4 3