Mailman 3 May 2016 - NumPy-Discussion

NumPy 1.11 docs
by Stephan Hoyer 31 May '16

31 May '16

These still are missing from the SciPy.org page, several months after the release. What do we need to do to keep these updated? Is there someone at Enthought we should ping? Or do we really just need to transition to different infrastructure?

2 2

Cross-correlation PR stuck in limbo
by Honi Sanders 27 May '16

27 May '16

Hello all, I have completed a pull request to add a “maxlag” functionality to numpy.correlate. See here: https://github.com/numpy/numpy/pull/5978 <https://github.com/numpy/numpy/pull/5978>. This pull request has passed all tests and has been ready to be merged for around six months. Several people have commented requesting for it to be included on stackoverflow, the listserve, and github. Can someone please let me know what needs to be done or can it be merged? Here is some background: What was 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! I have introduced this question as a numpy issue <https://github.com/numpy/numpy/issues/5954>, a scipy issue <https://github.com/scipy/scipy/issues/4940> and on the scipy-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. Previous discussion of this functionality can be found at another discussion on numpy correlate (and convolution) <http://numpy-discussion.10968.n7.nabble.com/another-discussion-on-numpy-cor…>. Other issues related to correlate functions include ENH: Fold fftconvolve into convolve/correlate functions as a parameter #2651 <https://github.com/scipy/scipy/issues/2651>, Use FFT in np.correlate/convolve? (Trac #1260) #1858 <https://github.com/numpy/numpy/issues/1858>, and normalized cross-correlation (Trac #1714) #2310 <https://github.com/numpy/numpy/issues/2310>. The new implementation allows new types of the “mode” argument, to include an int value, which defines the maximum lag for which cross-correlation should be calculated, or a tuple, which defines the minlag, maxlag, and lagstep to be used in the same format as the arguments to numpy.arange. Please let me know what should be done to move this pull request forward. Honi

4 3

Numpy 1.11.1rc1 release
by Charles R Harris 27 May '16

27 May '16

Hi All, I am pleased to announce the release of Numpy 1.11.1rc1. The sources may be found on sourceforge <https://sourceforge.net/projects/numpy/files/NumPy/1.11.1rc1/> and wheels for OS X, Windows, and Linux will be available on pypi sometime in the next few days. The pypi release is delayed due to the decision that the wheels should go up before the sources in order that people not get a source install when what they want are wheels. The Python versions supported are 2.6-2.7 and 3.2-3.5. This release has mostly small fixes and build enhancements and should be good out of the starting gate, but prudence requires a release candidate as there are a few bits not tested in master. The following fixes have been applied: - #7506 BUG: Make sure numpy imports on python 2.6 when nose is unavailable. - #7530 BUG: Floating exception with invalid axis in np.lexsort. - #7535 BUG: Extend glibc complex trig functions blacklist to glibc < 2.18. - #7551 BUG: Allow graceful recovery for no compiler. - #7558 BUG: Constant padding expected wrong type in constant_values. - #7578 BUG: Fix OverflowError in Python 3.x. in swig interface. - #7590 BLD: Fix configparser.InterpolationSyntaxError. - #7597 BUG: Make np.ma.take work on scalars. - #7608 BUG: linalg.norm(): Don't convert object arrays to float. - #7638 BLD: Correct C compiler customization in system_info.py. - #7654 BUG: ma.median of 1d array should return a scalar. - #7656 BLD: Remove hardcoded Intel compiler flag -xSSE4.2. - #7660 BUG: Temporary fix for str(mvoid) for object field types. - #7665 BUG: Fix incorrect printing of 1D masked arrays. - #7670 BUG: Correct initial index estimate in histogram. - #7671 BUG: Boolean assignment no GIL release when transfer needs API. - #7676 BUG: Fix handling of right edge of final histogram bin. - #7680 BUG: Fix np.clip bug NaN handling for Visual Studio 2015. The following people have contributed to this release - Allan Haldane - Amit Aronovitch - Charles Harris - Eric Wieser - Evgeni Burovski - Loïc Estève - Mathieu Lamarre - Matthew Brett - Matthias Geier - Nathaniel J. Smith - Nikola Forró - Ralf Gommers - Robert Kern - Sebastian Berg - Simon Conseil - Simon Gibbons - Sorin Sbarnea - chiffa Chuck

2 1

Fwd: ifft padding
by Allen Welkie 26 May '16

26 May '16

I'd like to get some feedback on my [pull request]( https://github.com/numpy/numpy/pull/7593). This pull request adds a function `ifftpad` which pads a spectrum by inserting zeros where the highest frequencies would be. This is necessary because the padding that `ifft` does simply inserts zeros at the end of the array. But because of the way the spectrum is laid out, this changes which bins represent which frequencies and in general messes up the result of `ifft`. If you pad with the proposed `ifftpad` function, the zeros will be inserted in the middle of the spectrum and the time signal that results from `ifft` will be an interpolated version of the unpadded time signal. See the discussion in [issue #1346]( https://github.com/numpy/numpy/issues/1346). The following is a script to demonstrate what I mean: ``` import numpy from numpy import concatenate, zeros from matplotlib import pyplot def correct_padding(a, n, scale=True): """ A copy of the proposed `ifftpad` function. """ spectrum = concatenate((a[:len(a) // 2], zeros(n - len(a)), a[len(a) // 2:])) if scale: spectrum *= (n / len(a)) return spectrum def plot_real(signal, label): time = numpy.linspace(0, 1, len(signal) + 1)[:-1] pyplot.plot(time, signal.real, label=label) def main(): spectrum = numpy.zeros(10, dtype=complex) spectrum[-1] = 1 + 1j signal = numpy.fft.ifft(spectrum) signal_bad_padding = numpy.fft.ifft(10 * spectrum, 100) signal_good_padding = numpy.fft.ifft(correct_padding(spectrum, 100)) plot_real(signal, 'No padding') plot_real(signal_bad_padding, 'Bad padding') plot_real(signal_good_padding, 'Good padding') pyplot.legend() pyplot.show() if __name__ == '__main__': main() ```

2 1

AUTO: Leon Woo is out of the office (returning 06/06/2016)
by Leon Woo 26 May '16

26 May '16

I am out of the office until 06/06/2016. For standard requests within the scope of EMG PWM Berlin, please write to EMG PWM Berlin@DBEMEA. For non-standard requests, please cc Hien Pham-Thu. Note: This is an automated response to your message "NumPy-Discussion Digest, Vol 116, Issue 39" sent on 26.05.2016 14:00:01. This is the only notification you will receive while this person is away.-- Informationen (einschließlich Pflichtangaben) zu einzelnen, innerhalb der EU tätigen Gesellschaften und Zweigniederlassungen des Konzerns Deutsche Bank finden Sie unter http://www.deutsche-bank.de/de/content/pflichtangaben.htm. Diese E-Mail enthält vertrauliche und/ oder rechtlich geschützte Informationen. Wenn Sie nicht der richtige Adressat sind oder diese E-Mail irrtümlich erhalten haben, informieren Sie bitte sofort den Absender und vernichten Sie diese E-Mail. Das unerlaubte Kopieren sowie die unbefugte Weitergabe dieser E-Mail ist nicht gestattet. Please refer to http://www.db.com/en/content/eu_disclosures.htm for information (including mandatory corporate particulars) on selected Deutsche Bank branches and group companies registered or incorporated in the European Union. This e-mail may contain confidential and/or privileged information. If you are not the intended recipient (or have received this e-mail in error) please notify the sender immediately and delete this e-mail. Any unauthorized copying, disclosure or distribution of the material in this e-mail is strictly forbidden.

1 0

Re: [Numpy-discussion] Numpy arrays shareable among related processes (PR #7533)
by Elliot Hallmark 24 May '16

24 May '16

Strula, this sounds brilliant! To be clear, you're talking about serializing the numpy array and reconstructing it in a way that's faster than pickle? Or using shared memory and signaling array creation around that shared memory rather than using pickle? For what it's worth, I have used shared memory with numpy arrays as IPC (no queue), with one process writing to it and one process reading from it, and liked it. Your point #5 did not apply because I was reusing the shared memory. Do you have a public repo where you are working on this? Thanks! Elliot On Wed, May 11, 2016 at 3:29 AM, Sturla Molden <sturla.molden(a)gmail.com> wrote: > I did some work on this some years ago. I have more or less concluded that > it was a waste of effort. But first let me explain what the suggested > approach do not work. As it uses memory mapping to create shared memory > (i.e. shared segments are not named), they must be created ahead of > spawning processes. But if you really want this to work smoothly, you want > named shared memory (Sys V IPC or posix shm_open), so that shared arrays > can be created in the spawned processes and passed back. > > Now for the reason I don't care about shared memory arrays anymore, and > what I am currently working on instead: > > 1. I have come across very few cases where threaded code cannot be used in > numerical computing. In fact, multithreading nearly always happens in the > code where I write pure C or Fortran anyway. Most often it happens in > library code that are already multithreaded (Intel MKL, Apple Accelerate > Framework, OpenBLAS, etc.), which means using it requires no extra effort > from my side. A multithreaded LAPACK library is not less multithreaded if I > call it from Python. > > 2. Getting shared memory right can be difficult because of hierarchical > memory and false sharing. You might not see it if you only have a multicore > CPU with a shared cache. But your code might not scale up on computers with > more than one physical processor. False sharing acts like the GIL, except > it happens in hardware and affects your C code invisibly without any > explicit locking you can pinpoint. This is also why MPI code tends to scale > much better than OpenMP code. If nothing is shared there will be no false > sharing. > > 3. Raw C level IPC is cheap – very, very cheap. Even if you use pipes or > sockets instead of shared memory it is cheap. There are very few cases > where the IPC tends to be a bottleneck. > > 4. The reason IPC appears expensive with NumPy is because multiprocessing > pickles the arrays. It is pickle that is slow, not the IPC. Some would say > that the pickle overhead is an integral part of the IPC ovearhead, but i > will argue that it is not. The slowness of pickle is a separate problem > alltogether. > > 5. Share memory does not improve on the pickle overhead because also NumPy > arrays with shared memory must be pickled. Multiprocessing can bypass > pickling the RawArray object, but the rest of the NumPy array is pickled. > Using shared memory arrays have no speed advantage over normal NumPy arrays > when we use multiprocessing. > > 6. It is much easier to write concurrent code that uses queues for message > passing than anything else. That is why using a Queue object has been the > popular Pythonic approach to both multitreading and multiprocessing. I > would like this to continue. > > I am therefore focusing my effort on the multiprocessing.Queue object. If > you understand the six points I listed you will see where this is going: > What we really need is a specialized queue that has knowledge about NumPy > arrays and can bypass pickle. I am therefore focusing my efforts on > creating a NumPy aware queue object. > > We are not doing the users a favor by encouraging the use of shared memory > arrays. They help with nothing. > > > Sturla Molden > >

5 8

Proposal: numpy.random.random_seed
by Stephan Hoyer 23 May '16

23 May '16

I have recently encountered several use cases for randomly generate random number seeds: 1. When writing a library of stochastic functions that take a seed as an input argument, and some of these functions call multiple other such stochastic functions. Dask is one such example [1]. 2. When a library needs to produce results that are reproducible after calling numpy.random.seed, but that do not want to use the functions in numpy.random directly. This came up recently in a pandas pull request [2], because we want to allow using RandomState objects as an alternative to global state in numpy.random. A major advantage of this approach is that it provides an obvious alternative to reusing the private numpy.random._mtrand [3]. The implementation of this function (and the corresponding method on RandomState) is almost trivial, and I've already written such a utility for my code: def random_seed(): # numpy.random uses uint32 seeds np.random.randint(2 ** 32) The advantage of adding a new method is that it avoids the need for explanation by making the intent of code using this pattern obvious. So I think it is a good candidate for inclusion in numpy.random. Any opinions? [1] https://github.com/dask/dask/blob/e0b246221957c4bd618e57246f3a7ccc8863c494/… [2] https://github.com/pydata/pandas/pull/13161 [3] On a side note, if there's no longer a good reason to keep this object private, perhaps we should expose it in our public API. It would certainly be useful -- scikit-learn is already using it (see links in the pandas PR above).

7 20

A numpy based Entity-Component-System
by Elliot Hallmark 23 May '16

23 May '16

I have a Data Oriented programing library I'm writing that uses the Entity-Component-System model. https://github.com/Permafacture/data-oriented-pyglet I have initially called it Numpy-ECS but I don't know if that name is okay. The numpy license says: Neither the name of the NumPy Developers nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. but doesn't say anything about using the name of NumPy. Is this okay? For anyone interested in what this project is, I suggest looking at the project README. Briefly, a Component is the Data Oriented version of what Object Oriented calls an attribute (a Numpy array of values that represents an attribute that Entities might have). An Entity is the Data Oriented version of an object instance (something subscribed to some subset of Components, with an ID that can be used to access each instance's values in those components). A System is like an object method (a function that operates on the Components of all Entities that have that Component). This is something that one would use in a game engine or simulation. It allows you to have instances, but where the logic can be bulk applied to all involved instances by operations on slices of Numpy arrays. In an example shown in this video <https://vimeo.com/66736654>, I went from 50 FPS in an object oriented approach to using ufuncs operating on shared memory in multiprocessing and getting 180 space steps per second on the physics and 500-600 FPS in the rendering. Thanks for reading, Elliot

3 2

Behavior of .reduceat()
by Jaime Fernández del Río 23 May '16

23 May '16

Two of the oldest issues in the tracker (#834 <https://github.com/numpy/numpy/issues/834> and #835 <https://github.com/numpy/numpy/issues/835>) are about how .reduceat() handles its indices parameter. I have been taking a look at the source code, and it would be relatively easy to modify, the hardest part being to figure out what the exact behavior should be. Current behavior is that np.ufunc.reduceat(x, ind) returns [np.ufunc.reduce(a[ind[i]:ind[i+1]] for i in range(len(ind))] with a couple of caveats: 1. if ind[i] >= ind[i+1], then a[ind[i]] is returned, rather than a reduction over an empty slice. 2. an index of len(ind) is appended to the indices argument, to be used as the endpoint of the last slice to reduce over. 3. aside from this last case, the indices are required to be strictly inbounds, 0 <= index < len(x), or an error is raised The proposed new behavior, with some optional behaviors, would be: 1. if ind[i] >= ind[i+1], then a reduction over an empty slice, i.e. the ufunc identity, is returned. This includes raising an error if the ufunc does not have an identity, e.g. np.minimum. 2. to fully support the "reduction over slices" idea, some form of out of bounds indices should be allowed. This could mean either that: 1. only index = len(x) is allowed without raising an error, to allow computing the full reduction anywhere, not just as the last entry of the return, or 2. allow any index in -len(x) <= index <= len(x), with the usual meaning given to negative values, or 3. any index is allowed, with reduction results clipped to existing values (and the usual meaning for negative values). 3. Regarding the appending of that last index of len(ind) to indices, we could: 1. keep appending it, or 2. never append it, since you can now request it without an error being raised, or 3. only append it if the last index is smaller than len(x). My thoughts on the options: - The minimal, more conservative approach would go with 2.1 and 3.1. And of course 1, if we don't implement that none of this makes sense. - I kind of think 2.2 or even 2.3 are a nice enhancement that shouldn't break too much stuff. - 3.2 I'm not sure about, probably hurts more than it helps at this point, although in a brand new design you probably would either not append the last index or also prepend a zero, as in np.split. - And 3.3 seems too magical, probably not a good idea, only listed it for completeness. Any other thoughts or votes on what, if anything, should we implement, and what the deprecation of current behavior should look like? Jaime -- (\__/) ( O.o) ( > <) Este es Conejo. Copia a Conejo en tu firma y ayúdale en sus planes de dominación mundial.

3 2

FFT and reconstruct
by Vasco Gervasi 21 May '16

21 May '16

Hi all, I am trying to understand how FFT work, so I wrote the attached script. The idea is to extract amplitude and phase of a signal and then reconstruct using amplitude and phase information. As you can see, I create some cosine curve on the interval t0-t1. Let's start with t0=1.0 and t1=3.0 and consider just y['1'] = cos(1.0*omega*t), the signal is. [image: Immagine incorporata 1] The amplitude and phase for each order are: [image: Immagine incorporata 2] But if I try to reconstruct the signal using amplitude and phase: [image: Immagine incorporata 3] So as you can see there is a shifting of 180 deg. Now let's consider another case, t0=2 and t1=3, the signal is y['Signal'] = 1.0*cos(1.0*omega*t) + 2.0*cos(2.0*omega*t) + 3.0*cos(3.0*omega*t + pi/4) + 4.0*cos(4.0*omega*t) + 5.0*cos(5.0*omega*t) + 1.0 The reconstructed signal is very similar to the initial one: [image: Immagine incorporata 4] but is not exactly the same. Any advice? Thanks

2 3