Mailman 3 February 2014 - NumPy-Discussion

argsort speed
by josef.pktd＠gmail.com Feb. 24, 2014

Feb. 24, 2014

currently using numpy 1.6.1 What's the fastest argsort for a 1d array with around 28 Million elements, roughly uniformly distributed, random order? Is there a reason that np.argsort is almost 3 times slower than np.sort? I'm doing semi-systematic timing for a stats(models) algorithm. Josef

7 16

Proposal: Chaining np.dot with mdot helper function
by Stefan Otte Feb. 24, 2014

Feb. 24, 2014

Hey guys, I wrote myself a little helper function `mdot` which chains np.dot for multiple arrays. So I can write mdot(A, B, C, D, E) instead of these A.dot(B).dot(C).dot(D).dot(E) np.dot(np.dot(np.dot(np.dot(A, B), C), D), E) I know you can use `numpy.matrix` to get nicer formulas. However, most numpy/scipy function return arrays instead of numpy.matrix. Therefore, sometimes you actually use array multiplication when you think you use matrix multiplication. `mdot` is a simple way to avoid using numpy.matrix but to improve the readability. What do you think? Is this useful and worthy to integrate in numpy? I already created an issuer for this: https://github.com/numpy/numpy/issues/4311 jaimefrio also suggested to do some reordering of the arrays to minimize computation: https://github.com/numpy/numpy/issues/4311#issuecomment-35295857 Best, Stefan

6 13

How exactly ought 'dot' to work?
by Nathaniel Smith Feb. 24, 2014

Feb. 24, 2014

Hi all, Currently numpy's 'dot' acts a bit weird for ndim>2 or ndim<1. In practice this doesn't usually matter much, because these are very rarely used. But, I would like to nail down the behaviour so we can say something precise in the matrix multiplication PEP. So here's one proposal. # CURRENT: dot(0d, any) -> scalar multiplication dot(any, 0d) -> scalar multiplication dot(1d, 1d) -> inner product dot(2d, 1d) -> treat 1d as column matrix, matrix-multiply, then discard added axis dot(1d, 2d) -> treat 1d as row matrix, matrix-multiply, then discard added axis dot(2d, 2d) -> matrix multiply dot(2-or-more d, 2-or-more d) -> a complicated outer product thing: Specifically, if the inputs have shapes (r, n, m), (s, m, k), then numpy returns an array with shape (r, s, n, k), created like: for i in range(r): for j in range(s): output[i, j, :, :] = np.dot(input1[i, :, :], input2[j, :, :]) # PROPOSED: General rule: given dot on shape1, shape2, we try to match these shapes against two templates like (..., n?, m) and (..., m, k?) where ... indicates zero or more dimensions, and ? indicates an optional axis. ? axes are always matched before ... axes, so for an input with ndim>=2, the ? axis is always matched. An unmatched ? axis is treated as having size 1. Next, the ... axes are broadcast against each other in the usual way (prepending 1s to make lengths the same, requiring corresponding entries to either match or have the value 1). And then the actual computations are performed using the usual broadcasting rules. Finally, we return an output with shape (..., n?, k?). Here "..." indicates the result of broadcasting the input ...'s against each other. And, n? and k? mean: "either the value taken from the input shape, if the corresponding entry was matched -- but if no match was made, then we leave this entry out." The idea is that just as a column vector on the right is "m x 1", a 1d vector on the right is treated as "m x <nothing>". For purposes of actually computing the product, <nothing> acts like 1, as mentioned above. But it makes a difference in what we return: in each of these cases we copy the input shape into the output, so we can get an output with shape (n, <nothing>), or (<nothing>, k), or (<nothing>, <nothing>), which work out to be (n,), (k,) and (), respectively. This gives a (somewhat) intuitive principle for why dot(1d, 1d), dot(1d, 2d), dot(2d, 1d) are handled the way they are, and a general template for extending this behaviour to other operations like gufunc 'solve'. Anyway, the end result of this is that the PROPOSED behaviour differs from the current behaviour in the following ways: - passing 0d arrays to 'dot' becomes an error. (This in particular is an important thing to know, because if core Python adds an operator for 'dot', then we must decide what it should do for Python scalars, which are logically 0d.) - ndim>2 arrays are now handled by aligning and broadcasting the extra axes, instead of taking an outer product. So dot((r, m, n), (r, n, k)) returns (r, m, k), not (r, r, m, k). Comments? -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org

6 9

np.savetxt() default format
by Daniele Nicolodi Feb. 24, 2014

Feb. 24, 2014

Hello, I've noticed that numpy default format for saving data in ascii representation with np.savetxt() is "%.18e". Given that the default data type for numpy is double and that the resolution of doubles is 15 decimal digits, what's the reason of the the additional three digits? The three additional digits are definitely not an issue, but I would like to understand the reason why they are there. Thanks. Best, Daniele

2 1

GSOC
by Charles R Harris Feb. 23, 2014

Feb. 23, 2014

Thought I'd forward this to the lists in case we need to do something. Hi everyone, > > Just a friendly reminder that applications for mentoring organizations > close in about 24 hours. Please get your applications in soon, we will > not accept late applications for any reason! > > Thanks, > Carol Chuck

6 10

OpenBLAS on Mac
by Jurgen Van Gael Feb. 23, 2014

Feb. 23, 2014

Hi All, I run Mac OS X 10.9.1 and was trying to get OpenBLAS working for numpy. I've downloaded the OpenBLAS source and compiled it (thanks to Olivier Grisel). I installed everything to /usr/local/lib (I believe): e.g. "ll /usr/local/lib/ | grep openblas" lrwxr-xr-x 1 37B 10 Feb 14:51 libopenblas.a@ -> libopenblas_sandybridgep-r0.2.9.rc1.a lrwxr-xr-x 1 56B 10 Feb 14:51 libopenblas.dylib@ -> /usr/local/lib/libopenblas_sandybridgep-r0.2.9.rc1.dylib -rw-r--r-- 1 18M 7 Feb 16:02 libopenblas_sandybridgep-r0.2.9.rc1.a -rwxr-xr-x 1 12M 10 Feb 14:51 libopenblas_sandybridgep-r0.2.9.rc1.dylib* Then I download the numpy sources and add a site.cfg with the only three lines uncommented being: [openblas] libraries = openblas library_dirs = /usr/local/lib include_dirs = /usr/local/include When I run python setup.py config I get message that say that openblas_info has been found. I then run setup build and setup install (into a virtualenv). In the virtualenv, when I then check what _dotblas.so is linked to, I keep getting that it is linked to Accelerate. E.g. otool -L .../numpy/core/_dotblas.so => /System/Library/Frameworks/Accelerate.framework/Versions/A/Accelerate Any suggestions on getting my numpy working with OpenBLAS? Thanks, Jurgen

5 11

Re: [Numpy-discussion] How exactly ought 'dot' to work?
by Jaime Fernández del Río Feb. 22, 2014

Feb. 22, 2014

On Feb 22, 2014 2:03 PM, "Nathaniel Smith" <njs(a)pobox.com> wrote: > > Hi all, > > Currently numpy's 'dot' acts a bit weird for ndim>2 or ndim<1. In > practice this doesn't usually matter much, because these are very > rarely used. But, I would like to nail down the behaviour so we can > say something precise in the matrix multiplication PEP. So here's one > proposal. > > # CURRENT: > > dot(0d, any) -> scalar multiplication > dot(any, 0d) -> scalar multiplication > dot(1d, 1d) -> inner product > dot(2d, 1d) -> treat 1d as column matrix, matrix-multiply, then > discard added axis > dot(1d, 2d) -> treat 1d as row matrix, matrix-multiply, then discard added axis > dot(2d, 2d) -> matrix multiply > dot(2-or-more d, 2-or-more d) -> a complicated outer product thing: > Specifically, if the inputs have shapes (r, n, m), (s, m, k), then > numpy returns an array with shape (r, s, n, k), created like: > for i in range(r): > for j in range(s): > output[i, j, :, :] = np.dot(input1[i, :, :], input2[j, :, :]) > > # PROPOSED: > > General rule: given dot on shape1, shape2, we try to match these > shapes against two templates like > (..., n?, m) and (..., m, k?) > where ... indicates zero or more dimensions, and ? indicates an > optional axis. ? axes are always matched before ... axes, so for an > input with ndim>=2, the ? axis is always matched. An unmatched ? axis > is treated as having size 1. > > Next, the ... axes are broadcast against each other in the usual way > (prepending 1s to make lengths the same, requiring corresponding > entries to either match or have the value 1). And then the actual > computations are performed using the usual broadcasting rules. > > Finally, we return an output with shape (..., n?, k?). Here "..." > indicates the result of broadcasting the input ...'s against each > other. And, n? and k? mean: "either the value taken from the input > shape, if the corresponding entry was matched -- but if no match was > made, then we leave this entry out." The idea is that just as a column > vector on the right is "m x 1", a 1d vector on the right is treated as > "m x <nothing>". For purposes of actually computing the product, > <nothing> acts like 1, as mentioned above. But it makes a difference > in what we return: in each of these cases we copy the input shape into > the output, so we can get an output with shape (n, <nothing>), or > (<nothing>, k), or (<nothing>, <nothing>), which work out to be (n,), > (k,) and (), respectively. This gives a (somewhat) intuitive principle > for why dot(1d, 1d), dot(1d, 2d), dot(2d, 1d) are handled the way they > are, and a general template for extending this behaviour to other > operations like gufunc 'solve'. > > Anyway, the end result of this is that the PROPOSED behaviour differs > from the current behaviour in the following ways: > - passing 0d arrays to 'dot' becomes an error. (This in particular is > an important thing to know, because if core Python adds an operator > for 'dot', then we must decide what it should do for Python scalars, > which are logically 0d.) > - ndim>2 arrays are now handled by aligning and broadcasting the extra > axes, instead of taking an outer product. So dot((r, m, n), (r, n, k)) > returns (r, m, k), not (r, r, m, k). > > Comments? The proposed behavior for ndim > 2 is what matrix_multiply (is it still in umath_tests?) does. The nice thing of the proposed new behavior is that the old behavior is easy to reproduce by fooling a little around with the shape of the first argument, while the opposite is not true. Jaime > > -- > Nathaniel J. Smith > Postdoctoral researcher - Informatics - University of Edinburgh > http://vorpus.org > _______________________________________________ > NumPy-Discussion mailing list > NumPy-Discussion(a)scipy.org > http://mail.scipy.org/mailman/listinfo/numpy-discussion

1 0

A PEP for adding infix matrix multiply to Python
by Nathaniel Smith Feb. 22, 2014

Feb. 22, 2014

[Apologies for wide distribution -- please direct followups to either the github PR linked below, or else numpy-discussion(a)scipy.org] After the numpy-discussion thread about np.matrix a week or so back, I got curious and read the old PEPs that attempted to add better matrix/elementwise operators to Python. http://legacy.python.org/dev/peps/pep-0211/ http://legacy.python.org/dev/peps/pep-0225/ And I was a bit surprised -- if I were BDFL I probably would have rejected these PEPs too. One is actually a proposal to make itertools.product into an infix operator, which no-one would consider seriously on its own merits. And the other adds a whole pile of weirdly spelled new operators with no clear idea about what they should do. But it seems to me that at this point, with the benefit of multiple years more experience, we know much better what we want -- basically, just a nice clean infix op for matrix multiplication. And that just asking for this directly, and explaining clearly why we want it, is something that hasn't been tried. So maybe we should try and see what happens. As a starting point for discussion, I wrote a draft. It can be read and commented on here: https://github.com/numpy/numpy/pull/4351 It's important that if we're going to do this at all, we do it right, and that means being able to honestly say that this document represents our consensus when going to python-dev. So if you think you might object please do so now :-) -n -- Nathaniel J. Smith Postdoctoral researcher - Informatics - University of Edinburgh http://vorpus.org

1 0

Header files in windows installers
by Matt Newell Feb. 21, 2014

Feb. 21, 2014

I have a small c++ extension used to feed a 1d numpy array into a QPainterPath. Very simple just using PyArray_[Check|FLAGS|SIZE|DATA]. I developed it on debian which was of course very straightforward, but now I need to deploy on windows, which is of course where the fun always begins. I was unpleasantly surprised to discover that the numpy installers for window do not have an option to install the include files, which appear to be all that is needed to compile extensions that use numpy's C-api. I have actually managed to copy my set of include files from linux and with a few modifications to _numpyconfig.h got my extension compiled and working. Is there any chance that the includes could be included in future releases? Thanks, Matt Newell

3 3

Document server error.
by Charles R Harris Feb. 21, 2014

Feb. 21, 2014

>From issue #1951 <https://github.com/numpy/numpy/issues/1951>: The following URL shows an 500 internal server error: > http://docs.scipy.org/doc/numpy/reference/generated/numpy.var.html > Can someone with access to the server take a look? TIA, Chuck

2 5