Mailman 3 May 2016 - NumPy-Discussion

mixed mode arithmetic
by Neal Becker 11 Jul '23

11 Jul '23

I've been browsing the numpy source. I'm wondering about mixed-mode arithmetic on arrays. I believe the way numpy handles this is that it never does mixed arithmetic, but instead converts arrays to a common type. Arguably, that might be efficient for a mix of say, double and float. Maybe not. But for a mix of complex and a scalar type (say, CDouble * Double), it's clearly suboptimal in efficiency. So, do I understand this correctly? If so, is that something we should improve?

4 6

Invalid value encoutered : how to prevent numpy.where to do this?
by Eric Emsellem 18 Feb '23

18 Feb '23

Dear all, I have a code using lots of "numpy.where" to make some constrained calculations as in: data = arange(10) result = np.where(data == 0, 0., 1./data) # or data1 = arange(10) data2 = arange(10)+1.0 result = np.where(data1 > data2, np.sqrt(data1-data2), np.sqrt(data2-data2)) which then produces warnings like: /usr/bin/ipython:1: RuntimeWarning: invalid value encountered in sqrt or for the first example: /usr/bin/ipython:1: RuntimeWarning: divide by zero encountered in divide How do I avoid these messages to appear? I know that I could in principle use numpy.seterr. However, I do NOT want to remove these warnings for other potential divide/multiply/sqrt etc errors. Only when I am using a "where", to in fact avoid such warnings! Note that the warnings only happen once, but since I am going to release that code, I would like to avoid the user to get such messages which are irrelevant here (because I am testing, with the where, when NOT to divide by zero or take a sqrt of a negative number). thanks! Eric

5 4

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

Re: [Numpy-discussion] axis parameter for count_nonzero
by G Young 17 Jun '16

17 Jun '16

Hi, I have had a PR <https://github.com/numpy/numpy/pull/7177> open (first draft can be found here <https://github.com/numpy/numpy/pull/7138>) for quite some time now that adds an 'axis' parameter to *count_nonzero*. While the functionality is fully in-place, very robust, and actually higher-performing than the original *count_nonzero* function, the obstacle at this point is the implementation, as most of the functionality is now surfaced at the Python level instead of at the C level. I have made several attempts to move the code into C to no avail and have not received much feedback from maintainers unfortunately to move this forward, so I'm opening this up to the mailing list to see what you guys think of the changes and whether or not it should be merged in as is or be tabled until a more C-friendly solution can be found. Thanks!

2 4

Re: [Numpy-discussion] broadcasting for randint
by G Young 17 Jun '16

17 Jun '16

Hi, I have had a PR <https://github.com/numpy/numpy/pull/6938> open for quite some time now that allows arguments to broadcast in *randint*. While the functionality is fully in-place and very robust, the obstacle at this point is the implementation. When the *dtype* parameter was added to *randint* (see here <https://github.com/numpy/numpy/pull/6910>), a big issue with the implementation was that it created so much duplicate code that it would be a huge maintenance nightmare. However, this was dismissed in the original PR message because it was believed that template-ing would be trivial, which seemed reasonable at the time. When I added broadcasting, I introduced a template system to the code that dramatically cut down on the duplication. However, the obstacle has been whether or not this template system is too *ad hoc* to be merged into the library. Implementing a template in Cython was not considered sufficient and is in fact very tricky to do, and unfortunately, I have not received any constructive suggestions from maintainers about how to proceed, so I'm opening this up to the mailing to see whether or not there are better alternatives to what I did, whether this should be merged as it, or whether this should be tabled until a better template can be found. Thanks!

1 3

EuroSciPy 2016
by Pierre de Buyl 07 Jun '16

07 Jun '16

Dear NumPy and SciPy communities, No annoucement was made here for EuroSciPy 2016 I believe. The call for contributions (talks, posters, sprints) is still open for a few days. EuroSciPy 2016 takes place in Erlangen, Germany, from the 23 to the 27 of August and consists of two days of tutorials (beginner and advanced tracks) and two days of conference representing many fields of science, with a focus on Python tools for science. A day of sprints follows (sprints TBA). The keynote speakers are Gaël Varoquaux and Abby Cabunoc Mayes and we can expect a rich tutorial and scientific program! Videos from previous years are available at https://www.youtube.com/playlist?list=PLYx7XA2nY5GeQCCugyvtnHMVLdhYlrRxH and https://www.youtube.com/playlist?list=PLYx7XA2nY5Gcpabmu61kKcToLz0FapmHu Visit us, register and submit an abstract on our website! https://www.euroscipy.org/2016/ SciPythonic regards, The EuroSciPy 2016 team

1 1

ENH: compute many inner products quickly
by Scott Sievert 07 Jun '16

07 Jun '16

I recently ran into an application where I had to compute many inner products quickly (roughy 50k inner products in less than a second). I wanted a vector of inner products over the 50k vectors, or `[x1.T @ A @ x1, …, xn.T @ A @ xn]` with A.shape = (1k, 1k). My first instinct was to look for a NumPy function to quickly compute this, such as np.inner. However, it looks like np.inner has some other behavior and I couldn’t get tensordot/einsum to work for me. Then a labmate pointed out that I can just do some slick matrix multiplication to compute the same quantity, `(X.T * A @ X.T).sum(axis=0)`. I opened [a PR] with this, and proposed that we define a new function called `inner_prods` for this. However, in the PR, @shoyer pointed out > The main challenge is to figure out how to transition the behavior of all these operations, while preserving backwards compatibility. Quite likely, we need to pick new names for these functions, though we should try to pick something that doesn't suggest that they are second class alternatives. Do we choose new function names? Do we add a keyword arg that changes what np.inner returns? [a PR]:https://github.com/numpy/numpy/pull/7690

8 9

Integers to integer powers
by Charles R Harris 04 Jun '16

04 Jun '16

Hi All, There are currently several pull requests apropos integer arrays/scalars to integer powers and, because the area is messy and involves tradeoffs, I'd like to see some discussion here on the list before proceeding. *Scalars in 1.10* In [1]: 1 ** -1 Out[1]: 1.0 In [2]: int16(1) ** -1 Out[2]: 1 In [3]: int32(1) ** -1 Out[3]: 1 In [4]: int64(1) ** -1 Out[4]: 1.0 In [5]: 2 ** -1 Out[5]: 0.5 In [6]: int16(2) ** -1 Out[6]: 0 In [7]: int32(2) ** -1 Out[7]: 0 In [8]: int64(2) ** -1 Out[8]: 0.5 In [9]: 0 ** -1 --------------------------------------------------------------------------- ZeroDivisionError Traceback (most recent call last) <ipython-input-9-fd405d6cf4bc> in <module>() ----> 1 0 ** -1 ZeroDivisionError: 0.0 cannot be raised to a negative power In [10]: int16(0) ** -1 /home/charris/.local/bin/ipython:1: RuntimeWarning: divide by zero encountered in power #!/usr/bin/python /home/charris/.local/bin/ipython:1: RuntimeWarning: invalid value encountered in power #!/usr/bin/python Out[10]: -9223372036854775808 In [11]: int32(0) ** -1 Out[11]: -9223372036854775808 In [12]: int64(0) ** -1 /home/charris/.local/bin/ipython:1: RuntimeWarning: divide by zero encountered in long_scalars #!/usr/bin/python Out[12]: inf Proposed - for non-zero numbers the return type should be float. - for zero numbers a zero division error should be raised. *Scalar Arrays in 1.10* In [1]: array(1, dtype=int16) ** -1 Out[1]: 1 In [2]: array(1, dtype=int32) ** -1 Out[2]: 1 In [3]: array(1, dtype=int64) ** -1 Out[3]: 1 In [4]: array(2, dtype=int16) ** -1 Out[4]: 0 In [5]: array(2, dtype=int32) ** -1 Out[5]: 0 In [6]: array(2, dtype=int64) ** -1 Out[6]: 0 In [7]: array(0, dtype=int16) ** -1 /home/charris/.local/bin/ipython:1: RuntimeWarning: divide by zero encountered in power #!/usr/bin/python /home/charris/.local/bin/ipython:1: RuntimeWarning: invalid value encountered in power #!/usr/bin/python Out[7]: -9223372036854775808 In [8]: array(0, dtype=int32) ** -1 Out[8]: -9223372036854775808 In [9]: array(0, dtype=int64) ** -1 Out[9]: -9223372036854775808 In [10]: type(array(1, dtype=int64) ** -1) Out[10]: numpy.int64 In [11]: type(array(1, dtype=int32) ** -1) Out[11]: numpy.int64 In [12]: type(array(1, dtype=int16) ** -1) Out[12]: numpy.int64 Note that the return type is always int64 in all these cases. However, type is preserved in non-scalar arrays, although the value of int16 is not compatible with int32 and int64 for zero division. In [22]: array([0]*2, dtype=int16) ** -1 Out[22]: array([0, 0], dtype=int16) In [23]: array([0]*2, dtype=int32) ** -1 Out[23]: array([-2147483648, -2147483648], dtype=int32) In [24]: array([0]*2, dtype=int64) ** -1 Out[24]: array([-9223372036854775808, -9223372036854775808]) Proposed: - Raise an ZeroDivisionError for zero division, that is, in the ufunc. - Scalar arrays to return scalar arrays Thoughts? Chuck

11 30

Changing FFT cache to a bounded LRU cache
by Lion Krischer 02 Jun '16

02 Jun '16

Hi all, I was told to take this to the mailing list. Relevant pull request: https://github.com/numpy/numpy/pull/7686 NumPy's FFT implementation caches some form of execution plan for each encountered input data length. This is currently implemented as a simple dictionary which can grow without bounds. Calculating lots of different FFTs thus cause a memory leak from the users' perspective. We encountered a real world situation where this is an issue. The PR linked above proposes to replace the simple dictionary with an LRU (least recently used) cache. It will remove the least recently used pieces of data if it grows beyond a specified size (currently an arbitrary limit of 100 MB per cache). Thus almost all users will still benefit from the caches but their total memory size is now limited. Things to consider: * This adds quite some additional complexity but I could not find a simple way to achieve the same result. * What is a good limit on cache size? I used 100 MB because it works for my uses cases. Cheers! Lion

8 12

Developer Meeting at EuroScipy?
by Sebastian Berg 01 Jun '16

01 Jun '16

Hi all, since we had decided to do a regular developers meeting last year, how would EuroScipy (Aug. 23.-27., Erlangen, Germany) look as a possible place and time to have one? I believe EuroScipy would include a few people who were not able to come to SciPy last year, and it seems SciPy itself already sees some of the devs quite busy anyway. Not having checked back for room availability at the conference or anything, one option may be: August 24. (parallel to the second/last day of Tutorials) Does this seem like a good plan or do you have alternative suggestions or scheduling difficulties with this? Regards, Sebastian

2 1