Mailman 3 June 2007 - NumPy-Discussion

Is this an indexing bug?
by Tom K. June 20, 2007

June 20, 2007

>>> h = zeros((1, 4, 100)) >>> h[0,:,arange(14)].shape (14, 4)

5 13

Re: [Numpy-discussion] Radix sort?
by Jon Wright June 20, 2007

June 20, 2007

> "Charles R Harris" <charlesr.harris(a)gmail.com> wrote: > > Straight radix sort might be an interesting option for some things. > However, its performance can depend on whether the input data is random or > not and it takes up more space than merge sort. Other potential drawbacks > arise from the bit twiddling needed for signed numbers and floats, the > former solved by converting to offset binary numbers (complement the sign > bit), and the latter in the way your links indicate, but both leading to a > proliferation of special cases. Maintaining byte order and byte addressing > portability between cpu architectures might also require masking and > shifting that will add computational expense and may lead to more special > cases for extended precision floats and so on. That said, I would be curious > to see how it works out if you want to give it a try. I agree completely about the proliferation of special cases, and mess that will make. This radix sort is bad for tiny arrays, but good for big random ones (there is no insertion sort either?). An intelligent algorithm chooser is needed, something like an "atlas"/"fftw" but for sorting, which has been invented already it seems. Platform and datatype and even the data themselves seem to be important. eg: http://www.spiral.net/software/sorting.html http://www.cgo.org/cgo2004/papers/09_cgo04.pdf Seems like a significant amount of work - and for the numpy case it should work with strided/sliced arrays without copying. Is there a list somewhere of the numpy numeric types, together with their binary representations on all of the numpy supported platforms? I'm guessing integers are almost always: [signed/unsigned] [8|16|32|64 bits] [big|little endian] ... and that many popular platforms only use ieee745 floats and doubles, which might be big or little endian. Is there an important case I miss, such as decimal hardware? The experimental observation is that the code from Michael Herf appears to win for float32 for random arrays >1024 elements or sorted arrays >2M elements, compared any of the 3 algorithms already in numpy (ymmv). The methods could also have a positive impact on the numpy.histogram function for the smaller data types, and also lead to other order statistics, like median, trimeans and n-th largest with a performance improvement. Since sorting is one of the most studied problems in all of computer science, I am reluctant to start writing a new library! Does anyone know of some good libraries we could try out? Best, Jon

2 1

Bugfix for numpy.info bug
by Fernando Perez June 20, 2007

June 20, 2007

Bug === In [8]: N.info(N.ones(3)) class: ndarray shape: (3,) strides: (8,) itemsize: 8 aligned: True contiguous: True fortran: True --------------------------------------------------------------------------- TypeError Traceback (most recent call last) /home/fperez/research/code/mwadap-stable/mwadap/test/<ipython console> in <module>() /home/fperez/usr/opt/lib/python2.5/site-packages/numpy/lib/utils.py in info(object, maxwidth, output, toplevel) 298 elif isinstance(object, ndarray): 299 import numpy.numarray as nn --> 300 nn.info(object, output=output, numpy=1) 301 elif isinstance(object, str): 302 if _namedict is None: /home/fperez/usr/opt/lib/python2.5/site-packages/numpy/numarray/functions.py in info(obj, output, numpy) 377 extra = "" 378 tic = "" --> 379 print >> output, "data pointer: %s%s" % (hex(obj.ctypes._as_parameter_), extra) 380 print >> output, "byteorder: ", 381 endian = obj.dtype.byteorder TypeError: hex() argument can't be converted to hex In [9]: debug > /home/fperez/usr/opt/lib/python2.5/site-packages/numpy/numarray/functions.py(379)info() 378 tic = "" --> 379 print >> output, "data pointer: %s%s" % (hex(obj.ctypes._as_parameter_), extra) 380 print >> output, "byteorder: ", Fix === planck[numpy]> svn diff Index: numpy/numarray/functions.py =================================================================== --- numpy/numarray/functions.py (revision 3874) +++ numpy/numarray/functions.py (working copy) @@ -376,7 +376,7 @@ else: extra = "" tic = "" - print >> output, "data pointer: %s%s" % (hex(obj.ctypes._as_parameter_), extra) + print >> output, "data pointer: %s%s" % (hex(obj.ctypes._as_parameter_.value), extra) print >> output, "byteorder: ", endian = obj.dtype.byteorder if endian in ['|','=']: Question ======== any objection if I commit this? Since I don't really touch the codebase often, I'd rather ask the real core people. I also don't know if it's really the right thing to do, I just tabbed into the object and picked what seemed to be the most reasonable answer. It's trivial, but I'd rather double check. Cheers, f

2 2

Radix sort?
by Jon Wright June 19, 2007

June 19, 2007

Dear numpy experts, I see from the docs that there seem to be 3 sorting algorithms for array data (quicksort, mergesort and heapsort). After hearing a rumour about radix sorts and floats I google'd and now I'm wondering about a radix sort for numpy (and Numeric) scalars? See: http://www.stereopsis.com/radix.html http://en.wikipedia.org/wiki/Radix_sort The algorithm is apparently a trick where no comparisons are used. A shockingly bad benchmark of a swig wrapped test function below suggests it really is quicker than the array.sort() numpy method for uint8. At least with >256 element uint8 test arrays (numpy1.0.3, mingw32, winXP), for me, today; of course ymmv. With large enough arrays of all of the integer numeric types and also ieee reals, appropriate versions of the radix sort might be able to: "... kick major ass in the speed department." [http://www.cs.berkeley.edu/~jrs/61b/lec/36] Have these sorts already been implemented in scipy? Can someone share some expertise in this area? There is a problem about the sorting not being done in-place (eg: better for argsort than sort). I see there is a lexsort, but it is not 100% obvious to me how to use that to sort scalars, especially ieee floats. If numpy is a library for handling big chunks of numbers with knowable binary representations, I guess there might be some general interest in having this radix sort compiled for the platforms where it works? Thanks for any opinions, Jon --- ==setup.py== from distutils.core import setup, Extension e = Extension("_sort_radix_uint8",sources=["radix_wrap.c"]) setup(ext_modules=[e]) ==radix.i== %module sort_radix_uint8 %{ #define SWIG_FILE_WITH_INIT void sort_radix_uint8(unsigned char * ar, int len){ int hits[256]; int i, j; for(i=0;i<256;i++) hits[i]=0; for(i=0;i<len;i++) hits[(int)ar[i]]+=1; i=0; j=0; while(i<256){ while(hits[i]>0){ /* shortcut for uint8 */ ar[j++] = (unsigned char) i; hits[i]--; } i++; } } %} %init %{ import_array(); %} %include "numpy.i" %apply (unsigned char* INPLACE_ARRAY1, int DIM1) { (unsigned char * ar, int len)} void sort_radix_uint8(unsigned char * ar, int len); == test.py == import numpy,sort_radix_uint8 from time import clock def clearcache(): t = numpy.random.random(1024*1024) t=t*t def test(i): t = numpy.random.random(i) a1 = (t*255).astype(numpy.uint8) a2 = (t*255).astype(numpy.uint8) clearcache() tick=clock() a1.sort() t1 = clock()-tick clearcache() tick = clock() sort_radix_uint8.sort_radix_uint8(a2) t2 = clock() - tick assert a1.all() == a2.all() clearcache() tick=clock() a1.sort() # already sorted... t3 = clock()-tick return t1,t2,t3 for j in range(25): r = test(pow(2,j)) print numpy.argmin(r), j,pow(2,j),r

2 1

question about numpy
by Xuemei Tang June 19, 2007

June 19, 2007

Dear Sir/Madam, I meet a problem when I installed numpy. I installed numpy by the command "python setup.py install". Then I tested it by "python -c 'import numpy; numpy.test()'". But it doesn't work. There is an error message: "Running from numpy source directory. Traceback (most recent call last): File "<string>", line 1, in ? AttributeError: 'module' object has no attribute 'test'" Could you help me to figure out what's wrong with it? Thanks! Best wishes! Xuemei Tang

5 5

Issues with the memmap object
by Sturla Molden June 18, 2007

June 18, 2007

After struggling with NumPy's memmap object, I examined the code and detected three severe problems. I suggest that memmap is removed from NumPy, at least on Windows, as it's shortcomings is severe and undocumented. Problem 1: I/O errors are never detected on Win32: On Windows, i/o errors are trapped using structured exception handling when using memory mapped objects. Neither NumPy nor Python use structured exception handling on Win32. This means that i/o errors (such as network or disk failure) will go undetected, and be a source of obscure bugs. The bugfix for this is to wrap any access attempt to an PyArrayObject's "data" pointer with __try and __except blocks, and using an MSVC compiler on Windows. GCC/MinGW cannot be used, as it does not support structured exception handling. In other words, PyArrayObject *memmap; __try { /* safe read/write access to memmap->data here */ } __except( GetExceptionCode() == EXCEPTION_IN_PAGE_ERROR ? EXCEPTION_EXECUTE_HANDLER : EXCEPTION_CONTINUE_SEARCH) { /* Windows signaled an I/O error, handle the problem here */ } Not only must NumPy itself be rewritten, but also any library getting a data pointer from a NumPy memmap array. Fixing this will be extremely difficult, if not impossible. The only safe way to access file data from NumPy is numpy.fromfile() and numpy.array.tofile(). Problem 2: Mapping always starts from the beginning of the file: Python's standard mmap object from the beginning of the file, regardless of the size. NumPy's memmap object depends on Python's mmap through the buffer protocol. Even though NumPy's memmap object takes an offset parameter, the actual memory mapping starts from the beginning of the file. Thus, virtual memory equal to the memmap object's offset parameter will be leaked until the memmap object is deleted. Problem 3: No 64 bit support on Windows or Linux: On Linux, large files must be memory mapped using mmap64 (or mmap2 if 4k boundaries are acceptable). On Windows, CreateFileMapping/MapViewOfFile has 64 bit support, but Python's mmap does not use it (the high offset DWORD is always zero). Only files smaller than 4 GB can be memory mapped. Regards, Sturla Molden

2 3

Anyone written an SQL-like interface to numpy/PyTables?
by Stephen Simmons June 18, 2007

June 18, 2007

Hi, Has anyone written a parser for SQL-like queries against PyTables HDF tables or numpy recarrays? I'm asking because I have written code for grouping then summing rows of source data, where the groups are defined by functions of the source data, or looking up a related field in a separate lookup tables. I use this for tracking the performance of customer segments by status, with 48 monthly files of product usage/customer status data on 5m customers. Each of these is a 5m row HDF file, with several other 5m row HDF files that are used to work out which segment a customer belongs to. This grouping and summing is equivalent to something like the following SQL code: SELECT grp_fn1(table1.*), grp_fn2(table1.*), grp_fn3(table2.*), grp_fn4(table3.*), count(table1.*), sum(table1.field1), sum(table1.field2), ..., sum(table1.fieldK) FROM table1 PARTITION(date) LEFT JOIN table2 ON table1.field0=table2.field0 LEFT JOIN table3 ON table2.field0=table3.field0 WHERE min_date<=date<=max_date GROUP BY grp_fn1(table1.*), grp_fn2(table1.*), grp_fn3(table2.*), grp_fn4(table3.*) I'm using numpy.bincount() function to do the grouping/summing, numpy.searchsorted() for fast lookup tables implementing the grouping functions grp_fn(), and some other C functions for a fast "zip" join of related tables whose primary keys are in the same order as the monthly date partitions. The Python code that specifies the grouping/summing fields looks like this: agg = HDFAggregator('table1.hdf') agg.add_group_function('MONTH', ...<min_date, max_date> ) agg.add_group_function('SEGMENT', ... <lookup on table 2> ) agg.add_group_function('STATUS', ... <lookup on table 3> ) agg.do_aggregation(groupby='MONTH SEGMENT STATUS', count='CUST_NO', sum='<list of fields to sum>) agg.add_calculated_field('PROFIT', 'VOLUME*(PRICE* (1-DISCOUNT)-COGS) - COSTOFSALES') agg.save('output.hdf') On my laptop, this zips over my data at a speed of 400k rows/sec, aggregating it into 230,000 groups (48 months x 120 customer segments/subsegments x 5 product groups x 8 statuses) with subtotals for 30 data fields in each group. This is essentially as fast as PyTables can read in the HDF files from disk; peak speeds with fewer groups (e.g. 48x5x1x4) are above 1Mrows/sec if the HDF files are already in the disk cache. One option I am considering now is bolting an SQL-like parser on the front to provide a more natural interface for those unfortunate people who prefer SQL to Python. I don't want to write an SQL parser from scratch, so it would be great to know if there are any existing projects to put an SQL-like interface on numpy or PyTables (other than numexpr). So has anyone looked at using an SQL-like syntax for querying numpy/PyTables data? Cheers Stephen

2 1

f2py and type construct
by lorenzo bolla June 18, 2007

June 18, 2007

hi all. I'm trying to compile an F90 source file with f2py, but it fails with the construct "type ... end type". here is an example: -------------------- ! file test19.f90 module basic implicit none save integer, parameter :: ciao = 17 end module basic module basic2 implicit none save type test_t integer :: x end type test_t type(test_t) :: ciao end module basic2 ---------------------- $>f2py -c test19.f90 -m test --fcompiler=intele --compiler=intel (I'm compiling on an SGI Altix) and this is the error message: ---------------------- running build running config_fc running build_src building extension "test" sources f2py options: [] f2py:> /tmp/tmprBrnf7/src.linux-ia64-2.5/testmodule.c creating /tmp/tmprBrnf7 creating /tmp/tmprBrnf7/src.linux-ia64-2.5 Reading fortran codes... Reading file 'test19.f90' (format:free) Post-processing... Block: test Block: basic Block: basic2 Block: test_t Post-processing (stage 2)... Block: test Block: unknown_interface Block: basic Block: basic2 Block: test_t Building modules... Building module "test"... Constructing F90 module support for "basic"... Variables: ciao Constructing F90 module support for "basic2"... Variables: ciao getctype: No C-type found in "{'typespec': 'type', 'typename': 'test_t'}", assuming void. Traceback (most recent call last): File "/xlv1/labsoi_devices/bollalo001/bin/f2py", line 26, in <module> main() File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/f2py/f2py2e.py", line 552, in main run_compile() File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/f2py/f2py2e.py", line 539, in run_compile setup(ext_modules = [ext]) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/distutils/core.py", line 174, in setup return old_setup(**new_attr) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/distutils/core.py", line 151, in setup dist.run_commands() File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/distutils/dist.py", line 974, in run_commands self.run_command(cmd) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/distutils/dist.py", line 994, in run_command cmd_obj.run() File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/distutils/command/build.py", line 112, in run self.run_command(cmd_name) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/distutils/cmd.py", line 333, in run_command self.distribution.run_command(command) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/distutils/dist.py", line 994, in run_command cmd_obj.run() File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/distutils/command/build_src.py", line 87, in run self.build_sources() File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/distutils/command/build_src.py", line 106, in build_sources self.build_extension_sources(ext) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/distutils/command/build_src.py", line 218, in build_extension_sources sources = self.f2py_sources(sources, ext) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/distutils/command/build_src.py", line 471, in f2py_sources ['-m',ext_name]+f_sources) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/f2py/f2py2e.py", line 362, in run_main ret=buildmodules(postlist) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/f2py/f2py2e.py", line 314, in buildmodules dict_append(ret[mnames[i]],rules.buildmodule(modules[i],um)) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/f2py/rules.py", line 1130, in buildmodule mr,wrap = f90mod_rules.buildhooks(m) File "/xlv1/labsoi_devices/bollalo001/lib/python2.5/site-packages/numpy/f2py/f90mod_rules.py", line 127, in buildhooks at = capi_maps.c2capi_map[ct] KeyError: 'void' Exit 1 ------------------------------------ module basic gives no problems, but module basic2 yes, because of the type construct. what am I doing wrong? thank you very much, lorenzo.

1 1

SharedLibrary builder (ticket 213)
by David Cournapeau June 18, 2007

June 18, 2007

Hi, following the discussion on using ctypes in scipy, the main problem seems the inability to build a usable dll usable by ctypes using distutils. Numpy ticket 213 tackles the issue; as I need it personnally, I am willing to work on the feature, but I would need a bit some advices on numpy.distutils: - First, is anyone working on the SharedLibrary right now ? - how to add a new argument to numpy.distutils.setup ? The way I see it, but I don't understand the whole distutils arch yet so I may be wrong, would be to add one argument to the setup function, like shared_libraries, which would contains a list of objects representing the libraries (sources, name, etc...), and which would be build by a special builder. But right now, I don't know how the whole calling sequence work such as my own builder is called for the items in a given argument of setup. cheers, David

1 0

Re: [Numpy-discussion] Problems with Numeric
by Ali Santacruz June 17, 2007

June 17, 2007

Thank you very much Chuck for your help. Regards, Ali S. ----- Mensaje original ---- De: Charles R Harris <charlesr.harris(a)gmail.com> Para: Discussion of Numerical Python <numpy-discussion(a)scipy.org> Enviado: domingo, 17 de junio, 2007 10:25:01 Asunto: Re: [Numpy-discussion] Problems with Numeric On 6/16/07, Ali Santacruz <amsd2013(a)yahoo.com> wrote: Hi dear list, I am designing an application that uses GDAL (which seems to bind to Numeric) to read an image and convert it to an array. The code runs without problem in the console, then I compile it (apparently without problems), but when I try to launch the application, it fails and the next error appears in the .log: Traceback (most recent call last): File "ViewGdal_0.0-1.py", line 10, in ? File "gdal\__init__.pyc", line 11, in ? File "gdal\gdalnumeric.pyc", line 85, in ? #from Numeric import * File "Numeric.pyc", line 93, in ? #from Precision import * File "Precision.pyc", line 26, in ? #_code_table = _fill_table(typecodes) File "Precision.pyc", line 23, in _fill_table #table[key] = _get_precisions(value) File "Precision.pyc", line 18, in _get_precisions TypeError: data type not understood This refers to the next in Precision.py in Numeric: typecodes = {'Character':'c', 'Integer':'1sil', 'UnsignedInteger':'bwu', 'Float':'fd', 'Complex':'FD'} def _get_precisions(typecodes): lst = [] for t in typecodes: lst.append( (zeros( (1,), t ).itemsize()*8, t) ) #this is line 18 return lst As far as I understand (I've been working with Python for just a couple of months), there is a problem with the typecodes. What can I do to solve this problem? I've not found a version of GDAL that binds to Numpy. I am using GDAL from hobutools (v. 1.75), Numeric v. 24.2, Python 2.4.3 on WinXP. Looks like GDAL supports numpy These facilities have evolved somewhat over time. In the past the package was known as "Numeric" and imported using "import Numeric". A new generation is imported using "import numpy". Currently the old generation bindings only support the older Numeric package, and the new generatio bindings only support the new generation numpy package. They are mostly compatible, and by importing gdalnumeric you will get whichever is appropriate to the current bindings type. I would look for a package with a more recent version of GDAL that supports numpy. Apart from that, I can't help you. Chuck _______________________________________________ Numpy-discussion mailing list Numpy-discussion(a)scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion __________________________________________________ Correo Yahoo! Espacio para todos tus mensajes, antivirus y antispam ¡gratis! Regístrate ya - http://correo.espanol.yahoo.com/

1 0