NumPy-Discussion July 2015

numpy-discussion@python.org

37 participants
32 discussions

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

Generalize hstack/vstack --> stack; Block matrices like in matlab
by Stefan Otte 09 Jan '16

09 Jan '16

Hey, quite often I work with block matrices. Matlab offers the convenient notation [ a b; c d ] to stack matrices. The numpy equivalent is kinda clumsy: vstack([hstack([a,b]), hstack([c,d])]) I wrote the little function `stack` that does exactly that: stack([[a, b], [c, d]]) In my case `stack` replaced `hstack` and `vstack` almost completely. If you're interested in including it in numpy I created a pull request [1]. I'm looking forward to getting some feedback! Best, Stefan [1] https://github.com/numpy/numpy/pull/5057

12 24

future of f2py and Fortran90+
by Eric Firing 04 Dec '15

04 Dec '15

F2py is a great tool, but my impression is that it is being left behind by the evolution of Fortran from F90 onward. This is unfortunate; it would be nice to be able to easily wrap new Fortran libraries. I'm curious: has anyone been looking into what it would take to enable f2py to handle modern Fortran in general? And into prospects for getting such an effort funded? Eric

5 6

Shared memory check on in-place modification.
by Anton Akhmerov 07 Aug '15

07 Aug '15

Hi everyone, I have encountered an initially rather confusing problem in a piece of code that attempted to symmetrize a matrix: `h += h.T` The problem of course appears due to `h.T` being a view of `h`, and some elements being overwritten during the __iadd__ call. What made it nasty is that it appeared in code that was already covered by tests, but due to the way __iadd__ works internally, the result only became wrong for matrix sizes > 90. Even worse, it was actually `h += h.T.conj()`, but h was real so conj() was also returning a view. I imagine this was already discussed before (also Charles Harris in https://github.com/numpy/numpy/issues/6119 indicated the issue is known), but I would like to still reopen the question since I was unable to find previous iterations of this discussion. First of all does anyone know of a simple way to test for this issue? Specifically I'm concerned about all in-place array modifications that rely on the data being modified. I considered monkey-patching just to figure out where the problem appears, but numpy.ndarray isn't written in Python so I'm not sure it can be monkey-patched. Secondly, is there no way to automatically check for this behavior in numpy itself? I'm not sure how reliable it is, but would adding np.may_share_memory be a reasonably reliable and fast check? This behavior very much c-like, and in my case it was only discovered by accident in production code that's been out for a couple of years. Best, Anton

4 5

Proposal: Deprecate np.int, np.float, etc.?
by Nathaniel Smith 04 Aug '15

04 Aug '15

Hi all, So one of the things exposed in the numpy namespace are objects called np.int np.float np.bool etc. These are commonly used -- in fact, just yesterday on another project I saw a senior person reviewing a pull request instruct a more junior person that they should use np.float instead of float or np.float64. But AFAICT everyone who is actually using them is doing this based on a very easy-to-fall-for misconception, i.e., that these objects have something to do with numpy. In fact they are just aliases for the regular builtin Python types: 'int', 'float', 'bool', etc. NumPy *does have* special numpy-specific types -- but these are called np.int_, np.float_, np.bool_, etc. Apparently they were set up this way back in numpy 0.something, as a backwards compatibility (!) hack: https://github.com/numpy/numpy/pull/6103#issuecomment-123801937 Now, 10+ years later, they continue to confuse people on a regular, ongoing basis, and new users are still being taught misleading "facts" about them. I suggest that we should deprecate them, with no fixed schedule for actually removing them. (I have no idea if/when people will actually stop using them to the point that we can get away with removing them entirely, but in the mean time we should at least be publicizing that any code which is using them is almost certainly based on a misunderstanding.) The technical challenge here is that historically it has simply been impossible to deprecate a global constant like this without using version-specific hacks or accepting unacceptable slowdowns on every attribute access. But, python 3.5 finally adds the necessary machinery to do this in a future-proof way, so now it can be done safely across all versions of Python that we care about, including future unreleased versions: https://github.com/njsmith/metamodule/ Hence: https://github.com/numpy/numpy/pull/6103 Thoughts? -n P.S.: using metamodule.py also gives us the option of making np.testing lazily imported, which last time this came up was benchmarked to improve numpy's import speed by ~35% [1] -- not too bad given that most production code will never touch np.testing. But this is just a teaser postscript; I'm not proposing that we actually do this at this time :-). [1] http://mail.scipy.org/pipermail/numpy-discussion/2012-July/063147.html -- Nathaniel J. Smith -- http://vorpus.org

10 27

ANN: PyWavelets 0.3.0 release
by Ralf Gommers 31 Jul '15

31 Jul '15

Dear all, On behalf of the PyWavelets development team I'm excited to announce the availability of PyWavelets 0.3.0. This is the first release of the package in 3 years. It is the result of a significant effort of a growing development team to modernize the package, to provide Python 3.x support and to make a start with providing new features as well as improved performance. A 0.4.0 release will follow shortly, and will contain more significant new features as well as changes/deprecations to streamline the API. This release requires Python 2.6, 2.7 or 3.3-3.5 and Numpy 1.6.2 or greater. Sources and release notes can be found on https://pypi.python.org/pypi/PyWavelets and https://github.com/PyWavelets/pywt/releases. Activity on the project is picking up quickly. If you're interested in wavelets in Python, you are welcome and invited to join us at https://github.com/PyWavelets/pywt Enjoy, Ralf ============================== PyWavelets 0.3.0 Release Notes ============================== PyWavelets 0.3.0 is the first release of the package in 3 years. It is the result of a significant effort of a growing development team to modernize the package, to provide Python 3.x support and to make a start with providing new features as well as improved performance. A 0.4.0 release will follow shortly, and will contain more significant new features as well as changes/deprecations to streamline the API. This release requires Python 2.6, 2.7 or 3.3-3.5 and NumPy 1.6.2 or greater. Highlights of this release include: - Support for Python 3.x (>=3.3) - Added a test suite (based on nose, coverage up to 61% so far) - Maintenance work: C style complying to the Numpy style guide, improved templating system, more complete docstrings, pep8/pyflakes compliance, and more. New features ============ Test suite ---------- The test suite can be run with ``nosetests pywt`` or with:: >>> import pywt >>> pywt.test() n-D Inverse Discrete Wavelet Transform -------------------------------------- The function ``pywt.idwtn``, which provides n-dimensional inverse DWT, has been added. It complements ``idwt``, ``idwt2`` and ``dwtn``. Thresholding ------------ The function `pywt.threshold` has been added. It unifies the four thresholding functions that are still provided in the ``pywt.thresholding`` namespace. Backwards incompatible changes ============================== None in this release. Other changes ============= Development has moved to `a new repo <https://github.com/PyWavelets/pywt>`_. Everyone with an interest in wavelets is welcome to contribute! Building wheels, building with ``python setup.py develop`` and many other standard ways to build and install PyWavelets are supported now. Authors ======= * Ankit Agrawal + * François Boulogne + * Ralf Gommers + * David Menéndez Hurtado + * Gregory R. Lee + * David McInnis + * Helder Oliveira + * Filip Wasilewski * Kai Wohlfahrt + A total of 9 people contributed to this release. People with a "+" by their names contributed a patch for the first time. This list of names is automatically generated, and may not be fully complete.

1 0

numpy where and dtype in 1.9
by Nathan Jensen 30 Jul '15

30 Jul '15

Hi, The numpy.where() function was rewritten in numpy 1.9 to speed it up. I traced it to this changeset. https://github.com/numpy/numpy/commit/593e3c30c24f0c61a271dc883c614724d7a57… The weird thing is the 1.9 behavior changed the resulting dtype in some situations when using scalar values as the second or third argument. To try and illustrate, I wrote a simple test script and ran it against both numpy 1.7 and 1.9. Here are the results: 2.7.9 (default, Jul 25 2015, 03:06:43) [GCC 4.4.7 20120313 (Red Hat 4.4.7-3)] ***** numpy version 1.7.2 ***** === testing numpy.where with NaNs === numpy.where([True], numpy.float32(1.0), numpy.NaN).dtype float64 numpy.where([True], [numpy.float32(1.0)], numpy.NaN).dtype float32 numpy.where([True], numpy.float32(1.0), [numpy.NaN]).dtype float64 numpy.where([True], [numpy.float32(1.0)], [numpy.NaN]).dtype float64 === testing numpy.where with integers === numpy.where([True], [numpy.float32(1.0)], 65535).dtype float32 numpy.where([True], [numpy.float32(1.0)], 65536).dtype float32 numpy.where([True], [numpy.float32(1.0)], -32768).dtype float32 numpy.where([True], [numpy.float32(1.0)], -32769).dtype float32 2.7.9 (default, Mar 10 2015, 09:26:44) [GCC 4.4.7 20120313 (Red Hat 4.4.7-3)] ***** numpy version 1.9.2 ***** === testing numpy.where with NaNs === numpy.where([True], numpy.float32(1.0), numpy.NaN).dtype float64 numpy.where([True], [numpy.float32(1.0)], numpy.NaN).dtype float32 numpy.where([True], numpy.float32(1.0), [numpy.NaN]).dtype float64 numpy.where([True], [numpy.float32(1.0)], [numpy.NaN]).dtype float64 === testing numpy.where with integers === numpy.where([True], [numpy.float32(1.0)], 65535).dtype float32 numpy.where([True], [numpy.float32(1.0)], 65536).dtype float64 numpy.where([True], [numpy.float32(1.0)], -32768).dtype float32 numpy.where([True], [numpy.float32(1.0)], -32769).dtype float64 Regarding the NaNs with where, the behavior does not differ between 1.7 and 1.9. But it's a little odd that the one scenario returns a dtype of float32 where the other three scenarios return dtype of float64. I'm not sure if that was intentional or a bug? Regarding using ints with where, in 1.7 the resulting dtype is consistent but then in 1.9 the resulting dtype is influenced by the value of the int. It appears it is somehow related to whether the value falls within the range of a short. I'm not sure if this was a side effect of the performance improvement or was intentional? At the very least I think this change in where() should probably be noted in the release notes for 1.9. Our project saw an increase in memory usage with 1.9 due to where(cond, array, scalar) returning arrays of dtype float64 when using scalars not within that limited range. I've attached my simple script if you're interested in running it.

2 2

stats functions with weight support
by Neal Becker 29 Jul '15

29 Jul '15

The histogram function supports a weights option, but most others (e.g., percentile) do not. For my problem, I have a trace of the amounts of time (floating point) that my machine under test is in each of N states. I'd like to produce histograms, kde, maybe nice pics with seaborn. I can use weights option to histogram, but cannot do any fancier ops, because they don't accept a weights= input. For that matter, it would even perhaps be nicer if they accepted a histogram input directly. I looked at seaborn distplot, but it doesn't accept weights and doesn't accept a histogram. I thought maybe I could modify it, but it uses functions like percentile, which also doesn't accept either.

1 0

strange casting rules
by Neal Becker 29 Jul '15

29 Jul '15

np.uint64(-1)+0 Out[36]: 1.8446744073709552e+19 I often work on signal processing requiring bit-exact integral arithmetic. Promoting to float is not helpful - I don't understand the logic of the above example.

2 1

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

2003

2002

2001

2000

NumPy-Discussion July 2015