Re: [Numpy-discussion] example reading binary Fortran file
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
Hi, I'm also not quite clear whether the optimized FFTW and UMFPACK libraries are being used or required in numpy at all as show_config() doesn't report it. I see that fftw and umfpack are being used for scipy. I have attached my site.cfg. Any help would be much appreciated. Cheers, Satra ---- snip site.cfg ---- [atlas] libraries = f77blas, cblas, atlas library_dirs = C:\MinGW\lib;C:\DOWNLOADS\BUILDS\lib include_dirs = C:\DOWNLOADS\BUILDS\include\ATLAS [lapack] libraries = flapack, f77blas, cblas, atlas library_dirs = C:\MinGW\lib;C:\DOWNLOADS\BUILDS\lib # UMFPACK # ------- # The UMFPACK library is used to factor large sparse matrices. It, in turn, # depends on the AMD library for reordering the matrices for better performance. # Note that the AMD library has nothing to do with AMD (Advanced Micro Devices), # the CPU company. # # http://www.cise.ufl.edu/research/sparse/umfpack/ # http://www.cise.ufl.edu/research/sparse/amd/ # [amd] library_dirs = C:\MinGW\lib;C:\DOWNLOADS\BUILDS\lib include_dirs = C:\DOWNLOADS\BUILDS\include amd_libs = amd [umfpack] library_dirs = C:\MinGW\lib;C:\DOWNLOADS\BUILDS\lib include_dirs = C:\DOWNLOADS\BUILDS\include\UMFPACK umfpack_libs = umfpack # FFT libraries # ------------- # There are two FFT libraries that we can configure here: FFTW (2 and 3) and djbfft. # # http://fftw.org/ # http://cr.yp.to/djbfft.html # # Given only this section, numpy.distutils will try to figure out which version # of FFTW you are using. [fftw] library_dirs = C:\MinGW\lib;C:\DOWNLOADS\BUILDS\lib include_dirs = C:\DOWNLOADS\BUILDS\include libraries = fftw3 ---- end snip
On Feb 11, 2007, at 22:51 , Satrajit Ghosh wrote:
Hi,
I'm also not quite clear whether the optimized FFTW and UMFPACK libraries are being used or required in numpy at all as show_config() doesn't report it.
I see that fftw and umfpack are being used for scipy.
I have attached my site.cfg. Any help would be much appreciated.
No, they're only in there for scipy (and for other packages that would like to use them). They're not required for Numpy. -- |>|\/|< /------------------------------------------------------------------\ |David M. Cooke http://arbutus.physics.mcmaster.ca/dmc/ |cookedm@physics.mcmaster.ca
Marjolaine, Solution: unique_index = [i for i,x in enumerate(l) if not or x != l[i-1]] Remember that enumerate gives the index,value pairs of the items in any iterable object. Try it for yourself. Here's the output from my IDLE session. In [1]: l = [(1,1), (2,3), (1, 1), (4,5), (2,3), (10,21)] In [2]: l.sort() In [3]: l Out[3]: [(1, 1), (1, 1), (2, 3), (2, 3), (4, 5), (10, 21)] In [4]: unique_index = [i for i, x in enumerate(l) if not i or x != l[i-1]] In [5]: unique_index Out[5]: [0, 2, 4, 5] BTW, I'm posting my response to the numpy-discussion group so others may benefit. It's best to address your questions to the group, as individuals are not always available to answer your question in a timely manner. And by posting your message to the group, you draw from a large body of very knowledgeable people who will gladly help you. Daran --
I saw your message on the numpy discussion board regarding the solution for the efficient removal of duplicates. I have the same problem but need to need to return the indices of the values as an input with associated z values. I was wondering if there was any ways to have the method you propoosed return the indices of the duplicate (or alternatively unique) values in a.
Here is the piece of code that you suggested at the time (Re: [Numpy-discussion] Efficient removal of duplicates, posted on Tue, 16 Dec 2008 01:10:00 -0800)
--------------------------------------------- import numpy as np
a = [(x0,y0), (x1,y1), ...] # A numpy array, but could be a list l = a.tolist() l.sort() unique = [x for i, x in enumerate(l) if not i or x != l[i-1]] # <---- a_unique = np.asarray(unique)
---------------------------------------------
Best regards, marjolaine
-- This message is subject to the CSIR's copyright terms and conditions, e-mail legal notice, and implemented Open Document Format (ODF) standard. The full disclaimer details can be found at http://www.csir.co.za/disclaimer.html.
This message has been scanned for viruses and dangerous content by MailScanner, and is believed to be clean. MailScanner thanks Transtec Computers for their support.
Hi All, I have spot a strange behavior of numpy.fromfunction(). The sample codes are as follows:
import numpy as np
def myOnes(i,j):
return 1.0
a = np.fromfunction(myOnes,(2000,2000))
a
1.0 Actually what I expected is that the function will return a 2000*2000 2d array with unit value. The returned single float value really confused me. Is this a known bug? The numpy version I used is 1.6.1. Regards, Cheng
On Thu, Jul 19, 2012 at 5:52 AM, Cheng Li <scrappedprince.li@gmail.com>wrote:
Hi All,****
** **
I have spot a strange behavior of numpy.fromfunction(). The sample codes are as follows:****
import numpy as np****
def myOnes(i,j):****
return 1.0****
a = np.fromfunction(myOnes,(2000,2000))****
a****
1.0****
** **
Actually what I expected is that the function will return a 2000*2000 2d array with unit value. The returned single float value really confused me. Is this a known bug? The numpy version I used is 1.6.1.****
**
Your function will be called *once*, with arguments that are *arrays* of coordinate values. It must handle these arrays when it computes the values of the array to be created. To see what is happening, print the values of i and j from within your function, e.g.: In [57]: def ijsum(i, j): ....: print "i =", i ....: print "j =", j ....: return i + j ....: In [58]: fromfunction(ijsum, (3, 4)) i = [[ 0. 0. 0. 0.] [ 1. 1. 1. 1.] [ 2. 2. 2. 2.]] j = [[ 0. 1. 2. 3.] [ 0. 1. 2. 3.] [ 0. 1. 2. 3.]] Out[58]: array([[ 0., 1., 2., 3.], [ 1., 2., 3., 4.], [ 2., 3., 4., 5.]]) Your `myOnes` function will work if you modify it something like this: In [59]: def myOnes(i, j): ....: return np.ones(i.shape) ....: In [60]: fromfunction(myOnes, (3, 4)) Out[60]: array([[ 1., 1., 1., 1.], [ 1., 1., 1., 1.], [ 1., 1., 1., 1.]]) The bug is in the docstring for fromfunction. In the description of the `function` argument, it says "`function` must be capable of operating on arrays, and should return a scalar value." But the function should *not* return a scalar value. It should return an array of values appropriate for the given arguments. Warren
**
Regards,****
Cheng****
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
Can you post the complete output from pip? On Jul 14, 2016 2:06 PM, "Brian M Belgodere" <bmbelgod@us.ibm.com> wrote:
While attempting to install numpy on a Power8 (ppc) machine Running RHEL 7.2 I'm encountering the following error below. Numpy 1.11.0 installs and works fine. pip install numpy==1.11.1 ERROR error: Command "gcc -pthread -fno-strict-aliasing -g -O2 -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes -fPIC -Inumpy/core/include -Ibuild/src.linux-ppc64le-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/src/npysort -I/opt/share/Python-2.7.12/ppc64le/include/python2.7 -Ibuild/src.linux-ppc64le-2.7/numpy/core/src/private -Ibuild/src.linux-ppc64le-2.7/numpy/core/src/private -Ibuild/src.linux-ppc64le-2.7/numpy/core/src/private -c build/src.linux-ppc64le-2.7/numpy/core/src/npymath/npy_math_complex.c -o build/temp.linux-ppc64le-2.7/build/src.linux-ppc64le-2.7/numpy/core/src/npymath/npy_math_complex.o" failed with exit status 1
---------------------------------------- Rolling back uninstall of numpy
Command "/opt/share/Python-2.7.12/ppc64le/bin/python -u -c "import setuptools, tokenize;__file__='/tmp/pip-build-LvkRpf/numpy/setup.py';exec(compile(getattr(tokenize, 'open', open)(__file__).read().replace('\r\n', '\n'), __file__, 'exec'))" install --record /tmp/pip-Bg_8Ko-record/install-record.txt --single-version-externally-managed --compile" failed with error code 1 in /tmp/pip-build-LvkRpf/numpy/ ENVIRONMENT DETAILS OS: Red Hat Enterprise Linux 7.2 Maipo Kernel: ppc64le Linux 3.10.0-327.el7.ppc64le CPU: IBM PowerPC G3 POWER8 (raw) @ 160x 3.491GHz $pip list Cython (0.24) numpy (1.11.0) pip (8.1.2) scipy (0.17.1) setuptools (24.0.3) virtualenv (15.0.2) $gcc -v Using built-in specs. COLLECT_GCC=gcc COLLECT_LTO_WRAPPER=/usr/libexec/gcc/ppc64le-redhat-linux/4.8.5/lto-wrapper Target: ppc64le-redhat-linux Configured with: ../configure --prefix=/usr --mandir=/usr/share/man --infodir=/usr/share/info --with-bugurl= http://bugzilla.redhat.com/bugzilla --enable-bootstrap --enable-shared --enable-threads=posix --enable-checking=release --with-system-zlib --enable-__cxa_atexit --disable-libunwind-exceptions --enable-gnu-unique-object --enable-linker-build-id --with-linker-hash-style=gnu --enable-languages=c,c++,objc,obj-c++,java,fortran,ada,go,lto --enable-plugin --enable-initfini-array --disable-libgcj --with-isl=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-ppc64le-redhat-linux/isl-install --with-cloog=/builddir/build/BUILD/gcc-4.8.5-20150702/obj-ppc64le-redhat-linux/cloog-install --enable-gnu-indirect-function --enable-secureplt --with-long-double-128 --enable-targets=powerpcle-linux --disable-multilib --with-cpu-64=power7 --with-tune-64=power8 --build=ppc64le-redhat-linux Thread model: posix gcc version 4.8.5 20150623 (Red Hat 4.8.5-4) (GCC)
*Brian Michael Belgodere*
Systems Management Engineering and Optimization Software Engineer IBM Research Phone: +1-724-510-5308 e-mail: bmbelgod@us.ibm.com
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org https://mail.scipy.org/mailman/listinfo/numpy-discussion
Hi, I'm working with Numpy in the context of supporting different memory types such as persistent memory and CXL attached. I would like to propose a minor change, but figured I would get some initial feedback from the developer community before submitting a PR. In multiarray/alloc.c the allocator (beneath the cache) using the POSIX malloc/calloc/realloc/free. I propose that these should be changed to PyMem_RawXXX equivalents. The reason for this is that by doing so, one can use the python custom allocator functions (e.g. PyMem_GetAllocator/PyMem_SetAllocator) to intercept the memory allocator for NumPy arrays. This will be useful as heterogeneous memories need supporting. I don't think this will drastically change performance but it is an extra function redirection (and it will only impact when the cache can't deliver). There are likely other places in NumPy that could do with a rinse and repeat - may be someone could advise? Thanks, Daniel Waddington IBM Research --- Example patch for 1.19.x (I'm building with Python3.6) diff --git a/numpy/core/src/multiarray/alloc.c b/numpy/core/src/multiarray/alloc.c index 795fc7315..e9e888478 100644 --- a/numpy/core/src/multiarray/alloc.c +++ b/numpy/core/src/multiarray/alloc.c @@ -248,7 +248,7 @@ PyDataMem_NEW(size_t size) void *result; assert(size != 0); - result = malloc(size); + result = PyMem_RawMalloc(size); if (_PyDataMem_eventhook != NULL) { NPY_ALLOW_C_API_DEF NPY_ALLOW_C_API @@ -270,7 +270,7 @@ PyDataMem_NEW_ZEROED(size_t size, size_t elsize) { void *result; - result = calloc(size, elsize); + result = PyMem_RawCalloc(size, elsize); if (_PyDataMem_eventhook != NULL) { NPY_ALLOW_C_API_DEF NPY_ALLOW_C_API @@ -291,7 +291,7 @@ NPY_NO_EXPORT void PyDataMem_FREE(void *ptr) { PyTraceMalloc_Untrack(NPY_TRACE_DOMAIN, (npy_uintp)ptr); - free(ptr); + PyMem_RawFree(ptr);
On Thu, 2021-07-22 at 16:48 +0000, Daniel Waddington wrote:
Hi, I'm working with Numpy in the context of supporting different memory types such as persistent memory and CXL attached. I would like to propose a minor change, but figured I would get some initial feedback from the developer community before submitting a PR.
Hi Daniel, you may want to have a look at Matti's NEP to allow custom allocation strategies: https://numpy.org/neps/nep-0049.html When implemented, this will allow to explicitly modify the behaviour here (which means you could make it use the Python version). In principle, once that work is done, we could also use the Python allocator as you are proposing. It may be a follow-up discussion. The difficulty is that the NumPy ABI is fully open: 1. A user can create an array with data they allocated 2. In theory, a user could `realloc` or even replace an arrays `data` In practice, hopefully nobody does the second one, but we can't be sure. The first means we have to wait for the NEP, because it will allow us to work around the problem: We can use different `free`/`realloc` if a user provided the data. The second means that we have to be careful when consider changing the default even after implementing the NEP. But it may be possible, at least if we do it slowly/gently. Cheers, Sebastian
In multiarray/alloc.c the allocator (beneath the cache) using the POSIX malloc/calloc/realloc/free. I propose that these should be changed to PyMem_RawXXX equivalents. The reason for this is that by doing so, one can use the python custom allocator functions (e.g. PyMem_GetAllocator/PyMem_SetAllocator) to intercept the memory allocator for NumPy arrays. This will be useful as heterogeneous memories need supporting. I don't think this will drastically change performance but it is an extra function redirection (and it will only impact when the cache can't deliver). There are likely other places in NumPy that could do with a rinse and repeat - may be someone could advise? Thanks, Daniel Waddington IBM Research --- Example patch for 1.19.x (I'm building with Python3.6) diff --git a/numpy/core/src/multiarray/alloc.c b/numpy/core/src/multiarray/alloc.c index 795fc7315..e9e888478 100644 --- a/numpy/core/src/multiarray/alloc.c +++ b/numpy/core/src/multiarray/alloc.c @@ -248,7 +248,7 @@ PyDataMem_NEW(size_t size) void *result; assert(size != 0); - result = malloc(size); + result = PyMem_RawMalloc(size); if (_PyDataMem_eventhook != NULL) { NPY_ALLOW_C_API_DEF NPY_ALLOW_C_API @@ -270,7 +270,7 @@ PyDataMem_NEW_ZEROED(size_t size, size_t elsize) { void *result; - result = calloc(size, elsize); + result = PyMem_RawCalloc(size, elsize); if (_PyDataMem_eventhook != NULL) { NPY_ALLOW_C_API_DEF NPY_ALLOW_C_API @@ -291,7 +291,7 @@ NPY_NO_EXPORT void PyDataMem_FREE(void *ptr) { PyTraceMalloc_Untrack(NPY_TRACE_DOMAIN, (npy_uintp)ptr); - free(ptr); + PyMem_RawFree(ptr);
_______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@python.org https://mail.python.org/mailman/listinfo/numpy-discussion
participants (10)
-
Brian M Belgodere -
Cheng Li -
Daniel Waddington -
Daran L. Rife -
David M. Cooke -
Nathaniel Smith -
Neil Martinsen-Burrell -
Satrajit Ghosh -
Sebastian Berg -
Warren Weckesser