Mailman 3 March 2015 - 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 23 Jul '21

23 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

C-API: multidimensional array indexing?
by Johann Bauer 31 Mar '16

31 Mar '16

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

Generalize hstack/vstack --> stack; Block matrices like in matlab
by Stefan Otte 10 Jan '16

10 Jan '16

Hey, quite often I work with block matrices. Matlab offers the convenient notation [ a b; c d ] to stack matrices. The numpy equivalent is kinda clumsy: vstack([hstack([a,b]), hstack([c,d])]) I wrote the little function `stack` that does exactly that: stack([[a, b], [c, d]]) In my case `stack` replaced `hstack` and `vstack` almost completely. If you're interested in including it in numpy I created a pull request [1]. I'm looking forward to getting some feedback! Best, Stefan [1] https://github.com/numpy/numpy/pull/5057

12 24

F2PY and multi-dimension arrays - unexpected array size error
by Backing Olof 01 Jul '15

01 Jul '15

Hi I am helping out with a Python and Fortran project. Let me give you some background: * Fortran source: C Bergstrom FCC C User subroutine VUMAT subroutine VUMAT( C Read only - * nblock, ndir, nshr, nstatev, nprops, * stepTime, dt, * props, * density, strainInc, * tempOld, * stressOld, stateOld, enerInternOld, enerInelasOld, C Write only - * stressNew, stateNew, enerInternNew, enerInelasNew ) C include 'vaba_param.inc' C dimension props(nprops), 1 density(nblock), strainInc(nblock,ndir+nshr), 2 tempOld(nblock), 5 stressOld(nblock,ndir+nshr), stateOld(nblock,nstatev), 6 enerInternOld(nblock), enerInelasOld(nblock), 2 stressNew(nblock,ndir+nshr), stateNew(nblock,nstatev), 3 enerInternNew(nblock), enerInelasNew(nblock) * Corresponding .pyf integer :: nblock integer :: ndir integer :: nshr integer :: nstatev integer :: nprops real :: steptime real :: dt real dimension(nprops) :: props real dimension(nblock) :: density real dimension(nblock,ndir+nshr) :: straininc real dimension(nblock) :: tempold real dimension(nblock,ndir+nshr) :: stressold real dimension(nblock,nstatev) :: stateold real dimension(nblock) :: enerinternold real dimension(nblock) :: enerinelasold real dimension(nblock,ndir+nshr),intent(out) :: stressnew real dimension(nblock,nstatev),intent(out) :: statenew real dimension(nblock),intent(out) :: enerinternnew real dimension(nblock),intent(out) :: enerinelasnew * Python source with call of Fortran routine: nblock = 1 ndir = 3 nshr = 3 nstatev = 3 nprops = 11 stepTime = 1 dt = 1 props = np.array([10, 0.5, 1e10, 5, 1e12, 3e-6, 8e-6, 27, 2], float) density = np.array([[7.8e3]], float) strainInc = np.array([[1,-0.5,-0.5,0,0,0]], float) tempOld = np.array([[1]], float) stressOld = np.array([[1,1,1,1,1,1]], float) stateOld = np.array([[1,1,1]], float) enerInternOld = np.array([1], float) enerInelasOld = np.array([1], float) stressNew = np.array([[]], float) stateNew = np.array([[]], float) enerInternNew = np.array([[]], float) enerInelasNew = np.array([[]], float) stressNew, stateNew, enerInternNew, enerInelasNew = vumat(nblock, ndir, nshr, nstatev, nprops, stepTime, dt, props, density, strainInc, tempOld, stressOld, stateOld, enerInternOld, enerInelasOld) When trying to run with Python 2.7 I get: olof@ubuntu:~$ ./demo.py unexpected array size: new_size=4, got array with arr_size=1 Traceback (most recent call last): File "./demo.py", line 33, in <module> main() File "./demo.py", line 30, in main stressNew, stateNew, enerInternNew, enerInelasNew = vumat(nblock, ndir, nshr, nstatev, nprops, stepTime, dt, props, density, strainInc, tempOld, stressOld, stateOld, enerInternOld, enerInelasOld) VUMAT_Bergstrom_FCC.error: failed in converting 9th argument `stressold' of VUMAT_Bergstrom_FCC.vumat to C/Fortran array Other stuff: * python 2.7.6 * numpy/f2py 1.8.2 * gcc/gfortran 4.8.2 * ubuntu 14.04 LTS 32-bit I have tried to google, read the f2py manual, fortran tutorials etc, but to no avail. I must also admit that my knowledge in python is so-so and fortran even less(!). What is the missing statement/syntax that I can’t get correct? Your humble programmer, Olof

1 1

PR added: frozen dimensions in gufunc signatures
by Jaime Fernández del Río 23 Jun '15

23 Jun '15

Hi, I have just sent a PR (https://github.com/numpy/numpy/pull/5015), adding the possibility of having frozen dimensions in gufunc signatures. As a proof of concept, I have added a `cross1d` gufunc to `numpy.core.umath_tests`: In [1]: import numpy as np In [2]: from numpy.core.umath_tests import cross1d In [3]: cross1d.signature Out[3]: '(3),(3)->(3)' In [4]: a = np.random.rand(1000, 3) In [5]: b = np.random.rand(1000, 3) In [6]: np.allclose(np.cross(a, b), cross1d(a, b)) Out[6]: True In [7]: %timeit np.cross(a, b) 10000 loops, best of 3: 76.1 us per loop In [8]: %timeit cross1d(a, b) 100000 loops, best of 3: 13.1 us per loop In [9]: c = np.random.rand(1000, 2) In [10]: d = np.random.rand(1000, 2) In [11]: cross1d(c, d) --------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-11-72c66212e40c> in <module>() ----> 1 cross1d(c, d) ValueError: cross1d: Operand 0 has a mismatch in its core dimension 0, with gufunc signature (3),(3)->(3) (size 2 is different from 3) The speed up over `np.cross` is nice, and while `np.cross` is not the best of examples, as it needs to handle more sizes, in many cases this will allow producing gufuncs that work without a Python wrapper redoing checks that are best left to the iterator, such as dimension sizes. It still needs tests, but before embarking on fully developing those, I wanted to make sure that there is an interest on this. I would also like to further enhance gufuncs providing computed dimensions, e.g. making it possible to e.g. define `pairwise_cross` with signature '(n, 3)->($m, 3)', where the $ indicates that m is a computed dimension, that would have to be calculated by a function passed to the gufunc constructor and stored in the gufunc object, based on the other core dimensions. In this case it would make $m be n*(n-1), so that all pairwise cross products between 3D vectors could be computed. The syntax with '$' is kind of crappy, so any suggestions on how to better express this in the signature are more than welcome, as well as any feedback on the merits (or lack of them) of implementing this. Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

4 4

Behavior of np.random.multivariate_normal with bad covariance matrices
by Blake Griffith 08 Apr '15

08 Apr '15

I have an open PR which lets users control the checks on the input covariance matrix. The matrix is required to be symmetric and positve semi-definite (PSD). The current behavior is that NumPy raises a warning if the matrix is not PSD, and does not even check for symmetry. I added a symmetry check, which raises a warning when the input is not symmetric. And added two keyword args which users can use to turn off the checks/warnings when the matrix is ill formed. So this would only cause another new warning to be raised in existing code. This is needed because sometimes the covariance matrix is only *almost* symmetric or PSD due to roundoff error. Thoughts? PR: https://github.com/numpy/numpy/pull/5726

2 2

GSoC students: please read
by Ralf Gommers 05 Apr '15

05 Apr '15

Hi all, It's great to see that this year there are a lot of students interested in doing a GSoC project with Numpy or Scipy. So far five proposals have been submitted, and it looks like several more are being prepared now. I'd like to give you a bit of advice as well as an idea of what's going to happen in the few weeks. The deadline for submitting applications is 27 March. Don't wait until the last day to submit your proposal! It has happened before that Melange was overloaded and unavailable - the Google program admins will not accept that as an excuse and allow you to submit later. So as soon as your proposal is in good shape, put it in. You can still continue revising it. >From 28 March until 13 April we will continue to interact with you, as we request slots from the PSF and rank the proposals. We don't know how many slots we will get this year, but to give you an impression: for the last two years we got 2 slots. Hopefully we can get more this year, but that's far from certain. Our ranking will be based on a combination of factors: the interaction you've had with potential mentors and the community until now (and continue to have), the quality of your submitted PRs, quality and projected impact of your proposal, your enthusiasm, match with potential mentors, etc. We will also organize a video call (Skype / Google Hangout / ...) with each of you during the first half of April to be able to exchange ideas with a higher communication bandwidth medium than email. Finally a note on mentoring: we will be able to mentor all proposals submitted or suggested until now. Due to the large interest and technical nature of a few topics it has in some cases taken a bit long to provide feedback on draft proposals, however there are no showstoppers in this regard. Please continue improving your proposals and working with your potential mentors. Cheers, Ralf

2 4

How to Force Storage Order
by Klemm, Michael 02 Apr '15

02 Apr '15

Dear all, I have found a bug in one of my codes and the way it passes a Numpy matrix to MKL's dgemm routine. Up to now I was assuming that the matrixes are using C order. I guess I have to correct this assumption :-). I have found that the numpy.linalg.svd algorithm creates the resulting U, sigma, and V matrixes with Fortran storage. Is there any way to force these kind of algorithms to not change the storage order? That would make passing the matrixes to the native dgemm operation much easier. Cheers, -michael Dr.-Ing. Michael Klemm Senior Application Engineer Software and Services Group Developer Relations Division Phone +49 89 9914 2340 Cell +49 174 2417583 Intel GmbH Dornacher Strasse 1 85622 Feldkirchen/Muenchen, Deutschland Sitz der Gesellschaft: Feldkirchen bei Muenchen Geschaeftsfuehrer: Christian Lamprechter, Hannes Schwaderer, Douglas Lusk Registergericht: Muenchen HRB 47456 Ust.-IdNr./VAT Registration No.: DE129385895 Citibank Frankfurt a.M. (BLZ 502 109 00) 600119052

4 4