Mailman 3 February 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

PyPy nditer context managers
by xoviat 27 Feb '18

27 Feb '18

The Pull Request related nditer has been dormant for a while, and there seems to be no consensus on moving forward at this point. Charles, are you opposed to requiring the .close() call on nditer and/or adding context managers? If so, then we might as well let mattip know that the Pull Request is not moving forward. If you aren’t opposed, perhaps we should clarify what the expectations for the Pull Request are and a path to having it merged.

1 0

improving arange()? introducing fma()?
by Benjamin Root 27 Feb '18

27 Feb '18

Note, the following is partly to document my explorations in computing lat/on grids in numpy lately, in case others come across these issues in the future. It is partly a plea for some development of numerically accurate functions for computing lat/lon grids from a combination of inputs: bounds, counts, and resolutions. I have been playing around with the decimal package a bit lately, and I discovered the concept of "fused multiply-add" operations for improved accuracy. I have come to realize that fma operations could be used to greatly improve the accuracy of linspace() and arange(). In particular, I have been needing improved results for computing latitude/longitude grids, which tend to be done in float32's to save memory (at least, this is true in data I come across). Since fma is not available yet in python, please consider the following C snippet: ``` $ cat test_fma.c #include<stdlib.h> #include<math.h> #include<stdio.h> int main(){ float res = 0.01; float x0 = -115.0; int cnt = 7001; float x1 = res * cnt + x0; float x2 = fma(res, cnt, x0); printf("x1 %.16f x2 %.16f\n", x1, x2); float err1 = fabs(x1 - -44.990000000000000); float err2 = fabs(x2 - -44.990000000000000); printf("err1 %.5e err2 %.5e\n", err1, err2); return 0; } $ gcc test_fma.c -lm -o test_fma $ ./test_fma x1 -44.9899978637695312 x2 -44.9900016784667969 err1 2.13623e-06 err2 1.67847e-06 ``` And if you do similar for double-precision, the fma still yields significant accuracy over double precision explicit multiply-add. Now, to the crux of my problem. It is next to impossible to generate a non-trivial numpy array of coordinates, even in double precision, without hitting significant numerical errors. Which has lead me down the path of using the decimal package (which doesn't play very nicely with numpy because of the lack of casting rules for it). Consider the following: ``` $ cat test_fma.py from __future__ import print_function import numpy as np res = np.float32(0.01) cnt = 7001 x0 = np.float32(-115.0) x1 = res * cnt + x0 print("res * cnt + x0 = %.16f" % x1) x = np.arange(-115.0, -44.99 + (res / 2), 0.01, dtype='float32') print("len(arange()): %d arange()[-1]: %16f" % (len(x), x[-1])) x = np.linspace(-115.0, -44.99, cnt, dtype='float32') print("linspace()[-1]: %.16f" % x[-1]) $ python test_fma.py res * cnt + x0 = -44.9900015648454428 len(arange()): 7002 arange()[-1]: -44.975044 linspace()[-1]: -44.9900016784667969 ``` arange just produces silly results (puts out an extra element... adding half of the resolution is typically mentioned as a solution on mailing lists to get around arange()'s limitations -- I personally don't do this). linspace() is a tad bit better than arange(), but still a little bit worse than just computing the value straight out. It also produces a result that is on par with fma(), so that is reassuring that it is about as accurate as can be at the moment. This also matches up almost exactly to what I would get if I used Decimal()s and then casted to float32's for final storage, so that's not too bad as well. So, does it make any sense to improve arange by utilizing fma() under the hood? Also, any plans for making fma() available as a ufunc? Notice that most of my examples required knowing the number of grid points ahead of time. But what if I didn't know that? What if I just have the bounds and the resolution? Then arange() is the natural fit, but as I showed, its accuracy is lacking, and you have to do some sort of hack to do a closed interval. Thoughts? Comments? Ben Root (donning flame suit)

9 16

Re: [Numpy-discussion] Prototype CorePRNG implementation feedback (Neal Becker)
by Kevin Sheppard 23 Feb '18

23 Feb '18

> > > What is the syntax to construct an initialized generator? > RandomGenerator(Xoroshiro128(my_seed))? > > Not 100% certain on this. There was talk in the earlier thread that seed should be killed, although I wasn't clear what mathod would be preferrable to get easily reproducible random numbers. Another option would be RandomGenerator(Xoroshiro128, seed=my_seed) Kevin

2 1

1.15 release manager
by Charles R Harris 22 Feb '18

22 Feb '18

Hi All, Would any of the other NumPy developers be interested in making the 1.15 release? It isn't difficult, and I don't mind doing it, but it would be a good thing if more people became familiar with the process. Chuck

1 0

Prototype CorePRNG implementation feedback
by Kevin Sheppard 22 Feb '18

22 Feb '18

I have implemented a working prototype of a pluggable RandomState. The repo is located at https://github.com/bashtage/core-prng. The main design has a small PRNG (currently only SplitMix64 and Xoroshiro128 are implemented) that does not have any public functions available to generate random numbers. These objects only expose methods to get and set state and to jump the PRNG. The random number generator is exposed via an opaque structure wrapped in a PyCapsule. This structure has two void pointers, one to the state and one to the function that generates the next uint 64. I think in the long-run it might make sense to expose a few additional interfaces, such as the next random double (makes sense for dSFMT) or the next random uint32 (makes sense for MT19937). For now, I have deliberately kept it simple. The user-facing object is called `RandomGenerator` and it can be initialized using the syntax `RandomGenerator` or `RandomGenerator(Xoroshiro128())`. The object exposes user-facing methods (currently only `random_integer` and `random_double`). A basic demo: from core_prng.generator import RandomGenerator # Splitmix 64 for now RandomGenerator().random_integer() from core_prng.xoroshiro128 import Xoroshiro128 RandomGenerator(Xoroshiro128()).random_integer() A few questions that have come up: 1. Should a passed core PRNG be initialized or not. The syntax RandomGenerator(Xoroshiro128) looks nicer than RandomGenerator(Xoroshiro128()) 2. Which low-level methods should be part of the core PRNG? These would all be scalar generators that would be consumed by RandomGenerator and exposed through void *. I could imagine some subset of the following: - uint64 - uint32 - double - float - uint53 (applicable to PRNGs that use doubles as default type) 3. Since RandomState is taken, a new name is needed for the user-facing object. RandomGenerator is a placeholder but ideally, something meaningful could be chosen. 4. I can't see how locking can be usefully implemented in the core PRNGs without a large penalty. This means that these are not threadsafe. This isn't a massive problem but they do access the state (get or set) although with GIL. 5. Should state be a property -- this isn't Java... Other feedback is of course also welcome. Kevin

2 1

NumPy 1.14.1 released
by Charles R Harris 21 Feb '18

21 Feb '18

Hi All, On behalf of the NumPy team, I am pleased to announce NumPy 1.14.1. This is a bugfix release for some problems reported following the 1.14.0 release. The major problems fixed are the following. - Problems with the new array printing, particularly the printing of complex values, Please report any additional problems that may turn up. - Problems with ``np.einsum`` due to the new ``optimized=True`` default. Some fixes for optimization have been applied and ``optimize=False`` is now the default. - The sort order in ``np.unique`` when ``axis=<some-number>`` will now always be lexicographic in the subarray elements. In previous NumPy versions there was an optimization that could result in sorting the subarrays as unsigned byte strings. - The change in 1.14.0 that multi-field indexing of structured arrays returns a view instead of a copy has been reverted but remains on track for NumPy 1.15. Affected users should read the 1.14.1 Numpy User Guide section "basics/structured arrays/accessing multiple fields" for advice on how to manage this transition. This release supports Python 2.7 and 3.4 - 3.6. Wheels for the release are available on PyPI. Source tarballs, zipfiles, release notes, and the changelog are available on github <https://github.com/numpy/numpy/releases/tag/v1.14.1>. *Contributors* A total of 14 people contributed to this release. People with a "+" by their names contributed a patch for the first time. * Allan Haldane * Charles Harris * Daniel Smith * Dennis Weyland + * Eric Larson * Eric Wieser * Jarrod Millman * Kenichi Maehashi + * Marten van Kerkwijk * Mathieu Lamarre * Sebastian Berg * Simon Conseil * Simon Gibbons * xoviat Cheers, Charles Harris

2 1

Using the C-API iterator with object arrays
by Eugen Wintersberger 17 Feb '18

17 Feb '18

Hi there, I have a question concerning the numpy iterator C-API. I want to create a numpy string array using NPY_OBJECT as a datatype for creating the array (the reason I am going for this approach is that I do not know the length of the individual strings at the time I construct the array, I only know its shape). The original strings are stored in a std::vector<char*> instance. The approach I took was something like this std::vector<char*> buffer = ....; NpyIter *iter = NpyIter_New(numpy_array, NPY_ITER_READWRITE | NPY_ITER_C_INDEX | NPY_ITER_REFS_OK, NPY_CORDER , NPY_NO_CASTING,nullptr); if(iter==NULL) { return; } NpyIter_IterNextFunc *iternext = NpyIter_GetIterNext(iter,nullptr); if(iternext == NULL) { std::cerr<<"Could not instantiate next iterator function"<<std::endl; return; } PyObject **dataptr = (PyObject**)NpyIter_GetDataPtrArray(iter); for(auto string: buffer) { dataptr[0] = PyString_FromSting(string); // this string construction seem to work iternext(iter); } NpyIter_Deallocate(iter); This code snippet is a bit stripped down with all the safety checks removed to make things more readable. However, the array I get back still contains only a None instance. Does anyone have an idea what I am doing wrong here? Thanks in advance. best regards Eugen

3 3