Mailman 3 July 2010 - NumPy-Discussion

np.asfortranarray: unnecessary copying?
by Kurt Smith July 31, 2010

July 31, 2010

What are the rules for when 'np.asarray' and 'np.asfortranarray' make a copy? This makes sense to me: In [3]: carr = np.arange(3) In [6]: carr2 = np.asarray(carr) In [8]: carr2[0] = 1 In [9]: carr Out[9]: array([1, 1, 2]) No copy is made. But doing the same with a fortran array makes a copy: In [10]: farr = np.arange(3).copy('F') In [12]: farr2 = np.asfortranarray(farr) In [13]: farr2[0] = 1 In [14]: farr Out[14]: array([0, 1, 2]) Could it be a 1D thing, since it's both C contiguous & F contiguous? Here's a 2D example: In [15]: f2D = np.arange(10).reshape((2,5), order='F') In [17]: f2D2 = np.asfortranarray(f2D) In [19]: f2D2[0,0] = 10 In [20]: f2D Out[20]: array([[10, 2, 4, 6, 8], [ 1, 3, 5, 7, 9]]) So it looks like np.asfortranarray makes an unnecessary copy if the array is simultaneously 1D, C contiguous and F contiguous. Coercing the array with np.atleast_2d() makes asfortranarry behave. Looking further, np.isfortran always returns false if the array is 1D, even if it's Fortran contiguous (and np.isfortran is documented as such). What is the rationale here? Is it a 'column' vs. 'row' thing? Kurt

2 2

Unable to get most basic numpy array creation from C++ working.
by rif July 31, 2010

July 31, 2010

Hi numpy-discuss . I've distilled my problem down to a very simple program that segfaults. I'm sure I'm doing something silly, but hopefully someone can point me in the right direction. I create a program called numpy_test.cc: --------------------- #include <Python.h> #include "numpy/arrayobject.h" int main(int argc, char **argv) { Py_Initialize(); npy_intp length[1]; length[0] = 10; PyObject* my_array = PyArray_SimpleNew(1, length, NPY_FLOAT); return 0; } --------------------- Standing in the directory above where I've done an in-place build of numpy-1.3.0, I compile with: g++ numpy_test.cc -I/usr/include/python2.5 -I/tmp/numpy-1.3.0/numpy/core/include -lpython2.5* *It compiles clean, but the resulting executable crashes in the PyArray_SimpleNew line [verified via debugger]. I have yet to get a call to PyArray_SimpleNew to not segfault. Any thoughts or advice are appreciated. Cheers, rif * *

2 2

test errors in the trunk
by Darren Dale July 31, 2010

July 31, 2010

I just upgraded my svn checkout and did a fresh install. When I try to run the test suite, I get a ton of errors: np.test() Running unit tests for numpy NumPy version 2.0.0.dev8550 NumPy is installed in /Users/darren/.local/lib/python2.6/site-packages/numpy Python version 2.6.5 (r265:79063, Jul 19 2010, 09:08:11) [GCC 4.2.1 (Apple Inc. build 5659)] nose version 0.11.3 ................................................................................................................Reloading numpy.lib Reloading numpy.lib.info Reloading numpy.lib.numpy Reloading numpy Reloading numpy.numpy Reloading numpy.show EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE ====================================================================== [...] File "/Users/darren/.local/lib/python2.6/site-packages/numpy/lib/__init__.py", line 23, in <module> __all__ += type_check.__all__ NameError: name 'type_check' is not defined I checked numpy/lib/__init__.py, and it does a bunch of imports like "from type_check import *" but not "import type_check", which are needed to append to __all__. Darren

3 3

Fwd: [Matplotlib-users] Vectorization
by George Nurser July 30, 2010

July 30, 2010

---------- Forwarded message ---------- From: George Nurser <gnurser(a)gmail.com> Date: 30 July 2010 22:37 Subject: Re: [Matplotlib-users] Vectorization To: Nicolas Bigaouette <nbigaouette(a)gmail.com>, Discussion of Numerical Python <numpy-discussion(a)scipy.org> > I want to do the same for the calculation of the kinetic energy: > <phi|p^2|phi>/2m. There is a laplacian in the volume integral which > complicates things: >> >> K = 0.0 >> for i in numpy.arange(len(dx)-1): >> for j in numpy.arange(len(dy)-1): >> for k in numpy.arange(len(dz)-1): >> K += -0.5 * m * phi[k,j,i] * ( >> (phi[k,j,i-1] - 2.0*phi[k,j,i] + phi[k,j,i+1]) / >> dx[i]**2 >> + (phi[k,j-1,i] - 2.0*phi[k,j,i] + phi[k,j+1,i]) / >> dy[j]**2 >> + (phi[k-1,j,i] - 2.0*phi[k,j,i] + phi[k+1,j,i]) / >> dz[k]**2 >> ) If dx ,dy, dz vary, I don't think this is quite correct. To get dphi/dz at the interface between k & k+1 you need the distance between the psi-points inside the kth and k+1th gridbox. If you assume the psi-points lie at the centres of the gridboxes (other assumptions are possible), this distance is .5*(dz[k]+dz[k+1]) So, writing rdzw[k] = 2./(dz[k]+dz[k+1]) d/dz(dphi/dz))[k] = {(phi[k+1]-phi[k])*rdzw[k] - (phi[k]-phi[k-1])*rdz[k-1]} / dz[k] = {phi[k+1]*rdz[k]- phi[k]*(rdz[k] + rdz[k-1]) + phi[k-1]*rdz[k-1]} / dz[k] & similarly for the d2phi/dx2 & d2phi/dy2 terms There are other metric terms which appear in the discretization of the Laplacian if dz depends on i or j, or dy depends on k or i, or dx depends on j or k, but these do not appear in your problem where dz depends only on k, dy only on j and dx only on i. HTH. George Nurser.

1 0

endian.h change
by Ralf Gommers July 30, 2010

July 30, 2010

Hi David, Commit r8541 broke building with numscons for me, does this fix look okay: http://github.com/rgommers/numpy/commit/1c88007ab00cf378ebe19fbe54e9e868212… I am puzzled though why my endian.h is not picked up in the build - I have a good collection of those on my system, at least in all OS X SDKs. Any idea? Cheers, Ralf

3 2

Re: [Numpy-discussion] [Matplotlib-users] Vectorization
by Nicolas Bigaouette July 30, 2010

July 30, 2010

2010/7/2 Charles R Harris <charlesr.harris(a)gmail.com> > > > On Fri, Jul 2, 2010 at 12:15 PM, Nicolas Bigaouette <nbigaouette(a)gmail.com > > wrote: > >> Hi all, >> >> I don't really know where to ask, so here it is. >> >> I was able to vectorize the normalization calculation in quantum >> mechanics: <phi|phi>. Basically it's a volume integral of a scalar field. >> Using: >> >>> norm = 0.0 >>> for i in numpy.arange(len(dx)-1): >>> for j in numpy.arange(len(dy)-1): >>> for k in numpy.arange(len(dz)-1): >>> norm += psi[k,j,i]**2 * dx[i] * dy[j] * dz[k] >>> >> if dead slow. I replaced that with: >> >>> norm = (psi**2 * >>> dx*dy[:,numpy.newaxis]*dz[:,numpy.newaxis,numpy.newaxis]).sum() >>> >> which is almost instantanious. >> >> I want to do the same for the calculation of the kinetic energy: >> <phi|p^2|phi>/2m. There is a laplacian in the volume integral which >> complicates things: >> >>> K = 0.0 >>> for i in numpy.arange(len(dx)-1): >>> for j in numpy.arange(len(dy)-1): >>> for k in numpy.arange(len(dz)-1): >>> K += -0.5 * m * phi[k,j,i] * ( >>> (phi[k,j,i-1] - 2.0*phi[k,j,i] + phi[k,j,i+1]) / >>> dx[i]**2 >>> + (phi[k,j-1,i] - 2.0*phi[k,j,i] + phi[k,j+1,i]) / >>> dy[j]**2 >>> + (phi[k-1,j,i] - 2.0*phi[k,j,i] + phi[k+1,j,i]) / >>> dz[k]**2 >>> ) >>> >> >> My question is, how would I vectorize such loops? I don't know how I would >> manage the "numpy.newaxis" code-foo with neighbours dependency... Any idea? >> >> Thanx! >> >> > Wouldn't the Laplacian look better in momentum space? That is, do a 3d > Fourier transform, multiply by the appropriate weight, and do the integral > in the transform space. > > Chuck > Thanx all for suggestions. The best approach was the convolution. I couldn't do a single convolution since the cell sizes are the same in x, y and z. I had to separate the laplacian into a sum over the three dimensions and convolve each individually. There is also terms appearing because of the non-uniform grid. Not that I work in atomic units, so m=1. I'm pasting the code I've finally written: > def Energy_Kinetic(phi, J1x, J1y, J1z, J2x, J2y, J2z): > > kernel_x1 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], > [[0,0,0], [1,-2,1], [0,0,0]], > [[0,0,0], [0,0,0], [0,0,0]]]) > kernel_y1 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], > [[0,1,0], [0,-2,0], [0,1,0]], > [[0,0,0], [0,0,0], [0,0,0]]]) > kernel_z1 = numpy.array([[[0,0,0], [0,1,0], [0,0,0]], > [[0,0,0], [0,-2,0], [0,0,0]], > [[0,0,0], [0,1,0], [0,0,0]]]) > > kernel_x2 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], > [[0,0,0], [1,0,-1], [0,0,0]], > [[0,0,0], [0,0,0], [0,0,0]]]) > kernel_y2 = numpy.array([[[0,0,0], [0,0,0], [0,0,0]], > [[0,1,0], [0,0,0], [0,-1,0]], > [[0,0,0], [0,0,0], [0,0,0]]]) > kernel_z2 = numpy.array([[[0,0,0], [0,1,0], [0,0,0]], > [[0,0,0], [0,0,0], [0,0,0]], > [[0,0,0], [0,-1,0], [0,0,0]]]) > > convolution = ( \ > scipy.signal.convolve(phi, kernel_x1, mode='same') / J1x**2 + scipy.signal.convolve(phi, kernel_y1, mode='same') / > J1y[:,numpy.newaxis]**2 > + scipy.signal.convolve(phi, kernel_z1, mode='same') / > J1z[:,numpy.newaxis,numpy.newaxis]**2 > + 0.5*(J2x / J1x**3)*scipy.signal.convolve(phi, kernel_x2, > mode='same') \ > + 0.5*(J2y[:,numpy.newaxis] / > J1y[:,numpy.newaxis]**3)*scipy.signal.convolve(phi, kernel_y2, mode='same') > \ > + 0.5*(J2z[:,numpy.newaxis,numpy.newaxis] / > J1z[:,numpy.newaxis,numpy.newaxis]**3)*scipy.signal.convolve(phi,kernel_z2, > mode='same') \ > ) > > K = -0.5 * (numpy.conjugate(phi) * convolution * J1x * > J1y[:,numpy.newaxis] * J1z[:,numpy.newaxis,numpy.newaxis]).real.sum() > > return K > # >

1 0

Compiling NumPy on 64 bit Solaris 10 sparc)
by Chris LeBlanc July 30, 2010

July 30, 2010

Hi, Firstly thanks to everyone that has helped bring NumPy to the point it is today. Its an excellent bit of software. I've recently managed to get NumPy to compile on a 64 bit Solaris 10 (sparc) machine. We've embedded 64 bit Python in some of our C code using the Python C-API, and therefore require a 64 bit version of NumPy. I had to use a bit of a hack to get NumPy to compile in 64 bit mode, and I'm just wondering if there's a more elegant way of doing it? If I set CFLAGS to "-m64 -xcode=pic32", and CC to "/usr/bin/cc", it will raise the following error: ... snipped most of backtrace ... File "/home/leblanc/source/numpy/numpy/distutils/command/build_src.py", line 385, in generate_sources source = func(extension, build_dir) File "numpy/core/setup.py", line 675, in get_mathlib_info raise RuntimeError("Broken toolchain: cannot link a simple C program") However, if I create a wrapper script and set CC to the path of the wrapper, then it will build normally and create a NumPy module I can use from the 64 bit version of Python installed on the system (/usr/bin/sparcv9/python). Here's the wrapper: #!/bin/sh /usr/bin/cc -m64 -xcode=pic32 $* I did a bit of googling and couldn't find much information on this issue (which is why I wanted to post this message, so others could find it). I'm satisfied with the results, but wondering if there's a better way to do it. Cheers, Chris Chris LeBlanc GNS Science

3 4

tofile problems
by Paul Probert July 29, 2010

July 29, 2010

All, I'm trying to write numpy arrays as binary data, to support a legacy file format. So I open a file and write to it: fp = open('somefile','w') ... oldpos = fp.tell() somenumpyarray.tofile(fp) newpos = fp.tell() diff = newpos - oldpos - somenumpyarray.nbytes if diff != 0: print 'ahhah! mismatch=',diff ... I'm observing that every once in while I get 'ahhah! mismatch=1', that is, the file position is advanced by one more byte than x.nbytes would give. I'm using 1.4.1 on windows, python 2.5. Any ideas? Thanks, Paul Probert University of Wisconsin

2 2

Re: [Numpy-discussion] Is there anyway to read raw binary file via pytable?
by 脑关(BrainGateway)生命科学仪器 July 29, 2010

July 29, 2010

2010/7/28 Ken Watford <kwatford+scipy(a)gmail.com>: > 2010/7/28 脑关生命科学仪器 <braingateway(a)gmail.com>: >> it seems like pytable only support HDF5. I had some 500GB numerical arrays >> to process. Pytable claims to have some advance feature to enhance >> processing speed and largely reduce physical memory requirement. However, I >> do not wanna touch the raw data I had. Simply because I do not have doubled >> diskspace to covert all 10TB data into HDF5. Is there any way to let pytable >> read raw binary files or alternatively to package raw files into HDF5 format >> without change the files themselves.? >> >> Thanks >> >> Brain Gateway > > HDF5 does support datasets with an external contiguous binary file as > the storage area. For documentation on this, see: > http://www.hdfgroup.org/HDF5/doc/UG/10_Datasets.html#Allocation Thanks a lot! I do found 3 API functions to handle this problem: H5P have 3 function to handle external raw data file. H5Pset_external http://davis.lbl.gov/Manuals/HDF5-1.6.1/RM_H5P.html#Property-SetExternal H5Pget_external_count http://davis.lbl.gov/Manuals/HDF5-1.6.1/RM_H5P.html#Property-GetExternalCou… H5Pget_external http://davis.lbl.gov/Manuals/HDF5-1.6.1/RM_H5P.html#Property-GetExternal However, right now h5py does not provide these functions, though existing datasets created with external storage should work fine. I need to use C, fortran, or matlab to generate datasets with external data. So far it works fine. LittleBigBrain

1 0

A bug in boolean indexing?
by Nadav Horesh July 29, 2010

July 29, 2010

The following does not raise an error: a = np.arange(5) a[a>0] = a although a.shape == (5,) while a[a>0].shape == (4,) I get in on python2.6.5, numpy 1.4.1 on win32, and python 2.6.5, numpy 2.0.0.dev8469 on linux64. Nadav.

2 2