Mailman 3 March 2007 - NumPy-Discussion

masked array with mask an instance of numpy.memmap
by Glen W. Mabey March 27, 2007

March 27, 2007

I've been using numpy.memmap to access large arrays very nicely, but now I desire to instantiate a masked array that uses memmap objects. Using a memmap instance for the data works just fine and as expected. However, using a memmap for the mask value seems to load the entire object into memory, even if it is already of type bool. How hard would it be to change that behavior? Also, it seems that median() is not implemented on masked arrays? Is there some method for a masked array like ravel() but that only includes the non-masked values? Thank you, Glen Mabey

2 6

A unique function...
by Pierre GM March 27, 2007

March 27, 2007

All, I could swear that I ran once into a numpy (or scipy) function that output the unique values of a 1d ndarray along with the number of occurences of each element. If I'm not completely mistaken, it 's a private function initially in Fortran. Does this ring a bell to anyone ? Where could I find this function ?

2 2

Tuning sparse stuff in NumPy
by David Koch March 27, 2007

March 27, 2007

Hi, I ran some tests on the very same matrices in Matlab/Numpy and it seems that for sparse matrix multipilcation to be faster than dense multiplication - the degree of sparsity has to be much higher in Numpy than in Matlab. Is there anything I can tune in the underlying routines? I need good performance at about 0.02 +- 0.01 nnz elements. Thanks, /David Matrix size: 10k * 1k, times are in seconds: sparsity: 0.01 nonzero elements Matlab dense 3.9, sparse 0.09 Numpy dense 2.9, sparse 15.24 sparsity: increase sparsity to 0.001 nnz Matlab dense 3.9, sparse 0.002 Numpy dense 2.9, sparse 0.27

2 7

defmatrix.py
by Nils Wagner March 27, 2007

March 27, 2007

Hi, Several tests didn't pass due to the recent changes in defmatrix.py. Nils ====================================================================== ERROR: check_matmat (scipy.sparse.tests.test_sparse.test_csc) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/sparse/tests/test_sparse.py", line 142, in check_matmat assert_array_almost_equal((a*bsp).todense(), a*b) File "/usr/local/lib64/python2.5/site-packages/numpy/core/defmatrix.py", line 162, in __mul__ return N.dot(self, other) ValueError: objects are not aligned ====================================================================== ERROR: check_rmatvec (scipy.sparse.tests.test_sparse.test_csc) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/sparse/tests/test_sparse.py", line 107, in check_rmatvec assert_array_almost_equal(row*M, row*M.todense()) File "/usr/local/lib64/python2.5/site-packages/numpy/core/defmatrix.py", line 162, in __mul__ return N.dot(self, other) ValueError: objects are not aligned ====================================================================== ERROR: check_matmat (scipy.sparse.tests.test_sparse.test_csr) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib64/python2.5/site-packages/scipy/sparse/tests/test_sparse.py", line 142, in check_matmat assert_array_almost_equal((a*bsp).todense(), a*b) File "/usr/local/lib64/python2.5/site-packages/numpy/core/defmatrix.py", line 162, in __mul__ return N.dot(self, other) ValueError: objects are not aligned ====================================================================== ERROR: check_rmatvec (scipy.sparse.tests.test_sparse.test_csr) ...

5 15

ATLAS problems
by Justin Bedo March 27, 2007

March 27, 2007

Hi All, I'm trying to compile NumPy using the ATLAS libraries but am having some difficulty with it. I've managed to build numpy such that lapack_lite.so is linked as follows: linux-gate.so.1 => (0xffffe000) libatlas.so => /usr/local/atlas/lib/libatlas.so (0xb79e9000) libptcblas.so => /usr/local/atlas/lib/libptcblas.so (0xb79c7000) libptf77blas.so => /usr/local/atlas/lib/libptf77blas.so (0xb79ac000) libpython2.5.so.1.0 => /usr/lib/libpython2.5.so.1.0 (0xb788f000) libgfortran.so.1 => /usr/lib/libgfortran.so.1 (0xb7817000) libm.so.6 => /lib/libm.so.6 (0xb77f1000) libgcc_s.so.1 => /lib/libgcc_s.so.1 (0xb77e5000) libc.so.6 => /lib/libc.so.6 (0xb76b7000) libpthread.so.0 => /lib/libpthread.so.0 (0xb769f000) libdl.so.2 => /lib/libdl.so.2 (0xb769a000) libutil.so.1 => /lib/libutil.so.1 (0xb7696000) /lib/ld-linux.so.2 (0x80000000) However despite this, the performance when multiplying matrices is extremely slow, and only a single thread is used even though lapack_lite.so is linked with the pt atlas libraries. Does anybody know how to fix this problem? Thanks Justin

2 1

two previously unresolved issues
by Stefan van der Walt March 26, 2007

March 26, 2007

Hi, I just went through my mail archive and found these two minor outstanding issues. Thought I'd ask for comments before the new release: """ From: "Charles R Harris" <charlesr.harris(a)gmail.com> Subject: Re: [Numpy-discussion] Assign NaN, get zero On 11/11/06, Lisandro Dalcin <dalcinl(a)gmail.com> wrote: > > On 11/11/06, Stefan van der Walt <stefan(a)sun.ac.za> wrote: > > NaN (or inf) is a floating point number, so seeing a zero in integer > > representation seems correct: > > > > In [2]: int(N.nan) > > Out[2]: 0L > > > > Just to learn myself: Why int(N.nan) should be 0? Is it C behavior? In [1]: int32(0)/int32(0) Warning: divide by zero encountered in long_scalars Out[1]: 0 In [2]: float32(0)/float32(0) Out[2]: nan In [3]: int(nan) Out[3]: 0L I think it was just a default for numpy. Hmmm, numpy now warns on integer division by zero, didn't used to. Looks like a warning should also be raised when casting nan to integer. It is probably a small bug not to. I also suspect int(nan) should return a normal python zero, not 0L. """ """ From: "Bill Baxter" <wbaxter(a)gmail.com> To: numpy-discussion(a)scipy.org Subject: [Numpy-discussion] linalg.lstsq for complex Is this code from linalg.lstsq for the complex case correct? lapack_routine = lapack_lite.zgelsd lwork = 1 rwork = zeros((lwork,), real_t) work = zeros((lwork,),t) results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, 0, work, -1, rwork, iwork, 0) lwork = int(abs(work[0])) rwork = zeros((lwork,),real_t) a_real = zeros((m,n),real_t) bstar_real = zeros((ldb,n_rhs,),real_t) results = lapack_lite.dgelsd(m, n, n_rhs, a_real, m, bstar_real, ldb, s, rcond, 0, rwork, -1, iwork, 0) lrwork = int(rwork[0]) work = zeros((lwork,), t) rwork = zeros((lrwork,), real_t) results = lapack_routine(m, n, n_rhs, a, m, bstar, ldb, s, rcond, The middle call to dgelsd looks unnecessary to me. At the very least, allocating astar_real and bstar_real shouldn't be necessary since they aren't referenced anywhere else in the lstsq function. The lapack documentation for zgelsd also doesn't mention any need to call dgelsd to compute the size of the work array. """ Cheers Stéfan

1 0

Fixed scalar coercion model in NumPy
by Travis Oliphant March 26, 2007

March 26, 2007

I've finally made the changes to fix the scalar coercion model problems in NumPy 1.0.1 Now, scalar coercion rules only apply when involved types are of the same basic "kind". Thus, array([1,2,3],int8)*10 returns an int8 array but array([1,2,3],int8)*10.0 returns a float64 array. If you want a float32 array you must use array([1,2,3],int8)*float32(10.0) This is actually a behavioral change which is why I asked earlier if we should change it. However, the previous behavior was not documented anywhere and any previous discussion on this point should have been interpreted as the behavior now in SVN. It's also a rare use-case and so should not create too many issues. In running tests with NumPy and SciPy there appear to be three tests in ndimage breaking now. I really do want to get 1.0.2 out the door soon. What still needs to be fixed before then? -Travis

3 3

Re: [Numpy-discussion] New Operators in Python
by jh＠physics.ucf.edu March 26, 2007

March 26, 2007

I strongly agree with the need for different infix operators for matrix and array multiplication. In IDL, I have used both in the same code. IDL uses * for array and # for matrix (hey, # even *looks* like a matrix, very intuitive). It uses ## for right-multiply (why people don't just M-t the variable names I don't know, though maybe people who do more linear algebra than I can say). There's an important justification that hasn't come up yet, and that is advocacy. Languages are good or not good depending on how readable and writeable they are for various purposes. LISP is a technically very powerful language. It was object-oriented a zillion years before C++. It had and has powerful proponents within CS, such as Stallman. But, whenever CS people try to push LISP outside their ivory tower, it gets trashed. The LISP dialect Guile was written as the GNU project's extension language and nobody used it. Now projects like GIMP have two extension languages, Guile and Python. Why? LISP is great for analyzing the flow and structure of computer code, but LISP math does not look like real math, which makes it very difficult and counterintuitive to use for math or math-heavy topics like graphics. It has little syntactic sugar to make common math tasks either quick to write or easy to understand or verify, until you've retrained your brain. So, even though there are tons of technical people exposed to LISP, nobody outside CS adopts it who doesn't have to. On the other hand, Python, IDL, and Matlab are attractive to us mainly because their syntaxes allow us to see the math, understand it on inspection, and verify its correctness. The math we write in these languages looks as much like the math we do on paper as ASCII will allow. (Obviously, we also choose them because we don't like writing loops and declaring variables.) So, whenever I hear someone suggest using a functional syntax for a concept that exists notationally, I cringe. We're alienating a class of users each time we do that. Those are people who will never come to Python. There are extremes to which this argument cannot go - a prime will never be a derivative because quotes are more important - but I think that matrix multiplication is a no-brainer here. We should let the Python developers know we want it now and then follow up with a syntax and implementation. --jh-- Prof. Joseph Harrington Department of Physics MAP 420 University of Central Florida Orlando, FL 32816-2385 (407) 823-3416 voice (407) 823-5112 fax (407) 823-2325 physics office jh(a)physics.ucf.edu (direct) jh(a)alum.mit.edu (permanent forwarding address, will work forever)

8 11

fastest way to do multiplication with diagonal matrices from left or right
by daniel.egloff＠zkb.ch March 26, 2007

March 26, 2007

Dear list what is the fastet way to multiply with a diagonal matrix from left or right and without to build a square matrix from the diagonal. Here it what I am looking for: import numpy as N def diagmult(X, Y): """ Matrix multiplication X*Y where either X or Y is a diagonal matrix. """ if X.ndim == 1 and Y.ndim == 2: R = Y.copy() for i, d in enumerate(X): R[i,:] *= d return R elif X.ndim == 2 and Y.ndim == 1: R = X.copy() for i, d in enumerate(Y): R[:,i] *= d return R elif X.ndim == 1 and Y.ndim == 1: return X*Y else raise ValueError('diagmult dimension mismatch X.ndim = %d, Y.ndim = %d' % (X.ndim, Y.ndim)) Freundliche Grüsse Daniel Egloff Zürcher Kantonalbank Leiter(in) Financial Computing, ZEF Josefstrasse 222, 8005 Zürich Telefon 044 292 45 33, Fax 044 292 45 95 Briefadresse: Postfach, 8010 Zürich, http://www.zkb.ch ___________________________________________________________________ Disclaimer: Diese Mitteilung ist nur fuer die Empfaengerin / den Empfaenger bestimmt. Fuer den Fall, dass sie von nichtberechtigten Personen empfangen wird, bitten wir diese hoeflich, die Mitteilung an die ZKB zurueckzusenden und anschliessend die Mitteilung mit allen Anhaengen sowie allfaellige Kopien zu vernichten bzw. zu loeschen. Der Gebrauch der Information ist verboten. This message is intended only for the named recipient and may contain confidential or privileged information. If you have received it in error, please advise the sender by return e-mail and delete this message and any attachments. Any unauthorised use or dissemination of this information is strictly prohibited.

3 2

header file location upon installation
by Daniel Wheeler March 26, 2007

March 26, 2007

Hi, Should the header files in .../lib/python2.4/site-packages/python/numpy/core/include/ be copied to .../include/python2.4/numpy/ upon installation of numpy? As far as I can tell this is not happening. Just wondering what the default behavior should be. Thanks -- Daniel Wheeler

2 2