NumPy-Discussion June 2015

numpy-discussion@python.org

35 participants
30 discussions

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

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 09 Jan '16

09 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

Video meeting this week
by Nathaniel Smith 12 Jul '15

12 Jul '15

Hi all, In a week and a half, this is happening: https://github.com/numpy/numpy/wiki/SciPy-2015-developer-meeting It's somewhat short notice (my bad :-/), but I think it would be good to have a short video meeting sometime this week as a kind of "pre-meeting" -- to at least briefly go over the main issues we see facing the project to prime the pump, get a better idea about what we want to accomplish at the meeting itself, and gather some early feedback from anyone who won't be able to make it to SciPy (we'll miss you). The obligatory doodle: http://doodle.com/6b4s6thqt9xt4vnh Depending on the interest level, I'm thinking we'll either use Google Hangouts or Bluejeans (https://bluejeans.com/ -- same as what Ralf used for the similar SciPy meeting a few months ago; needs a plugin installed but is available for Windows / OS X / 64-bit Linux / Android / iOS, or regular telephone, or h323 softphone). -n -- Nathaniel J. Smith -- http://vorpus.org

9 31

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

Python 3 and isinstance(np.int64(42), int)
by Jens Jørgen Mortensen 23 Jun '15

23 Jun '15

Hi! I just finished porting a large code-base to Python 3 (making it work on 2.6, 2.7 and 3.4). It wasn't that difficult, but one thing gave me a hard time and it was this: Python 2.7.9 (default, Apr 2 2015, 15:33:21) [GCC 4.9.2] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> a = np.zeros(7, int) >>> n = a[3] >>> type(n) <type 'numpy.int64'> >>> isinstance(n, int) True With Python 3.4 you get False. I think I understand why (np.int64 is no longer a subclass of int). So, I did this instead: import numbers isinstance(n, numbers.Integral) which works fine (with numpy-1.9). Is this the "correct" way or is there a better way to do it? I would imagine that a lot of code will break because of this - so it would be nice if isinstance(n, int) could be made to work the same way in 2 and 3, but I don't know if this is possible (or desirable). Jens Jørgen

6 5

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

Removal of Deprecated Keywords/functionality
by Charles R Harris 22 Jun '15

22 Jun '15

Hi All, There are three long ago deprecations that I am not sure how to handle. - keywords skiprows and missing in genfromtxt, deprecated in 1.5. - keyword old_behavior (default False) in correlate. added in 1.5 at least, but default value changed later. The documentation says they will be removed in numpy 2.0, but we might want to try ealier. The case of the correlation function is trickier, as we probabaly need to provide a function with the old behavior before removing the keyword. I've left these cases as is, but the more old stuff hanging about the greater our technical debt. Chuck

6 6

How to limit cross correlation window width in Numpy?
by Honi Sanders 21 Jun '15

21 Jun '15

I am learning numpy/scipy, coming from a MATLAB background. The xcorr function in Matlab has an optional argument "maxlag" that limits the lag range from –maxlag to maxlag. This is very useful if you are looking at the cross-correlation between two very long time series but are only interested in the correlation within a certain time range. The performance increases are enormous considering that cross-correlation is incredibly expensive to compute. What is troubling me is that numpy.correlate does not have a maxlag feature. This means that even if I only want to see correlations between two time series with lags between -100 and +100 ms, for example, it will still calculate the correlation for every lag between -20000 and +20000 ms (which is the length of the time series). This (theoretically) gives a 200x performance hit! Is it possible that I could contribute this feature? I have introduced this question as a scipy issue https://github.com/scipy/scipy/issues/4940 and on the spicy-dev list (http://mail.scipy.org/pipermail/scipy-dev/2015-June/020757.html). It seems the best place to start is with numpy.correlate, so that is what I am requesting. I have done a simple implementation (https://gist.github.com/bringingheavendown/b4ce18aa007118e4e084) which gives 50x speedup under my conditions (https://github.com/scipy/scipy/issues/4940#issuecomment-110187847). This is my first experience with contributing to open-source software, so any pointers are appreciated.

3 7

Clarification sought on Scipy Numpy version requirements.
by Charles R Harris 21 Jun '15

21 Jun '15

Hi All, I'm looking to change some numpy deprecations into errors as well as remove some deprecated functions. The problem I see is that SciPy claims to support Numpy >= 1.5 and Numpy 1.5 is really, really, old. So the question is, does "support" mean compiles with earlier versions of Numpy ? If that is the case there is very little that can be done about deprecation. OTOH, if it means Scipy can be compiled with more recent numpy versions but used with earlier Numpy versions (which is a good feat), I'd like to know. I'd also like to know what the interface requirements are, as I'd like to remove old_defines.h Chuck

4 14

2022

2021

2020

2019

2018

2017

2016

2015

2014

2013

2012

2011

2010

2009

2008

2007

2006

2005

2004

2003

2002

2001

2000

NumPy-Discussion June 2015