Mailman 3 February 2006 - NumPy-Discussion

A case for rank-0 arrays
by Sasha Feb. 24, 2006

Feb. 24, 2006

I have reviewed mailing list discussions of rank-0 arrays vs. scalars and I concluded that the current implementation that contains both is (almost) correct. I will address the "almost" part with a concrete proposal at the end of this post (search for PROPOSALS if you are only interested in the practical part). The main criticism of supporting both scalars and rank-0 arrays is that it is "unpythonic" in the sense that it provides two almost equivalent ways to achieve the same result. However, I am now convinced that this is the case where practicality beats purity. If you take the "one way" rule to it's logical conclusion, you will find that once your language has functions, it does not need numbers or any other data type because they all can be represented by functions (see http://en.wikipedia.org/wiki/Church_numeral). Another example of core python violating the "one way rule" is the presence of scalars and length-1 tuples. In S+, for example, scalars are represented by single element lists. The situation with ndarrays is somewhat similar. A rank-N array is very similar to a function with N arguments, where each argument has a finite domain (i-th domain of a is range(a.shape[i])). A rank-0 array is just a function with no arguments and as such it is quite different from a scalar. Just as a function with no arguments cannot be replaced by a constant in the case when a value returned may change during the run of the program, rank-0 array cannot be replaced by an array scalar because it is mutable. (See http://projects.scipy.org/scipy/numpy/wiki/ZeroRankArray for use cases). Rather than trying to hide rank-0 arrays from the end-user and treat it as an implementation artifact, I believe numpy should emphasize the difference between rank-0 arrays and scalars and have clear rules on when to use what. PROPOSALS ========== Here are three suggestions: 1. Probably the most controversial question is what getitem should return. I believe that most of the confusion comes from the fact that the same syntax implements two different operations: indexing and projection (for the lack of better name). Using the analogy between ndarrays and functions, indexing is just the application of the function to its arguments and projection is the function projection ((f, x) -> lambda (*args): f(x, *args)). The problem is that the same syntax results in different operations depending on the rank of the array. Let >>> x = ones((2,2)) >>> y = ones(2) then x[1] is projection and type(x[1]) is ndarray, but y[1] is indexing and type(y[1]) is int32. Similarly, y[1,...] is indexing, while x[1,...] is projection. I propose to change numpy rules so that if ellipsis is present inside [], the operation is always projection and both y[1,...] and x[1,1,...] return zero-rank arrays. Note that I have previously rejected Francesc's idea that x[...] and x[()] should have different meaning for zero-rank arrays. I was wrong. 2. Another source of ambiguity is the various "reduce" operations such as sum or max. Using the previous example, type(x.sum(axis=0)) is ndarray, but type(y.sum(axis=0)) is int32. I propose two changes: a. Make x.sum(axis) return ndarray unless axis is None, making type(y.sum(axis=0)) is ndarray true in the example. b. Allow axis to be a sequence of ints and make x.sum(axis=range(rank(x))) return rank-0 array according to the rule 2.a above. c. Make x.sum() raise an error for rank-0 arrays and scalars, but allow x.sum(axis=()) to return x. This will make numpy sum consistent with the built-in sum that does not work on scalars. 3. This is a really small change currently >>> empty(()) array(0) but >>> ndarray(()) Traceback (most recent call last): File "<stdin>", line 1, in ? ValueError: need to give a valid shape as the first argument I propose to make shape=() valid in ndarray constructor.

3 8

repmat equivalent?
by Albert Strasheim Feb. 24, 2006

Feb. 24, 2006

Hello all I recently started using NumPy and one function that I am really missing from MATLAB/Octave is repmat. This function is very useful for implementing algorithms as matrix multiplications instead of for loops. Here's my first attempt at repmat for 1d and 2d (with some optimization by Stefan van der Walt): def repmat(a, m, n): if a.ndim == 1: a = array([a]) (origrows, origcols) = a.shape rows = origrows * m cols = origcols * n b = a.reshape(1,a.size).repeat(m, 0).reshape(rows, origcols).repeat(n, 0) return b.reshape(rows, cols) print repmat(array([[1,2],[3,4]]), 2, 3) produces: [[1 2 1 2 1 2] [3 4 3 4 3 4] [1 2 1 2 1 2] [3 4 3 4 3 4]] which is the same as in MATLAB. There are various issues with my function that I don't quite know how to solve: - How to handle scalar inputs (probably need asarray here) - How to handle more than 2 dimensions More than 2 dimensions is tricky, since NumPy and MATLAB don't seem to agree on how more-dimensional data is organised? As such, I don't know what a NumPy user would expect repmat to do with more than 2 dimensions. Here are some test cases that the current repmat should pass, but doesn't: a = repmat(1, 1, 1) assert_equal(a, 1) a = repmat(array([1]), 1, 1) assert_array_equal(a, array([1])) a = repmat(array([1,2]), 2, 3) assert_array_equal(a, array([[1,2,1,2,1,2], [1,2,1,2,1,2]])) a = repmat(array([[1,2],[3,4]]), 2, 3) assert_array_equal(a, array([[1,2,1,2,1,2], [3,4,3,4,3,4], [1,2,1,2,1,2], [3,4,3,4,3,4]])) Any suggestions on how do repmat in NumPy would be appreciated. Regards Albert

6 11

Re: [Numpy-discussion] Matrix times scalar is wacky
by Travis Oliphant Feb. 24, 2006

Feb. 24, 2006

Bill Baxter wrote: > Excellent. > Is working on numpy your full time job? You certainly seem to be > putting a full-time effort into it at any rate. It is appreciated. No :-) But, success of NumPy is critical for me. This rate of effort will have to significantly decrease very soon. I'm very anxious to get NumPy stable, though, and so I try to respond quickly to serious errors when I can --- I can be somewhat of a perfectionist so that problems eat at me until they are fixed --- perhaps there is medicine I can take ;-) Fortunately, several people are becoming familar with the internals of the code which is essential so that it can carry forward when my time to spend on it wanes. Thanks for the appreciation. -Travis

1 0

Re: [Numpy-discussion] A case for rank-0 arrays
by Francesc Altet Feb. 24, 2006

Feb. 24, 2006

A Divendres 24 Febrer 2006 02:47, vàreu escriure: > Could these be considered as dimensionless, to avoid having to explain > to people that the word rank doesn't have the same meaning as the matrix > rank? Colin Williams was proposing calling arrays coming from array(5) as 'dimensionless'. So, for the moment, we have three ways to name such a beasts: - 'rank-0' - '0-dimensional' or '0-dim' for short - 'dimensionless' Perhaps the time has come to choose a number for them (if we have to live with such array for a long time, as it seems to be the case). However, as more I think on this the more I'm convinced that, similarly to their higher dimension counterparts, we will not arrive to any definite agreement and people will continue to call them in any way he would be more comfortable with. IMO, this should be not a problem at all because all three words express a 'lack of dimensionality'. Cheers, -- >0,0< Francesc Altet http://www.carabos.com/ V V Cárabos Coop. V. Enjoy Data "-"

2 1

Re: [Numpy-discussion] Re: indexing problem
by Tim Hochberg Feb. 23, 2006

Feb. 23, 2006

>> >>Ryan Krauss wrote: >> >> >> >>>This may only be a problem for ridiculously large numbers. I actually >>>meant to be dealing with these values: >>> >>>In [75]: d >>>Out[75]: >>>array([ 246.74011003, 986.96044011, 2220.66099025, 3947.84176044, >>> 6168.50275068, 8882.64396098, 12090.26539133, 15791.36704174, >>> 19985.94891221, 24674.01100272]) >>> >>>In [76]: s=d[-1]*1.0j >>> >>>In [77]: s >>>Out[77]: 24674.011002723393j >>> >>>In [78]: type(s) >>>Out[78]: <type 'complex128scalar'> >>> >>>In [79]: s**2 >>>Out[79]: (-608806818.96251547+7.4554869875188623e-08j) >>> >>>So perhaps the previous difference of 26 orders of magnitude really >>>did mean that the imaginary part was negligibly small, that just got >>>obscured by the fact that the real part was order 1e+135. >>> >>>On 2/13/06, Ryan Krauss <ryanlists(a)gmail.com> wrote: >>> >>> I got myself all tied up in a knot over this because I couldn't figure out how multiplying two purely complex numbers was going to result in something with a complex portion. Since I couldn't find the complex routines my imagination went wild: perhaps, I thought, numpy uses the complex multiplication routine that uses 3 multiplies instead of the more straightforward one that uses 4 multiples, etc, etc. None of these panned out, and of course they all evaporated when I got pointed to the code that implements this which is pure vanilla. All the time I was overlooking the obvious: Ryan is using s**2, not s*s. So the obvious answer, is that he's just seeing normal error in the function that is implementing pow. If this is inacuracy is problem, I'd just replace s**2 with s*s. It will probably be both faster and more accurate anyway Foolishly, -tim

7 33

help(xxx) vs print(xxx.__doc__)
by Bill Baxter Feb. 23, 2006

Feb. 23, 2006

Can someone explain why help(numpy.r_) doesn't contain all the information in print(numpy.r_.__doc__)? Namely you don't get the helpful example showing usage with 'help' that you get with '.__doc__'. I'd rather be able to use 'help' as the one-stop shop for built-in documentation. It's less typing and just looks nicer. Thanks, --Bill

2 1

array.argmax() question
by Ryan Krauss Feb. 23, 2006

Feb. 23, 2006

Is this the expected behavior for array.argmax(): ipdb> abs(real(disps)).max() Out[38]: 1.7373584411866401e-05 ipdb> abs(real(disps)).argmax() Out[38]: 32 ipdb> shape(disps) Out[38]: (11, 3) ipdb> disps[11,1] *** IndexError: invalid index ipdb> disps[10,1] Out[38]: 0j ipdb> disps[10,2] Out[38]: (-1.7373584411866401e-05+5.2046737124258386e-21j) Basically, I want to find the element with the largest absolute value in a matrix and use it to scale by. But I need to correct for the possibility that the largest abs value may be from a negative number. So, I need to get the corresponding element itself. My array is shape (11,3) and argmax without an axis argument returns 32, which would be the index if the matrix was reshaped into a (33,) vector. Is there a clean way to extract the element based on the output of argmax? (and in my case it is actually using the output of argmax to extract the element from the matrix without the abs). Or do I need to reshape the matrix into a vector first? Thanks, Ryan

1 1

A proposal to implement round in C Was: Rookie problems - Why is C-code much faster?
by Sasha Feb. 23, 2006

Feb. 23, 2006

[I am reposting this under a different subject because my original post got buried in a long thread that went on to discussing unrelated topics. Sorry if you had to read this post twice.] It turns out that around (or round_) is implemented in python: def round_(a, decimals=0): """Round 'a' to the given number of decimal places. Rounding behaviour is equivalent to Python. Return 'a' if the array is not floating point. Round both the real and imaginary parts separately if the array is complex. """ a = asarray(a) if not issubclass(a.dtype.type, _nx.inexact): return a if issubclass(a.dtype.type, _nx.complexfloating): return round_(a.real, decimals) + 1j*round_(a.imag, decimals) if decimals is not 0: decimals = asarray(decimals) s = sign(a) if decimals is not 0: a = absolute(multiply(a, 10.**decimals)) else: a = absolute(a) rem = a-asarray(a).astype(_nx.intp) a = _nx.where(_nx.less(rem, 0.5), _nx.floor(a), _nx.ceil(a)) # convert back if decimals is not 0: return multiply(a, s/(10.**decimals)) else: return multiply(a, s) I see many ways to improve the performance here. First, there is no need to check for "decimals is not 0" three times. This can be done once, maybe at the expense of some code duplication. Second, _nx.where(_nx.less(rem, 0.5), _nx.floor(a), _nx.ceil(a)) seems to be equivalent to _nx.floor(a+0.5). Finally, if rint is implemented as a ufunc as Mads originally suggested, "decimals is 0" branch can just call that. It is tempting to rewrite the whole thing in C, but before I do that I have a few questions about current implementation. 1. It is implemented in oldnumeric.py . Does this mean it is deprecated. If so, what is the recommended replacement? 2. Was it intended to support array and fractional values for decimals or is it an implementation artifact. Currently: >>> around(array([1.2345]*5),[1,2,3,4,5]) array([ 1.2 , 1.23 , 1.235 , 1.2345, 1.2345]) >>> around(1.2345,2.5) array(1.2332882874656679) 3. It does nothing to exact types, even if decimals<0 >>> around(1234, -2) array(1234) Is this a bug? Consider that >>> round(1234, -2) 1200.0 and >>> around(1234., -2) array(1200.0) Docstring is self-contradictory: "Rounding behaviour is equivalent to Python" is not consistent with "Return 'a' if the array is not floating point." I propose to deprecate around and implement a new "round" member function in C that will only accept scalar "decimals" and will behave like a properly vectorized builtin round. I will do the coding if there is interest. In any case, something has to be done here. I don't think the following timings are acceptable: > python -m timeit -s "from numpy import array; x = array([1.5]*1000)" "(x+0.5).astype(int).astype(float)" 100000 loops, best of 3: 18.8 usec per loop > python -m timeit -s "from numpy import array, around; x = array([1.5]*1000)" "around(x)" 10000 loops, best of 3: 155 usec per loop

5 6

inconsistent use of axis= keyword argument?
by Vidar Gundersen Feb. 23, 2006

Feb. 23, 2006

(i've been updating the cross reference of MATLAB synonymous commands in Numeric Python to NumPy. I've kept Numeric/numarray alternatives in the source XML, but omitted it in the PDF outputs. see, http://37mm.no/download/matlab-python-xref.pdf. feedback is highly appreciated.) as i was working on this, i started wondering why a.max(0), a.min(0), a.ptp(0), a.flatten(0), ... does not allow the axis=0 keyword argument used with the exact same meaning for: m.mean(axis=0), m.sum(axis=0), ... and i also wonder why concatenate can't be used to stack 1-d arrays on top of each other, returning a 2-d array? axis relates to the number of axes in the original array(s)? n [3]: v = arange(9) In [7]: concatenate((v,v)) Out[7]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8]) In [8]: concatenate((v,v),axis=0) Out[8]: array([0, 1, 2, 3, 4, 5, 6, 7, 8, 0, 1, 2, 3, 4, 5, 6, 7, 8]) In [15]: concatenate((v,v)).reshape(2,-1) Out[15]: array([[0, 1, 2, 3, 4, 5, 6, 7, 8], [0, 1, 2, 3, 4, 5, 6, 7, 8]]) In [5]: m = v.reshape(3,-1) In [10]: concatenate((m,m)) Out[10]: array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [0, 1, 2], [3, 4, 5], [6, 7, 8]]) In [11]: concatenate((m,m), axis=0) Out[11]: array([[0, 1, 2], [3, 4, 5], [6, 7, 8], [0, 1, 2], [3, 4, 5], [6, 7, 8]]) In [12]: concatenate((m,m), axis=1) Out[12]: array([[0, 1, 2, 0, 1, 2], [3, 4, 5, 3, 4, 5], [6, 7, 8, 6, 7, 8]])

4 5

RE: repmat equivalent
by Alex Liberzon Feb. 23, 2006

Feb. 23, 2006

Another thought: def repmat(a,m,n): from scipy import hstack, vstack a = eval('hstack(('+n*'a,'+'))') return eval('vstack(('+m*'a,'+'))') may hstack() and vstack() be better for 1d arrays? >>> from scipy import * numerix Numeric 24.2 >>> a = array([1,2]) >>> repmat(a,2,3) array([[1, 2, 1, 2, 1, 2], [1, 2, 1, 2, 1, 2]]) >>> equal(repmat(1,1,1),1) array([ [1]],'b') of course, scipy.linalg.kron(ones((m,n,p,...)),a) is more robust and works for higher dimensions. probably it's the best.

1 0