Re: [Numpy-discussion] example reading binary Fortran file
![](https://secure.gravatar.com/avatar/af440d5cd6bfcb9f69493548c1c72546.jpg?s=120&d=mm&r=g)
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
When I run the above program compiled with gfortran on my Intel Mac, I can read it back with::
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.]])
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
![](https://secure.gravatar.com/avatar/6146b81a7f666f75ad4e9a92c87af9eb.jpg?s=120&d=mm&r=g)
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
![](https://secure.gravatar.com/avatar/5c85708f2eed0869671a7d303ca55b85.jpg?s=120&d=mm&r=g)
On Feb 11, 2007, at 22:51 , Satrajit Ghosh wrote:
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
![](https://secure.gravatar.com/avatar/0030ea3910806505d6f0d956ffa0fd83.jpg?s=120&d=mm&r=g)
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 --
![](https://secure.gravatar.com/avatar/02fc64e61fb49e66f33c83dad1f57937.jpg?s=120&d=mm&r=g)
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
![](https://secure.gravatar.com/avatar/b0f62d137f9ea1d0b6cc4e7e6f61b119.jpg?s=120&d=mm&r=g)
On Thu, Jul 19, 2012 at 5:52 AM, Cheng Li <scrappedprince.li@gmail.com>wrote:
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
![](https://secure.gravatar.com/avatar/a3b4ac86a7da4964c3071091ede3d2a7.jpg?s=120&d=mm&r=g)
![](https://secure.gravatar.com/avatar/fe867604a9d80bcc1d9037834a1891b0.jpg?s=120&d=mm&r=g)
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);
![](https://secure.gravatar.com/avatar/b4f6d4f8b501cb05fd054944a166a121.jpg?s=120&d=mm&r=g)
On Thu, 2021-07-22 at 16:48 +0000, Daniel Waddington wrote:
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
![](https://secure.gravatar.com/avatar/6146b81a7f666f75ad4e9a92c87af9eb.jpg?s=120&d=mm&r=g)
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
![](https://secure.gravatar.com/avatar/5c85708f2eed0869671a7d303ca55b85.jpg?s=120&d=mm&r=g)
On Feb 11, 2007, at 22:51 , Satrajit Ghosh wrote:
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
![](https://secure.gravatar.com/avatar/0030ea3910806505d6f0d956ffa0fd83.jpg?s=120&d=mm&r=g)
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 --
![](https://secure.gravatar.com/avatar/02fc64e61fb49e66f33c83dad1f57937.jpg?s=120&d=mm&r=g)
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
![](https://secure.gravatar.com/avatar/b0f62d137f9ea1d0b6cc4e7e6f61b119.jpg?s=120&d=mm&r=g)
On Thu, Jul 19, 2012 at 5:52 AM, Cheng Li <scrappedprince.li@gmail.com>wrote:
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
![](https://secure.gravatar.com/avatar/a3b4ac86a7da4964c3071091ede3d2a7.jpg?s=120&d=mm&r=g)
![](https://secure.gravatar.com/avatar/fe867604a9d80bcc1d9037834a1891b0.jpg?s=120&d=mm&r=g)
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);
![](https://secure.gravatar.com/avatar/b4f6d4f8b501cb05fd054944a166a121.jpg?s=120&d=mm&r=g)
On Thu, 2021-07-22 at 16:48 +0000, Daniel Waddington wrote:
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
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