Mailman 3 October 2011 - NumPy-Discussion

speeding up operations on small vectors
by Christoph Groth 11 Oct '11

11 Oct '11

Dear numpy experts, I could not find a satisfying solution to the following problem, so I thought I would ask: In one part of a large program I have to deal a lot with small (2d or 3d) vectors and matrices, performing simple linear algebra operations with them (dot products and matrix multiplications). For several reasons I want the code to be in python. The most natural choice seemed to be to use numpy. However, the constant time cost when dealing with numpy arrays seems to be immense, as demonstrated by this toy program: **************************************************************** import numpy as np from time import time def points_tuple(radius): rr = radius**2 def inside(point): return point[0]**2 + point[1]**2 < rr M = ((1, 0), (0, 1)) for x in xrange(-radius, radius + 1): for y in xrange(-radius, radius + 1): point = (M[0][0] * x + M[0][1] * y, M[1][0] * x + M[1][1] * y) if inside(point): yield point def points_numpy(radius): rr = radius**2 def inside(point): return np.dot(point, point) < rr M = np.identity(2, dtype=int) for x in xrange(-radius, radius + 1): for y in xrange(-radius, radius + 1): point = np.dot(M, (x, y)) if inside(point): yield point def main(): r = 200 for func in [points_tuple, points_numpy]: t = time() for point in func(r): pass print func.__name__, time() - t, 'seconds' if __name__ == '__main__': main() **************************************************************** On my trusty old machine the output (python 2.6, numpy 1.5.1) is: points_tuple 0.36815404892 seconds points_numpy 6.20338392258 seconds I do not need C performance here, but the latter is definitely too slow. In the real program it's not so easy (but possible) to use tuples because the code is dimension-independent. Before considering writing an own "small vector" module, I'd like about other possible solutions. Other people must have stumbled across this before! Thanks, Christoph

6 12

genfromtxt
by Nils Wagner 11 Oct '11

11 Oct '11

Hi all, How do I use genfromtxt to read a file with the following lines 1 1 2.2592365264892578D+01 2 2 2.2592365264892578D+01 1 3 2.6666669845581055D+00 3 3 2.2592365264892578D+01 2 4 2.6666669845581055D+00 4 4 2.2592365264892578D+01 3 5 2.6666669845581055D+00 5 5 2.2592365264892578D+01 4 6 2.6666669845581055D+00 6 6 2.2592365264892578D+01 1 7 2.9814243316650391D+00 7 7 1.7259031295776367D+01 2 8 2.9814243316650391D+00 8 8 1.7259031295776367D+01 ... names =("i","j","v") A = np.genfromtxt('bmll.mtl',dtype=[('i','int'),('j','int'),('v','d')],names=names) V = A[:]['v'] >>> V array([ NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN, NaN]) yields NaN, while convertfunc = lambda x: x.replace('D','E') names =("i","j","v") A = np.genfromtxt('bmll.mtl',dtype=[('i','int'),('j','int'),('v','|S24')],names=names,converters={"v":convertfunc}) V = A[:]['v'].astype(float) >>> V array([ 22.59236526, 22.59236526, 2.66666698, 22.59236526, 2.66666698, 22.59236526, 2.66666698, 22.59236526, 2.66666698, 22.59236526, 2.98142433, 17.2590313 , 2.98142433, 17.2590313 , 2.98142433, 2.98142433, 2.66666698, 22.59236526, 2.98142433, 2.98142433, 2.66666698, 22.59236526, 2.98142433, 2.98142433, 2.66666698, 22.59236526, 2.98142433, 2.98142433, 2.66666698, 22.59236526, 2.98142433, 2.66666698, 17.2590313 , 2.98142433, 2.66666698, 17.2590313 ]) works fine. Nils

2 1

kind of a matrix multiplication
by Martin Raspaud 11 Oct '11

11 Oct '11

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Hi all, I have a stack of vectors: v1 = np.arange(3) v2 = np.arange(3) + 3 stack = np.vstack(v1, v2) (now stack is : array([[0, 1, 2], [3, 4, 5]])) and a 3d matrix: mat = np.dstack((np.eye(3), np.eye(3) * 2)) (mat is now array([[[ 1., 2.], [ 0., 0.], [ 0., 0.]], [[ 0., 0.], [ 1., 2.], [ 0., 0.]], [[ 0., 0.], [ 0., 0.], [ 1., 2.]]])) I'm looking for the operation needed to get the two (stacked) vectors array([[0, 1, 2], [6, 8, 10]])) or its transpose. I tried various combinations of tensor products, but I always get a result in 3 dimensions, while I just want two. Any suggestions ? Thanks, Martin -----BEGIN PGP SIGNATURE----- Version: GnuPG v2.0.14 (GNU/Linux) Comment: Using GnuPG with Red Hat - http://enigmail.mozdev.org/ iQEcBAEBAgAGBQJOk9jWAAoJEBdvyODiyJI4y30IAJu6YIHK+ED8pN5M2TFrEKj8 k/K22MjitlQ8wTFDxwc5xBRI+yoniqgAfpzWjdU3pc5MxzXRgbZrRZagYWjepZyI CtN/CHy+BfM8EPJulFeVcInAgo1pgfAhH4xwEakbu88XhKSgat1Y9xlNRcrohTUQ oBVd+DNmBYpEUAa0pDjkMYXM8vaJqzePZZGaviZxY0AY2MBDrbZN/z6t4u2Unajn 8X1vjCg/XfDbm9v7FK/52MUorAJinZRdHiWBTE9rOmAqjJxTBoFKkN+0FMTUk6Sj acJNjr5KFjl6o3JPxqU4jRfw1zFRO9BEouzosKfYcs/kLozNjTBmfZztg0np/dg= =Ykij -----END PGP SIGNATURE-----

2 2

Lookup array
by Andrey N. Sobolev 11 Oct '11

11 Oct '11

Hi all, I have 2 arrays - A with the dimensions of 1000x4 and B with the dimensions of 5000x4. B doesn't (hopefully) contain any rows that are not in A. I need to create a lookup array C, the i-th value of which will be the index of B[i] in A. In the (very rare) case when B[i] is not in A C[i] should be equal to -1. Now I'm dealing with it using which: ... import numpy as np for i, row in enumerate(B): try: C[i] = np.which(np.all(A == row, axis = 1))[0][0] except (IndexError, ): C[i] = -1 ... but that's very slow (it consumes 70% of cpu time needed by the whole program). I guess that it's because of a slow pythonic loop, but I just can't get how to get rid of it. Any suggestions would be appreciated. Thanks in advance, Andrey -- Researcher, General and theoretical physics dept., South Ural State University 454080, Pr. Lenina, 76, Chelyabinsk, Russia Tel: +7 351 265-47-13 andrey(a)physics.susu.ac.ru

3 4

np.where and broadcasting
by Jesper Larsen 10 Oct '11

10 Oct '11

Hi numpy-users I have a 2d array of shape (ny, nx). I want to broadcast (copy) this array to a target array of shape (nt, nz, ny, nz) or (nt, ny, nx) so that the 2d array is repeated for each t and z index (corresponding to nt and nz). I am not sure how to do this (generic solution, for different target array shapes). Can anyone help? Best regards, Jesper

2 1

failure to build
by Amos Anderson 10 Oct '11

10 Oct '11

hello -- i'm trying to build numpy-1.6.1. so far, it works on several machines that i've tried, but now it's failing. I put some details on the machine below. I tried to change the compiler, but setup.py objected to the things I tried (--compiler=g++4, --compiler=icc). thanks for any advice! Amos. building 'numpy.core.multiarray' extension compiling C sources C compiler: gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC creating build/temp.linux-x86_64-2.7/numpy/core/src/multiarray compile options: '-Inumpy/core/include -Ibuild/src.linux-x86_64-2.7/numpy/core/include/numpy -Inumpy/core/src/private -Inumpy/core/src -Inumpy/core -Inumpy/core/src/npymath -Inumpy/core/src/multiarray -Inumpy/core/src/umath -Inumpy/core/include -I/home/amosa/triad/tools/python/include/python2.7 -Ibuild/src.linux-x86_64-2.7/numpy/core/src/multiarray -Ibuild/src.linux-x86_64-2.7/numpy/core/src/umath -c' gcc: numpy/core/src/multiarray/multiarraymodule_onefile.c numpy/core/src/multiarray/descriptor.c: In function `_convert_divisor_to_multiple': numpy/core/src/multiarray/descriptor.c:606: warning: 'q' might be used uninitialized in this function numpy/core/src/multiarray/einsum.c.src: In function `float_sum_of_products_contig_outstride0_one': numpy/core/src/multiarray/einsum.c.src:852: error: unrecognizable insn: (insn:HI 440 228 481 14 /usr/lib/gcc/x86_64-redhat-linux/3.4.6/include/xmmintrin.h:915 (set (reg:SF 148) (vec_select:SF (reg/v:V4SF 67 [ accum_sse ]) (parallel [ (const_int 0 [0x0]) ]))) -1 (insn_list 213 (nil)) (nil)) numpy/core/src/multiarray/einsum.c.src:852: internal compiler error: in extract_insn, at recog.c:2083 Please submit a full bug report, with preprocessed source if appropriate. See <URL:http://bugzilla.redhat.com/bugzilla> for instructions. Preprocessed source stored into /tmp/ccWyirX3.out file, please attach this to your bugreport. numpy/core/src/multiarray/descriptor.c: In function `_convert_divisor_to_multiple': numpy/core/src/multiarray/descriptor.c:606: warning: 'q' might be used uninitialized in this function numpy/core/src/multiarray/einsum.c.src: In function `float_sum_of_products_contig_outstride0_one': numpy/core/src/multiarray/einsum.c.src:852: error: unrecognizable insn: (insn:HI 440 228 481 14 /usr/lib/gcc/x86_64-redhat-linux/3.4.6/include/xmmintrin.h:915 (set (reg:SF 148) (vec_select:SF (reg/v:V4SF 67 [ accum_sse ]) (parallel [ (const_int 0 [0x0]) ]))) -1 (insn_list 213 (nil)) (nil)) numpy/core/src/multiarray/einsum.c.src:852: internal compiler error: in extract_insn, at recog.c:2083 Please submit a full bug report, with preprocessed source if appropriate. See <URL:http://bugzilla.redhat.com/bugzilla> for instructions. Preprocessed source stored into /tmp/ccWyirX3.out file, please attach this to your bugreport. error: Command "gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -Inumpy/core/include -Ibuild/src.linux-x86_64-2.7/numpy/core/include/numpy -Inumpy/core/src/private -Inumpy/core/src -Inumpy/core -Inumpy/core/src/npymath -Inumpy/core/src/multiarray -Inumpy/core/src/umath -Inumpy/core/include -I/home/amosa/triad/tools/python/include/python2.7 -Ibuild/src.linux-x86_64-2.7/numpy/core/src/multiarray -Ibuild/src.linux-x86_64-2.7/numpy/core/src/umath -c numpy/core/src/multiarray/multiarraymodule_onefile.c -o build/temp.linux-x86_64-2.7/numpy/core/src/multiarray/multiarraymodule_onefile.o" failed with exit status 1 cat /etc/redhat-release Red Hat Enterprise Linux AS release 4 (Nahant Update 4) $ svnversion . /repos/svn/trunk 4168 gcc: _configtest.c gcc -pthread _configtest.o -L/usr/lib64 -lf77blas -lcblas -latlas -o _configtest ATLAS version 3.7.11 built by root on Mon Jun 5 10:14:12 EDT 2006: UNAME : Linux intel1.lsf.platform.com 2.6.9-34.ELsmp #1 SMP Fri Feb 24 16:56:28 EST 2006 x86_64 x86_64 x86_64 GNU/Linux INSTFLG : MMDEF : /export/madison/src/roll/hpc/BUILD/ATLAS/CONFIG/ARCHS/P4E64SSE3/gcc/gemm ARCHDEF : /export/madison/src/roll/hpc/BUILD/ATLAS/CONFIG/ARCHS/P4E64SSE3/gcc/misc F2CDEFS : -DAdd__ -DStringSunStyle CACHEEDGE: 393216 F77 : /usr/bin/g77, version GNU Fortran (GCC) 3.4.5 20051201 (Red Hat 3.4.5-2) F77FLAGS : -fomit-frame-pointer -O -m64 CC : /usr/bin/gcc, version gcc (GCC) 3.4.5 20051201 (Red Hat 3.4.5-2) CC FLAGS : -fomit-frame-pointer -O3 -funroll-all-loops -m64 MCC : /usr/bin/gcc, version gcc (GCC) 3.4.5 20051201 (Red Hat 3.4.5-2) MCCFLAGS : -fomit-frame-pointer -O -m64 success!

1 0

scikits.image 0.3 release
by Stéfan van der Walt 10 Oct '11

10 Oct '11

Announcement: scikits.image 0.3 =============================== After a brief (!) absence, we're back with a new and shiny version of scikits.image, the image processing toolbox for SciPy. This release runs under all major operating systems where Python (>=2.6 or 3.x), NumPy and SciPy can be installed. For more information, visit our website http://scikits-image.org or the examples gallery at http://scikits-image.org/docs/0.3/auto_examples/ New Features ------------ - Shortest paths - Total variation denoising - Hough and probabilistic Hough transforms - Radon transform with reconstruction - Histogram of gradients - Morphology, including watershed, connected components - Faster homography transformations (rotations, zoom, etc.) - Image dtype conversion routines - Line drawing - Better image collection handling - Constant time median filter - Edge detection (canny, sobel, etc.) - IO: Freeimage, FITS, Qt and other image loaders; video support. - SIFT feature loader - Example data-sets ... as well as many bug fixes and minor updates. Contributors for this release ----------------------------- Martin Bergtholdt Luis Pedro Coelho Chris Colbert Damian Eads Dan Farmer Emmanuelle Gouillart Brian Holt Pieter Holtzhausen Thouis (Ray) Jones Lee Kamentsky Almar Klein Kyle Mandli Andreas Mueller Neil Muller Zachary Pincus James Turner Stefan van der Walt Gael Varoquaux Tony Yu

1 0

numpy.lib.npyio.load
by Jean-Louis Durrieu 10 Oct '11

10 Oct '11

Hi everyone, I was just wondering something: lately, I had to use the load function, to load arrays stored in npz files. During one session, I need to read quite a few times several files (or even the same files), for some model training. I however just found out that the batch processing I ran failed because of a "too many open files" problem. After checking, with lsof, it seems that the use of np.load(filename), where filename is a string (= path to the file), worked an unexpected way. When I do the following, in a ipython 0.11 session, with the --pylab option : In [1]: np.__version__ Out[1]: '1.6.1' In [2]: np.load Out[2]: <function numpy.lib.npyio.load> In [3]: struc = np.load('path/to/file.npz') In [4]: ar1 = struc['ar1'] I would expect to have opened a file, read the array in it, and closed it. However, 'lsof' proved me wrong, and I found out that I need to explicitly do 'struc.close()' in order to close the file. While this is not a big issue, since one has to do that anyway when opening/closing files, I was a bit surprised that I had to do it when using the load function. Maybe this behaviour should be made explicit in the documentation? (at least, in the help(load), nothing is said about that, as far as I could read). I just wanted to share that, just in case someone was crossing the same type of issues. I hope this is going to be useful for others! Best regards, Jean-Louis Durrieu

2 1

what python module to modify NetCDF data?
by Chao YUE 09 Oct '11

09 Oct '11

Dear all, I want to change some variable values in a series of NetCDF file. Did anybody else did this before using python? Now I use pupynere for reading data from NetCDF files and making plots. but the document of pupynere for writing data to NetCDF file is quite simple and I still feel difficult to do this with pupynere. the NetCDF file I want to change is a global data (0.5X0.5d resolution, 360X720grid with 12 time steps) and have approx. 10 variables. I just want to change some points for a specific variable for all 12 time steps. I know it's possible use NCO ncap2 utility to do the job. but now I have some problem in using ncap2 within a shell script. I guess there is some easy way to use some python module to do the job? like mainly altering the data that need to change while let the others remaining intact? Any idea will be greatly appreciated. I will all a good weekend, Chao -- *********************************************************************************** Chao YUE Laboratoire des Sciences du Climat et de l'Environnement (LSCE-IPSL) UMR 1572 CEA-CNRS-UVSQ Batiment 712 - Pe 119 91191 GIF Sur YVETTE Cedex Tel: (33) 01 69 08 29 02; Fax:01.69.08.77.16 ************************************************************************************

6 6

Non-rectangular neighborhoods
by T J 09 Oct '11

09 Oct '11

While reading the documentation for the neighborhood iterator, it seems that it can only handle rectangular neighborhoods. Have I understood this correctly? If it is possible to do non-rectangular regions, could someone post an example/sketch of how to do this?

1 0