Mailman 3 March 2016 - 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

Behavior of .reduceat()
by Jaime Fernández del Río May 23, 2016

May 23, 2016

Two of the oldest issues in the tracker (#834 <https://github.com/numpy/numpy/issues/834> and #835 <https://github.com/numpy/numpy/issues/835>) are about how .reduceat() handles its indices parameter. I have been taking a look at the source code, and it would be relatively easy to modify, the hardest part being to figure out what the exact behavior should be. Current behavior is that np.ufunc.reduceat(x, ind) returns [np.ufunc.reduce(a[ind[i]:ind[i+1]] for i in range(len(ind))] with a couple of caveats: 1. if ind[i] >= ind[i+1], then a[ind[i]] is returned, rather than a reduction over an empty slice. 2. an index of len(ind) is appended to the indices argument, to be used as the endpoint of the last slice to reduce over. 3. aside from this last case, the indices are required to be strictly inbounds, 0 <= index < len(x), or an error is raised The proposed new behavior, with some optional behaviors, would be: 1. if ind[i] >= ind[i+1], then a reduction over an empty slice, i.e. the ufunc identity, is returned. This includes raising an error if the ufunc does not have an identity, e.g. np.minimum. 2. to fully support the "reduction over slices" idea, some form of out of bounds indices should be allowed. This could mean either that: 1. only index = len(x) is allowed without raising an error, to allow computing the full reduction anywhere, not just as the last entry of the return, or 2. allow any index in -len(x) <= index <= len(x), with the usual meaning given to negative values, or 3. any index is allowed, with reduction results clipped to existing values (and the usual meaning for negative values). 3. Regarding the appending of that last index of len(ind) to indices, we could: 1. keep appending it, or 2. never append it, since you can now request it without an error being raised, or 3. only append it if the last index is smaller than len(x). My thoughts on the options: - The minimal, more conservative approach would go with 2.1 and 3.1. And of course 1, if we don't implement that none of this makes sense. - I kind of think 2.2 or even 2.3 are a nice enhancement that shouldn't break too much stuff. - 3.2 I'm not sure about, probably hurts more than it helps at this point, although in a brand new design you probably would either not append the last index or also prepend a zero, as in np.split. - And 3.3 seems too magical, probably not a good idea, only listed it for completeness. Any other thoughts or votes on what, if anything, should we implement, and what the deprecation of current behavior should look like? Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

3 2

linux wheels coming soon
by Nathaniel Smith April 22, 2016

April 22, 2016

Hi all, Just a heads-up that we're planning to upload Linux wheels for numpy to PyPI soon. Unless there's some objection, these will be using ATLAS, just like the current Windows wheels, for the same reasons -- moving to something faster like OpenBLAS would be good, but given the concerns about OpenBLAS's reliability we want to get something working first and then worry about making it fast. (Plus it doesn't make sense to ship different BLAS libraries on Windows versus Linux -- that just multiplies our support burden for no reason.) -n -- Nathaniel J. Smith -- https://vorpus.org

13 60

Using OpenBLAS for manylinux wheels
by Matthew Brett April 13, 2016

April 13, 2016

Hi, Olivier Grisel and I are working on building and testing manylinux wheels for numpy and scipy. We first thought that we should use ATLAS BLAS, but Olivier found that my build of these could be very slow [1]. I set up a testing grid [2] which found test errors for numpy and scipy using ATLAS wheels. On the other hand, the same testing grid finds no errors or failures [3] using latest OpenBLAS (0.2.17) and running tests for: numpy scipy scikit-learn numexpr pandas statsmodels This is on the travis-ci ubuntu VMs. Please do test on your own machines with something like this script [4]: source test_manylinux.sh We have worried in the past about the reliability of OpenBLAS, but I find these tests reassuring. Are there any other tests of OpenBLAS that we should run to assure ourselves that it is safe to use? Matthew [1] https://github.com/matthew-brett/manylinux-builds/issues/4#issue-143530908 [2] https://travis-ci.org/matthew-brett/manylinux-testing/builds/118780781 [3] I disabled a few pandas tests which were failing for reasons not related to BLAS. Some of the statsmodels test runs time out. [4] https://gist.github.com/matthew-brett/2fd9d9a29e022c297634

6 15

Re: [Numpy-discussion] Starting work on ufunc rewrite
by Irwin Zaid April 3, 2016

April 3, 2016

Hey guys, I figured I'd just chime in here. Over in DyND-town, we've spent a lot of time figuring out how to structure DyND callables, which are actually more general than NumPy gufuncs. We've just recently got them to a place where we are very happy, and are able to represent a wide range of computations. Our callables use a two-fold approach to evaluation. The first pass is a resolution pass, where a callable can specialize what it is doing based on the input types. It is able to deduce the return type, multidispatch, or even perform some sort of recursive analysis in the form of computations that call themselves. The second pass is construction of a kernel object that is exactly specialized to the metadata (e.g., strides, contiguity, ...) of the array. The callable itself can store arbitrary data, as can each pass of the evaluation. Either (or both) of these passes can be done ahead of time, allowing one to have a callable exactly specialized for your array. If NumPy is looking to change it's ufunc design, we'd be happy to share our experiences with this. Irwin On Thu, Mar 31, 2016 at 4:00 PM, Jaime Fernández del Río <jaime.frio at gmail.com <https://mail.scipy.org/mailman/listinfo/numpy-discussion>> wrote: >* I have started discussing with Nathaniel the implementation of the ufunc ABI *>* break that he proposed in a draft NEP a few months ago: *>>* http://thread.gmane.org/gmane.comp.python.numeric.general/61270 <http://thread.gmane.org/gmane.comp.python.numeric.general/61270> *>>* His original proposal was to make the public portion of PyUFuncObject be: *>>* typedef struct { *>* PyObject_HEAD *>* int nin, nout, nargs; *>* } PyUFuncObject; *>>* Of course the idea is that internally we would use a much larger struct that *>* we could change at will, as long as its first few entries matched those of *>* PyUFuncObject. My problem with this, and I may very well be missing *>* something, is that in PyUFunc_Type we need to set the tp_basicsize to the *>* size of the extended struct, so we would end up having to expose its *>* contents. This is somewhat similar to what now happens with PyArrayObject: *>* anyone can #include "ndarraytypes.h", cast PyArrayObject* to *>* PyArrayObjectFields*, and access the guts of the struct without using the *>* supplied API inline functions. Not the end of the world, but if you want to *>* make something private, you might as well make it truly private. *>>* I think it would be to have something similar to what NpyIter does:: *>>* typedef struct { *>* PyObject_HEAD *>* NpyUFunc *ufunc; *>* } PyUFuncObject; *>>* where NpyUFunc would, at this level, be an opaque type of which nothing *>* would be known. We could have some of the NpyUFunc attributes cached on the *>* PyUFuncObject struct for easier access, as is done in NewNpyArrayIterObject. *>* This would also give us more liberty in making NpyUFunc be whatever we want *>* it to be, including a variable-sized memory chunk that we could use and *>* access at will. NpyIter is again a good example, where rather than storing *>* pointers to strides and dimensions arrays, these are made part of the *>* NpyIter memory chunk, effectively being equivalent to having variable sized *>* arrays as part of the struct. And I think we will probably no longer trigger *>* the Cython warnings about size changes either. *>>* Any thoughts on this approach? Is there anything fundamentally wrong with *>* what I'm proposing here? *>>* Also, this is probably going to end up being a rewrite of a pretty large and *>* complex codebase. I am not sure that working on this on my own and *>* eventually sending a humongous PR is the best approach. Any thoughts on how *>* best to handle turning this into a collaborative, incremental effort? Anyone *>* who would like to join in the fun? *>>* Jaime *>>* -- *>* (\__/) *>* ( O.o) *>* ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de *>* dominación mundial. *>>* _______________________________________________ *>* NumPy-Discussion mailing list *>* NumPy-Discussion at scipy.org <https://mail.scipy.org/mailman/listinfo/numpy-discussion> *>* https://mail.scipy.org/mailman/listinfo/numpy-discussion <https://mail.scipy.org/mailman/listinfo/numpy-discussion> *>

2 1

Starting work on ufunc rewrite
by Jaime Fernández del Río April 3, 2016

April 3, 2016

I have started discussing with Nathaniel the implementation of the ufunc ABI break that he proposed in a draft NEP a few months ago: http://thread.gmane.org/gmane.comp.python.numeric.general/61270 His original proposal was to make the public portion of PyUFuncObject be: typedef struct { PyObject_HEAD int nin, nout, nargs; } PyUFuncObject; Of course the idea is that internally we would use a much larger struct that we could change at will, as long as its first few entries matched those of PyUFuncObject. My problem with this, and I may very well be missing something, is that in PyUFunc_Type we need to set the tp_basicsize to the size of the extended struct, so we would end up having to expose its contents. This is somewhat similar to what now happens with PyArrayObject: anyone can #include "ndarraytypes.h", cast PyArrayObject* to PyArrayObjectFields*, and access the guts of the struct without using the supplied API inline functions. Not the end of the world, but if you want to make something private, you might as well make it truly private. I think it would be to have something similar to what NpyIter does:: typedef struct { PyObject_HEAD NpyUFunc *ufunc; } PyUFuncObject; where NpyUFunc would, at this level, be an opaque type of which nothing would be known. We could have some of the NpyUFunc attributes cached on the PyUFuncObject struct for easier access, as is done in NewNpyArrayIterObject. This would also give us more liberty in making NpyUFunc be whatever we want it to be, including a variable-sized memory chunk that we could use and access at will. NpyIter is again a good example, where rather than storing pointers to strides and dimensions arrays, these are made part of the NpyIter memory chunk, effectively being equivalent to having variable sized arrays as part of the struct. And I think we will probably no longer trigger the Cython warnings about size changes either. Any thoughts on this approach? Is there anything fundamentally wrong with what I'm proposing here? Also, this is probably going to end up being a rewrite of a pretty large and complex codebase. I am not sure that working on this on my own and eventually sending a humongous PR is the best approach. Any thoughts on how best to handle turning this into a collaborative, incremental effort? Anyone who would like to join in the fun? Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

3 3

ANN: bcolz 1.0.0 RC2 is out!
by Francesc Alted March 31, 2016

March 31, 2016

========================== Announcing bcolz 1.0.0 RC2 ========================== What's new ========== Yeah, 1.0.0 is finally here. We are not introducing any exciting new feature (just some optimizations and bug fixes), but bcolz is already 6 years old and it implements most of the capabilities that it was designed for, so I decided to release a 1.0.0 meaning that the format is declared stable and that people can be assured that future bcolz releases will be able to read bcolz 1.0 data files (and probably much earlier ones too) for a long while. Such a format is fully described at: https://github.com/Blosc/bcolz/blob/master/DISK_FORMAT_v1.rst Also, a 1.0.0 release means that bcolz 1.x series will be based on C-Blosc 1.x series (https://github.com/Blosc/c-blosc). After C-Blosc 2.x (https://github.com/Blosc/c-blosc2) would be out, a new bcolz 2.x is expected taking advantage of shiny new features of C-Blosc2 (more compressors, more filters, native variable length support and the concept of super-chunks), which should be very beneficial for next bcolz generation. Important: this is a Release Candidate, so please test it as much as you can. If no issues would appear in a week or so, I will proceed to tag and release 1.0.0 final. Enjoy! For a more detailed change log, see: https://github.com/Blosc/bcolz/blob/master/RELEASE_NOTES.rst What it is ========== *bcolz* provides columnar and compressed data containers that can live either on-disk or in-memory. Column storage allows for efficiently querying tables with a large number of columns. It also allows for cheap addition and removal of column. In addition, bcolz objects are compressed by default for reducing memory/disk I/O needs. The compression process is carried out internally by Blosc, an extremely fast meta-compressor that is optimized for binary data. Lastly, high-performance iterators (like ``iter()``, ``where()``) for querying the objects are provided. bcolz can use numexpr internally so as to accelerate many vector and query operations (although it can use pure NumPy for doing so too). numexpr optimizes the memory usage and use several cores for doing the computations, so it is blazing fast. Moreover, since the carray/ctable containers can be disk-based, and it is possible to use them for seamlessly performing out-of-memory computations. bcolz has minimal dependencies (NumPy), comes with an exhaustive test suite and fully supports both 32-bit and 64-bit platforms. Also, it is typically tested on both UNIX and Windows operating systems. Together, bcolz and the Blosc compressor, are finally fulfilling the promise of accelerating memory I/O, at least for some real scenarios: http://nbviewer.ipython.org/github/Blosc/movielens-bench/blob/master/queryi… Other users of bcolz are Visualfabriq (http://www.visualfabriq.com/) the Blaze project (http://blaze.pydata.org/), Quantopian (https://www.quantopian.com/) and Scikit-Allel (https://github.com/cggh/scikit-allel) which you can read more about by pointing your browser at the links below. * Visualfabriq: * *bquery*, A query and aggregation framework for Bcolz: * https://github.com/visualfabriq/bquery * Blaze: * Notebooks showing Blaze + Pandas + BColz interaction: * http://nbviewer.ipython.org/url/blaze.pydata.org/notebooks/timings-csv.ipynb * http://nbviewer.ipython.org/url/blaze.pydata.org/notebooks/timings-bcolz.ip… * Quantopian: * Using compressed data containers for faster backtesting at scale: * https://quantopian.github.io/talks/NeedForSpeed/slides.html * Scikit-Allel * Provides an alternative backend to work with compressed arrays * https://scikit-allel.readthedocs.org/en/latest/model/bcolz.html Installing ========== bcolz is in the PyPI repository, so installing it is easy:: $ pip install -U bcolz Resources ========= Visit the main bcolz site repository at: http://github.com/Blosc/bcolz Manual: http://bcolz.blosc.org Home of Blosc compressor: http://blosc.org User's mail list: bcolz(a)googlegroups.com http://groups.google.com/group/bcolz License is the new BSD: https://github.com/Blosc/bcolz/blob/master/LICENSES/BCOLZ.txt Release notes can be found in the Git repository: https://github.com/Blosc/bcolz/blob/master/RELEASE_NOTES.rst ---- **Enjoy data!** -- Francesc Alted

1 0

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