Mailman 3 June 2017 - 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

Scipy 2017 NumPy sprint
by Charles R Harris July 8, 2017

July 8, 2017

Hi All, I will be running the NumPy sprint at Scipy 2017 and I'm trying to put together a suitable list of things to sprint on. In my experience, sprinting on NumPy is hard, enhancements generally need lengthy review and even finding and doing simple bug fixes can take time. What I have in mind at this point, apart from what might be a getting started tutorial, could mostly be classified as janitorial work. 1. Triage issues and close those that no longer apply. This is mind numbing work, but it has been almost three years since the last pass. 2. Move the contents of `numpy/doc` into `doc/source` and make them normal *.rst files. 3. Convert the doctest in `numpy/lib/tests/test_polynomial.py` to regular tests. Might be tricky as it mostly checks print formatting.

17 31

proposed changes to array printing in 1.14
by Allan Haldane July 1, 2017

July 1, 2017

Hello all, There are various updates to array printing in preparation for numpy 1.14. See https://github.com/numpy/numpy/pull/9139/ Some are quite likely to break other projects' doc-tests which expect a particular str or repr of arrays, so I'd like to warn the list in case anyone has opinions. The current proposed changes, from most to least painful by my reckoning, are: 1. For float arrays, an extra space previously used for the sign position will now be omitted in many cases. Eg, `repr(arange(4.))` will now return 'array([0., 1., 2., 3.])' instead of 'array([ 0., 1., 2., 3.])'. 2. The printing of 0d arrays is overhauled. This is a bit finicky to describe, please see the release note in the PR. As an example of the effect of this, the `repr(np.array(0.))` now prints as 'array(0.)` instead of 'array(0.0)'. Also the repr of 0d datetime arrays is now like "array('2005-04-04', dtype='datetime64[D]')" instead of "array(datetime.date(2005, 4, 4), dtype='datetime64[D]')". 3. User-defined dtypes which did not properly implement their `repr` (and `str`) should do so now. Otherwise it now falls back to `object.__repr__`, which will return something ugly like `<mytype object at 0x7f37f1b4e918>`. (Previously you could depend on only implementing the `item` method and the repr of that would be printed. But no longer, because this risks infinite recursions.). 4. Bool arrays of size 1 with a 'True' value will now omit a space, so that `repr(array([True]))` is now 'array([True])' instead of 'array([ True])'. Allan

9 14

Array blitting (pasting one array into another)
by Mikhail V June 30, 2017

June 30, 2017

Hello all I often need to copy one array into another array, given an offset. This is how the "blit" function can be understood, i.e. in every graphical lib there is such a function. The common definition is like: blit ( dest, src, offset ): where dest is destination array, src is source array and offset is coordinates in destination where the src should pe blitted. Main feature of such function is that it never gives an error, so if the source does not fit into the destination array, it is simply trimmed. And respectively if there is no intersection area then nothing happens. Hope this is clear. So to make it work with Numpy arrays one need to calculate the slices before copying the data. I cannot find any Numpy or Python method to help with that so probably it does not exist yet. If so, my proposal is to add a Numpy method which helps with that. Namely the proposal will be to add a method which returns the slices for the intersection areas of two arbitrary arrays, given an offset, so then one can "blit" the array into another with simple assignment =. Here is a Python function I use for 2d arrays now: def interslice ( dest, src, offset ): y,x = offset H,W = dest.shape h,w = src.shape dest_starty = max (y,0) dest_endy = min (y+h,H) dest_startx = max (x,0) dest_endx = min (x+w,W) src_starty = 0 src_endy = h if y<0: src_starty = -y by = y+h - H # Y bleed if by>0: src_endy = h - by src_startx = 0 src_endx = w if x<0: src_startx = -x bx = x+w - W # X bleed if bx>0: src_endx = w - bx dest_sliceY = slice(dest_starty,dest_endy) dest_sliceX = slice(dest_startx,dest_endx) src_sliceY = slice(src_starty, src_endy) src_sliceX = slice(src_startx, src_endx) if dest_endy <= dest_starty: print "No Y intersection !" dest_sliceY = ( slice(0, 0) ) src_sliceY = ( slice(0, 0) ) if dest_endx <= dest_startx: print "No X intersection !" dest_sliceX = ( slice(0, 0) ) src_sliceX = ( slice(0, 0) ) dest_slice = ( dest_sliceY, dest_sliceX ) src_slice = ( src_sliceY, src_sliceX ) return ( dest_slice, src_slice ) ------ I have intentionally made it expanded and without contractions so that it is better understandable. It returns the intersection area of two arrays given an offset. First returned tuple element is the slice for DEST array and the second element is the slice for SRC array. If there is no intersection along one of the axis at all it returns the corresponding slice as (0,0) With this helper function one can blit arrays easily e.g. example code: W = 8; H = 8 DEST = numpy.ones([H,W], dtype = "uint8") w = 4; h = 1 SRC = numpy.zeros([h,w], dtype = "uint8") SRC[:]=8 offset = (0,9) ds, ss = interslice (DEST, SRC, offset ) # blit SRC into DEST DEST[ds] = SRC[ss] So changing the offset one can observe how the SRC array is trimmed if it is crosses the DEST boundaries. I think it is very useful function in general and it has well defined behaviour. It has usage not only for graphics, but actually any data copying-pasting between arrays. So I am looking forward to comments on this proposal. Mikhail

3 4

Boolean binary '-' operator
by Charles R Harris June 28, 2017

June 28, 2017

Hi All, The boolean binary '-' operator was deprecated back in NumPy 1.9 and changed to an error in 1.13. This caused a number of failures in downstream projects. The choices now are to continue the deprecation for another couple of releases, or simply give up on the change. For booleans, `a - b` was implemented as `a xor b`, which leads to the somewhat unexpected identity `a - b == b - a`, but it is a handy operator that allows simplification of some functions, `numpy.diff` among therm. At this point I'm inclined to give up on the deprecation and retain the old behavior. It is a bit impure but perhaps we can consider it a feature rather than a bug. The unary `-` operator for booleans, now an error, was also deprecated in 1.9 and changed to an error in 1.13. There have been no complaints about that (yet), and it seems like a reasonable thing to do, so I am inclined to leave that error in place. What do others think the correct way forward is? Chuck

13 16

why a[0][0].__mul__(a[0][0]) where a is np.array, gives 'missing 1 required positional argument'?
by Larry Evans June 27, 2017

June 27, 2017

The following code: --{--cut here-- #Purpose: # Demonstrate class with __add__ and __mult__ methods # to demonstrate shortestPath calculations. #OriginalSource: # ./shortestPath.py #References: # [1] http://www.marinamele.com/2014/04/modifying-add-method-of-python-class.html # [2] https://docs.scipy.org/doc/numpy/reference/c-api.types-and-structures.html#… # [3] https://www.abebooks.com/book-search/title/graphs-networks/author/carre-ber… # Table 3.1, path algebra P2. # [4] https://docs.python.org/3/library/abc.html ########################################### import numpy as np from abc import ABCMeta, abstractmethod class PathAbstract(metaclass=ABCMeta): @abstractmethod def __repr__(self): pass @abstractmethod def __add__(self, other): pass @abstractmethod def __radd__(self, other): pass @abstractmethod def __mul__(self, other): pass @abstractmethod def __rmul__(self, other): pass class ShortestNull(PathAbstract): '''Null object in ShortestPath algebra[3]''' def __str__(self): return "<0>" def __repr__(self): return self.__str__() def __addreal__(self, other): return other def __add__(self, other): return other def __radd__(self, other): return other def __mulreal__(self, other): return self def __mul__(self, other): return self def __rmul__(self, other): return self class ShortestUnit(PathAbstract): '''Unit object in ShortestPath algebra[3]''' def __str__(self): return "<1>" def __repr__(self): return self.__str__() def __addreal__(self, other): return self def __add__(self, other): return self def __radd__(self, other): return self def __mulreal__(self, other): return other def __mul__(self, other): return other def __rmul__(self, other): return other class ShortestPath(PathAbstract): def __init__(self, distance): self.distance = distance def __str__(self): return "{0:4.1f}".format(self.distance) def __repr__(self): return self.__str__() def __addreal__(self, other): #paths composed of following shortest of self path or other path. min_dist=min(self.distance,other.distance) return ShortestPath(min_dist) def __add__(self, other): #paths composed of following shortest of self path or other path. if(other == 0): n=ShortestNull() return n.__addreal__(self) else: return other.__addreal__(self) def __radd__(self, other): return self.__add__(other) def __mulreal__(self, other): #path composed of self path followed by other path. min_dist=self.distance+other.distance return ShortestPath(min_dist) def __mul__(self, other): return other.__mulreal__(self) def __rmul__(self, other): return other.__mulreal__(self) d0=ShortestNull() d1=ShortestUnit() d2=ShortestPath(3) d3=ShortestPath(6) print(' d0=',d0,' d1=',d1,' d2=',d2,' d3=',d3) print(' type(d0)=',type(d0).__name__) d4=d0+d2+d3 print('d4=',d4) d5=d2*d3 print('d5=',d5) d6=d0*d2 print('d6=',d6) d7=d1*d2 print('d7=',d7) d8=d2*d0 print('d8=',d8) d9=d2*d1 print('d9=',d9) a=np.array([[ShortestNull,ShortestPath(12)],[ShortestPath(12),ShortestNull()]],dtype=object) print('a=',a,sep='\n') print('a[0]=',a[0],sep='\n') print('a[0][0]=',a[0][0],sep='\n') #a00mul=a[0][0]*a[0][0] a00mul=a[0][0].__mul__(a[0][0]) print('a00mul=',a00mul,sep='\n') #at=a.transpose() #print('at=',at,sep='\n') #aat=a.dot(at) print('aat=',aat,sep='\n') #mp=np.matmul(a,at)#this gives error about Object arrays not supported. #print('mp=',mp) --}--cut here-- When run by make, gives this result: --{--cut here-- make -k python3 shortestPathABC.py d0= <0> d1= <1> d2= 3.0 d3= 6.0 type(d0)= ShortestNull d4= 3.0 d5= 9.0 d6= <0> d7= 3.0 d8= <0> d9= 3.0 a= [[<class '__main__.ShortestNull'> 12.0] [12.0 <0>]] a[0]= [<class '__main__.ShortestNull'> 12.0] a[0][0]= <class '__main__.ShortestNull'> Traceback (most recent call last): File "shortestPathABC.py", line 123, in <module> a00mul=a[0][0].__mul__(a[0][0]) TypeError: __mul__() missing 1 required positional argument: 'other' Makefile:7: recipe for target 'all' failed make: *** [all] Error 1 --}--cut here-- I don't understand why. Apparently, a[0][0] is not a ShortestNull because otherwise, the .__mul__ would have the positional argument, 'other' equal to a[0][0]. What am I missing? TIA. -regards, Larry

3 2

EuroSciPy 2017 call for contributions - extension of deadline
by Pierre de Buyl June 26, 2017

June 26, 2017

(Apologies if you receive multiple copies of this message) 10th European Conference on Python in Science August 28 - September 1, 2017 in Erlangen, Germany The Call for Papers is extended to July 02, 2017 23:00 CEST Description: 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. Erlangen is one of Germany's major science hubs and located north of Munich (90 minutes by train). The Call for Papers is extended to July 02, 2017 23:00 CEST Regards, The EuroSciPy team https://www.euroscipy.org/2017/

1 0

PyRSB: Python interface to librsb sparse matrices library
by Michele Martone June 25, 2017

June 25, 2017

Hi. I'm the author of the high performance multithreaded sparse matrix library `librsb' (mostly C, LGPLv3): http://librsb.sourceforge.net/ I'm *not* a user of SciPy/NumPy/Python, but using Cython I have written a proof-of-concept interface to librsb, named `PyRSB': https://github.com/michelemartone/pyrsb PyRSB is in a prototypal state; e.g. still lacks good error handling. Its interface is trivial, as it mimicks that of SciPy's 'csr_matrix'. Advantages over csr_matrix are in fast multithreaded multiplication of huge sparse matrices. Intended application area is iterative solution of linear systems; particularly fast if with symmetric matrices and many rhs. With this email I am looking for prospective: - users/testers - developers (any interest to collaborate/adopt/include the project?) Looking forward for your feedback, Michele

6 10