Mailman 3 October 2014 - NumPy-Discussion

numpy.mean slicing in a netCDF file
by Fadzil Mnor 16 Oct '14

16 Oct '14

Hi all, I wrote a script and plot monthly mean zonal wind (from a netcdf file names uwnd.mon.mean.nc) and I'm not sure if I'm doing it correctly. What I have: ******************************************************************************** *#this part calculates mean values for january only, from 1980-2010; thus the index looks like this 396:757:12* def ujan(): f = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc') u10_1 = f.variables['uwnd'] u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0) return u10_2 uJan = ujan()* #calling function* *#this part is only to define lon, lat and level * q = nc.Dataset('~/data/ncep/uwnd.mon.mean.nc') lon=q.variables['lon'] lat=q.variables['lat'] lev=q.variables['level'] *#for some reason I need to define this unless it gave error "length of x must be number of column in z"* lon=lon[39:43] *#begin plotting* clevs=np.arange(-10.,10.,0.5) fig = plt.figure(figsize=(11, 8)) fig.clf() ax = fig.add_subplot(111) ax.axis([97.5, 105., 1000., 10.]) ax.tick_params(direction='out', which='both') ax.set_xlabel('Lon (degrees)') ax.set_ylabel('Pressure (mb)') ax.set_xticks(np.arange(97.5, 105., .5)) ax.set_yticks([1000, 700, 500, 300, 100, 10]) cs=ax.contourf(lon, lev, uJan, clevs, extend='both',cmap='seismic') plt.title('Zonal winds average (Jan, 1981-2010)') cax = fig.add_axes([0.99, 0.1, 0.03, 0.8]) aa=fig.colorbar(cs,cax=cax,orientation='vertical') aa.set_label('m/s') plt.savefig('~/uwind-crossection-test.png', bbox_inches='tight') ******************************************************************************* the result is attached. I have no idea how to confirm the result (at least until this email is written) , but I believe the lower altitude values should be mostly negative, because over this region, the zonal wind are usually easterly (thus,negative values), but I got positive values. Put the information above aside, *I just want to know if my slicing in the ujan() function is correct*. If it is, then, there must be nothing wrong(except my above mentioned assumption). The data file dimension is: *[time,level,latitude,longitude]* This part: *u10_2 = np.mean(u10_1[396:757:12,:,38,39:43],axis=0)* The line above will calculate the mean of zonal wind (uwnd) in a range of time index 396 to 757 for each year (january only), for all vertical level, at latitude index 38 (5 N) and in between longitude index 39 to 43 (97.5E-105E). I assume it will calculate a 30-year average of zonal wind for january only. Is this correct? Thank you. Fadzil

4 5

Best way to expose std::vector to be used with numpy
by Daniele Nicolodi 16 Oct '14

16 Oct '14

Hello, I have a C++ application that collects float, int or complex data in a possibly quite large std::vector. The application has some SWIG generated python wrappers that expose this vector to python. However, the standard way in which SWIG exposes the data is to create a touple and pass this to python, where it is very often converted to a numpy array for processing. Of course this is not efficient. I would like therefore to generate a smarter python interface. In python 3 I would without doubt go for implementing the buffer protocol, which enables seamless interfacing with numpy. However, I need to support also python 2 and there the buffer protocol is not as nice. What is the best thing to do to expose the data buffer to python 2, without introducing a dependency on numpy in the application? Bonus points if the proposed solution comes with a pointer on how to implement in in SWIG (which I don't know much and feel a bit too magic for my taste). Thank you! Cheers, Daniele

10 12

Endianness not detected in OS X Lion
by Nicola 16 Oct '14

16 Oct '14

Hi, I am trying to install numpy 1.9.0 on OS X Lion 10.7.5 and I get the same error as reported in this old thread: http://thread.gmane.org/gmane.comp.python.numeric.general/47711/ that is, numpy/core/src/npymath/npy_math_private.h:78:3: error: conflicting types for 'ieee_double_shape_type' } ieee_double_shape_type; ^ numpy/core/src/npymath/npy_math_private.h:64:3: note: previous declaration of 'ieee_double_shape_type' was here } ieee_double_shape_type; ^ In file included from numpy/core/src/npymath/npy_math.c.src:56:0: numpy/core/src/npymath/npy_math_private.h:78:3: error: conflicting types for 'ieee_double_shape_type' } ieee_double_shape_type; ^ numpy/core/src/npymath/npy_math_private.h:64:3: note: previous declaration of 'ieee_double_shape_type' was here } ieee_double_shape_type; ^ no matter whether I use the system's clang or gcc 4.9 (installed with Homebrew). Is there a workaround for this? Nicola

2 2

array indexing question
by Neal Becker 15 Oct '14

15 Oct '14

I'm using np.nonzero to construct the tuple: (array([0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2]), array([1, 3, 5, 7, 2, 3, 6, 7, 4, 5, 6, 7])) Now what I want is the 2-D index array: [1,3,5,7, 2,3,6,7, 4,5,6,7] Any ideas? -- -- Those who don't understand recursion are doomed to repeat it

2 1

Any interest in a generalized piecewise function?
by Per.Brodtkorb＠ffi.no 14 Oct '14

14 Oct '14

I have worked on a generalized piecewise function (genpiecewise) that are simpler and more general than the current numpy.piecewise implementation. The new generalized piecewise function allows functions of the type f(x0, x1,.. , xn) i.e. to have arbitrary number of input arguments that are evaluated conditionally. The generalized piecewise function passes all the tests for numpy.piecewise function except the undocumented features of numpy.piecewise which allows condlist to be a single bool list/array or a single int array. A new keyword "fillvalue" opens up the possibility to specify "out of bounds values" to other values than 0 eg. Np.nan. Examples: Example 1) >>> x = numpy.linspace(-2,2,5) >>> numpy.piecewise(x, x<1, [1]) # = numpy.where(x<0, 1, 0) array([ 1., 1., 1., 0., 0.]) # can be written as >>> genpiecewise([x<1],[1]) array([1, 1, 1, 0, 0]) # or >>> genpiecewise([x<1],[1], x) array([ 1., 1., 1., 0., 0.]) # or >>> genpiecewise([x<1],[1.0]) array([ 1., 1., 1., 0., 0.]) Example 2) >>> numpy.piecewise(x, [x < 0, x >= 0], [lambda x: -x, lambda x: x]) array([ 2., 1., 0., 1., 2.]) # can be written as >>> genpiecewise([x < 0, x >= 0], [lambda x: -x, lambda x: x], (x,)) array([ 2., 1., 0., 1., 2.]) # or genpiecewise([x < 0,], [lambda x: -x, lambda x: x], x) array([ 2., 1., 0., 1., 2.]) Example 3) # New functionality >>> X,Y = numpy.meshgrid(x,x) >>> genpiecewise([X*Y<0,], [lambda x,y: -x*y, lambda x,y: x*y], xi=(X,Y)) array([[ 4., 2., -0., 2., 4.], [ 2., 1., -0., 1., 2.], [-0., -0., 0., 0., 0.], [ 2., 1., 0., 1., 2.], [ 4., 2., 0., 2., 4.]]) >>> genpiecewise([X*Y<-0.5, X*Y>0.5], [lambda x,y: -x*y, lambda x,y: x*y], xi=(X,Y), fillvalue=numpy.nan) array([[ 4., 2., nan, 2., 4.], [ 2., 1., nan, 1., 2.], [ nan, nan, nan, nan, nan], [ 2., 1., nan, 1., 2.], [ 4., 2., nan, 2., 4.]]) >>> genpiecewise([X*Y<-0.5, X*Y>0.5], [lambda x,y: -x*y, lambda x,y: x*y, numpy.nan], (X,Y)) array([[ 4., 2., nan, 2., 4.], [ 2., 1., nan, 1., 2.], [ nan, nan, nan, nan, nan], [ 2., 1., nan, 1., 2.], [ 4., 2., nan, 2., 4.]]) My question is: are there any interest in the community for such a function? Could some or all of this functionality replace the current numpy.piecewise? (This function could replace the function _lazywhere (in scipy.stats._distn_infrastructure) which is heavily used in scipy.stats) Per A. Brodtkorb The code looks like this: def valarray(shape, value=np.nan, typecode=None): """Return an array of all value. """ out = ones(shape, dtype=bool) * value if typecode is not None: out = out.astype(typecode) if not isinstance(out, np.ndarray): out = asarray(out) return out def genpiecewise(condlist, funclist, xi=None, fillvalue=0, args=(), **kw): """ Evaluate a piecewise-defined function. Given a set of conditions and corresponding functions, evaluate each function on the input data wherever its condition is true. Parameters ---------- condlist : list of bool arrays Each boolean array corresponds to a function in `funclist`. Wherever `condlist[i]` is True, `funclist[i](x0,x1,...,xn)` is used as the output value. Each boolean array in `condlist` selects a piece of `xi`, and should therefore be of the same shape as `xi`. The length of `condlist` must correspond to that of `funclist`. If one extra function is given, i.e. if ``len(funclist) - len(condlist) == 1``, then that extra function is the default value, used wherever all conditions are false. funclist : list of callables, f(*(xi + args), **kw), or scalars Each function is evaluated over `x` wherever its corresponding condition is True. It should take an array as input and give an array or a scalar value as output. If, instead of a callable, a scalar is provided then a constant function (``lambda x: scalar``) is assumed. xi : tuple input arguments to the functions in funclist, i.e., (x0, x1,...., xn) fillvalue : scalar fillvalue for out of range values. Default 0. args : tuple, optional Any further arguments given here passed to the functions upon execution, i.e., if called ``piecewise(..., ..., args=(1, 'a'))``, then each function is called as ``f(x0, x1,..., xn, 1, 'a')``. kw : dict, optional Keyword arguments used in calling `piecewise` are passed to the functions upon execution, i.e., if called ``piecewise(..., ..., lambda=1)``, then each function is called as ``f(x0, x1,..., xn, lambda=1)``. Returns ------- out : ndarray The output is the same shape and type as x and is found by calling the functions in `funclist` on the appropriate portions of `x`, as defined by the boolean arrays in `condlist`. Portions not covered by any condition have undefined values. See Also -------- choose, select, where Notes ----- This is similar to choose or select, except that functions are evaluated on elements of `xi` that satisfy the corresponding condition from `condlist`. The result is:: |-- |funclist[0](x0[condlist[0]],x1[condlist[0]],...,xn[condlist[0]]) out = |funclist[1](x0[condlist[1]],x1[condlist[1]],...,xn[condlist[1]]) |... |funclist[n2](x0[condlist[n2]],x1[condlist[n2]],...,xn[condlist[n2]]) |-- Examples -------- Define the sigma function, which is -1 for ``x < 0`` and +1 for ``x >= 0``. >>> x = np.linspace(-2.5, 2.5, 6) >>> genpiecewise([x < 0, x >= 0], [-1, 1]) array([-1., -1., -1., 1., 1., 1.]) Define the absolute value, which is ``-x`` for ``x <0`` and ``x`` for ``x >= 0``. >>> genpiecewise([x < 0, x >= 0], [lambda x: -x, lambda x: x], (x,)) array([ 2.5, 1.5, 0.5, 0.5, 1.5, 2.5]) """ def otherwise_condition(condlist): return ~np.logical_or.reduce(condlist, axis=0) def check_shapes(condlist, funclist): nc, nf = len(condlist), len(funclist) if nc not in [nf-1, nf]: raise ValueError("function list and condition list" + " must be the same length") check_shapes(condlist, funclist) if xi is not None and not isinstance(xi, tuple): xi = (xi,) condlist = np.broadcast_arrays(*condlist) if len(condlist) == len(funclist)-1: condlist.append(otherwise_condition(condlist)) arrays = dtype = None if xi is not None: arrays = np.broadcast_arrays(*xi) dtype = np.result_type(*arrays) else: # funclist is a list of scalars only dtype = np.result_type(*funclist) out = valarray(condlist[0].shape, fillvalue, dtype) for cond, func in zip(condlist, funclist): if isinstance(func, collections.Callable): temp = tuple(np.extract(cond, arr) for arr in arrays) + args np.place(out, cond, func(*temp, **kw)) else: # func is a scalar value np.place(out, cond, func) return out

2 3

non-linear rebin function recipe?
by RayS 14 Oct '14

14 Oct '14

Most of my work has used the Fourrier based method for "linear" rebin of evenly sampled time data of length m (say 1500) to a new number of samples n (say 2048); the delta time change per sample is a constant over the array. I'd like to test the effect a non-constant delta t, ie, stretching some regions of data while leaving others static, and transitioning in a smooth way. - Ray Schumacher

1 0

numpy 1.9.1 rc1 probably next weekend
by Julian Taylor 14 Oct '14

14 Oct '14

hi, we have collected a couple bugfixes in the 1.9 branch that should get released soon. The most important bugs fixes will be: - fix matplotlibs build (#5067) - avoid triggering a compiler bug on redhat 5, gcc 4.1.2 (#5163) - fix compatibility of npy files between py2 and py3 on win32/64 (#5170) There is still a bug in nanmedian if you have infinities in the array (#5138) which I hope will be fixed soon. After that I'd like to cut a 1.9.1rc1 with final release the week after. If there is something missing in the 1.9.x branch please report or fix it soon. Cheers, Julian

1 0

Detect if array has been transposed
by Mads Ipsen 13 Oct '14

13 Oct '14

Hi, In part of my C++ code, I often do void foo(PyObject * matrix) { do stuff } where matrix is a numpy mxn matrix created on the Python side, where foo() eg. is invoked as a = numpy.array([[1,2],[3,5]]) foo(a) However, if you call transpose() on a, some care should be taken, since numpy's internal matrix data first gets transposed on demand. In that case I must do a = numpy.array([[1,2],[3,5]]) a.transpose() foo(a.copy()) to make sure the correct data of the array gets transferred to the C++ side. Is there any way for me to detect (on the Python side) that transpose() has been invoked on the matrix, and thereby only do the copy operation when it really is needed? For example if a_has_transposed_data: foo(a.copy()) else: foo(a) Best regards, Mads -- +---------------------------------------------------------+ | Mads Ipsen | +----------------------+----------------------------------+ | Gåsebæksvej 7, 4. tv | phone: +45-29716388 | | DK-2500 Valby | email: mads.ipsen(a)gmail.com | | Denmark | map : www.tinyurl.com/ns52fpa | +----------------------+----------------------------------+

4 5

[ANN] bcolz 0.7.2
by Valentin Haenel 13 Oct '14

13 Oct '14

====================== Announcing bcolz 0.7.2 ====================== What's new ========== This is a maintenance release that fixes various bits and pieces. Importantly, compatibility with Numpy 1.9 and Cython 0.21 has been fixed and the test suit no longer segfaults on 32 bit UNIX. Feature-wise a new ``carray.view()`` method has been introduced which allows carrays to share the same raw data. ``bcolz`` is a renaming of the ``carray`` project. The new goals for the project are to create simple, yet flexible compressed containers, that can live either on-disk or in-memory, and with some high-performance iterators (like `iter()`, `where()`) for querying them. 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… For more detailed info, see the release notes in: https://github.com/Blosc/bcolz/wiki/Release-Notes What it is ========== bcolz provides columnar and compressed data containers. 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, a high-performance compressor that is optimized for binary data. 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, 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. 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 ---- **Enjoy data!**

1 0

parallel compilation with numpy.distutils in numpy 1.10
by Julian Taylor 11 Oct '14

11 Oct '14

hi, To speed up compilation of extensions I have made a PR to compile extension files in parallel: https://github.com/numpy/numpy/pull/5161 It adds the --jobs/-j flags to the build command of setup.py which defines the number of parallel compile processes. E.g. python setup.py build --jobs 4 install --prefix /tmp/local Additionally it adds the environment variable NPY_NUM_BUILD_JOBS which is used if no commandline is set. This helps e.g. with pip installations, travis builds (which give you 1.5 cpus) or to put in your .bashrc. The parallelization is only with the files of an extension so it is not super efficient but an uncached numpy build goes down from 1m40s to 1m00s with 3 cores on my machine which is quite decent. Building scipy from scratch decreased from 10minutes to 6m30s on my machine. Unfortunately projects using cython will not profit as cython tends to build an extension from one single file. (You may want to look into gccs internal parallelization for that, -flto=jobs)- Does some see issues with the interface I have currently set? Please speak up soon. There is still one problem in regards to parallelizing fortran 90. The ccompiler.py contains following comment: # build any sources in same order as they were originally specified # especially important for fortran .f90 files using modules This indicates the f90 builds cannot be trivially parallelized. I do not know much fortran, can someone explain to me when ordering of single file compiles is an issue in f90? Cheers, Julian Taylor

2 5