Mailman 3 August 2011 - NumPy-Discussion

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

SVD errors
by mtrumpis＠berkeley.edu 04 Jun '13

04 Jun '13

Hello list.. I've run into two SVD errors over the last few days. Both errors are identical in numpy/scipy. I've submitted a ticket for the 1st problem (numpy ticket #990). Summary is: some builds of the lapack_lite module linking against system LAPACK (not the bundled dlapack_lite.o, etc) give a "LinAlgError: SVD did not converge" exception on my matrix. This error does occur using Mac's Accelerate framework LAPACK, and a coworker's Ubuntu LAPACK version. It does not seem to happen using ATLAS LAPACK (nor using Octave/Matlab on said Ubuntu) Just today I've come across a negative singular value cropping up in an SVD of a different matrix. This error does occur on my ATLAS LAPACK based numpy, as well as on the Ubuntu setup. And once again, it does not happen in Octave/Matlab. I'm using numpy 1.3.0.dev6336 -- don't know what the Ubuntu box is running. Here are some npy files for the two different cases: https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/noconverge_operator.npy https://cirl.berkeley.edu/twiki/pub/User/MikeTrumpis/negsval_operator.npy Mike

6 7

example reading binary Fortran file
by David Froger 13 Feb '13

13 Feb '13

Hy, My question is about reading Fortran binary file (oh no this question again...) Until now, I was using the unpack module like that : def lread(f,fourBeginning,fourEnd,*tuple): from struct import unpack """Reading a Fortran binary file in litte-endian""" if fourBeginning: f.seek(4,1) for array in tuple: for elt in xrange(array.size): transpose(array).flat[elt] = unpack(array.dtype.char,f.read(array.itemsize))[0] if fourEnd: f.seek(4,1) After googling, I read that fopen and npfille was deprecated, and we should use numpy.fromfile and ndarray.tofile, but despite of the documentaion, the cookbook, the mailling list and google, I don't succed in making a simple example. Considering the simple Fortran code below what is the Python script to read the four arrrays? What about if my pc is little endian and the file big endian? I think it will be a good idea to put the Fortran writting-arrays code and the Python reading-array script in the cookbook and maybe a page to help people comming from Fortran to start with Python ? Best, David Froger 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

8 22

dtype '|S0' not understood
by David Warde-Farley 07 Oct '12

07 Oct '12

Howdy, It seems it's possible using e.g. In [25]: dtype([('foo', str)])Out[25]: dtype([('foo', '|S0')]) to get yourself a zero-length string. However dtype('|S0') results in a TypeError: data type not understood. I understand the stupidity of creating a 0-length string field but it's conceivable that it's accidental. For example, it could lead to a situation where you've created that field, are missing all the data you had meant to put in it, serialize with np.save, and upon np.load aren't able to get _any_ of your data back because the dtype descriptor is considered bogus (can you guess why I thought of this scenario?). It seems that either dtype(str) should do something more sensible than zero-length string, or it should be possible to create it with dtype('| S0'). Which should it be? David

3 3

Conditional random variables
by Ted To 29 Oct '11

29 Oct '11

Hi, Is there an easy way to make random draws from a conditional random variable? E.g., draw a random variable, x conditional on x>=\bar x. Thank you, Ted To

3 6

Weird upcast behavior with 1.6.x, working as intended?
by Olivier Delalleau 30 Sep '11

30 Sep '11

Hi, This is with numpy 1.6.1 under Linux x86_64, testing the upcast mechanism of "scalar + array": >>> import numpy; print (numpy.array(3, dtype=numpy.complex128) + numpy.ones(3, dtype=numpy.float32)).dtype complex64 Since it has to upcast my array (float32 is not "compatible enough" with complex128), why does it upcast it to complex64 instead of complex128? As far as I can tell 1.4.x and 1.5.x versions of numpy are indeed upcasting to complex128. Thanks, -=- Olivier

3 5

datetimes with date vs time units, local time, and time zones
by Mark Wiebe 30 Sep '11

30 Sep '11

Towards a reasonable behavior with regard to local times, I've made the default repr for datetimes use the C standard library to print them in a local ISO format. Combined with the ISO8601-prescribed behavior of interpreting datetime strings with no timezone specifier to be in local times, this allows the following cases to behave reasonably: >>> np.datetime64('now') numpy.datetime64('2011-06-15T15:16:51-0500','s') >>> np.datetime64('2011-06-15T18:00') numpy.datetime64('2011-06-15T18:00-0500','m') As noted in another thread, there can be some extremely surprising behavior as a consequence: >>> np.array(['now', '2011-06-15'], dtype='M') array(['2011-06-15T15:18:26-0500', '2011-06-14T19:00:00-0500'], dtype='datetime64[s]') Having the 15th of June print out as 7pm on the 14th of June is probably not what one would generally expect, so I've come up with an approach which hopefully deals with this in a good way. One firm principal of the datetime in NumPy is that it is always stored as a POSIX time (referencing UTC), or a TAI time. There are two categories of units that can be used, which I will call *date units* and *time units*. The date units are 'Y', 'M', 'W', and 'D', while the time units are 'h', 'm', 's', ..., 'as'. Time zones are only applied to datetimes stored in time units, so there's a qualitative difference between date and time units with respect to string conversions and calendar operations. I would like to place an 'unsafe' casting barrier between the date units and the time units, so that the above conversion from a date into a datetime will raise an error instead of producing a confusing result. This only applies to datetimes and not timedeltas, because for timedeltas the day <-> hour case is fine, it is just the year/month <-> other units which has issues, and that is already treated with an 'unsafe' casting barrier. Two new functions will facilitate the conversions between datetimes with date units and time units: date_as_datetime(datearray, hour, minute, second, microsecond, timezone='local', unit=None, out=None), which converts the provided dates into datetimes at the specified time, according to the specified timezone. If 'unit' is specified, it controls the output unit, otherwise it is the units in 'out' or the amount of precision specified in the function. datetime_as_date(datetimearray, timezone='local', out=None), which converts the provided datetimes into dates according to the specified timezone. In both functions, timezone can be any of 'UTC', 'TAI', 'local', '+/-####', or a datetime.tzinfo object. The latter will allow NumPy datetimes to work with the pytz library for flexible time zone support. I would also like to extend the 'today' input string parsing to accept strings like 'today 12:30' to allow a convenient way to express different local times occurring today, mostly useful for interactive usage. I welcome any comments on this design, particularly if you can find a case where this doesn't produce a reasonable behavior. Cheers, Mark

3 4

replacing the mechanism for dispatching ufuncs
by Mark Wiebe 27 Sep '11

27 Sep '11

NumPy has a mechanism built in to allow subclasses to adjust or override aspects of the ufunc behavior. While this goal is important, this mechanism only allows for very limited customization, making for instance the masked arrays unable to work with the native ufuncs in a full and proper way. I would like to deprecate the current mechanism, in particular __array_prepare__ and __array_wrap__, and introduce a new method I will describe below. If you've ever used these mechanisms, please review this design to see if it meets your needs. Any class type which would like to override its behavior in ufuncs would define a method called _numpy_ufunc_, and optionally an attribute __array_priority__ as can already be done. The class which wins the priority battle gets its _numpy_ufunc_ function called as follows: return arr._numpy_ufunc_(current_ufunc, *args, **kwargs) To support this overloading, the ufunc would get a new support method, result_type, and there would be a new global function, broadcast_empty_like. The function ufunc.empty_like behaves like the global np.result_type, but produces the output type or a tuple of output types specific to the ufunc, which may follow a different convention than regular arithmetic type promotion. This allows for a class to create an output array of the correct type to pass to the ufunc if it needs to be different than the default. The function broadcast_empty_like is just like empty_like, but takes a list or tuple of arrays which are to be broadcast together for producing the output, instead of just one. Thanks, Mark A simple class which overrides the ufuncs might look as follows: def sin(ufunc, *args, **kwargs): # Convert degrees to radians args[0] = np.deg2rad(args[0]) # Return a regular array, since the result is not in degrees return ufunc(*args, **kwargs) class MyDegreesClass: """Array-like object with a degrees unit""" def __init__(arr): self.arr = arr def _numpy_ufunc_(ufunc, *args, **kwargs): override = globals().get(ufunc.name) if override: return override(ufunc, *args, **kwargs) else: raise TypeError, 'ufunc %s incompatible with MyDegreesClass' % ufunc.name A more complex example will be something like this: def my_general_ufunc(ufunc, *args, **kwargs): # Extract the 'out' argument. This only supports ufuncs with # one output, currently. out = kwargs.get('out') if len(args) > ufunc.nin: if out is None: out = args[ufunc.nin] else: raise ValueError, "'out' given as both a position and keyword argument" # Just want the inputs from here on args = args[:ufunc.nin] # Strip out MyArrayClass, but allow operations with regular ndarrays raw_in = [] for a in args: if isinstance(a, MyArrayClass): raw_in.append(a.arr) else: raw_in.append(a) # Allocate the output array if not out is None: if isinstance(out, MyArrayClass): raise TypeError, "'out' must have type MyArrayClass" else: # Create the output array, obeying the 'order' parameter, # but disallowing subclasses out = np.broadcast_empty_like([args, order=kwargs.get('order'), dtype=ufunc.result_type(args), subok=False) # Override the output argument kwargs['out'] = out.arr # Call the ufunc ufunc(*args, **kwargs) # Return the output return out class MyArrayClass: def __init__(arr): self.arr = arr def _numpy_ufunc_(ufunc, *args, **kwargs): override = globals().get(ufunc.name) if override: return override(ufunc, *args, **kwargs) else: return my_general_ufunc(ufunc, *args, **kwargs)

5 15

load from text files Pull Request Review
by Christopher Jordan-Squire 14 Sep '11

14 Sep '11

Hi-- I've submitted a pull request for a new method for loading data from text files into a record array/masked record array. https://github.com/numpy/numpy/pull/143 Click on the link for more info, but the general idea is to create a regular expression for what entries should look like and loop over the file, updating the regular expression if it's wrong. Once the types are determined the file is loaded line by line into a pre-allocated numpy array. Compared to genfromtxt this function has several advantages/potential advantages. *More modular (genfromtxt is a rather large, nearly 500 line, monolithic function. In my pull request no individual method is longer than around 80 lines, and they're fairly self-contained.) *delimiters can be specified via regex's *missing data can be specified via regex's *it's bit simpler and has sensible defaults *it actually works on some (unfortunately proprietary) data that genfromtxt doesn't seem robust enough for *it supports datetimes *fairly extensible for the power user *makes two passes through the file, the first to determine types/sizes for strings and the second to read in the data, and pre-allocates the array for the second pass. So no giant memory bloating for reading large text files *fairly fast, though I think there is plenty of room for optimizations All that said, it's entirely possible that the innards which determine the type should be ripped out and submitted as a function on their own. I'd love suggestions for improvements, as well as suggestions for a better name. (Currently it's called loadtable, which I don't really like. It was just a working name.) -Chris Jordan-Squire

7 24

NumPy-Discussion August 2011