Mailman 3 June 2021 - NumPy-Discussion

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

Documentation Team meeting - Monday June 8th
by Melissa Mendonça 04 Dec '22

04 Dec '22

Hi all! A reminder that on Monday, June 8, we have another documentation team meeting at 3PM UTC**. If you wish to join on Zoom, you need to use this link https://zoom.us/j/420005230 Here's the permanent hackmd document with the meeting notes: https://hackmd.io/oB_boakvRqKR-_2jRV-Qjg <https://www.google.com/url?q=https%3A%2F%2Fhackmd.io%2FoB_boakvRqKR-_2jRV-Q…> Hope to see you around (especially if you want to introduce yourself or discuss ideas for Google Season of Docs). ** You can click this link to get the correct time at your timezone: https://www.timeanddate.com/worldclock/fixedtime.html?msg=NumPy+Documentati… - Melissa

6 74

copy="never" discussion and no deprecation cycle?
by Sebastian Berg 22 Jul '22

22 Jul '22

Hi all, (sorry for the length, details/discussion below) On the triage call, there seemed a preference to just try to skip the deprecation and introduce `copy="never"`, `copy="if_needed"`, and `copy="always"` (i.e. string options for the `copy` keyword argument). Strictly speaking, this is against the typical policy (one year of warning/errors). But nobody could think of a reasonable chance that anyone actually uses it. (For me just "policy" will be enough of an argument to just take it slow.) BUT: If nobody has *any* concerns at all, I think we may just end up introducing the change right away. The PR is: https://github.com/numpy/numpy/pull/19173 ## The Feature There is the idea to add `copy=never` (or similar). This would modify the existing `copy` argument to make it a 3-way decision: * `copy=always` or `copy=True` to force a copy * `copy=if_needed` or `copy=False` to prefer no-copy behavior * `copy=never` to error when no-copy behavior is not possible (this ensures that a view is returned) this would affect the functions: * np.array(object, copy=...) * arr.astype(new_dtype, copy=...) * np.reshape(arr, new_shape, copy=...), and the method arr.reshape() * np.meshgrid and possibly Where `reshape` currently does not have the option and would benefit by allowing for `arr.reshape(-1, copy=never)`, which would guarantee a view. ## The Options We have three options that are currently being discussed: 1. We introduce a new `np.CopyMode` or `np.<something>.Copy` Enum with values `np.CopyMode.NEVER`, `np.CopyMode.IF_NEEDED`, and `np.CopyMode.ALWAYS` * Plus: No compatibility concerns * Downside(?): This would be a first in NumPy, and is untypical API due to that. 2. We introduce `copy="never"`, `copy="if_needed"` and `copy="always"` as strings (all other strings will be a `TypeError`): * Problem: `copy="never"` currently means `copy=True` (the opposite) Which means new code has to take care when it may run on older NumPy versions. And in theory could make old code return the wrong thing. * Plus: Strings are the typical for options in NumPy currently. 3. Same as 2. But we take it very slow: Make strings an error right now and only introduce the new options after two releases as per typical deprecation policy. ## Discussion We discussed it briefly today in the triage call and we were leaning towards strings. I was honestly expecting to converge to option 3 to avoid compatibility issues (mainly surprises with `copy="never"` on older versions). But considering how weird it is to currently pass `copy="never"`, the question was whether we should not change it with a release note. The probability of someone currently passing exactly one of those three (and no other) strings seems exceedingly small. Personally, I don't have a much of an opinion. But if *nobody* voices any concern about just changing the meaning of the string inputs, I think the current default may be to just do it. Cheers, Sebastian

12 28

[Feature Request] Add alias of np.concatenate as np.concat
by Iordanis Fostiropoulos 11 May '22

11 May '22

In regard to Feature Request: https://github.com/numpy/numpy/issues/16469 It was suggested to sent to the mailing list. I think I can make a strong point as to why the support for this naming convention would make sense. Such as it would follow other frameworks that often work alongside numpy such as tensorflow. For backward compatibility, it can simply be an alias to np.concatenate I often convert portions of code from tf to np, it is as simple as changing the base module from tf to np. e.g. np.expand_dims -> tf.expand_dims. This is done either in debugging (e.g. converting tf to np without eager execution to debug portion of the code), or during prototyping, e.g. develop in numpy and convert in tf. I find myself more than at one occasion to getting syntax errors because of this particular function np.concatenate. It is unnecessarily long. I imagine there are more people that also run into the same problems. Pandas uses concat (torch on the other extreme uses simply cat, which I don't think is as descriptive).

7 6

request to remove the numpy-aarch64 package from PyPI
by Ralf Gommers 14 Mar '22

14 Mar '22

Hi all, FYI, I noticed this package that claimed to be maintained by us: https://pypi.org/project/numpy-aarch64/. That's not ours, so I tried to contact the author (no email provided, but guessed the same username on GitHub) and asked to remove it: https://github.com/tomasriv/DNA_Sequence/issues/1. There are a very large number of packages with "numpy" in the name on PyPI, and there's no way we can audit/police that effectively, but if it's a rebuild that pretends like it's official then I think it's worth doing something about. It could contain malicious code for all we know. Cheers, Ralf

5 7

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

Add smallest_normal and smallest_subnormal attributes to finfo
by Stephannie Jiménez Gacha 20 Jul '21

20 Jul '21

Good afternoon, Given the discussions happened in the Data API consortium when looking into the attributes of `finfo` used in the wild, we found that `tiny` is used regularly but in a good amount of cases not for its intended purpose but rather as "just give me a small number". Following this we are proposing the addition of `smallest_normal` and `smallest_subnormal` attributes. Personally, I think that the `tiny` name is a little bit odd and misleading, so it will be great to leave that as an alias but have a clear name in this class. Right now the PR: https://github.com/numpy/numpy/pull/18536 has all the changes and all the values added were checked against IEEE-754 standard. One of the main concerns is the support of subnormal numbers in certain architectures, where the values can't be calculated accurately. Given the state of the discussion, we don't know if the best alternative is to not add the `smallest_subnormal` attribute and just add the `smallest_number` attribute as an alias to `tiny`. We open this to discussion to see what way we can go in order to get this PR merged. *Stephannie Jimenez Gacha*Software developer *Quansight* | Your Data Experts w: www.quansight.com e: sgacha(a)quansight.com <https://www.linkedin.com/company/quansight> <https://twitter.com/quansightai>

2 1

Proposal to accept NEP 49: Data allocation strategies
by Matti Picus 18 Jul '21

18 Jul '21

Here is the current rendering of the NEP:https://numpy.org/neps/nep-0049.html The mailing list discussion, started on April 20 did not bring up any objections to the proposal, nor were there objections in the discussion around the text of the NEP. There were questions around details of the implementation, thank you reviewers for carefully looking at them and suggesting improvements. If there are no substantive objections within 7 days from this email, then the NEP will be accepted; see NEP 0 for more details. Matti

5 13

is_triangular, is_diagonal, is_symmetric et al. in NumPy or SciPy linalg
by Ilhan Polat 03 Jul '21

03 Jul '21

Dear all, I'm writing some helper Cythpm functions for scipy.linalg which is kinda performant and usable. And there is still quite some wiggle room for more. In many linalg routines there is a lot of performance benefit if the structure can be discovered in a cheap and reliable way at the outset. For example if symmetric then eig can delegate to eigh or if triangular then triangular solvers can be used in linalg.solve and lstsq so forth Here is the Cythonized version for Jupyter notebook to paste to discover the lower/upper bandwidth of square array A that competes well with A != 0 just to use some low level function (note the latter returns an array hence more cost is involved) There is a higher level supervisor function that checks C-contiguousness otherwise specializes to different versions of it Initial cell %load_ext Cython %load_ext line_profiler import cython import line_profiler Then another cell %%cython # cython: language_level=3 # cython: linetrace=True # cython: binding = True # distutils: define_macros=CYTHON_TRACE=1 # distutils: define_macros=CYTHON_TRACE_NOGIL=1 cimport cython cimport numpy as cnp import numpy as np import line_profiler ctypedef fused np_numeric_t: cnp.int8_t cnp.int16_t cnp.int32_t cnp.int64_t cnp.uint8_t cnp.uint16_t cnp.uint32_t cnp.uint64_t cnp.float32_t cnp.float64_t cnp.complex64_t cnp.complex128_t cnp.int_t cnp.long_t cnp.longlong_t cnp.uint_t cnp.ulong_t cnp.ulonglong_t cnp.intp_t cnp.uintp_t cnp.float_t cnp.double_t cnp.longdouble_t @cython.linetrace(True) @cython.initializedcheck(False) @cython.boundscheck(False) @cython.wraparound(False) cpdef inline (int, int) band_check_internal(np_numeric_t[:, ::1]A): cdef Py_ssize_t n = A.shape[0], lower_band = 0, upper_band = 0, r, c cdef np_numeric_t zero = 0 for r in xrange(n): # Only bother if outside the existing band: for c in xrange(r-lower_band): if A[r, c] != zero: lower_band = r - c break for c in xrange(n - 1, r + upper_band, -1): if A[r, c] != zero: upper_band = c - r break return lower_band, upper_band Final cell for use-case --------------- # Make arbitrary lower-banded array n = 50 # array size k = 3 # k'th subdiagonal R = np.zeros([n, n], dtype=np.float32) R[[x for x in range(n)], [x for x in range(n)]] = 1 R[[x for x in range(n-1)], [x for x in range(1,n)]] = 1 R[[x for x in range(1,n)], [x for x in range(n-1)]] = 1 R[[x for x in range(k,n)], [x for x in range(n-k)]] = 2 Some very haphazardly put together metrics %timeit band_check_internal(R) 2.59 µs ± 84.7 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each) %timeit np.linalg.solve(R, zzz) 824 µs ± 6.24 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each) %timeit R != 0. 1.65 µs ± 43.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each) So the worst case cost is negligible in general (note that the given code is slower as it uses the fused type however if I go with tempita standalone version is faster) Two questions: 1) This is missing np.half/float16 functionality since any arithmetic with float16 is might not be reliable including nonzero check. IS it safe to view it as np.uint16 and use that specialization? I'm not sure about the sign bit hence the question. I can leave this out since almost all linalg suite rejects this datatype due to well-known lack of supprt. 2) Should this be in NumPy or SciPy linalg? It is quite relevant to be on SciPy but then again this stuff is purely about array structures. But if the opinion is for NumPy then I would need a volunteer because NumPy codebase flies way above my head. All feedback welcome Best ilhan

3 4

Interface for wrapping random distribution functions (__array_function__, __ufunc__, ?)
by Yoann Mocquin 01 Jul '21

01 Jul '21

Hello there, here is a feature request for the possibility to wrap numpy random distributions to be wrapped by a mechanism like __array_function__ or __array_ufunc__. You can find the GH issue at : https://github.com/numpy/numpy/issues/19382 <https://github.com/numpy/numpy/issues/19382>. This post follows [this one](https://github.com/numpy/numpy/issues/18902 <https://github.com/numpy/numpy/issues/18902>). I figured I should open another issue for `np.random.normal`, but this applies for all distributions I guess. ## Feature Basicaly, I would like that `numpy.random.*` distributions could trigger an interface when passed instances from custom classes, "à la" `__array_ufunc__` or `__array_function__`. Here is an example concept : ```python arr = np.arange(10) # Reference result print(np.mean(arr)) print(np.random.normal(arr)) custom_obj = MyArrayLike(arr) print(np.mean(custom_obj)) # OK : np.mean will trigger __array_function__ interface print(np.random.normal(custom_obj)) # KO : np.random.normal will "only" try to cast the object to float ``` And here is a MWE : ```python import numpy as np np.random.seed(1234) HANDLED_FUNCTIONS = {} class NumericalLabeled(): def __init__(self, value, label=""): self.value = value self.label = label def __repr__(self): return "NumericalLabelled<"+str(self.value) + "," + self.label+">" def __array_function__(self, func, types, args, kwargs): if func not in HANDLED_FUNCTIONS: return NotImplemented return HANDLED_FUNCTIONS[func](*args, **kwargs) def make_numericallabelled(x, label=""): """ Helper function to cast anything into a NumericalLabelled object. """ if isinstance(x, NumericalLabeled): return x else: return NumericalLabeled(x, label=label) # Numpy functions # Override functions - used with __array_function__ def implements(np_function): def decorator(func): HANDLED_FUNCTIONS[np_function] = func return func return decorator @implements(np.random.normal) def np_random_normal(loc=0.0, scale=1.0, **kwargs): # cast both loc and scale into Numericallabelled loc = make_numericallabelled(loc) scale = make_numericallabelled(scale) # check their label is "compatible" if not loc.label == scale.label: raise ValueError return NumericalLabeled(np.random.rand(loc=loc.value, scale=scale.value, **kwargs), loc.label+scale.label) @implements(np.mean) def np_mean(a, *args, **kwargs): return NumericalLabeled(np.mean(a.value, *args, **kwargs), a.label) def main(): # reference result for standard array arr = np.arange(10) print(np.mean(arr)) print(np.random.normal(arr)) # array-like object num_labeled = NumericalLabeled(arr, "toto") print(np.mean(num_labeled)) try: print(np.random.normal(num_labeled)) except Exception as e: print(e) main() ``` which results in ``` 4.5 [ 0.47143516 -0.19097569 3.43270697 2.6873481 3.27941127 5.88716294 6.85958841 6.3634765 8.01569637 6.75731505] NumericalLabelled<4.5,toto> float() argument must be a string or a number, not 'NumericalLabeled' ``` Since the distribution functions accept array as input, I would expect them to be wrappable like many other array functions. ### Versions ``` Python : 3.8.5 (default, Sep 4 2020, 02:22:02) [Clang 10.0.0 ] Numpy : 1.21.0 ``` Cheers

3 2

NumPy-Discussion June 2021