Mailman 3 October 2007 - NumPy-Discussion

Memory leak in ndarray, more info
by Robert Crida 29 Oct '07

29 Oct '07

Hi all I recently posted about a memory leak in numpy and failed to mention the version. The leak manifests itself in numpy-1.0.3.1 but is not present in numpy-1.0.2 The following code reproduces the bug: import numpy as np a = np.array([1.0, 2.0, 3.0]) while True: b = str(a) What happens above is that is repeatedly converted to a string. The process size grow quite rapidly. Has anyone else come across this? Where do I look to try to correct it? Thanks Robert

5 7

Array struct
by Sean Ross-Ross 28 Oct '07

28 Oct '07

Hi, I am looking to create an __array_struct__ member for a c++ vector<vector<cppObj>> i have wrapped in python, where cppObj contains an array of floats. Currently I have created an __array_struct__ for the cppObj and in python I iterate over the vector of vectors and put all of the data into a 1-d array. Is it possible to do this with an __array_struct__ within c++ for the whole object at once? Perhaps creating a float **. Thanks for your help ~Sean

2 1

fortran array storage question
by Georg Holzmann 26 Oct '07

26 Oct '07

Hallo! I have the following problem: I get a data array in column major storage order and want to use it as numpy array without copying data. Therefore I do the following (2D example): obj = PyArray_FromDimsAndData(2, dim0, PyArray_DOUBLE, (char*)data); PyArrayObject *tmp = (PyArrayObject*)obj; tmp->flags = NPY_FARRAY; In python I can see that the flag F_CONTIGUOUS is true, but the python numpy array is still displayed in the wrong (C-style) order, same as if I don't set the NPY_FARRAY flag. Do I have to consider something more - or are there any other problems ? (sorry - I'm new to numpy ;) Thanks for any hint, LG Georg

4 6

Bug in scalartypes.inc.src
by Nils Wagner 26 Oct '07

26 Oct '07

Hi all, I cannot install numpy from recent svn creating build/temp.linux-x86_64-2.5 creating build/temp.linux-x86_64-2.5/numpy creating build/temp.linux-x86_64-2.5/numpy/core creating build/temp.linux-x86_64-2.5/numpy/core/src compile options: '-Ibuild/src.linux-x86_64-2.5/numpy/core/src -Inumpy/core/include -Ibuild/src.linux-x86_64-2.5/numpy/core -Inumpy/core/src -Inumpy/core/include -I/usr/include/python2.5 -c' gcc: numpy/core/src/multiarraymodule.c In file included from numpy/core/src/arrayobject.c:509, from numpy/core/src/multiarraymodule.c:97: numpy/core/src/scalartypes.inc.src: In function gentype_reduce: numpy/core/src/scalartypes.inc.src:1312: error: label at end of compound statement In file included from numpy/core/src/arrayobject.c:509, from numpy/core/src/multiarraymodule.c:97: numpy/core/src/scalartypes.inc.src: In function gentype_reduce: numpy/core/src/scalartypes.inc.src:1312: error: label at end of compound statement error: Command "gcc -pthread -fno-strict-aliasing -DNDEBUG -O2 -fmessage-length=0 -Wall -D_FORTIFY_SOURCE=2 -g -fPIC -Ibuild/src.linux-x86_64-2.5/numpy/core/src -Inumpy/core/include -Ibuild/src.linux-x86_64-2.5/numpy/core -Inumpy/core/src -Inumpy/core/include -I/usr/include/python2.5 -c numpy/core/src/multiarraymodule.c -o build/temp.linux-x86_64-2.5/numpy/core/src/multiarraymodule.o" failed with exit status 1 Nils

1 0

C-API for non-contiguous arrays
by Oliver Kranz 26 Oct '07

26 Oct '07

Hi, I am working on a Python extension module using of the NumPy C-API. The extension module is an interface to an image processing and analysis library written in C++. The C++ functions are exported with boos::python. Currently I am implementing the support of three-dimensional data sets which can consume a huge amount of memory. The 3D data is stored in a numpy.ndarray. This array is passed to C++ functions which do the calculations. In general, multi-dimensional arrays can be organized in memory in four different ways: 1. C order contiguous 2. Fortran order contiguous 3. C order non-contiguous 4. Fortran order non-contiguous Am I right that the NumPy C-API can only distinguish between three ways the array is organized in memory? These are: 1. C order contiguous e.g. with PyArray_ISCONTIGUOUS(arr) 2. Fortran order contiguous e.g. with PyArray_ISFORTRAN(arr) 3. non-contiguous e.g. with !PyArray_ISCONTIGUOUS(arr) && !PyArray_ISFORTRAN(arr) So there is no way to find out if a non-contiguous array has C order or Fortran order. The same holds for Python code e.g. by use of the flags. a.flags.contiguous a.flags.fortran This is very important for me because I just want to avoid to copy every non-contiguous array into a contiguous array. This would be very inefficient. But I can't find an other solution than copying the array. Also the iterator provided by the C-API only loops over the array in C order. Even if the array is in Fortran non-contiguous order. Or are there just no Fortran order non-contiguous arrays? I think I can construct one. a = numpy.ndarray((3,4,5), order="F") b = a[:,1:2,:] Now, I think b's elements are organized in memory in Fortran non-contiguous order. But the flags only tell me that it is non-contiguous and not if it is in Fortran order or in C order. And if b would be passed to a C++ function it would not be possible to find out with the C-API if it is in Fortran order or in C order, too. Any ideas? Or do I always have to create contiguous arrays? Cheers, Oliver

3 4

[Announce] numpy.scons , ALPHA version
by David Cournapeau 26 Oct '07

26 Oct '07

Hi there, I've finally managed to implement most of the things I wanted for numpy.scons, hence a first alpha. Compared to the 2d milestone from a few days ago, a few optimized libraries are supported (ATLAS, Sun sunperf, Apple Accelerate and vecLib, Intel MKL). Who === Outside people interested in numpy.scons, people who have trouble building numpy may want to try this. In particular, using MKL or sunperf should be easier (for sunperf, it should work out of the box if sun compilers are used). Also, compiling with intel compilers on linux should be easier. DO NOT USE THIS FOR PRODUCTION USE ! I am interested in success (that is build and all the test pass) as well as failure. Numpy developers: I would be interested in comments about the code: - the package are built with the setupscons.py/SConstruct - most of the new code is in numpy/distutils/scons How to test =========== Grab the sources: svn co http://svn.scipy.org/svn/numpy/branches/numpy.scons And then just do: python setup.py scons; python setup.py install. If you have multiple CPU/Cores, you may want to use the parallel build options: python setup.py scons --jobs=N with N an integer (number of CPU/CORES + 1 is a rule of thumb). For non default compilers, use the --compiler and --fcompiler options (distutils names). For the MKL, you should add a mkl option in the site.cfg. For ATLAS or sunperf, if it is not in standard places, it won't work (yet). What works: =========== The tested platforms have been tested: - linux (ubuntu x86) with gcc, g77 and ATLAS. - linux (x86 and x64) with gcc, without BLAS or LAPACK. - linux (ubuntu x86) with gcc, g77, and mkl (9.1.023). - linux (ubuntu x86) with intel compilers. - windows (32 bits) with mingw. - windows (32 bits) with VS 2003. - solaris (solaris studio express 12) with sunperf on x86. - Mac Os X Intel with Accelerate framework. - People reported basic success on IRIX, too. What does NOT work ================== One important target missing is windows 64, but this should not be too difficult to solve. There are still many corner cases not yet solved (in particular some windows things, most libraries path cannot yet be overriden in site.cfg); also, I do not attempt to be 100 % backward compatible with the current build system (that for a given environment, you won't necessarily have the same build options). Some things are not yet supported either (overriding libraries with env variables, for example, is not supported). There is no real doc yet (I will do this once the basic public API is stabilized). Details ======= Library checkers ---------------- The main difference compared to the second milestone is the support for libraries checkers (think BLAS, Atlas, etc...). There is now support to check for different libraries, with different compilers. This effectively replaces numpy.distutils.system_info. The way I've implemented those checks are fundamentally different compared to the distutils way: instead of looking for files, I look for capabilities at runtime. For example, for mkl, instead of looking for libmkl.so, etc... I try to compile, link and run a small mkl program: if any of the step failed, then the mkl is not used. This means that I can detect config errors much sooner; also, since I use rpath, there is no need to set things like LD_LIBRARY_PATH and so on. If the path is set right in site.cfg, then everything is taken care of. Also, as before, a general library checker for new libraries is available: look at basic examples in numpy/scons_fake/. In particular, there is a NumpyCheckLib which can be used to check a library with some symbols, whose paths can be overriden with a site.cfg. Fortran support --------------- The other big work has been on fortran support. Here too, the way I implemented it is different than numpy.distutils. Instead of harcoding flags depending on versions, I use code snipets, to get link options, mangling, and so on. For example, if you call the following in a Sconscript: config.CheckF77Mangling() This will try to get the mangling of the F77 compiler, and if sucessfull, will set a python function mangle, which returns the fortran name from the 'common' name. For example, for g77, mangle('dgemm') will returns dgemm_. Since the mangling is found at runtime, this can be used for not yet known compilers. Next ==== This was quite a lot of work, and I wanted to get something working as fast as possible. As such, the code is not always clean. Once I am satisfied with the global 'feel' of the new build, I will sanitize the API and implement the doc. Still, I believe the code to be much easier, simpler and robust than distutils (everything, from fortran support to libraries checkers is less than 2000 lines of python code, including tests; of course, this is largely thanks to scons). The thing I am the less sure about is how to combine checkers in meta-checkers: for example, I have checker for mkl, ATLAS, sunperf, and we sometimes want a BLAS checker, which would check all those as long as nothing is found. API-wise, I am not yet satisfied with the current way. I also want to play nice with f2py and scipy (totally untested for now).

2 1

Memory leak in ndarray
by Robert Crida 26 Oct '07

26 Oct '07

Hi all I seem to have tracked down a memory leak in the string conversion mechanism of numpy. It is demonstrated using the following code: import numpy as np a = np.array([1.0, 2.0, 3.0]) while True: b = str(a) What happens above is that is repeatedly converted to a string. The process size grow quite rapidly. Has anyone else come across this? Where do I look to try to correct it? Thanks Robert

5 8

Best way to handle "wrapped" array accesses?
by Victor Putz 26 Oct '07

26 Oct '07

I've seen a few references on this, but hadn't found a proper solution... I'm doing Lattice-Boltzmann simulations with periodic boundary conditions, which always necessarily involve either padding the edges and doing additional steps, or making a "wrapped" array (for example, if I have an array that's N rows by M columns, and cells are moving right, I then make a "wrapped" array a la "D = hstack( A[:,0:-1], A[:,-1:] )" (creating a temporary) and adding them. I'd like to avoid the temporary and still make the operation as fast as possible. A really elegant solution would be a "wrapped array" subclass that would automatically handle the indexing (C = A[:,1:M+1] + B), but I'm not certain how friendly the array class is with subclassing. What's the best way to handle this? Am I pretty much stuck with temporaries or hacked functions ("add_wrapped(A, offset, B, C)")? (It's probably a moot point in the long term, as once things get big I will have to use the cluster, which is probably not SciPy-friendly and it'll likely be a bit before I can convince the sys admin to let me play outside the "installed software" boundaries. But now I just want to know). Thanks, -->VPutz

1 0

convolution and wiener khinchin
by John Hunter 25 Oct '07

25 Oct '07

I am working on an example to illustrate convolution in the temporal and spectral domains, using the property that a convolution in the time domain is a multiplication in the fourier domain. I am using numpy.fft and numpy.convolve to compute the solution two ways, and comparing them. I am getting an error for small times in my fourier solution. At first I thought this was because of edge affects, but I see it even when I apply a windowing function. Can anyone advise me about what I am doing wrong here? """ In signal processing, the output of a linear system to an arbitrary input is given by the convolution of the impulse response function (the system response to a Dirac-delta impulse) and the input signal. Mathematically: y(t) = \int_0^\t x(\tau)r(t-\tau)d\tau where x(t) is the input signal at time t, y(t) is the output, and r(t) is the impulse response function. In this exercise, we will compute investigate the convolution of a white noise process with a double exponential impulse response function, and compute the results 3 ways: * using numpy.convolve * in Fourier space using the property that a convolution in the temporal domain is a multiplication in the fourier domain """ import numpy as npy import matplotlib.mlab as mlab from pylab import figure, show # build the time, input, output and response arrays dt = 0.01 t = npy.arange(0.0, 10.0, dt) # the time vector from 0..5 Nt = len(t) def impulse_response(t): 'double exponential response function' return (npy.exp(-t) - npy.exp(-5*t))*dt win = npy.hanning(Nt) x = npy.random.randn(Nt) # gaussian white noise x = win*x r = impulse_response(t)*dt # evaluate the impulse function r = win*r y = npy.convolve(x, r, mode='full') # convultion of x with r y = y[:Nt] # plot t vs x, t vs y and t vs r in three subplots fig = figure() ax1 = fig.add_subplot(311) ax1.plot(t, x) ax1.set_ylabel('input x') ax2 = fig.add_subplot(312) ax2.plot(t, y, label='convolve') ax2.set_ylabel('output y') ax3 = fig.add_subplot(313) ax3.plot(t, r) ax3.set_ylabel('input response') ax3.set_xlabel('time (s)') # compute y via numerical integration of the convolution equation F = npy.fft.fft(r) X = npy.fft.fft(x) Y = F*X yi = npy.fft.ifft(Y).real ax2.plot(t, yi, label='fft') ax2.legend(loc='best') show()

2 1

Re: [Numpy-discussion] how to select indices from boolean array
by Pauli Virtanen 25 Oct '07

25 Oct '07

wBernhard Voigt <bernhard.voigt(a)gmail.com> kirjoitti: > is nonzero the appropriate function to get the indexes of a boolean array? > > foo = numpy.random.random(10000) > cut = foo > .5 > indexes = cut.nonzero()[0] I think writing indices = where(foo > 5)[0] might be more clear, although it does the same thing as .nonzero(). > Another question: Why is the the type returned by nonzero a tuple and not an > ndarray of the same shape as the input (or self if used on an object)? Because indexing with booleans in numpy works differently from indexing with numbers, especially for multidimensional arrays. You want to have ( foo[foo > 5] == foo[where(foo > 5)] ).all() so where/nonzero returns a tuple that contains 1D-arrays, one per dimension, which together enumerate the relevant entries. (IIRC, in Matlab, find returns indices in the flattened array.) -- Pauli Virtanen --

1 0