Mailman 3 October 2017 - 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

numpy grant update
by Nathaniel Smith 13 Nov '17

13 Nov '17

Hi all, I wanted to give everyone an update on what's going on with the NumPy grant [1]. As you may have noticed, things have been moving a bit slower than originally hoped -- unfortunately my health is improving but has continued to be rocky [2]. Fortunately, I have awesome co-workers, and BIDS has an institutional interest/mandate for figuring out how to make these things happen, so after thinking it over we've decided to reorganize how we're doing things internally and split up the work to let me focus on the core technical/community aspects without getting overloaded. Specifically, Fernando Pérez and Jonathan Dugan [3] are taking on PI/administration duties, Stéfan van der Walt will focus on handling day-to-day management of the incoming hires, and Nelle Varoquaux & Jarrod Millman will also be joining the team (exact details TBD). This shouldn't really affect any of you, except that you might see some familiar faces with @berkeley.edu emails becoming more engaged. I'm still leading the Berkeley effort, and in any case it's still ultimately the community and NumPy steering council who will be making decisions about the project – this is just some internal details about how we're planning to manage our contributions. But in the interest of full transparency I figured I'd let you know what's happening. In other news, the job ad to start the official hiring process has now been submitted for HR review, so it should hopefully be up soon -- depending on how efficient the bureaucracy is. I'll definitely let everyone know as soon as its posted. I'll also be giving a lunch talk at BIDS tomorrow to let folks locally know about what's going on, which I think will be recorded – I'll send around a link after in case others are interested. -n [1] https://mail.python.org/pipermail/numpy-discussion/2017-May/076818.html [2] https://vorpus.org/blog/emerging-from-the-underworld/ [3] https://bids.berkeley.edu/people/jonathan-dugan -- Nathaniel J. Smith -- https://vorpus.org

6 11

Re: [Numpy-discussion] is __array_ufunc__ ready for prime-time?
by William Sheffler 07 Nov '17

07 Nov '17

Thank you all kindly for your responses! Based on your encouragement, I will pursue an ndarray subclass / __array_ufunc__ implementation. I had been toying with np.set_numeric_ops, which is less than ideal (for example, np.ndarray.around segfaults if I use set_numeric_ops in any way). A second question: very broadly speaking, how much 'pain' can I expect trying to use an np.ndarray subclass in the broader python scientific computing ecosystem, and is there general consensus that projects 'should' support ndarray subclasses? Will > We spent a *long time* sorting out the messy details of __array_ufunc__ [1], especially for handling interactions between different types, e.g., between numpy arrays, non-numpy array-like objects, builtin Python objects, objects that override arithmetic to act in non-numpy-like ways, and of course subclasses of all the above. > We hope that we have it right this time, but as we wrote in the NumPy 1.13 release notes "The API is provisional, we do not yet guarantee backward compatibility as modifications may be made pending feedback." That said, let's give it a try! > If any changes are necessary, I expect it would likely relate to how we handle interactions between different types. That's where we spent the majority of the design effort, but debate is a poor substitute for experience. I would be very surprised if the basic cases (one argument or two arguments of the same type) need any changes. Best, Stephan

11 38

Re: [Numpy-discussion] numpy grant update
by Peter Creasey 31 Oct '17

31 Oct '17

> Date: Fri, 27 Oct 2017 17:52:23 -0400 > From: Marten van Kerkwijk > > Hi Peter, > > When using units, if `a` is not angular (or dimensionless), I don't > see how one could write code in which your example wouldn't fail... > But I may be missing something, since for your example one would just > realize that cos(ka)+i sin(ka) = exp(ika), in which case the log is > just ika and one can the whole complexity... > Hi Marten, Sorry, I thought I replied to you but somehow it didn’t go through. Yes, that example was a bit contrived, but it was just an example where something like sin(x) can be meaningful even if x is dimensional (though you much more typically see these things with log or exp). Best, Peter

1 0

is __array_ufunc__ ready for prime-time?
by William Sheffler 27 Oct '17

27 Oct '17

Right before 1.12, I arranged an API around an np.ndarray subclass, making use of __array_ufunc__ to customize behavior based on structured dtype (we come from c++ and really like operator overloading). Having seen __array_ufunc__ featured in Travis Oliphant's Guide to NumPy: 2nd Edition, I assumed this was the way to go. But it was removed in 1.12. Now that 1.13 has reintroduced __array_ufunc__, can I now rely on its continued availability? I am considering using it as a base-level component in several libraries... is this a dangerous idea? Thanks! Will -- William H. Sheffler Ph.D. Principal Engineer Institute for Protein Design University of Washington

4 4

MATLAB to Numpy
by Andrei Berceanu 27 Oct '17

27 Oct '17

Hi, I am new to Numpy, and would like to start by translating a (badly written?) piece of MATLAB code. What I have come up with so far is this: px = np.zeros_like(tmp_px); py = np.zeros_like(tmp_py); pz = np.zeros_like(tmp_pz) w = np.zeros_like(tmp_w) x = np.zeros_like(tmp_x); y = np.zeros_like(tmp_y); z = np.zeros_like(tmp_z) j=-1 for i in range(tmp_px.size): if tmp_px[i] > 2: j += 1 px[j] = tmp_px[i] py[j] = tmp_py[i] pz[j] = tmp_pz[i] w[j] = tmp_w[i] x[j] = tmp_x[i] y[j] = tmp_y[i] z[j] = tmp_z[i] px=px[:j+1]; py=py[:j+1]; pz=pz[:j+1] w=w[:j+1] x=x[:j+1]; y=y[:j+1]; z=z[:j+1] It works, but I'm sure it's probably the most inefficient way of doing it. What would be a decent rewrite? Thank you so much, Best regards, Andrei

7 11

Re: [Numpy-discussion] numpy grant update
by Peter Creasey 27 Oct '17

27 Oct '17

> Date: Thu, 26 Oct 2017 17:27:33 -0400 > From: Marten van Kerkwijk > > That sounds somewhat puzzling as units cannot really propagate without > them somehow telling how they would change! (e.g., the outcome of > sin(a) is possible only for angular units and then depends on that > unit). But in any case, the mailing list is probably not the best case > to discuss this - rather, I look forward to -- and will most happily > give feedback on -- a NEP or other more detailed explanation! > So whilst it’s true that trigonometric functions only make sense for dimensionless quantities, you might still want to compute them for dimensional quantities for reasons of computational efficiency. Taking your example of sin(a) in a spectral density identity: log(cos(ka) + i sin(ka)) = k log(cos(a) + i sin(a)) so if you are computing the LHS for many k and a single a (i.e k the wavenumber and ka dimensionless) then you might prefer the RHS, which actually uses sin(a). Peter

2 1

Re: [Numpy-discussion] Is there a better way to write a stacked matrix
by Peter Creasey 26 Oct '17

26 Oct '17

> > On Thu, Oct 26, 2017 at 12:11 PM, Daniele Nicolodi <daniele(a)grinta.net> > > wrote: > > > >> is there a better way to write the dot product between a stack of > >> matrices? In my case I need to compute > >> > >> y = A.T @ inv(B) @ A > >> > >> with A a 3x1 matrix and B a 3x3 matrix, N times, with N in the few > >> hundred thousands range. I thus "vectorize" the thing using stack of > >> matrices, so that A is a Nx3x1 matrix and B is Nx3x3 and I can write: > >> > >> y = np.matmul(np.transpose(A, (0, 2, 1)), np.matmul(inv(B), A)) > >> If you only ever multiply your matrix inverse by a single vector then you may also wish to consider np.linalg.solve(B,A) which usually has a better prefactor (although for 3x3 it's pretty marginal, your hardware may vary). Peter

1 0