Mailman 3 August 2014 - NumPy-Discussion

mixed mode arithmetic
by Neal Becker July 11, 2023

July 11, 2023

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 Feb. 18, 2023

Feb. 18, 2023

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 July 22, 2021

July 22, 2021

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 March 31, 2016

March 31, 2016

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

PR added: frozen dimensions in gufunc signatures
by Jaime Fernández del Río June 23, 2015

June 23, 2015

Hi, I have just sent a PR (https://github.com/numpy/numpy/pull/5015), adding the possibility of having frozen dimensions in gufunc signatures. As a proof of concept, I have added a `cross1d` gufunc to `numpy.core.umath_tests`: In [1]: import numpy as np In [2]: from numpy.core.umath_tests import cross1d In [3]: cross1d.signature Out[3]: '(3),(3)->(3)' In [4]: a = np.random.rand(1000, 3) In [5]: b = np.random.rand(1000, 3) In [6]: np.allclose(np.cross(a, b), cross1d(a, b)) Out[6]: True In [7]: %timeit np.cross(a, b) 10000 loops, best of 3: 76.1 us per loop In [8]: %timeit cross1d(a, b) 100000 loops, best of 3: 13.1 us per loop In [9]: c = np.random.rand(1000, 2) In [10]: d = np.random.rand(1000, 2) In [11]: cross1d(c, d) --------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-11-72c66212e40c> in <module>() ----> 1 cross1d(c, d) ValueError: cross1d: Operand 0 has a mismatch in its core dimension 0, with gufunc signature (3),(3)->(3) (size 2 is different from 3) The speed up over `np.cross` is nice, and while `np.cross` is not the best of examples, as it needs to handle more sizes, in many cases this will allow producing gufuncs that work without a Python wrapper redoing checks that are best left to the iterator, such as dimension sizes. It still needs tests, but before embarking on fully developing those, I wanted to make sure that there is an interest on this. I would also like to further enhance gufuncs providing computed dimensions, e.g. making it possible to e.g. define `pairwise_cross` with signature '(n, 3)->($m, 3)', where the $ indicates that m is a computed dimension, that would have to be calculated by a function passed to the gufunc constructor and stored in the gufunc object, based on the other core dimensions. In this case it would make $m be n*(n-1), so that all pairwise cross products between 3D vectors could be computed. The syntax with '$' is kind of crappy, so any suggestions on how to better express this in the signature are more than welcome, as well as any feedback on the merits (or lack of them) of implementing this. Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

4 4

New function `count_unique` to generate contingency tables.
by Warren Weckesser Jan. 25, 2015

Jan. 25, 2015

I created a pull request (https://github.com/numpy/numpy/pull/4958) that defines the function `count_unique`. `count_unique` generates a contingency table from a collection of sequences. For example, In [7]: x = [1, 1, 1, 1, 2, 2, 2, 2, 2] In [8]: y = [3, 4, 3, 3, 3, 4, 5, 5, 5] In [9]: (xvals, yvals), counts = count_unique(x, y) In [10]: xvals Out[10]: array([1, 2]) In [11]: yvals Out[11]: array([3, 4, 5]) In [12]: counts Out[12]: array([[3, 1, 0], [1, 1, 3]]) It can be interpreted as a multi-argument generalization of `np.unique(x, return_counts=True)`. It overlaps with Pandas' `crosstab`, but I think this is a pretty fundamental counting operation that fits in numpy. Matlab's `crosstab` (http://www.mathworks.com/help/stats/crosstab.html) and R's `table` perform the same calculation (with a few more bells and whistles). For comparison, here's Pandas' `crosstab` (same `x` and `y` as above): In [28]: import pandas as pd In [29]: xs = pd.Series(x) In [30]: ys = pd.Series(y) In [31]: pd.crosstab(xs, ys) Out[31]: col_0 3 4 5 row_0 1 3 1 0 2 1 1 3 And here is R's `table`: > x <- c(1,1,1,1,2,2,2,2,2) > y <- c(3,4,3,3,3,4,5,5,5) > table(x, y) y x 3 4 5 1 3 1 0 2 1 1 3 Is there any interest in adding this (or some variation of it) to numpy? Warren

5 11

taking a 2D uneven surface slice
by Bryan Woods Dec. 3, 2014

Dec. 3, 2014

I'm trying to do something that at first glance I think should be simple but I can't quite figure out how to do it. The problem is as follows: I have a 3D grid Values[Nx, Ny, Nz] I want to slice Values at a 2D surface in the Z dimension specified by Z_index[Nx, Ny] and return a 2D slice[Nx, Ny]. It is not as simple as Values[:,:,Z_index]. I tried this: >>> values.shape (4, 5, 6) >>> coords.shape (4, 5) >>> slice = values[:,:,coords] >>> slice.shape (4, 5, 4, 5) >>> slice = np.take(values, coords, axis=2) >>> slice.shape (4, 5, 4, 5) >>> Obviously I could create an empty 2D slice and then fill it by using np.ndenumerate to fill it point by point by selecting values[i, j, Z_index[i, j]]. This just seems too inefficient and not very pythonic.

3 3

Weighted Covariance/correlation
by Sebastian Berg Sept. 8, 2014

Sept. 8, 2014

Hi all, Tom Poole has opened pull request https://github.com/numpy/numpy/pull/4960 to implement weights into np.cov (correlation can be added), somewhat picking up the effort started by Noel Dawe in https://github.com/numpy/numpy/pull/3864. The pull request would currently implement an accuracy type `weights` keyword argument as default, but have a switch `repeat_weights` to use repeat type weights instead (frequency type are a special case of this I think). As far as I can see, the code is in a state that it can be tested. But since it is a new feature, the names/defaults are up for discussion, so maybe someone who might use such a feature has a preference. I know we had a short discussion about this before, but it was a while ago. For example another option would be to have the two weights as two keyword arguments, instead of a boolean switch. Regards, Sebastian

3 2

Does a `mergesorted` function make sense?
by Jaime Fernández del Río Sept. 4, 2014

Sept. 4, 2014

A request was open in github to add a `merge` function to numpy that would merge two sorted 1d arrays into a single sorted 1d array. I have been playing around with that idea for a while, and have a branch in my numpy fork that adds a `mergesorted` function to `numpy.lib`: https://github.com/jaimefrio/numpy/commit/ce5d480afecc989a36e5d2bf4ea1d1ba5… I drew inspiration from C++ STL algorithms, and merged into a single function what merge, set_union, set_intersection, set_difference and set_symmetric_difference do there. My first thought when implementing this was to not make it a public function, but use it under the hood to speed-up some of the functions of `arraysetops.py`, which are now merging two already sorted functions by doing `np.sort(np.concatenate((a, b)))`. I would need to revisit my testing, but the speed-ups weren't that great. One other thing I saw value in for some of the `arraysetops.py` functions, but couldn't fully figure out, was in providing extra output aside from the merged arrays, either in the form of indices, or of boolean masks, indicating which items of the original arrays made it into the merged one, and/or where did they end up in it. Since there is at least one other person out there that likes it, is there any more interest in such a function? If yes, any comments on what the proper interface for extra output should be? Although perhaps the best is to leave that out for starters and see what use people make of it, if any. Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

7 29

inplace unary operations?
by josef.pktd＠gmail.com Aug. 31, 2014

Aug. 31, 2014

Is there a way to negate a boolean, or to change the sign of a float inplace ? Josef random thoughts

4 5