Mailman 3 January 2010 - NumPy-Discussion

Module Index for numpy?
by Wayne Watson 25 Jan '10

25 Jan '10

I was just looking at the (Win) Python documentation via the Help on IDLE, and a Global Module Index. Does anything like that exist for numpy, matplotlib, scipy? -- Wayne Watson (Watson Adventures, Prop., Nevada City, CA) (121.015 Deg. W, 39.262 Deg. N) GMT-8 hr std. time) Obz Site: 39° 15' 7" N, 121° 2' 32" W, 2700 feet "I was thinking about how people seem to read the Bible a whole lot more as they get older; then it dawned on me . . they're cramming for their final exam." -- George Carlin Web Page: <www.speckledwithstars.net/>

5 5

Proposed fix for MKL and dynamic loading
by Dag Sverre Seljebotn 25 Jan '10

25 Jan '10

(Apologies if this has been fixed in trunk; I base this on 1.4.0 and no related comments of MKL on the mailing list) I finally got the latest version of MKL working. What appears to have changed is that the MKL shared libraries will themselves dynamically load different other libraries, depending on the detected CPU. This is in some ways great news for me, because it means I can avoid worrying about miscompiles when compiling one single version of NumPy/SciPy to use for our heterogenous cluster. So I'd rather *not* link statically [1]. Anyway, after modifying site.cfg [2], things almost work, but not quite. The problem is that Python by default imports shared libs using RTLD_LOCAL. With this patch to NumPy it does: Change in numpy/linalg/linalg.py: from numpy.linalg import lapack_lite to: try: import sys import ctypes _old_rtld = sys.getdlopenflags() sys.setdlopenflags(_old_rtld|ctypes.RTLD_GLOBAL) from numpy.linalg import lapack_lite finally: sys.setdlopenflags(_old_rtld) del sys; del ctypes; del _old_rtld Questions: a) Should I submit a patch? b) Negative consequences? Perhaps another Python module can now not load a different BLAS implementation? (That still seems better than not being able to use MKL IMO). c) Should this only be enabled by a flag somewhere? Where? Or can one just do it regardless of BLAS? d) Do I need a "if hasattr" for Windows, or will Windows just ignore it, or does this apply to Windows too? [1] BTW, I could not figure out how to link statically if I wanted -- is "search_static_first = 1" supposed to work? Perhaps MKL will insist on loading some parts dynamically even then *shrug*. Dag Sverre

3 9

Interact with matplotlib in Sage
by Gökhan Sever 24 Jan '10

24 Jan '10

Hello, I have thought of this might interesting to share. Register at www.sagenb.org or try on your local Sage-notebook and using the following code: # Simple example demonstrating how to interact with matplotlib directly. # Comment plt.clf() to get the plots overlay in each update. # Gokhan Sever & Harald Schilly (2010-01-24) from scipy import stats import numpy as np import matplotlib.pyplot as plt @interact def plot_norm(loc=(0,(0,10)), scale=(1,(1,10))): rv = stats.norm(loc, scale) x = np.linspace(-10,10,1000) plt.plot(x,rv.pdf(x)) plt.grid(True) plt.savefig('plt.png') plt.clf() A very easy to use example, also well-suited for learning and demonstration purposes. Posted at: http://wiki.sagemath.org/interact/graphics#Interactwithmatplotlib Have fun ;) -- Gökhan

1 0

fast duplicate of array
by Alan G Isaac 24 Jan '10

24 Jan '10

Suppose x and y are conformable 2d arrays. I now want x to become a duplicate of y. I could create a new array: x = y.copy() or I could assign the values of y to x: x[:,:] = y As expected the latter is faster (no array creation). Are there better ways? Thanks, Alan Isaac

4 7

random v. random_state in RandomState
by Alan G Isaac 23 Jan '10

23 Jan '10

As I understand it, numpy.random provides the function ``random`` as an alias for ``random_state``. Might this be moved into numpy.random.mtrand.RandomState, for interface consistency? Right now if I start with ``prng = np.random.RandomState(seed=myseed)`` I cannot use ``prng.random`` as it does not exist. Note that I can use the convenience function ``prng.rand``, and its documentation makes reference to ``random``! Alan Isaac

1 0

fromfile() for reading text (one more time!)
by Christopher Barker 23 Jan '10

23 Jan '10

Hi folks, I'm taking a look once again at fromfile() for reading text files. I often have the need to read a LOT of numbers form a text file, and it can actually be pretty darn slow do i the normal python way: for line in file: data = map(float, line.strip().split()) or various other versions that are similar. It really does take longer to read the text, split it up, convert to a number, then put that number into a numpy array, than it does to simply read it straight into the array. However, as it stands, fromfile() turn out to be next to useless for anything but whitespace separated text. Full set of ideas here: http://projects.scipy.org/numpy/ticket/909 However, for the moment, I'm digging into the code to address a particular problem -- reading files like this: 123, 65.6, 789 23, 3.2, 34 ... That is comma (or whatever) separated text -- pretty common stuff. The problem with the current code is that you can't read more than one line at time with fromfile: a = np.fromfile(infile, sep=",") will read until it doesn't find a comma, and thus only one line, as there is no comma after each line. As this is a really typical case, I think it should be supported. Here is the question: The work of finding the separator is done in: multiarray/ctors.c: fromfile_skip_separator() It looks like it wouldn't be too hard to add some code in there to look for a newline, and consider that a valid separator. However, that would break backward compatibility. So maybe a flag could be passed in, saying you wanted to support newlines. The problem is that flag would have to get passed all the way through to this function (and also for fromstring). I also notice that it supports separators of arbitrary length, which I wonder how useful that is. But it also does odd things with spaces embedded in the separator: ", $ #" matches all of: ",$#" ", $#" ",$ #" Is it worth trying to fix that? In the longer term, it would be really nice to support comments as well, tough that would require more of a re-factoring of the code, I think (though maybe not -- I suppose a call to fromfile_skip_separator() could look for a comment character, then if it found one, skip to where the comment ends -- hmmm. thanks for any feedback, -Chris -- Christopher Barker, Ph.D. Oceanographer Emergency Response Division NOAA/NOS/OR&R (206) 526-6959 voice 7600 Sand Point Way NE (206) 526-6329 fax Seattle, WA 98115 (206) 526-6317 main reception Chris.Barker(a)noaa.gov

9 30

generalized ufunc problem
by David Warde-Farley 23 Jan '10

23 Jan '10

I decided to take a crack at adding a generalized ufunc for logsumexp, i.e. collapsed an array along the last dimension by subtracting the maximum element E along that dimension, taking the exponential, adding, and then adding back E. Functionally the same logaddexp.reduce() but presumably faster and less prone to error accumulation. I added the following to umath_tests.c.src and got everything to compile, but for some reason it doesn't give me the behaviour I'm looking for. I was expecting a (500, 50) array to be collapsed down to a (500,) array. Is that not what the signature calls for? Thanks, David char *logsumexp_signature = "(i)->()"; /**begin repeat #TYPE=LONG,DOUBLE# #typ=npy_long, npy_double# #EXPFUN=expl, exp# #LOGFUN=logl, log# */ /* * This implements the function * out[n] = sum_i { in1[n, i] * in2[n, i] }. */ static void @TYPE@_logsumexp(char **args, intp *dimensions, intp *steps, void *NPY_UNUSED(func)) { INIT_OUTER_LOOP_3 intp di = dimensions[0]; intp i; intp is1=steps[0]; BEGIN_OUTER_LOOP_3 char *ip1=args[0], *op=args[1]; @typ@ max = (*(@typ@ *)ip1); @typ@ sum = 0; for (i = 0; i < di; i++) { max = max < (*(@typ@ *)ip1) ? (*(@typ@ *)ip1) : max; ip1 += is1; } ip1 = args[0]; for (i = 0; i < di; i++) { sum += @EXPFUN@((*(@typ@ *)ip1) - max); ip1 += is1; } *(@typ@ *)op = @LOGFUN@(sum + max); END_OUTER_LOOP } /**end repeat**/ static PyUFuncGenericFunction logsumexp_functions[] = { LONG_logsumexp, DOUBLE_logsumexp }; static void * logsumexp_data[] = { (void *)NULL, (void *)NULL }; static char logsumexp_signatures[] = { PyArray_LONG, PyArray_LONG, PyArray_DOUBLE, PyArray_DOUBLE }; /* and inside addUfuncs() */ ... f = PyUFunc_FromFuncAndData(logsumexp_functions, logsumexp_data, logsumexp_signatures, 2, 1, 1, PyUFunc_None, "logsumexp", "inner1d with a weight argument \\n""\" \"\n"" \\ \"(i)->()\\\" \\n""", 0); PyDict_SetItemString(dictionary, "logsumexp", f); Py_DECREF(f); ....

2 2

void * abuse
by Charles R Harris 22 Jan '10

22 Jan '10

Example: typedef void (PyArray_VectorUnaryFunc)(void *, void *, npy_intp, void *, void *) The actual functions that implement this type don't match the prototype, an unfortunate fact that is papered over by casting said functions to the type to avoid warnings about initialization from incompatible pointer types. Theoretically real problems could result from this, in practice the main effect is that type checking by the compiler doesn't take place leading to undisciplined implementions. The current prototype isn't much better than a simple typedef void (PyArray_VectorUnaryFunc)() So... this can be fixed but the fix could cause some current code in the wild to result in some warnings next time it is compiled. I doubt there will be any binary incompatibility with existing compiled code. Should I fix it? It is really ugly as it stands. Chuck

1 0

np.testing.assert_equal
by Keith Goodman 22 Jan '10

22 Jan '10

Should np.testing.assert_equal(np.array(1), 1) raise an AssertionError? (It currently doesn't.) The use case I have in mind is that scipy.stats.nanamedian incorrectly returns np.array(1.0) for the median of a 1d array while np.median correctly returns 1.0. It would be handy if the assert statement caught the difference.

2 1

Broadcasting and indexing
by Thomas Robitaille 22 Jan '10

22 Jan '10

Hello, I'm trying to understand how array broadcasting can be used for indexing. In the following, I use the term 'row' to refer to the first dimension of a 2D array, and 'column' to the second, just because that's how numpy prints them out. If I consider the following example: >>> a = np.random.random((4,5)) >>> b = np.random.random((5,)) >>> a + b array([[ 1.45499556, 0.60633959, 0.48236157, 1.55357393, 1.4339261 ], [ 1.28614593, 1.11265001, 0.63308615, 1.28904227, 1.34070499], [ 1.26988279, 0.84683018, 0.98959466, 0.76388223, 0.79273084], [ 1.27859505, 0.9721984 , 1.02725009, 1.38852061, 1.56065028]]) I understand how this works, because it works as expected as described in http://docs.scipy.org/doc/numpy/reference/ufuncs.html#broadcasting So b gets broadcast to shape (1,5), then because the first dimension is 1, the operation is applied to all rows. Now I am trying to apply this to array indexing. So for example, I want to set specific columns, indicated by a boolean array, to zero, but the following fails: >>> c = np.array([1,0,1,0,1], dtype=bool) >>> a[c] = 0 Traceback (most recent call last): File "<stdin>", line 1, in <module> IndexError: index (4) out of range (0<=index<3) in dimension 0 However, if I try reducing the size of c to 4, then it works, and sets rows, not columns, equal to zero >>> c = np.array([1,0,1,0], dtype=bool) >>> a[c] = 0 >>> a array([[ 0. , 0. , 0. , 0. , 0. ], [ 0.41526315, 0.7425491 , 0.39872546, 0.56141914, 0.69795153], [ 0. , 0. , 0. , 0. , 0. ], [ 0.40771227, 0.60209749, 0.7928894 , 0.66089748, 0.91789682]]) But I would have thought that the indexing array would have been broadcast in the same way as for a sum, i.e. c would be broadcast to have dimensions (1,5) and then would have been able to set certain columns in all rows to zero. Why is it that for indexing, the broadcasting seems to happen in a different way than when performing operations like additions or multiplications? For background info, I'm trying to write a routine which performs a set of operations on an n-d array, where n is not known in advance, with a 1D array, so I can use broadcasting rules for most operations without knowing the dimensionality of the n-d array, but now that I need to perform indexing, and the convention seems to change, this is a real issue. Thanks in advance for any advice, Thomas

3 2