Mailman 3 May 2018 - 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

Allowing broadcasting of code dimensions in generalized ufuncs
by Marten van Kerkwijk 03 Jul '18

03 Jul '18

Hi All, Following on a PR combining the ability to provide fixed and flexible dimensions [1] (useful for, e.g., 3-vector input with a signature like `(3),(3)->(3)`, and for `matmul`, resp.; based on earlier PRs by Jaime [2] and Matt (Picus) [3]), I've now made a PR with a further enhancement, which allows one can indicate that a core dimension can be broadcast [4]. A particular use case is `all_equal`, a new function suggested in a stalled PR by Matt (Harrigan) [5], which compares two arrays axis-by-axis, but short-circuits if a non-equality is found (unlike what is the case if one does `(a==b).all(axis)`). One thing that would be obviously useful for a routine like `all_equal` is to be able to provide an array as one argument and a constant as another, i.e., if the core dimensions can be broadcast if needed, just like they are in `(a==b).all(axis)`. This is currently not possible: with its signature of `(n),(n)->()`, the two arrays have to have the same trailing size. My PR provides the ability to indicate in the signature that a core dimension can be broadcast, by using a suffix of "|1". Thus, the signature of `all_equal` would become: ``` (n|1),(n|1)->() ``` Comments most welcome (yes, even on the notation - though I think it is fairly self-explanatory)! Marten p.s. There are some similarities to the new "flexible" dimensions implemented for `matmul` [1], but also differences. In particular, for a signature of `(n?),(n?)->()`, one could also pass in an array of trailing size n and a constant, but it would not be possible to pass in an array with trailing size 1: the dimensions with the same name have to be either present and the same or absent. In contrast, for broadcasting, dimensions with the same name can have trailing size n, size 1, or be absent (in which case they count as 1). For broadcasting, any output dimensions with the same name are never affected, while for flexible dimensions those are removed. [1] https://github.com/numpy/numpy/pull/11175 [2] https://github.com/numpy/numpy/pull/5015 [3] https://github.com/numpy/numpy/pull/11132 [4] https://github.com/numpy/numpy/pull/11179 [5] https://github.com/numpy/numpy/pull/8528

10 42

A roadmap for NumPy - longer term planning
by Matti Picus 04 Jun '18

04 Jun '18

At the recent NumPy sprint at BIDS (thanks to those who made the trip) we spent some time brainstorming about a roadmap for NumPy, in the spirit of similar work that was done for Jupyter. The idea is that a document with wide community acceptance can guide the work of the full-time developer(s), and be a source of ideas for expanding development efforts. I put the document up at https://github.com/numpy/numpy/wiki/NumPy-Roadmap, and hope to discuss it at a BOF session during SciPy in the middle of July in Austin. Eventually it could become a NEP or formalized in another way. Matti

12 14

Where to discuss NEPs (was: Re: new NEP: np.AbstractArray and np.asabstractarray)
by Nathaniel Smith 31 May '18

31 May '18

On Thu, Mar 8, 2018 at 7:06 AM, Marten van Kerkwijk <m.h.vankerkwijk(a)gmail.com> wrote: > Hi Nathaniel, > > Overall, hugely in favour! For detailed comments, it would be good to > have a link to a PR; could you put that up? Well, there's a PR here: https://github.com/numpy/numpy/pull/10706 But, this raises a question :-). (One which also came up here: https://github.com/numpy/numpy/pull/10704#issuecomment-371684170) There are sensible two workflows we could use (or at least, two that I can think of): 1. We merge updates to the NEPs as we go, so that whatever's in the repo is the current draft. Anyone can go to the NEP webpage at http://numpy.org/neps (WIP, see #10702) to see the latest version of all NEPs, whether accepted, rejected, or in progress. Discussion happens on the mailing list, and line-by-line feedback can be done by quote-replying and commenting on individual lines. From time to time, the NEP author takes all the accumulated feedback, updates the document, and makes a new post to the list to let people know about the updated version. This is how python-dev handles PEPs. 2. We use Github itself to manage the review. The repo only contains "accepted" NEPs; draft NEPs are represented by open PRs, and rejected NEPs are represented by PRs that were closed-without-merging. Discussion uses Github's commenting/review tools, and happens in the PR itself. This is roughly how Rust handles their RFC process, for example: https://github.com/rust-lang/rfcs Trying to do some hybrid version of these seems like it would be pretty painful, so we should pick one. Given that historically we've tried to use the mailing list for substantive features/planning discussions, and that our NEP process has been much closer to workflow 1 than workflow 2 (e.g., there are already a bunch of old NEPs already in the repo that are effectively rejected/withdrawn), I think we should maybe continue that way, and keep discussions here? So my suggestion is discussion should happen on the list, and NEP updates should be merged promptly, or just self-merged. Sound good? -n -- Nathaniel J. Smith -- https://vorpus.org

7 15

Inconsistent results for the covariance matrix between scipy.optimize.curve_fit and numpy.polyfit
by Jonathan Tammo Siebert 29 May '18

29 May '18

Hi, I hope this is the appropriate place to ask something like this, otherwise please let me know (or feel free to ignore this). Also I hope that I do not misunderstood something or did some silly mistake. If so, please let me know as well! TLDR: When scaling the covariance matrix based on the residuals, scipy.optimize.curve_fit uses a factor of chisq(popt)/(M-N) (with M=number of point, N=number of parameters) and numpy.polyfit uses chisq(popt)/(M-N-2). I am wondering which is correct. I am somewhat confused about different results I am getting for the covariance matrix of a simple linear fit, when comparing `scipy.optimize.curve_fit` and `numpy.polyfit`. I am aware, that `curve_fit` solves the more general non-linear problem numerically, while `polyfit` finds an analytical solution to the linear problem. However, both converge to the same solution, so I suspect that this difference is not important here. The difference, I am curious about is not in the returned parameters but in the estimate of the corresponding covariance matrix. As I understand, there are two different ways to estimate it, based either on the absolute values of the provided uncertainties or by interpreting those only as weights and then scaling the matrix to produce an appropriate reduced chisq. To that end, curve_fit has the parameter: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.cur ve_fit.html: "absolute_sigma : bool, optional If True, sigma is used in an absolute sense and the estimated parameter covariance pcov reflects these absolute values. If False, only the relative magnitudes of the sigma values matter. The returned parameter covariance matrix pcov is based on scaling sigma by a constant factor. This constant is set by demanding that the reduced chisq for the optimal parameters popt when using the scaled sigma equals unity. In other words, sigma is scaled to match the sample variance of the residuals after the fit. Mathematically, pcov(absolute_sigma=False) = pcov(absolute_sigma=True) * chisq(popt)/(M-N)" https://docs.scipy.org/doc/numpy/reference/generated/numpy.polyfit.html on the other hand, does not say anything about how the covariance matrix is estimated. To my understanding, its default should correspond to `absolute_sigma=False` for `curve_fit`. As `polyfit` has a weight parameter instead of an uncertainty parameter, I guess the difference in default behavior is not that surprising. However, even when specifying `absolute_sigma=False`, `curve_fit` and `polyfit` produce different covariance matrices as the applied scaling factors are chisq(popt)/(M-N-2) for `polyfit` (https://github.com/numpy/numpy/blob/6a58e25703cbecb6786faa09a04ae2ec82 21348b/numpy/lib/polynomial.py#L598-L605) and chisq(popt)/(M-N) for `curve_fit` (https://github.com/scipy/scipy/blob/607a21e07dad234f8e63fcf03b7994137a 3ccd5b/scipy/optimize/minpack.py#L781-L782). The argument given in a comment to the scaling `polyfit` is: "Some literature ignores the extra -2.0 factor in the denominator, but it is included here because the covariance of Multivariate Student-T (which is implied by a Bayesian uncertainty analysis) includes it. Plus, it gives a slightly more conservative estimate of uncertainty.", but honestly, in a quick search, I was not able to find any literature not ignoring the extra "factor". But obviously, I could very well be misunderstanding something. Nonetheless, as `curve_fit` ignores it as well, I was wondering whether those two shouldn't give consistent results and if so, which would be the correct solution. Best, Jonathan

2 3

matmul as a ufunc
by Matti Picus 29 May '18

29 May '18

I have made progress with resolving the issue that matmul, the operation which implements `a @ b`, is not a ufunc [2]. Discussion on the issue, which prevents the __array_ufunc__ mechanism for overriding matmul on subclasses of ndarray, yeilded two approaches: - create a wrapper that can convince the ufunc mechanism to call __array_ufunc__ even on functions that are not true ufuncs - expand the semantics of core signatures so that a single matmul ufunc can implement matrix-matrix, vector-matrix, matrix-vector, and vector-vector multiplication. I have put up prototypes of both approaches as pr 11061 [0] and 11133 [1], they are WIP to prove the concept and are a bit rough. Either approach can be made to work. Which is preferable? What are the criteria we should use to judge the relative merits (backward compatibility, efficiency, code clarity, enabling further enhancements, ...) of the approaches? Matti [0] https://github.com/numpy/numpy/pull/11061 [1] https://github.com/numpy/numpy/pull/11133 [2] https://github.com/numpy/numpy/issues/9028

7 9

ANN: EuroScipy 2018
by Valerio Maggio 29 May '18

29 May '18

*** Apologies if you receive multiple copies *** Dear Colleagues, We are delighted to invite you to join us for the **11th European Conference on Python in Science**. The EuroSciPy 2018 (https://www.euroscipy.org/2018/) Conference will be organised by Fondazione Bruno Kessler (FBK) and will take place from August 28 to September 1 in **Trento, Italy**. The EuroSciPy meeting is a cross-disciplinary gathering focused on the use and development of the Python language in scientific research. This event strives to bring together both users and developers of scientific tools, as well as academic research and state of the art industry. The conference will be structured as it follows: Aug, 28-29 : Tutorials and Hands-on Aug, 30-31 : Main Conference Sep, 1 : Sprint TOPICS OF INTEREST: =================== Presentations of scientific tools and libraries using the Python language, including but not limited to: - Algorithms implemented or exposed in Python - Astronomy - Data Visualisation - Deep Learning & AI - Earth, Ocean and Geo Science - General-purpose Python tools that can be of special interest to the scientific community. - Image Processing - Materials Science - Parallel computing - Political and Social Sciences - Project Jupyter - Reports on the use of Python in scientific achievements or ongoing projects. - Robotics & IoT - Scientific data flow and persistence - Scientific visualization - Simulation - Statistics - Vector and array manipulation - Web applications and portals for science and engineering - 3D Printing CALL FOR PROPOSALS: =================== EuroScipy will accept three different kinds of contributions: * Regular Talks: standard talks for oral presentations, allocated in time slots of 15, or 30 minutes, depending on your preference and scheduling constraints. Each time slot considers a Q&A session at the end of the talk (at least, 5 mins). * Hands-on Tutorials: These are beginner or advanced training sessions to dive into the subject with all details. These sessions are 90 minutes long, and the audience will be strongly encouraged to bring a laptop to experiment. * Poster: EuroScipy will host two poster sessions during the two days of Main Conference. So attendees and students are highly encourage to present their work and/or preliminary results as posters. Proposals should be submitted using the EuroScipy submission system at https://pretalx.com/euroscipy18. **Submission deadline is May, 31st 2018**. REGISTRATION & FEES: ==================== To register to EuroScipy 2018, please go to http://euroscipy2018.eventbrite.co.uk or to http://www.euroscipy.org/2018/ Fees: ----- | Tutorials | Student* | Academic/Individual | Industry| |-----------------------------|----------|---------------------|----------| | Early Bird (till July, 1st) | € 50,00 | € 70,00 | € 125,00 | | Regular (till Aug, 5th | € 100,00 | € 110,00 | € 250,00 | | Late (till Aug, 22nd) | € 135,00 | € 135,00 | € 300,00 | | Main Conference | Student* | Academic/Individual | Industry| |-----------------------------|----------|---------------------|----------| | Early Bird (till July, 1st) | € 50,00 | € 70,00 | € 125,00 | | Regular (till Aug, 5th | € 100,00 | € 110,00 | € 250,00 | | Late (till Aug, 22nd) | € 135,00 | € 135,00 | € 300,00 | * A proof of student status will be required at time of the registration. kind regards, Valerio EuroScipy 2018 Organising Committee, Email: info(a)euroscipy.org Website: http://www.euroscipy.org/2018 twitter: @euroscipy

1 0

Adding take_along_axis and put_along_axis functions
by Eric Wieser 28 May '18

28 May '18

These functions provide a vectorized way of using one array to look up items in another. In particular, they extend the 1d: a = np.array([4, 5, 6, 1, 2, 3]) b = np.array(["four", "five", "six", "one", "two", "three"]) i = a.argsort() b_sorted = b[i] To work for higher-dimensions: a = np.array([[4, 1], [5, 2], [6, 3]]) b = np.array([["four", "one"], ["five", "two"], ["six", "three"]]) i = a.argsort(axis=1) b_sorted = np.take_along_axis(b, i, axis=1) put_along_axis is the obvious but less useful dual to this operation, inserting elements rather than extracting them. (Unlike put and take which are not obvious duals). These have been merged in gh-11105 <https://github.com/numpy/numpy/pull/11105>, but as a new addition this probably should have gone by the mailing list first. There was a lack of consensus in gh-8714 <https://github.com/numpy/numpy/pull/8714> about how best to generalize to differing dimensions, so only the non-controversial case where the indices and array have the same dimensions was implemented. These names were chosen to mirror apply_along_axis, which behaves similarly. Do they seem reasonable? Eric

2 1