Mailman 3 January 2015 - NumPy-Discussion

Aligned array allocations
by Antoine Pitrou 22 Jan '15

22 Jan '15

Hello, In https://github.com/numpy/numpy/issues/5312 there's a request for an aligned allocator in Numpy (more than the default alignment of the platform's memory allocator). The reason is that on modern vectorization instruction sets, a certain alignment is required for optimal performance (even though unaligned data still works: it's just that performance is degraded... by how much will depend on the CPU micro-architecture). For example Intel recommends a 32-byte alignment for AVX loads and stores. In https://github.com/numpy/numpy/pull/5457 I have proposed a patch to wrap the system allocator in an aligned allocator. The proposed scheme makes the alignment configurable at runtime (through a Python API), because different platforms may have different desirable alignments, and it is not reasonable for Numpy to know about them all, nor for users to recompile Numpy each time they have a different CPU. By always using an aligned allocator there is some overhead: - all arrays occupy a bit more memory by a small average amount (probably 16 bytes average on a 64-bit machine, for a 16 byte guaranteed alignment) - array resizes can be more expensive in CPU time, when the physical start changes and its alignment changes too There is also a limitation: while the physical start of an array will always be aligned, this can be defeated when taking a view starting at a non-zero index. (note that to take advantage of certain instruction set features such as AVX, Numpy may need to be compiled with specific compiler flags... but Numpy's allocations also affect other packages such as Numba which is able to generate code at runtime) I would like to know if people are interested in this feature, and if the proposed approach is acceptable. Regards Antoine.

2 1

EDF+ specification
by Eric Firing 20 Jan '15

20 Jan '15

http://www.edfplus.info/specs/edfplus.html#additionalspecs Io, Is this the file format you have? Eric

2 2

structured arrays, recarrays, and record arrays
by Allan Haldane 18 Jan '15

18 Jan '15

Hello all, Documentation of recarrays is poor and I'd like to improve it. In order to do this I've been looking at core/records.py, and I would appreciate some feedback on my plan. Let me start by describing what I see. In the docs there is some confusion about 'structured arrays' vs 'record arrays' vs 'recarrays' - the docs use them often interchangeably. They also refer to structured dtypes alternately as 'struct data types', 'record data types' or simply 'records' (eg, see the reference/arrays.dtypes and reference/arrays.indexing doc pages). But by my reading of the code there are really three (or four) distinct types of arrays with structure. Here's a possible nomenclature: * "Structured arrays" are simply ndarrays with structured dtypes. That is, the data type is subdivided into fields of different type. * "recarrays" are a subclass of ndarrays that allow access to the fields by attribute. * "Record arrays" are recarrays where the elements have additionally been converted to 'numpy.core.records.record' type such that each data element is an object with field attributes. * (it is also possible to create arrays with dtype.dtype of numpy.core.records.record, but which are not recarrays. However I have never seen this done.) Here's code demonstrating the creation of the different types of array (in order: structured array, recarray, ???, record array). >>> arr = np.array([(1,'a'), (2,'b')], dtype=[('foo', int), ('bar', 'S1')]) >>> recarr = arr.view(type=np.recarray) >>> noname = arr.view(dtype=dtype(np.record, arr.dtype)) >>> recordarr = arr.view(dtype=dtype((np.record, arr.dtype)), type=np.recarray) >>> type(arr), arr.dtype.type (numpy.ndarray, numpy.void) >>> type(recarr), recarr.dtype.type (numpy.core.records.recarray, numpy.void) >>> type(recordarr), recordarr.dtype.type (numpy.core.records.recarray, numpy.core.records.record) Note that the functions numpy.rec.array, numpy.rec.fromrecords, numpy.rec.fromarrays, and np.recarray.__new__ create record arrays. However, in the docs you can see examples of the creation of recarrays, eg in the recarray and ndarray.view doctrings and in http://www.scipy.org/Cookbook/Recarray. The files numpy/lib/recfunctions.py and numpy/lib/npyio.py (and possibly masked arrays, but I haven't looked yet) make extensive use of recarrays (but not record arrays). The main functional difference between recarrays and record arrays is field access on individual elements: >>> recordarr[0].foo 1 >>> recarr[0].foo Traceback (most recent call last): File "<stdin>", line 1, in <module> AttributeError: 'numpy.void' object has no attribute 'foo' Also, note that recarrays have a small performance penalty relative to structured arrays, and record arrays have another one relative to recarrays because of the additional python logic. So my first goal in updating the docs is to use the right terms in the right place. In almost all cases, references to 'records' (eg 'record types') should be replaced with 'structured' (eg 'structured types'), with the exception of docs that deal specifically with record arrays. It's my guess that in the distant past structured datatypes were intended to always be of type numpy.core.records.record (thus the description in reference/arrays.dtypes) but that numpy.core.records.record became generally obsolete without updates to the docs. doc/records.rst.txt seems to document the transition. I've made a preliminary pass of the docs, which you can see here https://github.com/ahaldane/numpy/commit/d87633b228dabee2ddfe75d1ee9e41ba70… Mostly I renamed 'record type' to 'structured type', and added a very rough draft to numpy/doc/structured_arrays.py. I would love to hear from those more knowledgeable than myself on whether this works! Cheers, Allan

1 1

ANN: Scipy 0.15.1
by Pauli Virtanen 18 Jan '15

18 Jan '15

-----BEGIN PGP SIGNED MESSAGE----- Hash: SHA1 Dear all, We are pleased to announce the Scipy 0.15.1 release. Scipy 0.15.1 contains only bugfixes. The module ``scipy.linalg.calc_lwork`` removed in Scipy 0.15.0 is restored. This module is not a part of Scipy's public API, and although it is available again in Scipy 0.15.1, using it is deprecated and it may be removed again in a future Scipy release. Source tarballs, binaries, and full release notes are available at https://sourceforge.net/projects/scipy/files/scipy/0.15.1/ Best regards, Pauli Virtanen ========================== SciPy 0.15.1 Release Notes ========================== SciPy 0.15.1 is a bug-fix release with no new features compared to 0.15.0. Issues fixed - ------------ * `#4413 <https://github.com/scipy/scipy/pull/4413>`__: BUG: Tests too strict, f2py doesn't have to overwrite this array * `#4417 <https://github.com/scipy/scipy/pull/4417>`__: BLD: avoid using NPY_API_VERSION to check not using deprecated... * `#4418 <https://github.com/scipy/scipy/pull/4418>`__: Restore and deprecate scipy.linalg.calc_work -----BEGIN PGP SIGNATURE----- Version: GnuPG v1 iEYEARECAAYFAlS8CA4ACgkQ6BQxb7O0pWCmOQCgzg9AXDaqRaK5/QBWopIrv2OA WkEAn0ltDfDHFpw0zMzB9mUscAAb2xnE =JrGj -----END PGP SIGNATURE-----

1 0

using f2py with a module containing a derived type?
by Greg Fischer 17 Jan '15

17 Jan '15

Hello, I would like to use f2py to wrap a Fortran module that contains a derived data type. I don't necessarily need to access the data that is inside the derived type from Python, but I would really like to be able to call the subroutines that are contained inside the module. When I attempt to use f2py on this module, it appears to choke when it gets to the derived data type. Is there any way around this? Thanks, Greg

2 1

Re: [Numpy-discussion] Sorting refactor
by Lars Buitinck 16 Jan '15

16 Jan '15

2015-01-16 11:55 GMT+01:00 <numpy-discussion-request(a)scipy.org>: > Message: 2 > Date: Thu, 15 Jan 2015 21:24:00 -0800 > From: Jaime Fern?ndez del R?o <jaime.frio(a)gmail.com> > Subject: [Numpy-discussion] Sorting refactor > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: > <CAPOWHWkF6RnWcrGmcwsmq_LO3hShjgBVLsrN19z-MDPe25E2Aw(a)mail.gmail.com> > Content-Type: text/plain; charset="utf-8" > > This changes will make it easier for me to add a Timsort generic type > function to numpy's arsenal of sorting routines. And I think they have > value by cleaning the source code on their own. Yes, they do. I've been looking at the sorting functions as well and I've found the following: * The code is generally hard to read because it prefers pointer over indices. I'm wondering if it would get slower using indices. The closer these algorithms are to the textbook, the easier to insert fancy optimizations. * The heap sort exploits undefined behavior by using a pointer that points before the start of the array. However, rewriting it to always point within the array made it slower. I haven't tried rewriting it using indices. * Quicksort has a quadratic time worst case. I think it should be turned into an introsort [1] for O(n log n) worst case; we have the heapsort needed to do that. * Quicksort is robust to repeated elements, but doesn't exploit them. It can be made to run in linear time if the input array has only O(1) distinct elements [2]. This may come at the expense of some performance on arrays with no repeated elements. * Using optimal sorting networks instead of insertion sort as the base case can speed up quicksort on float arrays by 5-10%, but only if NaNs are moved out of the way first so that comparisons become cheaper [3]. This has consequences for the selection algorithms that I haven't fully worked out yet. * Using Cilk Plus to parallelize merge sort can make it significantly faster than quicksort at the expense of only a few lines of code, but I haven't checked whether Cilk Plus plays nicely with multiprocessing and other parallelism options (remember the trouble with OpenMP-ified OpenBLAS?). This isn't really an answer to your questions, more like a brain dump from someone who's stared at the same code for a while and did some experiments. I'm not saying we should implement all of this, but keep in mind that there are some interesting options besides implementing timsort. [1] https://en.wikipedia.org/wiki/Introsort [2] http://www.sorting-algorithms.com/static/QuicksortIsOptimal.pdf [3] https://github.com/larsmans/numpy/tree/sorting-nets

6 6

Re: [Numpy-discussion] NumPy-Discussion Digest, Vol 100, Issue 28
by Lars Buitinck 16 Jan '15

16 Jan '15

2015-01-16 13:29 GMT+01:00 <numpy-discussion-request(a)scipy.org>: > Date: Fri, 16 Jan 2015 12:43:43 +0100 > From: Julian Taylor <jtaylor.debian(a)googlemail.com> > Subject: Re: [Numpy-discussion] Sorting refactor > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: <54B8F96F.7090903(a)googlemail.com> > Content-Type: text/plain; charset=windows-1252 > > On 16.01.2015 12:33, Lars Buitinck wrote: >> * Quicksort has a quadratic time worst case. I think it should be >> turned into an introsort [1] for O(n log n) worst case; we have the >> heapsort needed to do that. > > This probably rarely happens in numeric data, and we do have guaranteed > nlog runtime algorithms available. It's no more likely or unlikely than in any other type of data (AFAIK), but it's possible for an adversary to DOS any (web server) code that uses np.sort. > I was also thinking about this, an advantage of a sorting network is > that it can be vectorized to be significantly faster than an insertion > sort. Handling NaN's should also be directly possible. Tried that, and it didn't give any speedup, at least not without explicit vector instructions. Just moving the NaNs aside didn't cost anything in my preliminary benchmarks (without sorting nets), the cost of the operation was almost exactly compensated by simpler comparisons. > The issue is that its probably too much complicated code for only a very > small gain. Maybe. The thing is that the selection algorithms are optimized for NaNs and seem to depend on the current comparison code. We'd need distinct <TYPE>_LT and <TYPE>_LT_NONAN for each <TYPE>. The sorting nets themselves aren't complicated, just lengthy. My branch has the length-optimal (not optimally parallel) ones for n <= 16. > But I assume you will get the same trouble as openmp but that needs > testing, also cilk+ in gcc is not really production ready yet, I got > lots of crashes when I last tried it (it will be in 5.0 though). The data parallel constructs tend to crash the compiler, but task spawning seems to be stable in 4.9.2. I've still to see how it handles multiprocessing/fork. What do you mean by will be in 5.0, did they do a big push? > Date: Fri, 16 Jan 2015 13:28:56 +0100 > From: Da?id <davidmenhur(a)gmail.com> > Subject: Re: [Numpy-discussion] Sorting refactor > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > Message-ID: > <CAJhcF=1O5Own_5ydzu+To8HHbm3e66k=iuNQrEiaSDY23DNW6w(a)mail.gmail.com> > Content-Type: text/plain; charset=UTF-8 > > On 16 January 2015 at 13:15, Eelco Hoogendoorn > <hoogendoorn.eelco(a)gmail.com> wrote: >> Perhaps an algorithm can be made faster, but often these multicore >> algorithms are also less efficient, and a less data-dependent way of putting >> my cores to good use would have been preferable. Potentially, other code >> could slow down due to cache trashing if too many parallel tasks run in >> parallel. Id rather be in charge of such matters myself; but I imagine >> adding a keyword arg for these matters would not be much trouble? > > As I understand it, that is where the strength of Cilk+ lies. It does > not force parallelisation, just suggests it. The decision to actually > spawn parallel is decided at runtime depending on the load of the > other cores. cilk+ guarantees that the amount of space used by a pool of P threads is at most P times the stack space used by the sequential version (+ a constant). The idea is that you can say for (i = 0; i < 1000000; i++) { cilk_spawn f(a[i]); } and it will never create more than P work items in memory, rather than 1e6, even if each f() spawns a bunch itself. Of course, it won't guarantee that OpenMP will not also spawn P threads and/or check that you're one of P processes cooperating on a task using multiprocessing. Personally I'd rather have an np.setnumthreads() to turn this on or off for a process and have the runtime distribute work for me instead of having to do it myself.

3 2

Re: [Numpy-discussion] Sorting refactor
by Lars Buitinck 16 Jan '15

16 Jan '15

2015-01-16 15:14 GMT+01:00 <numpy-discussion-request(a)scipy.org>: > Date: Fri, 16 Jan 2015 06:11:29 -0800 > From: Jaime Fern?ndez del R?o <jaime.frio(a)gmail.com> > Subject: Re: [Numpy-discussion] Sorting refactor > To: Discussion of Numerical Python <numpy-discussion(a)scipy.org> > > Most of my proposed original changes do not affect the core sorting > functionality, just the infrastructure around it. But if we agree that > sorting has potential for being an actively developed part of the code > base, then cleaning up its surroundings for clarity makes sense, so I'm > taking your brain dump as an aye for my proposal. ;-) It is!

1 0

Sorting refactor
by Jaime Fernández del Río 16 Jan '15

16 Jan '15

Hi all, I have been taking a deep look at the sorting functionality in numpy, and I think it could use a face lift in the form of a big code refactor, to get rid of some of the ugliness in the code and make it easier to maintain. What I have in mind basically amounts to: 1. Refactor _new_argsortlike to get rid of code duplication (there are two branches, one with buffering, one without, virtually identical, that could be merged into a single one). 2. Modify _new_argsortlike so that it can properly handle byte-swapped inputs of any dtype, see gh-5441. Add proper handling of types with references, in preparation for the rest of changes. 3. Add three functions to the npy_sort library: npy_aquicksort, npy_aheapsort, npy_amergesort, with a signature compatible with PyArray_ArgSortFunc , i.e. (char* start, npy_intp* result, npy_intp length, PyArrayObject *arr). These turn out to be almost identical to the string and unicode sort functions, but using the dtype's compare function to handle comparisons. 4. Modify PyArray_ArgSort (and PyArray_ArgPartition) to always call _new_argsortlike, even when there is no type specific argsort function, by using the newly added npy_axxx functions. This simplifies PyArray_ArgSort a lot, and gets rid of some of the global variable ugliness in the current code. And makes argsorting over non-contiguous axis more memory efficient. 5. Refactor _new_sortlike similarly to _new_argsortlike 6. Modify the npy_quicksort, npy_mergesort and npy_heapsort functions in npy_sort to have a signature compatible with PyArray_SortFunc, i.e. (char* start, npy_intp length, PyArrayObject *arr). npy_quicksort will no longer rely on libc's qsort, but be very similar to the string or unicode quicksort functions 7. Modify PyArray_Sort (and PyArray_Partition) to always call _new_sortlike, similarly to what was done with PyArray_ArgSort. This allows completing the removal of the remaining global variable ugliness, as well as similar benefits as for argsort before. This changes will make it easier for me to add a Timsort generic type function to numpy's arsenal of sorting routines. And I think they have value by cleaning the source code on their own. So my questions, mostly to the poor souls that will have to code review changes to several hundred lines of code: 1. Does this make sense, or is it better left alone? A subset of 1, 2 and 5 are a must to address the issues in gh-5441, the rest could arguably be left as is. 2. Would you rather see it submitted as one ginormous PR? Or split into 4 or 5 incremental ones? Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

2 2

Optimizing numpy's einsum expression (again)
by Daniel Smith 16 Jan '15

16 Jan '15

Hello everyone, I originally brought an optimized einsum routine forward a few weeks back that attempts to contract numpy arrays together in an optimal way. This can greatly reduce the scaling and overall cost of the einsum expression for the cost of a few intermediate arrays. The current version (and more details) can be found here: https://github.com/dgasmith/opt_einsum <https://github.com/dgasmith/opt_einsum> I think this routine is close to a finalized version, but there are two problems which I would like the community to weigh in on: Memory usage- Currently the routine only considers the maximum size of intermediates created rather than cumulative memory usage. If we only use np.einsum it is straightforward to calculate cumulative memory usage as einsum does not require any copies of inputs to be made; however, if we attempt to use a vendor BLAS the memory usage can becomes quite complex. For example, the np.tensordot routine always forces a copy for ndarrays because it uses ndarray.transpose(…).reshape(…). A more memory-conscious way to do this is to first try and do ndarray.reshape(…).T which does not copy the data and numpy can just pass a transpose flag to the vendor BLAS. The caveat here is that the summed indices must be in the correct order- if not a copy is required. Maximum array size is usually close to the total overhead of the opt_einsum function, but can occasionally be 2-5 times this size. I see the following ways forward: Ignore cumulative memory and stick with maximum array size as the limiting memory factor. Implement logic to figure out if the input arrays needs to be copied to use BLAS, compute the extra memory required, and add an extra dimension to the optimization algorithms (use BLAS or do not use BLAS at each step). Some of this is already done, but may get fairly complex. Build an in-place nd-transpose algorithm. Cut out BLAS entirely. Keeping in mind that vendor BLAS can be orders of magnitude faster than a pure einsum implementation, especially if the BLAS threading is used. Path algorithms- There are two algorithms “optimal” (a brute force algorithm, scales like N!) and “opportunistic” (a greedy algorithm, scales like N^3). The optimal path can take seconds to minutes to calculate for a 7-10 term expression while the opportunistic path takes microseconds even for 20+ term expressions. The opportunistic algorithm works surprisingly well and appears to obtain the correct scaling in all test cases that I can think of. Both algorithms use the maximum array size as a sieve- this is very beneficial from several aspects. The problem occurs when a much needed intermediate cannot be created due to this limit- on occasions not making this intermediate can have slowdowns of orders of magnitude even for small systems. This leads to certain (and sometimes unexpected) performance walls. Possible solutions: Warn the user if the ratio between an unlimited memory solution and a limited memory solution becomes large. Do not worry about it. Thank you for your time, -Daniel Smith

3 2