Mailman 3 September 2015 - 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

Generalize hstack/vstack --> stack; Block matrices like in matlab
by Stefan Otte Jan. 9, 2016

Jan. 9, 2016

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 Dec. 4, 2015

Dec. 4, 2015

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

Cython-based OpenMP-accelerated quartic polynomial solver
by Juha Jeronen Oct. 6, 2015

Oct. 6, 2015

Hi all, I recently developed a Cython-based, OpenMP-accelerated quartic (and cubic, quadratic) polynomial solver to address a personal research need for quickly solving a very large number of independent low-degree polynomial equations for both real and complex coefficients. For example, on an Intel i7-4710MQ, running with 8 threads this code solves approximately 1.2e7 quartic polynomial equations per second. (With OMP_NUM_THREADS=4 the solver gives approximately 8.9e6 solutions per second, so it seems HyperThreading in this case helps about 30%.) The algorithms for cubics and quadratics come from Numerical Recipes (3rd ed.), and the quartic problem is internally reduced to a cubic and two quadratics, using well-known standard tricks. Since to my understanding, for solving polynomial equations NumPy only provides roots(), which works one problem at a time, I'm writing to inquire whether there would be interest to include this functionality to NumPy, if I contribute the code (and clean it up for integration)? I have some previous open source contribution experience. I have authored the IVTC and Phosphor deinterlacers for VLC, and a modular postprocessing filter framework for the Panda3D game engine (posted at the forums on panda3d.org, projected for version 1.10). Currently the polynomial solver is available in a git repository hosted by our university, the relevant files being: https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworod… https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworod… https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworod… https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworod… https://yousource.it.jyu.fi/jjrandom2/miniprojects/blobs/master/misc/tworod… I'm working on a research grant, which allows me to hold the copyright; license is BSD. Polysolve2 is the fast Cython+OpenMP version, while polysolve is an earlier (slower) experimental version using only NumPy array operations. The example module contains a usage example, and setup (in the standard manner) instructs distutils and Cython as to how to build polysolve2 (including setting the relevant math flags and enabling OpenMP). (The setup script includes also some stuff specific to my solver for the Ziegler problem; basically, the "tworods" module can be ignored. I apologize for the inconvenience.) Thoughts? Best regards, -J ------------------------------------------------- Juha Jeronen, Ph.D. juha.jeronen(a)jyu.fi University of Jyväskylä Department of Mathematical Information Technology -------------------------------------------------

12 39

Sign of NaN
by Charles R Harris Sept. 30, 2015

Sept. 30, 2015

Hi All, Due to a recent commit, Numpy master now raises an error when applying the sign function to an object array containing NaN. Other options may be preferable, returning NaN for instance, so I would like to open the topic for discussion on the list. Thoughts? Chuck

12 22

composition of the steering council (was Re: Governance model request)
by Nathaniel Smith Sept. 30, 2015

Sept. 30, 2015

[splitting this out from the thread-o-doom] On Wed, Sep 23, 2015 at 8:47 AM, Chris Barker - NOAA Federal <chris.barker(a)noaa.gov> wrote: [snip] > So here is what I think is on the table: > > We have the steering council. Which leaves two questions: > > -How big should it be? > -Who will be on the original, "seed" council? > > Note that as I understand the current draft of the governance doc, > once established, the council itself decides who is on the council. So > at this point we really are ONLY deciding how it's going to start. It > has to be bootstrapped somehow. > > However, that had been contentious enough that it would probably be a > good idea to hash out some guidelines about the council membership. We actually do have some of those already -- dunno if they're perfect, but they exist :-). To make sure everyone's on the same page, here's a condensed summary of what the draft currently says: Joining the council: contributor must have produced "substantial" and "sustained" contributions over at least one year + be voted on by the current council + be interested and willing to serve. (Then there's some language emphasising that "contributions" should *not* be read narrowly as a euphemism for "lines of code".) Leaving the council: Happens by choice of member, or if inactive for one year and can't be contacted, or if inactive for two years. Former members are listed as "emeritus" to recognize their past service. Rejoining the council: aside from their entry on the list of emeritus members, a former-member and a never-member are treated identically in general, and the rules for re-joining are the same as the rules for joining. Proposal for seed council: "everyone who's merged a pull request since Jan 1, 2014". (Actual text is here: http://thread.gmane.org/gmane.comp.python.numeric.general/61106 , see section "Council membership".) We didn't talk much about these -- I think mostly on the theory that the exact details really aren't going to matter much in the end. These specific rules are exactly the rules that Jupyter/IPython use, stolen without changes. My interpretation is that these rules were designed to produce a council consisting of a broad spectrum of contributors who are actively engaged, in tune and up-to-date with the issues currently facing the project, and broadly respected by the broader community. The rationale for doing things this way (if we keep it) would be that the steering council's "primary responsibility is to facilitate the ordinary community-based decision making procedure", so you need people who are actively engaged in community discussions; and, if things break down then the people best positioned to resolve it are the ones who have the best view of what went wrong, understand the personalities involved, and so forth. In practice, I'm sure any council interventions would involve most members deferring to whoever they judge has the most expertise, whether or not that person is on the council -- it's not like they'll ever be making a decision in a vacuum. Regarding the seed council, I just tried to pick an objective criterion and an arbitrary date that seemed generally in keeping with idea of "should be active in the last 1-to-2-years-ish". Fiddling with the exact date in particular makes very little difference -- between pushing it back to 2 years ago today or forward to 1 year ago today, the only thing that changes is whether Pauli makes the list or not. (And Pauli is obviously a great council candidate, though I don't know whether he even wants to be on it.) > Personally, I have no idea how big the council should be. Too big, and > there is no point, consensus is harder to reach the larger the group, > and the main (only?) role of the council is to resolve issues where > consensus has not been reached in the larger community. But what is > too big? I'm a little wary of the idea of capping the council size. Assuming you're pre-filtering for basic competence and good faith (as we are), then the only way making the council smaller helps with decision making is that it arbitrarily throws away the opinion of some dedicated and valuable contributors. Plus then we have to start making judgements like "well, person A has been around for a while but pretty inactive recently, and person B is doing awesome stuff, should we kick person A off the council to let person B on or...?" Judging whether someone is or isn't a "substantial contributor" is fine, we can do that. Having to make a relative judgement of which of two people is the *more* "substantial contributor", though -- that sounds awful. And given how conflict-adverse groups can be, I suspect that capping the council size would in practice just have the effect that it never takes in new blood. (The old effect where "science advances one retirement at a time".) I'd be interested to hear what Jupyter/IPython's experience with this is, though: they currently have a 10 (!) person steering council, I assume they'd love to have more. Having lots of contributors who are active and engaged enough to meet the steering council qualifications is a good problem to have :-). Technically their situation is slightly different because their council runs on regular voting rather than Apache-style consensus voting (a luxury they can afford because they have a BDFL to step in if regular voting ends up disenfranchising the minority), but I sort of get the impression that in practice they just don't vote unless they know it will be unanimous, and it's worked out for them so far. (To understate the case.) > As for make-up of the council, I think we need to expand beyond people > who have recently contributed core code. > > Yes, the council does need to have expertise to make technical > decisions, but if you think about the likely contentious issues like > ABI breakage, a core-code focused view is incomplete. So there should > be representation by: > > Someone(s) with a long history of working with the code -- that > institutional memory of why decisions were made the way they were > could be key. Sure -- though I can't really imagine any way of framing a rule like this that *wouldn't* be satisfied by Chuck + Ralf + Pauli, so my guess is that such a rule would not actually have any effect on the council membership in practice. > Someone(s) that may not have worked on the core code, but is a major > player in the community, perhaps as the leader of a Numpy-dependent > project. This would provide representation for the broad community. Pointing out features of the current draft again for reference: In the current text, the "NumFOCUS subcommittee" does have an external member to provide some oversight. (So mathematically speaking, this means that the subcommittee is not a subset -- go figure. I blame IPython.) But, this oversight is specifically for financial matters only, not technical ones: "This Subcommittee shall NOT make decisions about the direction, scope or technical direction of the Project." Thomas Caswell (one of the leaders of matplotlib) volunteered to be our external member to start. We certainly could ask him to sit on the steering council as well, but honestly my guess is that this would have no effect, either positive or negative. (I know if someone asked me to serve on a hypothetical matplotlib steering council, then I would just... never do anything, because who am I to second-guess matplotlib's actual developers.) It's not like we don't already hear from downstream projects on a regular basis. And if things have gone *so* pear-shaped that we have a situation where the steering council feels they need to issue a ruling, *and simultaneously* the members of the council are so out-of-touch that they don't know or care about the needs of downstream projects, then the situation is unsalvageable and we should fork and start over. But I mean, it probably wouldn't hurt either. I'm not super-wedded to the current text. I just think we should limit how much effort we spend trying to cover every possible situation ahead of time. If we're clever enough to solve a problem now hypothetically, then we're probably also clever enough to solve it later when it becomes concrete. -n -- Nathaniel J. Smith -- http://vorpus.org

9 23

Cython-based OpenMP-accelerated quartic, polynomial solver
by Matti Picus Sept. 29, 2015

Sept. 29, 2015

2 1