Mailman 3 November 2012 - NumPy-Discussion

Re: [Numpy-discussion] example reading binary Fortran file
by Neil Martinsen-Burrell 22 Jul '21

22 Jul '21

David Froger <david.froger.info <at> gmail.com> writes: > Hy,My question is about reading Fortran binary file (oh no this question > again...) I've posted this before, but I finally got it cleaned up for the Cookbook. For this purpose I use a subclass of file that has methods for reading unformatted Fortran data. See http://www.scipy.org/Cookbook/FortranIO/FortranFile. I'd gladly see this in numpy or scipy somewhere, but I'm not sure where it belongs. > program makeArray > implicit none > integer,parameter:: nx=10,ny=20 > real(4),dimension(nx,ny):: ux,uy,p > integer :: i,j > open(11,file='uxuyp.bin',form='unformatted') > do i = 1,nx > do j = 1,ny > ux(i,j) = real(i*j) > uy(i,j) = real(i)/real(j) > p (i,j) = real(i) + real(j) > enddo > enddo > write(11) ux,uy > write(11) p > close(11) > end program makeArray When I run the above program compiled with gfortran on my Intel Mac, I can read it back with:: >>> import numpy as np >>> from fortranfile import FortranFile >>> f=FortranFile('uxuyp.bin', endian='<') >>> uxuy = f.readReals(prec='f') # 'f' for default reals >>> len(uxuy) 400 >>> ux = np.array(uxuy[:200]).reshape((20,10)).T >>> uy = np.array(uxuy[200:]).reshape((20,10)).T >>> p = f.readReals('f').reshape((20,10)).T >>> ux array([[ 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20.], [ 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 32., 34., 36., 38., 40.], [ 3., 6., 9., 12., 15., 18., 21., 24., 27., 30., 33., 36., 39., 42., 45., 48., 51., 54., 57., 60.], [ 4., 8., 12., 16., 20., 24., 28., 32., 36., 40., 44., 48., 52., 56., 60., 64., 68., 72., 76., 80.], [ 5., 10., 15., 20., 25., 30., 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 90., 95., 100.], [ 6., 12., 18., 24., 30., 36., 42., 48., 54., 60., 66., 72., 78., 84., 90., 96., 102., 108., 114., 120.], [ 7.,Proxy-Connection: keep-alive Cache-Control: max-age=0 14., 21., 28., 35., 42., 49., 56., 63., 70., 77., 84., 91., 98., 105., 112., 119., 126., 133., 140.], [ 8., 16., 24., 32., 40., 48., 56., 64., 72., 80., 88., 96., 104., 112., 120., 128., 136., 144., 152., 160.], [ 9., 18., 27., 36., 45., 54., 63., 72., 81., 90., 99., 108., 117., 126., 135., 144., 153., 162., 171., 180.], [ 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 110., 120., 130., 140., 150., 160., 170., 180., 190., 200.]]) >>> uy array([[ 1. , 0.5 , 0.33333334, 0.25 , 0.2 , 0.16666667, 0.14285715, 0.125 , 0.11111111, 0.1 , 0.09090909, 0.08333334, 0.07692308, 0.07142857, 0.06666667, 0.0625 , 0.05882353, 0.05555556, 0.05263158, 0.05 ], [ 2. , 1. , 0.66666669, 0.5 , 0.40000001, 0.33333334, 0.2857143 , 0.25 , 0.22222222, 0.2 , 0.18181819, 0.16666667, 0.15384616, 0.14285715, 0.13333334, 0.125 , 0.11764706, 0.11111111, 0.10526316, 0.1 ], [ 3. , 1.5 , 1. , 0.75 , 0.60000002, 0.5 , 0.42857143, 0.375 , 0.33333334, 0.30000001, 0.27272728, 0.25 , 0.23076923, 0.21428572, 0.2 , 0.1875 , 0.17647059, 0.16666667, 0.15789473, 0.15000001], [ 4. , 2. , 1.33333337, 1. , 0.80000001, 0.66666669, 0.5714286 , 0.5 , 0.44444445, 0.40000001, 0.36363637, 0.33333334, 0.30769232, 0.2857143 , 0.26666668, 0.25 , 0.23529412, 0.22222222, 0.21052632, 0.2 ], [ 5. , 2.5 , 1.66666663, 1.25 , 1. , 0.83333331, 0.71428573, 0.625 , 0.55555558, 0.5 , 0.45454547, 0.41666666, 0.38461539, 0.35714287, 0.33333334, 0.3125 , 0.29411766, 0.27777779, 0.2631579 , 0.25 ], [ 6. , 3. , 2. , 1.5 , 1.20000005, 1. , 0.85714287, 0.75 , 0.66666669, 0.60000002, 0.54545456, 0.5 , 0.46153846, 0.42857143, 0.40000001, 0.375 , 0.35294119, 0.33333334, 0.31578946, 0.30000001], [ 7. , 3.5 , 2.33333325, 1.75 , 1.39999998, 1.16666663, 1. , 0.875 , 0.77777779, 0.69999999, 0.63636363, 0.58333331, 0.53846157, 0.5 , 0.46666667, 0.4375 , 0.41176471, 0.3888889 , 0.36842105, 0.34999999], [ 8. , 4. , 2.66666675, 2. , 1.60000002, 1.33333337, 1.14285719, 1. , 0.8888889 , 0.80000001, 0.72727275, 0.66666669, 0.61538464, 0.5714286 , 0.53333336, 0.5 , 0.47058824, 0.44444445, 0.42105263, 0.40000001], [ 9. , 4.5 , 3. , 2.25 , 1.79999995, 1.5 , 1.28571427, 1.125 , 1. , 0.89999998, 0.81818181, 0.75 , 0.69230771, 0.64285713, 0.60000002, 0.5625 , 0.52941179, 0.5 , 0.47368422, 0.44999999], [ 10. , 5. , 3.33333325, 2.5 , 2. , 1.66666663, 1.42857146, 1.25 , 1.11111116, 1. , 0.90909094, 0.83333331, 0.76923078, 0.71428573, 0.66666669, 0.625 , 0.58823532, 0.55555558, 0.52631581, 0.5 ]]) >>> p array([[ 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21.], [ 3., 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22.], [ 4., 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23.], [ 5., 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24.], [ 6., 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25.], [ 7., 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26.], [ 8., 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27.], [ 9., 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28.], [ 10., 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29.], [ 11., 12., 13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24., 25., 26., 27., 28., 29., 30.]]) Note that you have to provide the shape information for ux and uy because fortran writes them together as a stream of 400 numbers. -Neil

10 10

C-API: multidimensional array indexing?
by Johann Bauer 31 Mar '16

31 Mar '16

Dear experts, is there a C-API function for numpy which implements Python's multidimensional indexing? Say, I have a 2d-array PyArrayObject * M; and an index int i; how do I extract the i-th row or column M[i,:] respectively M[:,i]? I am looking for a function which gives again a PyArrayObject * and which is a view to M (no copied data; the result should be another PyArrayObject whose data and strides points to the correct memory portion of M). I searched the API documentation, Google and mailing lists for quite a long time but didn't find anything. Can you help me? Thanks, Johann

3 3

SVD errors
by mtrumpis＠berkeley.edu 04 Jun '13

04 Jun '13

Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy. I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu) Just today I've come across a negative singular value cropping up in an SVD of a different matrix. This error does occur on my ATLAS LAPACK based numpy, as well as on the Ubuntu setup. And once again, it does not happen in Octave/Matlab. I'm using numpy 1.3.0.dev6336 -- don't know what the Ubuntu box is running. Here are some npy files for the two different cases: https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/noconverge_operator.npy https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/negsval_operator.npy Mike

6 7

pyhdf packaging
by Andreas Hilboll 23 May '13

23 May '13

Hi, I would like to package pyhdf for Ubuntu and make the package publicly available. Since the license is not totally clear to me (I cannot find any information in the sources, and the cheeseshop says "public", which doesn't mean anything to me), I tried to contact the maintainer, Andre Gosselin, however the email bounces, so I guess he's gone. Can anyone point me to how to proceed from here? Cheers, Andreas.

2 2

swig interface file (numpy.i) warning
by Tom Krauss 27 Feb '13

27 Feb '13

Hi, I've been happy to use numpy.i for generating SWIG interfaces to C++. For a while, I've noticed this warning while compiling: /Users/tkrauss/projects/dev_env/lib/python2.7/site-packages/numpy-1.8.0.dev_f2f0ac0_20120725-py2.7-macosx-10.8-x86_64.egg/numpy/core/include/numpy/npy_deprecated_api.h:11:2: warning: #warning "Using deprecated NumPy API, disable it by #defining NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" and today tried to get rid of the warning. So, in numpy.i, I followed the warning's advice. I added the # def here: %{ #ifndef SWIG_FILE_WITH_INIT # define NO_IMPORT_ARRAY #endif #include "stdio.h" *#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION* #include <numpy/arrayobject.h> %} SWIG was happy, but when compiling the C++ wrapper, there were many warnings followed by many errors. The warnings were for redefinition of NPY_MIN_BYTE and similar. The errors were for all kinds of stuff, excerpt here: native_wrap.cpp:3632: error: ‘PyArray_NOTYPE’ was not declared in this scope native_wrap.cpp:3633: error: cannot convert ‘PyObject*’ to ‘const PyArrayObject*’ for argument ‘1’ to ‘int PyArray_TYPE(const PyArrayObject*)’ native_wrap.cpp: At global scope: native_wrap.cpp:3877: error: ‘intp’ has not been declared native_wrap.cpp: In function ‘int require_fortran(PyArrayObject*)’: native_wrap.cpp:3929: error: ‘struct tagPyArrayObject’ has no member named ‘nd’ native_wrap.cpp:3933: error: ‘struct tagPyArrayObject’ has no member named ‘flags’ native_wrap.cpp:3933: error: ‘FARRAY’ was not declared in this scope native_wrap.cpp:20411: error: ‘struct tagPyArrayObject’ has no member named ‘data’ It looks like there is a new C API for numpy, and the version of numpy.i that I have doesn't use it. Is there a new version of numpy.i available (or in development) that works with the new API? Short term it will just get rid of a warning but I am interested in a good long term solution in case I need to upgrade numpy. Thanks, Tom

4 8

Fwd: Package: scipy-0.11.0-0.1.rc2.fc18 Tag: f18-updates-candidate Status: failed Built by: orion
by Orion Poplawski 13 Feb '13

13 Feb '13

This is a plea for some help. We've been having trouble getting scipy to pass all of the tests in the Fedora 18 build with python 3.3 (although it seems to build okay in Fedora 19). Below are the logs of the build. There appears to be some kind of memory corruption that manifests itself a little differently on 32-bit vs. 64-bit. I really have no idea myself how to pursue debugging this, though I'm happy to provide any more needed information. Thanks, Orion -------- Original Message -------- Subject: Package: scipy-0.11.0-0.1.rc2.fc18 Tag: f18-updates-candidate Status: failed Built by: orion Date: Thu, 20 Sep 2012 21:01:54 +0000 (UTC) From: Fedora Koji Build System <buildsys(a)fedoraproject.org> To: ausil(a)fedoraproject.org, jspaleta(a)fedoraproject.org, voronov(a)fedoraproject.org, torwangjl(a)fedoraproject.org, alagunambi(a)fedoraproject.org, urkle(a)fedoraproject.org, orion(a)fedoraproject.org Package: scipy-0.11.0-0.1.rc2.fc18 Tag: f18-updates-candidate Status: failed Built by: orion ID: 350761 Started: Thu, 20 Sep 2012 20:39:45 UTC Finished: Thu, 20 Sep 2012 21:01:32 UTC scipy-0.11.0-0.1.rc2.fc18 (350761) failed on buildvm-33.phx2.fedoraproject.org (x86_64), buildvm-35.phx2.fedoraproject.org (i386), buildvm-34.phx2.fedoraproject.org (noarch): BuildError: error building package (arch x86_64), mock exited with status 1; see build.log for more information SRPMS: scipy-0.11.0-0.1.rc2.fc18.src.rpm Failed tasks: ------------- Task 4509076 on buildvm-33.phx2.fedoraproject.org Task Type: buildArch (scipy-0.11.0-0.1.rc2.fc18.src.rpm, x86_64) logs: http://koji.fedoraproject.org/koji/getfile?taskID=4509076&name=build.log http://koji.fedoraproject.org/koji/getfile?taskID=4509076&name=mock_output.… http://koji.fedoraproject.org/koji/getfile?taskID=4509076&name=root.log http://koji.fedoraproject.org/koji/getfile?taskID=4509076&name=state.log Task 4509077 on buildvm-35.phx2.fedoraproject.org Task Type: buildArch (scipy-0.11.0-0.1.rc2.fc18.src.rpm, i686) logs: http://koji.fedoraproject.org/koji/getfile?taskID=4509077&name=build.log http://koji.fedoraproject.org/koji/getfile?taskID=4509077&name=mock_output.… http://koji.fedoraproject.org/koji/getfile?taskID=4509077&name=root.log http://koji.fedoraproject.org/koji/getfile?taskID=4509077&name=state.log Task 4509063 on buildvm-34.phx2.fedoraproject.org Task Type: build (f18-candidate, /scipy:cb69bd06f0d930fbe8840d89b918b617e28af63f) Closed tasks: ------------- Task 4509064 on buildvm-30.phx2.fedoraproject.org Task Type: buildSRPMFromSCM (/scipy:cb69bd06f0d930fbe8840d89b918b617e28af63f) logs: http://koji.fedoraproject.org/koji/getfile?taskID=4509064&name=build.log http://koji.fedoraproject.org/koji/getfile?taskID=4509064&name=checkout.log http://koji.fedoraproject.org/koji/getfile?taskID=4509064&name=mock_output.… http://koji.fedoraproject.org/koji/getfile?taskID=4509064&name=root.log http://koji.fedoraproject.org/koji/getfile?taskID=4509064&name=state.log Task Info: http://koji.fedoraproject.org/koji/taskinfo?taskID=4509063 Build Info: http://koji.fedoraproject.org/koji/buildinfo?buildID=350761

4 8

example reading binary Fortran file
by David Froger 13 Feb '13

13 Feb '13

Hy, My question is about reading Fortran binary file (oh no this question again...) Until now, I was using the unpack module like that : def lread(f,fourBeginning,fourEnd,*tuple): from struct import unpack """Reading a Fortran binary file in litte-endian""" if fourBeginning: f.seek(4,1) for array in tuple: for elt in xrange(array.size): transpose(array).flat[elt] = unpack(array.dtype.char,f.read(array.itemsize))[0] if fourEnd: f.seek(4,1) After googling, I read that fopen and npfille was deprecated, and we should use numpy.fromfile and ndarray.tofile, but despite of the documentaion, the cookbook, the mailling list and google, I don't succed in making a simple example. Considering the simple Fortran code below what is the Python script to read the four arrrays? What about if my pc is little endian and the file big endian? I think it will be a good idea to put the Fortran writting-arrays code and the Python reading-array script in the cookbook and maybe a page to help people comming from Fortran to start with Python ? Best, David Froger program makeArray implicit none integer,parameter:: nx=10,ny=20 real(4),dimension(nx,ny):: ux,uy,p integer :: i,j open(11,file='uxuyp.bin',form='unformatted') do i = 1,nx do j = 1,ny ux(i,j) = real(i*j) uy(i,j) = real(i)/real(j) p (i,j) = real(i) + real(j) enddo enddo write(11) ux,uy write(11) p close(11) end program makeArray

8 22

Do we want scalar casting to behave as it does at the moment?
by Matthew Brett 21 Jan '13

21 Jan '13

Hi, I wanted to check that everyone knows about and is happy with the scalar casting changes from 1.6.0. Specifically, the rules for (array, scalar) casting have changed such that the resulting dtype depends on the _value_ of the scalar. Mark W has documented these changes here: http://docs.scipy.org/doc/numpy/reference/ufuncs.html#casting-rules http://docs.scipy.org/doc/numpy/reference/generated/numpy.result_type.html http://docs.scipy.org/doc/numpy/reference/generated/numpy.promote_types.html Specifically, as of 1.6.0: In [19]: arr = np.array([1.], dtype=np.float32) In [20]: (arr + (2**16-1)).dtype Out[20]: dtype('float32') In [21]: (arr + (2**16)).dtype Out[21]: dtype('float64') In [25]: arr = np.array([1.], dtype=np.int8) In [26]: (arr + 127).dtype Out[26]: dtype('int8') In [27]: (arr + 128).dtype Out[27]: dtype('int16') There's discussion about the changes here: http://mail.scipy.org/pipermail/numpy-discussion/2011-September/058563.html http://mail.scipy.org/pipermail/numpy-discussion/2011-March/055156.html http://mail.scipy.org/pipermail/numpy-discussion/2012-February/060381.html It seems to me that this change is hard to explain, and does what you want only some of the time, making it a false friend. Is it the right behavior for numpy 2.0? Cheers, Matthew

12 87

site.cfg: Custom BLAS / LAPACK configuration
by Bradley M. Froehle 13 Dec '12

13 Dec '12

I recently installed NumPy 1.6.2 on a new computer and wanted to use ACML as the BLAS/LAPACK library. [I'm aware that ACML doesn't provide CBLAS, but that is easy to work around by compiling it yourself to produce libcblas.a or libcblas.so]. I experienced a great bit of difficulty in getting NumPy to use ACML (-lcblas -lacml), primarily stemming from the fact that there was a working ATLAS installation already in /usr/lib64. As far as I can tell, it's IMPOSSIBLE to create a site.cfg which will link to ACML when a system installed ATLAS is present. The detection routine for blas_opt (and similarly for lapack_opt) seem to operate as: * Is MKL present? If so, use it. * Is ATLAS present? If so, use it. * Use [blas] section from site.cfg. Instead I would have expected the detection routine to be more like: * Is [blas_opt] present in site.cfg? If so, use it. * Is MKL present? ... * Is ATLAS present? ... * Use [blas] section from site.cfg. This is not just a problem with ACML. I've also experienced this when using NumPy on some cray supercomputers where the default C compiler automatically links a preferred BLAS/LAPACK. I created a GitHub issue for this: https://github.com/numpy/numpy/issues/2728. In addition, I created a pull request with a "works for me" solution, but which should have needs some wider visibility https://github.com/numpy/numpy/pull/2751. I'd appreciate any reviews, workarounds, or other general feedback. If you want to test our the library detection mechanism you can run the following from within the NumPy source directory:: import __builtin__ as builtins builtins.__NUMPY_SETUP__ = True import numpy.distutils.system_info as si print si.get_info('blas_opt') Thanks, Brad

4 10

Hierarchical vs non-hierarchical ndarray.base and __array_interface__
by Gamblin, Todd 08 Dec '12

08 Dec '12

Hi all, I posted on the change in semantics of ndarray.base here: https://github.com/numpy/numpy/commit/6c0ad59#commitcomment-2153047 And some folks asked me to post my question to the numpy mailing list. I've implemented a tool for mapping processes in parallel applications to nodes in cartesian networks. It uses hierarchies of numpy arrays to represent the domain decomposition of the application, as well as corresponding groups of processes on the network. You can "map" an application to the network using assignment of through views. The tool is here if anyone is curious: https://github.com/tgamblin/rubik. I used numpy to implement this because I wanted to be able to do mappings for arbitrary-dimensional networks. Blue Gene/Q, for example, has a 5-D network. The reason I bring this up is because I rely on the ndarray.base pointer and some of the semantics in __array_interface__ to translate indices within my hierarchy of views. e.g., if a value is at (0,0) in a view I want to know that it's actually at (4,4) in its immediate parent array. After looking over the commit I linked to above, I realized I'm actually relying on a lot of stuff that's not guaranteed by numpy. I rely on .base pointing to its closest parent, and I rely on __array_interface__.data containing the address of the array's memory and its strides. None of these is guaranteed by the API docs: http://docs.scipy.org/doc/numpy/reference/arrays.interface.html So I guess I have a few questions: 1. Is translating indices between base arrays and views something that would be useful to other people? 2. Is there some better way to do this than using ndarray.base and __array_interface__? 3. What's the numpy philosophy on this? Should views know about their parents or not? They obviously have to know a little bit about their memory, but whether or not they know how they were derived from their owning array is a different question. There was some discussion on the vagueness of .base here: http://thread.gmane.org/gmane.comp.python.numeric.general/51688/focus=51703 But it doesn't look like you're deprecating .base in 1.7, only changing its behavior, which I tend to agree is worse than deprecating it. After thinking about all this, I'm not sure what I would like to happen. I can see the value of not keeping extra references around within numpy, and my domain is pretty different from the ways that I imagine people use numpy. I wouldn't have to change my code much to make it work without .base, but I do rely on __array_interface__. If that doesn't include the address and strides, t think I'm screwed as far as translating indices go. Any suggestions? Thanks! -Todd ______________________________________________________________________ Todd Gamblin, tgamblin(a)llnl.gov, http://people.llnl.gov/gamblin2 CASC @ Lawrence Livermore National Laboratory, Livermore, CA, USA

2 1

NumPy-Discussion November 2012