On Wed, Apr 23, 2008 at 6:21 AM, Stéfan van der Walt stefan@sun.ac.za wrote:

The question is: what blocks the release of 1.1? The following tickets deserve attention:

http://projects.scipy.org/scipy/numpy/ticket/750 http://projects.scipy.org/scipy/numpy/ticket/605 http://projects.scipy.org/scipy/numpy/ticket/551 http://projects.scipy.org/scipy/numpy/ticket/738 http://projects.scipy.org/scipy/numpy/ticket/743 http://projects.scipy.org/scipy/numpy/ticket/748 http://projects.scipy.org/scipy/numpy/ticket/751 http://projects.scipy.org/scipy/numpy/ticket/736

Further issues:

• Alan's matrix indexing: x[0][0] for matrices currently yield x[0]

I was still hoping to get 1.1.0 out ASAP. Unfortunately, I have been too busy for the last week to pay attention to the release blockers. The last time I looked the only issues I felt absolutely needed to be resolved before the release were: 1. Testing the Windows binary, which is now done. 2. Testing the Mac binary, which is now done. 3. Testing or removing fromregex, which is now done: http://projects.scipy.org/scipy/numpy/ticket/719

I haven't looked carefully at the above tickets, but at a cursory glance it didn't seem like there were any new regressions. Stefan could you confirm that? If so, I would like to branch 1.1.x and tag 1.1.0 on Friday night. Unless there is some concrete reason that we need to delay this release any longer, I will send out an email later today with an 'official' timeline for the release, which will probably be on Monday if I tag on Friday night.

Thanks,

2008/4/23, Stéfan van der Walt stefan@sun.ac.za:

Hi Jarrod

Of those tickets, the following are serious:

http://projects.scipy.org/scipy/numpy/ticket/605 (a patch is available?, David Huard) Fixing of histogram.

I haven't found a way to fix histogram reliably without breaking the current behavior. There is a patch attached to the ticket, if the decision is to break histogram.

David

I think a long term strategy needs to be adopted for histogram. Right now there is a great confusion in what the "bins" keyword does. Right now it is defined as the lower edge of each bin, meaning that the last bin is open ended and [inf,bin0> does not exist. While this may not be the right thing to fix in 1.1.0, I would really like to see it fixed somewhere down the line.

On Apr 24, 2008, at 10:28 AM, Pauli Virtanen wrote:

Wed, 23 Apr 2008 16:20:41 -0400, David Huard wrote:

...

I haven't found a way to fix histogram reliably without breaking the current behavior. There is a patch attached to the ticket, if the decision is to break histogram.

Summary of the facts (again...):

a) histogram's docstring does not match its behavior wrt discarding data

This is an easy fix and should definitively go into 1.1.0 :)

b) given variable-width bins, histogram(..., normed=True) the results are wrong

Also a quick fix that should be part of 1.1.0

c) it might make more sense to handle discarding data in some other way than what histogram does now

I would like to see this, but it does not have to happen in 1.1.0 :)

I think there are now a couple of choices what to do with this:

A) Change the semantics of histogram function. Old code using histogram will just simply break, maybe in mysterious ways

Not really a satisfactory approach. I really don't mind, even though it would break some code of mine. I would rather see a better function and have to do some code changes, than the current confusion. Other people will likely disagree.

B) Rename the bins parameter to bin_edges or something else, so that any old code using histogram immediately raises an exception that is easily understood.

Given this approach bin_edges would contain one more value than bins given that the right edge of the last bin has to be defined.

C) Create a new parameter with more sensible behavior and a name different from "bins", and deprecate (at least giving sequences to) the "bins" parameter: put up a DeprecationWarning if the user does this, but still produce the same results as the old histogram. This way the user can forward-port her code at leisure.

I think this is probably the best approach to accommodate everyone.

So which one (or something else) do we choose for 1.1.0?

-- Pauli Virtanen

Cheers Tommy

The problem I see with C is that it will break compatibility with the other histogram functions, which also use bins.

So here is suggestion E:

The most common use case ( I think) is the following: h, b = histogram(r, number_of_bins, normed=True/False) for which the function behaves correctly.

Assuming we want the next version to : ignore values outside of range and accept and return the bin edges instead of the left edges, here could be the new signature for 1.1: h, edges = histogram(a, bins=10, normed=False, range=None, normed=False, new=False)

If new=False, return the histogram and the left edges, with a warning that in the next version, the edges will be returned. If new=True, return the histogram and the edges. If range is given explicitly , raise a warning saying that in the next version, the outliers will be ignored. To ignore outliers, use new=True. If bins is a sequence, raise an error saying that bins should be an integer. To use explicit edges, use new=True.

In 1.2, set new=True as the default, and in 2.3, remove new altogether.

David

2008/4/24 Travis E. Oliphant oliphant@enthought.com:

Pauli Virtanen wrote:

...

C) Create a new parameter with more sensible behavior and a name different from "bins", and deprecate (at least giving sequences to) the "bins" parameter: put up a DeprecationWarning if the user does this, but still produce the same results as the old histogram. This way the user can forward-port her code at leisure.

...

D) Or, retain the old behavior (values below lowest bin ignored) and just fix the docstring and the normed=True bug? (I have a patch doing this.)

So which one (or something else) do we choose for 1.1.0?

I like either C or D, but prefer C if it can be done before 1.1.0.

-Travis

---

Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

2008/4/24 Jarrod Millman millman@berkeley.edu:

On Thu, Apr 24, 2008 at 1:22 PM, David Huard wrote:

...

Assuming we want the next version to : ignore values outside of range

and

...

accept and return the bin edges instead of the left edges, here could be

the

...

new signature for 1.1: h, edges = histogram(a, bins=10, normed=False, range=None,

normed=False,

...

new=False)

If new=False, return the histogram and the left edges, with a warning

that

...

in the next version, the edges will be returned. If new=True, return the histogram and the edges. If range is given explicitly , raise a warning saying that in the next version, the outliers will be ignored. To ignore outliers, use

new=True.

...

If bins is a sequence, raise an error saying that bins should be an integer. To use explicit edges, use new=True.

In 1.2, set new=True as the default, and in 2.3, remove new altogether.

+1 That sounds fine to me assuming 2.3 is 1.3.

Indeed.

Done in r5085. I added a bunch of tests, but I'd appreciate if someone could double check before the release. This is not the time to introduce new bugs.

Hopefully this is the end of the histogram saga.

David

-- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/ _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

On Fri, Apr 25, 2008 at 12:55 PM, David Huard david.huard@gmail.com wrote:

...

Done in r5085. I added a bunch of tests, but I'd appreciate if someone

could double check before the release. This is not the time to introduce new bugs.

...

Hopefully this is the end of the histogram saga.

Well, it's not... there is still an issue... give me a couple of minutes to fix it.

Hey David,

Excellent! I have one more request: Could you either send me a little blurb for the release notes or go ahead and edit it yourself: http://projects.scipy.org/scipy/numpy/milestone/1.1.0

Thanks,

Thanks Chuck,

I didn't know there were other tests for histogram outside of test_function_base.

The error is now raised only if bins are passed explicitly and normed=True.

David

2008/4/25 Charles R Harris charlesr.harris@gmail.com:

On Fri, Apr 25, 2008 at 12:55 PM, Jarrod Millman millman@berkeley.edu wrote:

...

On Fri, Apr 25, 2008 at 12:55 PM, David Huard david.huard@gmail.com wrote:

...

...

Done in r5085. I added a bunch of tests, but I'd appreciate if

someone

...

could double check before the release. This is not the time to

introduce new

...

bugs.

...

Hopefully this is the end of the histogram saga.

This one?

ERROR: Ticket #632

Traceback (most recent call last): File "/usr/lib/python2.5/site-packages/numpy/core/tests/test_regression.py", line 812, in check_hist_bins_as_list hist,edges = np.histogram([1,2,3,4],[1,2]) File "/usr/lib/python2.5/site-packages/numpy/lib/function_base.py", line 184, in histogram raise ValueError, 'Use new=True to pass bin edges explicitly.' ValueError: Use new=True to pass bin edges explicitly.

Chuck

---

Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

On Wed, Apr 23, 2008 at 4:56 PM, Sebastian Haase haase@msg.ucsf.edu wrote:

On Wed, Apr 23, 2008 at 2:32 PM, Ondrej Certik ondrej@certik.cz wrote:

...

On Mon, Apr 14, 2008 at 12:06 PM, Ondrej Certik ondrej@certik.cz wrote:

...

Hi Jarrod,

any news with the 1.0.5? If you have same prerelease, I'd like to test it. Debian has just moved from python2.4 to python2.5 yesterday, so I'd like to test numpy in advance, I am sure there will be some issues to fix.

What is the plan with the release? There are some minor problems in the Debian package, some of which are fixed by the new release. I didn't fix those in Debian as I thought the new release is coming out. But if it's going to take let's say month or two, I'll fix the current Debian package now.

The problem is, that it was noticed that non-backwards compatible changes were committed into svn. Most noticeably a complete rewrite of the "ma" (Masked Array) package was done. As a consequence it was decided, they no one would be interested in "temporarily taking those changes out" "just to release 1.0.5".

Instead the next release will be called 1.1 -- which is (only; rather still) "mostly" compatible with 1.0.x What used to be referred to a the 1.1 version, that can break more stuff, to allow for a cleaner design, will now be 1.2

OK, so I think we are going to try to fix the package in Debian right now and not wait for the release.

I think this uncertainty about releases and about backward compatibility is not good. The same about uncertainty where f2py is going to end up. But you know the drill "talk is cheap, show me the code", so I know I could step up and help Jarrod with the release, but unfortunately, I don't have time for it. :)

Ondrej

On 23/04/2008, Alan Isaac aisaac@american.edu wrote:

On Wed, 23 Apr 2008, Sebastian Haase wrote:

...

What used to be referred to a the 1.1 version, that can break more stuff, to allow for a cleaner design, will now be 1.2

So ... fixing x[0][0] for matrices should wait until 1.2. Is that correct?

It seems to me that everyone agrees that the current situation is broken, but there is considerable disagreement on what the correct fix is. My inclination is to apply the minimal fix you suggest, with tests and documentation, as a stopgap, and let the different factions sort it out by discussion on the way to 1.2.

Objections?

Anne

2008/4/23 Jarrod Millman millman@berkeley.edu:

On Wed, Apr 23, 2008 at 11:32 AM, Anne Archibald peridot.faceted@gmail.com wrote:

...

...

So ... fixing x[0][0] for matrices should wait until 1.2. Is that correct?

It seems to me that everyone agrees that the current situation is broken, but there is considerable disagreement on what the correct fix is. My inclination is to apply the minimal fix you suggest, with tests and documentation, as a stopgap, and let the different factions sort it out by discussion on the way to 1.2.

+1 That sound reasonable to me.

Done in r5072.

Cheers Stéfan

Alan Isaac wrote:

I have updated URL:http://www.scipy.org/MatrixIndexing to reflect this change (and its provisional status).

Thanks for writing this up -- it really clarifies what's being proposed. A few comments on that write up:

For matrix x, should x[0].A[0] == x.A[0][0]? That is, should the array attribute of a Row Vector or Column Vector be 1d?

Absolutely -- a vector is a 1-d array. The only difference here is that we can preserve the concept of a row vs. a column vector -- key to linear algebra, which is the whole point of the Matrix class. I'd be just as happy if a RowVector an ColumnVector were a more unique object -- a 1-d array with the operation overloaded, rather than a (1,n) or (n,1) matrix, but I understand that would be more work.

Deviations from the behavior of ndarrays. Specifically, iteration over a matrix will not yield 1d NumPy arrays.

Why should they? indexing a nested list gives a different result than indexing an ndarray, etc, etc, etc. The POINT of matrixes is that they are different! This proposal makes them a little more similar to ndarrays than the old behavior (always returning a matrix).

Inconsistency in the indexing style for producing row vectors and column vectors

not really:

row vector: M[i,:] column vector: M[:,i]

The fact that there is another way to get a row vector: M[i] is a bonus. I once proposed making M[i] raise an exception for the sake of enforcing consistency, but no one liked that!

Loss of a standard way to extract a row or a column as a submatrix (rather than a vector)

What's wrong with:

M[i:i+1, :] and M[:, i:i+1]

This is totally consistent with arrays (and other python sequences). Yes, it's a bit more awkward, but with the new vectors, it's also going to be needed far less.

No clear gain in functionality over the simpler proposal 2.

Well, I won't judge what's "clear" but I think this adds functionality. The concept of row and column vectors is a key part of linear algebra. Sure, you can model them as matrixes that happen to have only one row or one column, but that's what we had, but there seem to be real issues with use of that. Option (2), is "Let a Matrix be a container of 1d arrays.". Sounds simple, but do those arrays represent row or column vectors? Either would make sense. The proposal is for them to be row vectors, which is fine, but then if you want them to act like a row vectors, do you need to convert them back into matrices? Even if not, how do you get a column vector -- by indexing the current way, and getting a (n,1) matrix, so that now we have inconsistency between how row and column vectors are treated -- not a clean API.

And this is about a clean API for linear algebra -- otherwise, we'd just all use arrays (which many do!)

Aside from the fact that someone needs to write the code -- why don't people like the row/column vector idea? It just feels so natural to me:

A matrix is a 2-d array with some functionality overloaded to make it convenient to do linear algebra.

A vector is a 1-d array with some functionality overloaded to make it convenient to do linear algebra.

A scalar is the same in both contexts, so no need to change anything there.

And yes, maybe this will lead to tensors some day!

Another couple questions: -- would there be constructors for vectors?

np.RowVector((1,2,3,4,5)) np.ColumnVector((1,2,3,4,5,6))

and would: A_RowVector.T == A_ColumnVector

I would think so.

One more note:

All the example code I've seen during this discussion have been examples of iterating and indexing -- not one has actually been linear algebra -- perhaps we should see some of those before making any decisions.

-Chris

[CHOP]

The proposals thus far don't address two of the major issues I have with the matrix class:

1. The matrices and arrays should become more alike if possible and should share more of the same code base. From what I've seen, the people who write the code (for numpy) don't actually use matrices, they use the arrays. This is one of the main reasons matrices still have so many problems IMO; if one of the main developers were using them, they'd never have put up with it. This may be changing to some extent now, but the more we can make matrices and arrays share code and interface, the happier we'll be. None of the new rules strike me as helping in this arena and may in fact hurt. Again, IMO. 2. This doesn't support the concept of arrays of matrices/vectors, which is one thing I've always wanted out of the matrix class. For example it would be nice to be able to spell: 'A' is a 1D array of matrices (3D) overall, 'B' is a 1D array of vectors (3D) overall, matrix multiply them together, yielding a 1D array of matrices (3D) overall. I have frequently run into various permutations of this problem. You can fake it with loops over dot, but the efficiency is very bad for small arrays, plus it's ugly.

I have lot's of vague ideas on this subject that have been floating around my head for the past couple of years. Basically, I think the way to go is to probably add some meta information to arrays that make them look a little bit more like tensors. Then, if one is careful (and things work out as I hope) matrices could be implemented as a thin layer on top of arrays. Or perhaps, if we are very lucky, done away with altogether.

Anyway, that's all very vague and hand-wavey; if I can spring free some time, I'll try to write something up in the not too distant future. I'm glad to see that any major changes are being put off for a while: I think we should be able to do better than just patch up the old matrix class.

-tim

2008/4/24 Nadav Horesh nadavh@visionsense.com:

+1

+1 to what? I'm glad that Tim chimed in -- his opinion is always valued -- but I didn't see any concrete suggestions on how to improve the situation. I'll gladly wait for his code, though :)

Regards Stéfan

Hi, I would like to use the matrix functionality but I, like others, get by without it.

+1 to Tim's and Nadav's comments. As Tim said, there should be seamless integration between concepts of vectors and matrices - after all there really no distinction between them in linear algebra.

-2 for having rowvector and columnvector - both numpy and I should know the orientation, if not, I deserve garbage or an error. 0 for vector - I don't see the need for it as a unique entity but understand the simplicity of saying vector.

Yes, I agree that talk is cheap but it avoids wasting time on code that will not be used and allows time to understand what people really want.

Bruce

Nadav Horesh wrote:

Some, including me, stated that the matrix class suffer some deficiencies, so, after playing a bit with matrices, we always go back to the good old ndarray class. There were some arguments on this list what/if should be dome to correct it. I welcome any attempt to build a code that demonstrates an alternative.

Nadav.

-----הודעה מקורית----- מאת: numpy-discussion-bounces@scipy.org בשם St?fan van der Walt נשלח: ה 24-אפריל-08 13:13 אל: Discussion of Numerical Python נושא: Re: [Numpy-discussion] numpy release

2008/4/24 Nadav Horesh nadavh@visionsense.com:

...

+1

+1 to what? I'm glad that Tim chimed in -- his opinion is always valued -- but I didn't see any concrete suggestions on how to improve the situation. I'll gladly wait for his code, though :)

Regards St?fan _______________________________________________ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

---

---

Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion

On Thu, Apr 24, 2008 at 11:16 AM, Timothy Hochberg tim.hochberg@ieee.org wrote:

[CHOP]

The proposals thus far don't address two of the major issues I have with the matrix class:

The thing that seems missing to me is support for LAPACK's banded and packed (triangular) storage formats. I don't have an urgent need for it, but it seems like a serious front end for doing linear algebra should at least provide a convenient way to access all of what LAPACK offers.

Perhaps these other formats could/should be supported by different classes, though, the way it is done for sparse matrices. But the fact that more matrix types could exist some day is all the more reason to make sure the interface makes sense for more than just strided memory. My hunch is if the API is workable for both ndarray matrices and sparse matrices, then it should be ok for these other formats too.

--bb

Alan G Isaac wrote:

the cost of complexity should be justified by a gain in functionality.

I don't think functionality is the right word here. the Matrix class(es) is all about clean, convenient API, i.e. style, not functionality -- we have all the functionality already, indeed we have it with plain old arrays, so I think that's really beside the point.

I just like the "feel" of Row and Column Vectors.

But this brings me back to:

We need examples of doing linear algebra operations, not just iterating operations. Frankly, if you have a matrix, and you want to iterate through all the rows, and you want to index into them as 1-d arrays, why they heck don't you just:

for r in MyMatrix.A: ....

The answer (I think), is that folks have simply been confused and annoyed by the fact that they need that extra ".A", which brings us back to style.

One extra functionally: you can index a ColumnVector with a scalar:

col = M[:,i] col[i]

You can't do that now.

This brings me to:

the people who write the code (for numpy) don't actually use matrices, they use the arrays

I don't really write numpy code, but I don't really use arrays either. Which, I suppose, means I should just shut up. But maybe I would if the API was better.

So we're back to: who does want matrices? is it only newbies coming from Matlab, where they were used to them?

For once, we do have someone that is writing the code that wants improvements, so I guess that's why this discussion has been revived.

Timothy Hochberg wrote:

1. The matrices and arrays should become more alike if possible

I'm not sure I agree -- the more alike they are, the less point there is to having them.

  should share more of the same code base. 

Don't they share almost all the same code already? My understanding is that they are arrays with a couple operators overloaded -- how much more code could they share?

Nevertheless, I'm interested to see what Tim comes up with.

2. This doesn't support the concept of arrays of matrices/vectors, which is one thing I've always wanted out of the matrix class. For example it would be nice to be able to spell: 'A' is a 1D array of matrices (3D) overall, 'B' is a 1D array of vectors (3D) overall,

You could do this now with a 1-d object array, with a bunch of matrices in it.

  matrix multiply them together, yielding a 1D array of matrices
  (3D) overall.

But this, or course, would fail. However, if we expand this scenario to A and B being 1-d arrays of other arrays that may not all be the same size (and thus you couldn't just use a 3-d array), and be able to run various ufuncs on them, that would be cool!

Bill Spotz wrote:

...

The fact that there is another way to get a row vector: M[i] is a bonus.

Except that the behavior of M[i] is one of the driving issues of the conversation.

well, yes, but a bit of a red herring, I think -- it's syntactic sugar. The issue, I think, is that even if you do:

r = M[:,i]

you can't then index that with a single index -- it's a matrix, NOT a 1-d array. Again, I proposed disallowing M[i] altogether, as it is ambiguous, but that really is gratuitous consistency.

...

A vector is a 1-d array with some functionality overloaded to make it convenient to do linear algebra.

Again, I would argue for Vectors inheriting from Matrix. I would make the case based on code reuse and elegance, although these might be overridden by some other concern.

I think that's a bit of an implementation detail -- we need to define behavior that we want, and decide how best to implement that. i.e.:

RowVector[i] -> scalar RowVector.A -> 1-d array (same for column)

How do we want it pretty-printed?

We need to focus on column vectors here -- It's easy to say: a row vector is a 1-d array -- that's easy and clean. It's column vectors that are tricky -- indeed they are a pain with regular old arrays.

What about

np.Matrix(<container of RowVector>) np.Matrix(<container of ColumnVector>)

I'd say you get a matrix from each of those, but they would be transposes of each-other -- another good reason for the Vector classes. though maybe those should be spelled:

np.row_stack and np.column_stack

I have generally thought about this in the context of, say, a Krylov-space iterative method, and what that type of interface would lead to the most readable code.

Can you whip up a small example, starting with the current implementation?

Bruce Southey wrote:

Hi, I would like to use the matrix functionality but I, like others, get by without it.

What does it lack that keeps you from using it? That's the key question?

+1 to Tim's and Nadav's comments. As Tim said, there should be seamless integration between concepts of vectors and matrices - after all there really no distinction between them in linear algebra.

This is all about being able to index a "vector" with a single index -- that's it, I think. Matlab handles this by special casing matrices that have a single dimension of one. It also has an easy way to spell a literal for column vector, though I suppose:

np.Matrix((1,2,3,4,5)).T

isn't so bad.

-2 for having rowvector and columnvector - both numpy and I should know the orientation, if not, I deserve garbage or an error.

But HOW can you know the orientation? a 1-d array has no orientation -- which is why the current version always returns a matrix when indexing.

0 for vector - I don't see the need for it as a unique entity but understand the simplicity of saying vector.

I don't think we can have vector without ColumnVector, though I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.

By the way, maybe a column vector could be handy for doing array broadcasting, as an cleaner way to do:

x = np.arange(10) y = np.arange(20).reshape((-1,1))

z = FunctionOfXandY(x,y)

Final note:

If no one has an linear algebra examples, then we're wasting out time with even this cheap talk.

-Chris

On Thu, 24 Apr 2008, Christopher Barker apparently wrote:

I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.

This bothers me.

So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.

Or are you still having separate classes and not just changing ``__getitem__``? OK, I see that is most likely what you meant.

Cheers, Alan

2008/4/25 Alan G Isaac :

...

So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.

On Fri, 25 Apr 2008, Stefan van der Walt apparently wrote:

Yup. That's the current state of things.

I do not understand. The released "state of things" is that for matrix ``x`` we have that ``x[0]`` is a **matrix**.

I'm not working with SVN NumPy, but my understanding was that as of revision r5072 ``x[0]`` always returns a 1d array. The core requirement was that ``x[0][0]`` produce the first element of the matrix.

I do not have time to look at the revision right now, but if a matrix ``x`` we have that ``x[0]`` can return a scalar, that is very undesirable.

Cheers, Alan

On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:

In current SVN: In [6]: x[0] Out[6]: matrix([[0, 1, 2]])

I must have misunderstood: I thought the agreement was to provisionally return a 1d array for x[0], while we hashed through the other proposals.

Cheers, Alan

On 25/04/2008, Stéfan van der Walt stefan@sun.ac.za wrote:

2008/4/25 Alan G Isaac aisaac@american.edu:

...

I must have misunderstood: I thought the agreement was to provisionally return a 1d array for x[0], while we hashed through the other proposals.

The agreement was:

a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)

This breaks as little functionality as possible, at the cost of one (instead of two) API changes.

Hold on. There has definitely been some confusion here. This is not what I thought I was suggesting, or what Alan thought he was suggesting. I do not think special-casing matrices for which one dimension happens to be one is a good idea at all, even temporarily. This is the kind of thing that drives users crazy.

My suggested stopgap fix was to make x[0] return a 1D *array*; I feel that this will result in less special-casing. In fact I wasn't aware that anyone had proposed the fix you implemented. Can we change the stopgap solution?

We should now discuss the proposals on the table, choose a good one, and implement all the API changes necessary for 1.2 or 2. It's a pity we have to change the API again, but the current situation is not tenable.

Yes, well, it really looks unlikely we will be able to agree on what the correct solution is before 1.1, so I would like to have something non-broken for that release.

Anne

2008/4/25 Anne Archibald peridot.faceted@gmail.com:

...

The agreement was:

Clearly that was only the agreement in my mind. I apologise if I jumped the gun.

...

a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)

This breaks as little functionality as possible, at the cost of one (instead of two) API changes.

Hold on. There has definitely been some confusion here. This is not what I thought I was suggesting, or what Alan thought he was suggesting. I do not think special-casing matrices for which one dimension happens to be one is a good idea at all, even temporarily. This is the kind of thing that drives users crazy.

Apparently, not having x[0][0] == x[0,0] drives them even crazier :) I was trying to satisfy Alan's request without breaking the API in more than one place. The whole idea of matrices was that you always work with 2-dimensional arrays, even when you just extract a row. Unless you have a proper hierarchical container (as illustrated in the other patch I sent), returning a 1D array breaks at least some expectation. If that's what the matrix users want, though, then we should change it.

My suggested stopgap fix was to make x[0] return a 1D *array*; I feel that this will result in less special-casing. In fact I wasn't aware that anyone had proposed the fix you implemented. Can we change the stopgap solution?

Yes, of course. NumPy is a community project, after all.

...

We should now discuss the proposals on the table, choose a good one, and implement all the API changes necessary for 1.2 or 2. It's a pity we have to change the API again, but the current situation is not tenable.

Yes, well, it really looks unlikely we will be able to agree on what the correct solution is before 1.1, so I would like to have something non-broken for that release.

Unless we agree on something, and soon, that won't happen. Your workaround would break x[0] == x[0,:], so we're just swapping one set of broken functionality for another.

I'm starting to see Chris Barker's point; allowing x[0] is causing more problems than it is worth. On the other hand, how would you index into a vector (as in http://en.wikipedia.org/wiki/Vector_(spatial)) without it?

Regards Stéfan

2008/4/25 Alan G Isaac aisaac@american.edu:

...

I must have misunderstood: I thought the agreement was to provisionally return a 1d array for x[0], while we hashed through the other proposals.

On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:

The agreement was: a) That x[0][0] should be equal to x[0,0] and b) That x[0,:] should be equal to x[0] (as for ndarrays)

1. This is **not** what I understood as the agreement (and I think the current solution is bad). I certainly did not mean to support the change that as implemented, and it is not clear to me that others did either.

2. These two goals are substantially in conflict for matrices, as we are seeing.

3. The important goal, (goal a., which everyone agrees on), has NOT been accomplished by the current change: x[0][0] raises a TypeError when x is a 1 by N matrix.

This breaks as little functionality as possible, at the cost of one (instead of two) API changes.

I cannot agree with this assessment.

the current situation is not tenable.

I agree with this! Better than the SVN behavior would be to raise an error in response to scalar indexing of matrices. I strongly hope that this behavior will not show up in a release!

Please return 1d arrays in response to scalar indexing as the provisional arrangement. This is certainly better than the new behavior.

Alan Isaac

On Fri, 25 Apr 2008, Stéfan van der Walt wrote:

workaround would break x[0] == x[0,:]

But there is not universal agreement that x[0] == x[0,:] is desirable. In contrast, there *is* universal agreement that x[0][0]==x[0,0] is desirable. Or so I've understood the discussion.

As you know, my view is that nonscalar indexing should return submatrices, while scalar indexing should return 1d arrays. Thus I do not see x[0] != x[0,:] as breakage: it is desirable. Indeed, without it, we will have an untenable situation, under all the different proposals. So it is a goner one way or another.

Cheers, Alan Isaac

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

I don't see where you get your universality.

Lists, tuples, arrays ... Where is the exception? Only matrices are the only 2d container I can think of that are broken this way.

Cheers, Alan Isaac

...

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

...

I don't see where you get your universality.

On Sat, Apr 26, 2008 at 11:15:14AM -0400, Alan G Isaac wrote:

...

Lists, tuples, arrays ... Where is the exception? Only matrices are the only 2d container I can think of that are broken this way.

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

List and Tuples are not 2D containers.

A list of lists is naturally conceived as 2d, which is why you can initialize a 2d matrix with a list of lists.

Try a Google search on "python two dimensional list" to confirm that I am not promoting an anomalous view.

So I do not buy the way you are trying to parse language here. (Although I'm open to a fuller argument.)

Cheers, Alan

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

We may desagree on what is Right, but the breaks backward compatibility, and thus is a breakage and should be given a lot of thought.

I agree with this. I am a matrix user, and I have given it a lot of thought. I have been making this case for a *long* time. It has come to a head because of the announced desire to severely constrain API changes moving forward.

As I am best able to understand you abstract view, the right thing to do would be for matrices to raise an error in response to a scalar index. However, I deduce, you would just leave things alone to avoid backward incompatibility.

I weight the future more heavily. We are approaching a last chance to do things better, and we should seize it.

The right questions looking forward:

- what behavior allows the most generic code? - what behavior breaks fewest expectations? - what behavior is most useful?

Cheers, Alan Isaac

PS I cannot recall: do you use matrices?

oops, typo!

Christopher Barker wrote:

Gael Varoquaux wrote:

...

I don't know why people are indexing matrices with A[x][y], but they shouldn't.

I think there has been a misunderstanding here. I don't think anyone is suggesting that if a coder wants an element of a matrix, that s/he should write that:

element = M[i,j]

I meant -- no one is suggesting:

element = M[i][j]

because, of course, M[i,j] is exactly what you should do do get an element.

Which, by the way, brings up another inconsistency in the current behavior:

...

...

M

matrix([[1, 2, 3], [4, 5, 6], [7, 8, 9]])

...

...

M[1,1]

5

...

...

M[1,:]

matrix([[4, 5, 6]])

so if you pull out a single element, you get a scalar(0-d), but if you pull out a row, you get a matrix(2-d). It makes me want vectors more...

-Chris

On Tue, Apr 29, 2008 at 12:20 AM, Alan G Isaac aisaac@american.edu wrote:

On Mon, 28 Apr 2008, Christopher Barker apparently wrote:

...

I'm going to try to put down what I think are the key, very simple, questions:

For useful reference, these are now here: URL:http://www.scipy.org/MatrixIndexing#guidelines

Cheers, Alan

May I add that if I edit defmatrix.py to act like an array for scalar indexing, then the following works.

In [1]: a = matrix(eye(2))

In [2]: array([a,a]) Out[2]: array([[[ 1., 0.], [ 0., 1.]],

[[ 1., 0.], [ 0., 1.]]])

This generates an error with the current version of matrix and, frankly, I am not going to be bothered going all through the numpy c sources to special case matrices to fix that. Someone else can do it if they wish. There are recursive routines that expect the dimensions to decrease on each call.

Chuck

On Tue, 29 Apr 2008, Stéfan van der Walt apparently wrote:

We were thinking of raising a warning on scalar indexing for 1.1, but given the above, would that be sensical?

There seem to be three basic proposals for scalar indexing:

- raise an error - return a 1d array - return a new type of 1d container (apparently to be called a "vector" for some reason ...)

I do not think anyone is still arguing in favor of the keeping the current behavior (where scalar indexing returns a matrix)?

So the question seems to be: how best to warn users that the behavior of scalar indexing is very likely to change.

Cheers, Alan Isaac

Charles R Harris wrote:

On Tue, Apr 29, 2008 at 7:24 AM, Stéfan van der Walt <stefan@sun.ac.za mailto:stefan@sun.ac.za> wrote:

Hi Charles

2008/4/29 Charles R Harris <charlesr.harris@gmail.com
<mailto:charlesr.harris@gmail.com>>:
> May I add that if I edit defmatrix.py to act like an array for
scalar
> indexing, then the following works.
>
> In [1]: a = matrix(eye(2))
>
> In [2]: array([a,a])
> Out[2]:
> array([[[ 1.,  0.],
>          [ 0.,  1.]],
>
>        [[ 1.,  0.],
>         [ 0.,  1.]]])
>
> This generates an error with the current version of matrix and,
frankly, I
> am not going to be bothered going all through the numpy c
sources to special
> case matrices to fix that. Someone else can do it if they wish.
There are
> recursive routines that expect the dimensions to decrease on
each call.

I'd also like to see matrices become proper hierarchical containers --
the question is just how to do that.  Thus far, I'm most convinced by
the arguments for RowVectors/Columns, which leaves us with a sane
model for doing linear algebra, while providing the enhancements you
mentioned here and in comments to another ticket.

We were thinking of raising a warning on scalar indexing for 1.1, but
given the above, would that be sensical?

Hi Stefan,

The numpy c routines call PySequence_GetItem(s,i) as well as ask for the length (first index), so I think these should be left as they are for arrays in order to guarantee that matrices are compatible with all the normal array operations. This means either returning special row objects that index as expected or returning 1D arrays. I don't think the '*' operator has these problems, but in any case that is a well documented feature of matrices.

Thanks for looking in to this Chuck,

I'm quite persuaded now that a[i] should return a 1-d object for arrays. In addition to the places Chuck identified, there are at least 2 other places where special code was written to work-around the expectation that item selection returns an object with reduced dimensionality (a special-cased .tolist for matrices and a special-cased getitem in the fancy indexing code).

As the number of special-case work-arounds grows the more I'm convinced the conceptualization is wrong. So, I now believe we should change the a[i] for matrices to return a 1-d array.

The only down-side I see is that a[i] != a[i,:] for matrices. However, matrix(a[i]) == a[i,:], and so I'm not sure there is really a problem, there. I also don't think that whatever problem may actually be buried in the fact that type(a[i]) != type(a[i,:]) is worse than the problem that several pieces of NumPy code actually expect hierarchical container behavior of multi-dimensional sequences.

I don't think making the small change to have a[i] return 1-d arrays precludes us from that 1-d array being a bit more formalized in 1.2 as a RowVector should we choose that direction. It does, however, fix several real bugs right now.

So, I'm

+1 on making the following change:

def __getitem__(self, index): + if isscalar(index): + return self.__array__()[index]

-Travis

On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:

I really don't like a[i] != a[i,:].

Tim H's proposal avoids that problem. What do you think of it?

I also don't like loosing the information that you are doing linear algebra.

Hmmm. Is it a given that asking for ``a[i]`` is "doing linear algebra"? Is ``a[i:i+1,:]`` so bad if you want a matrix? I think some more use cases, like Bill's, might prove clarifying. URL:http://www.scipy.org/ConjugateGradientExample

Alan

On Tue, Apr 29, 2008 at 02:41:45PM -0400, Alan G Isaac wrote:

On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:

...

I really don't like a[i] != a[i,:].

Tim H's proposal avoids that problem. What do you think of it?

Sorry, I guess I am lost. Could you remind me which one it is? Given how well you are summing up the discussion on http://www.scipy.org/MatrixIndexing , I suspect it is already in the list of proposal.

...

I also don't like loosing the information that you are doing linear algebra.

Hmmm. Is it a given that asking for ``a[i]`` is "doing linear algebra"? Is ``a[i:i+1,:]`` so bad if you want a matrix?

Very, very bad. This is completly unreadable. Actually to the beginner, this looks like a 2*N matrix, not a 1*N matrix.

To me you are really replacing a problem by another. We should fix the problem, not replace it.

Gaël

On Tue, Apr 29, 2008 at 03:32:17PM -0400, Alan G Isaac wrote:

On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:

...

Could you remind me which one it is?

Proposal #6 at URL:http://www.scipy.org/MatrixIndexing#list-of-proposals is my effort to summarize what Tim was proposing. It does not include the oriented 1d vectors that I think you need (if I have understood you).

I prefer proposal number 1.

PS For beginners, perhaps ``x[[i],:]`` would be better than ``x[i:i+1,:]``. It has other advantages too...

I think this will get beginners lost. From my experience, people take a long time to understand indexing with an iterable. It is actually a source of confusion, but I am not even remotely thinking of removing it, because it is soooo nice, when you know how to use it.

Gaël

Gael Varoquaux wrote:

On Tue, Apr 29, 2008 at 10:49:50AM -0500, Travis E. Oliphant wrote:

...

As the number of special-case work-arounds grows the more I'm convinced the conceptualization is wrong. So, I now believe we should change the a[i] for matrices to return a 1-d array.

...

The only down-side I see is that a[i] != a[i,:] for matrices. However, matrix(a[i]) == a[i,:], and so I'm not sure there is really a problem, there.

I think we have simply replaced one problem by another one. I think this is a bad route to go.

The problem is that I know exactly what the problem we are replacing, but I have no idea only vagueries about the problem we are getting. So, in the short term, it seems like the right thing to do, because we aren't preventing a change to RowVector/ColumnVector in the future.

-Travis

On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:

I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?

It breaks A*B[i, :] where A and B are matrices.

Gaël

On Tue, Apr 29, 2008 at 11:16:15PM +0200, Anne Archibald wrote:

On 29/04/2008, Gael Varoquaux gael.varoquaux@normalesup.org wrote:

...

On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:

...

I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?

...

It breaks A*B[i, :] where A and B are matrices.

Really? How?

In [26]: A = np.matrix([[1,0],[0,1],[1,1]])

In [28]: A*np.ones(2) Out[28]: matrix([[ 1., 1., 2.]])

In [29]: np.ones(3)*A Out[29]: matrix([[ 2., 2.]])

Yes, sorry, I am shameful. I should have thought a bit more before posting.

There is no big problem with x[i,:] not being a 1-d object. The problem is for x[:, i]. However I would find it nice that, for linear algebra, x[i, :] == x[:, i].T

This is the kind of behavior I expect, and we won't be getting it with 1D arrays.

Cheers,

Gaël

On Tue, Apr 29, 2008 at 03:03:51PM -0700, Christopher Barker wrote:

Gael Varoquaux wrote:

...

However I would find it nice that, for linear algebra, x[i, :] == x[:, i].T

...

This is the kind of behavior I expect, and we won't be getting it with 1D arrays.

But you WOULD get it with 1-d row/column objects.

Yes. I agree 100% with you that this is the best solution. I fully agree with your complete message.

Gaël

On Tue, Apr 29, 2008 at 3:07 PM, Gael Varoquaux < gael.varoquaux@normalesup.org> wrote:

On Tue, Apr 29, 2008 at 11:03:58PM +0200, Anne Archibald wrote:

...

I am puzzled by this. What is the rationale for x[i,:] not being a 1-d object?

It breaks A*B[i, :] where A and B are matrices.

A*B[:,i], surely.

So it's a notational problem, what is the easiest way to get row/col vectors out of a matrix. And it will be helpful if the notation is consistent and easy to understand. In any case, dot doesn't care, so the fault above comes from the matrix class itself.

Chuck

On Tue, 29 Apr 2008, Christopher Barker apparently wrote:

a really good argument for: A * B.col[i]

Also see the syntax discussed in Proposal 5. (I am not expressing an opinion.) URL:http://www.scipy.org/MatrixIndexing#proposal-5

One possibility is to let the rows and cols methods take an argument with a default of None. The default is to return an iterator over the rows/columns. An iterator over subsets of rows/columns could be specified with a sequence. A single row/column could be specified with a scalar argument.

Cheers, Alan

On Tue, 29 Apr 2008, "Travis E. Oliphant" apparently wrote:

x[i,:] * A * x[:,i]

This need is also fully met by Proposal 5. It is just a matter of syntax and sticking with matrices. URL:http://www.scipy.org/MatrixIndexing#proposal-5

E.g., x.rows(i) * A * x.cols(i)

Cheers, Alan Isaac

Alan G Isaac wrote:

For useful reference, these are now here: URL:http://www.scipy.org/MatrixIndexing#guidelines

Thanks Alan, Here are a few comments on that text (I know, I really should go dig up my scipy password...)

1. Do we want to be able to extract a single row or column from a

matrix, and have it be indexable like a 1-d object? (Answer:

certainly > yes for rows, but columns could be gotten from the transpose.

I don't get this -- transpose of what? if you get a 1-d ndarray from: M[i,:], that can't be transposed -- it's 1-d. You could do M.T[0] but then you've lost the concept of it being a column. If we stick with M[i,:] being a (n,1) matrix, then it doesn't index like a 1-d object, and it doesn't need transposing anyway.

This is most useful, allows the most generic code, and breaks fewest expectations.)

2. Do we want to be able to do (1) without calling the Array

attribute: M.A[i]? (Yes: this is not less useful, allows the most generic code, and breaks fewest expectations.)

3. Do we want to have those 1-d objects to have an orientation so

that they act like Row or Column vectors in linear algebra operations (and broadcasting)? (Unclear: the resulting functionality could be provided by rows and columns attributes that yield corresponding matrices.)

no, it couldn't -- we can already use M[:,i] (or M[:,i:i+1] to get a matrix that is essentially a column -- however that can't be indexed with a scalar (see (1). That's kind of the point of this discussion. Or did you mean that you get a column by: np.matrix(M[:,i]).T ? yow!

Maybe we need to separate the row and column questions -- particularly with Charles' discovery of how other parts of numpy expect rank reduction with indexing, I think we probably have a consensus that we need to be able to extract 1-d objects from matrixes with indexing -- at least rows.

That leaves the column question: -- do we need an object that represents a column, that otherwise acts like a 1-d array (returns a scalar with scalar indexing)?

I guess my key point is that the fact that seeing a Matrix as a collection of rows rather than columns is an accident of C memory arrangement -- if numpy (and python) had been written in Fortran, we'd have the opposite. I'd rather that not be built into the API more than required. We can only have one definition for M[i], it might as well be a row, but other than that, we can have symmetry:

M[i,:] is row M[:,i] is a column

and rows and columns behave the same way everywhere they "should": indexing, converting to ndarrays, ???, but not in linear algebra operations.

Then, or course there are implementation details...

-Chris

Alan G Isaac wrote:

Perhaps I should note that I have come around to what I take to be Tim H's view that matrix indexing should should work just like ndarray indexing, except that when it produces a 2d result then a matrix is returned. (This means giving up the current behavior of ``x[0,:]``.) I also agree with those who want ``rows`` and ``cols`` attributes for iteration.

+1 on all of that.

A small difference: I would have these yield matrices. (I.e., no need for a new "vector" type.)

So the only thing we "give up" (which we don't have now anyway) is something that acts like column and can be indexed like a 1-d object. I still think that would be nice, but maybe we're coming close to a consensus about the rest of it.

-Chris

On Tue, 29 Apr 2008, Keith Goodman apparently wrote:

I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.

1. Could you give an example of the circumstances of this use? 2. Would any or all of the following be just as good a result? a. 1d matrix b. row and column vectors (1d but "oriented" for linear algebra purposes) 3. If 2a. and 2b. are no good for you, a. would e.g. ``x[i:i+1,:]`` be too burdensome? b. would ``rows`` and ``columns`` attributes that yielded matrics be an offsetting convenience?

Thank you, Alan Isaac

PS It is nice to have another matrix user in the conversation!

Let me throw out a couple of more thoughts:

First, there seems to be disagreement about what a row_vector and column_vector are (and even if they are sensible concepts, but let's leave that aside for moment). One school of thought is that they are one-dimensional objects that have some orientation (hence row/column). They correspond, more or less, to covariant and contravariant tensors, although I can never recall which is which. The second view, which I suspect is influenced by MatLab and its ilk, is that they are 2-dimensional 1xN and Nx1 arrays. It's my view that the pseudo tensor approach is more powerful, but it does require some metainformation be added to the array. This metadata can either take the form of making the different objects different classes, which leads to the matrix/row/column formulation, or adding some sort of tag to the array object (proposal #5, which so far lacks any detail).

Second, most of the stuff that we have been discussing so far is primarily about notational convenience. However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors). As a concrete example, I at times need to work with stacks of small matrices. If I do the operations one by one, the overhead is prohibitive, however, most of that overhead can be avoided. For example, I rewrote some of the linalg routines to work on stacks of matrices and inverse is seven times faster for a 100x10x10 array (a stack of 100 10x10 matrices) when operating on a stack than when operating on the matrices one at a time. This is a result of sharing the setup overhead, the C routines that called are the same in either case.

Hi,

On Tue, Apr 29, 2008 at 2:22 PM, Timothy Hochberg tim.hochberg@ieee.org wrote:

Let me throw out a couple of more thoughts:

First, there seems to be disagreement about what a row_vector and column_vector are (and even if they are sensible concepts, but let's leave that aside for moment). One school of thought is that they are one-dimensional objects that have some orientation (hence row/column). They correspond, more or less, to covariant and contravariant tensors, although I can never recall which is which. The second view, which I suspect is influenced by MatLab and its ilk, is that they are 2-dimensional 1xN and Nx1 arrays. It's my view that the pseudo tensor approach is more powerful, but it does require some metainformation be added to the array. This metadata can either take the form of making the different objects different classes, which leads to the matrix/row/column formulation, or adding some sort of tag to the array object (proposal #5, which so far lacks any detail).

Actually I think that this simply stems from the different variants of linear (and more specifically matrix) algebra. My background is very heavy in the statistical favor of this (what the heck are tensors... http://mathworld.wolfram.com/Tensor.html). It is also clear (from that definition) that others are from more physics and engineering backgrounds. Hence all the back and forth because people have slightly different usage, terminology and expectations.

The ability to treat vectors as matrices would be sufficient for my needs because these are almost always used in the context of vector-matrix multiplication. There is no additional benefit from having row or column shapes or metadata because the row/column nature is usually predetermined and would be represented by the shape of the corresponding matrix. It really would be annoying to find that for an n by 1 vector/matrix, the product of X.T*X is a scalar or an error rather than an n by n matrix. Or requires additional code to get the desired effect.

Regards Bruce

The problem is that ``x[0]`` being 2d has produced a variety of anomalies, and the natural fix is for ``x[0]`` to be 1d.

Gael has argued strongly that she should be able to use the following notation: ``x[0,:]*A*x[:,0]``. But this will work only if ``x[0,:]`` is 2d or if it is 1d but has an "orientation".

So *if* you think ``x[0] == x[0,:]`` is desirable, *and* you want to satisfy Gael, *then* it seems you must introduce 1d "oriented" vectors.

I believe Travis is also suggesting that we travel that road, taking a first step as follows: for now let ``x[0]`` be a 1d array to quickly fix the anomalies, but let ``x[0,:]`` continue to be a matrix until the vector code is written, at which point ``x[0]`` and ``x[0,:]`` we be the same "row vector".

Or so I have understood things.

You've characterized my current thinking pretty well. I'm less concerned that x[0] != x[0,:] than I think Gael is.

-Travis

On Tue, 29 Apr 2008, Keith Goodman apparently wrote:

I hope that changing x[0,:] is considered a major change since it will break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.

What if the extant matrix class would continue to be available for awhile as "oldmatrix", which you could then import as "matrix" if you wished. Would that meet your needs? (I am simply assuming this would be feasible, without being sure how a lot of the special casing for matrices has been done.)

Cheers, Alan Isaac

On Apr 29, 2008, at 6:01 PM, Keith Goodman wrote:

I hope that changing x[0,:] is considered a major change since it will break most any package based on matrices (mine). And so I hope that such a change wouldn't show up, if at all, until 2.0.

The only code that should break would be indexing the extracted row/ column with two indexes. Multiplication with the extracted row/column would still work as expected (and that sounded to me like what the application typically was . . . maybe I'm wrong).

** Bill Spotz ** ** Sandia National Laboratories Voice: (505)845-0170 ** ** P.O. Box 5800 Fax: (505)284-0154 ** ** Albuquerque, NM 87185-0370 Email: wfspotz@sandia.gov **

On Tue, 29 Apr 2008, Alan G Isaac apparently wrote:

Gael has argued strongly that he should be able to use the following notation: ``x[0,:]*A*x[:,0]``.

By the way Gael, is x.rows(0) * A * x.cols(0) a good replacement in your view? I find it nicely explicit and, to meet another of your concerns, usable by newbies.

Cheers, Alan

Bruce Southey wrote:

The ability to treat vectors as matrices would be sufficient for my needs because these are almost always used in the context of vector-matrix multiplication. There is no additional benefit from having row or column shapes or metadata because the row/column nature is usually predetermined and would be represented by the shape of the corresponding matrix.

The benefit is that they can be indexed with a scalar and converted to a 1-d array with r.A, and no reshaping. Also that indexing a matrix reduces its rank, which is expected in a lot of places.

It really would be annoying to find that for an n by 1 vector/matrix, the product of X.T*X is a scalar or an error rather than an n by n matrix.

yes, it would, which is the whole point of the matrix object in the first place, and the point of the proposed row/column objects. They would provide 1-d object that behave as row an column vectors with * and **.

-Chris

On Tue, Apr 29, 2008 at 2:43 PM, Christopher Barker Chris.Barker@noaa.gov wrote:

Timothy Hochberg wrote:

...

However, there is matrix related stuff that is at best poorly supported now, namely operations on stacks of arrays (or vectors).

Tim, this is important, but also appears to be an orthogonal issue to me -- whatever we do with matrices, rows, columns, whatever, we still need to solve this problem, and any of the schemes being proposed strikes me as amenable to having a "array or matrices", and the LA functions that act on them. What am I missing?

Ah, that's a very good question. I'll do my best to answer it, it's been a long time since I dug into this seriously, so I may make some missteps: holler if you see something odd. Basically, I believe that you need two things:

1. The ability to make matrix/row/column objects of arbitrary dimensionality, where the extra dimensions contain the matrices/rows/columns. Current matrix objects are always 2D. 2. A way to tell whether a given object represents a row, column or matrix. For instance, given a 2D object is it matrix or a stack of vectors. You might think that there's enough information implicit in the shape of the arrays in conjunction with the nature of the function called that this wouldn't be necessary, but it turns out not to be so. As I recall, I had to fudge things when I was trying to containerize linalg.solve, because there wasn't enough information available to figure out unambiguously the kinds of the arguments.

Any of the solutions that attach information to the array, such as your matrix/row/col approach, or the quasi-tensor approach that I (and someone else) have mentioned but never detailed, should work for this.

However, any of the solutions that involve trying to encode the information into the array/matrix shape will probably not work. At least not well. Fortunately, they seem to be losing traction.

2008/4/30 Charles R Harris charlesr.harris@gmail.com:

Some operations on stacks of small matrices are easy to get, for instance, +,-,*,/, and matrix multiply. The last is the interesting one. If A and B are stacks of matrices with the same number of dimensions with the matrices stored in the last two indices, then

sum(A[...,:,:,newaxis]*B[...,newaxis,:,:], axis=-2)

is the matrix-wise multiplication of the two stacks. If B is replaced by a stack of 1D vectors, x, it is even simpler:

sum(A[...,:,:]*x[...,newaxis,:], axis=-1)

This doesn't go through BLAS, but for large stacks of small matrices it might be even faster than BLAS because BLAS is kinda slow for small matrices.

Yes and no. For the first operation, you have to create a temporary that is larger than either of the two input arrays. These invisible (potentially) gigantic temporaries are the sort of thing that puzzle users when as their problem size grows they suddenly find they hit a massive slowdown because it starts swapping to disk, and then a failure because the temporary can't be allocated. This is one reason we have dot() and tensordot() even though they can be expressed like this. (The other is of course that it lets us use optimized BLAS.)

This rather misses the point of Timothy Hochberg's suggestion (as I understood it), though: yes, you can write the basic operations in numpy, in a more or less efficient fashion. But it would be very valuable for arrays to have some kind of metadata that let them keep track of which dimensions represented simple array storage and which represented components of a linear algebra object. Such metadata could make it possible to use, say, dot() as if it were a binary ufunc taking two matrices. That is, you could feed it two arrays of matrices, which it would broadcast to the same shape if necessary, and then it would compute the elementwise matrix product.

The question I have is, what is the right mathematical model for describing these arrays-some-of-whose-dimensions-represent-linear-algebra-objects?

One idea is for each dimension to be flagged as one of "replication", "vector", or "covector". A column vector might then be a rank-1 vector array, a row vector might be a rank-1 covector array, a linear operator might be a rank-2 object with one covector and one vector dimension, a bilinear form might be a rank-2 object with two covector dimensions. Dimensions designed for holding repetitions would be flagged as such, so that (for example) an image might be an array of shape (N,M,3) of types ("replication","replication","vector"); then to apply a color-conversion matrix one would simply use dot() (or "*" I suppose). without too much concern for which index was which. The problem is, while this formalism sounds good to me, with a background in differential geometry, if you only ever work in spaces with a canonical metric, the distinction between vector and covector may seem peculiar and be unhelpful.

Implementing such a thing need not be too difficult: start with a new subclass of ndarray which keeps a tuple of dimension types. Come up with an adequate set of operations on them, and implement them in terms of numpy's functions, taking advantage of the extra information about each axis. A few operations that spring to mind:

* Addition: it doesn't make sense to add vectors and covectors; raise an exception. Otherwise addition is always elementwise anyway. (How hard should addition work to match up corresponding dimensions?) * Multiplication: elementwise across "repetition" axes, it combines vector axes with corresponding covector axes to get some kind of generalized matrix product. (How is "corresponding" defined?) * Division: mostly doesn't make sense unless you have an array of scalars (I suppose it could compute matrix inverses?) * Exponentiation: very limited (though I suppose matrix powers could be implemented if the shapes are right) * Change of basis: this one is tricky because not all dimensions need come from the same vector space * Broadcasting: the rules may become a bit byzantine... * Dimension metadata fiddling

Is this a useful abstraction? It seems like one might run into trouble when dealing with arrays whose dimensions represent vectors from unrelated spaces.

Anne

On Wed, Apr 30, 2008 at 9:16 PM, Anne Archibald peridot.faceted@gmail.com wrote:

2008/4/30 Charles R Harris charlesr.harris@gmail.com:

...

Some operations on stacks of small matrices are easy to get, for

instance,

...

+,-,*,/, and matrix multiply. The last is the interesting one. If A and

B

...

are stacks of matrices with the same number of dimensions with the

matrices

...

stored in the last two indices, then

sum(A[...,:,:,newaxis]*B[...,newaxis,:,:], axis=-2)

is the matrix-wise multiplication of the two stacks. If B is replaced by

a

...

stack of 1D vectors, x, it is even simpler:

sum(A[...,:,:]*x[...,newaxis,:], axis=-1)

This doesn't go through BLAS, but for large stacks of small matrices it might be even faster than BLAS because BLAS is kinda slow for small matrices.

Yes and no. For the first operation, you have to create a temporary that is larger than either of the two input arrays. These invisible (potentially) gigantic temporaries are the sort of thing that puzzle users when as their problem size grows they suddenly find they hit a massive slowdown because it starts swapping to disk, and then a failure because the temporary can't be allocated. This is one reason we have dot() and tensordot() even though they can be expressed like this. (The other is of course that it lets us use optimized BLAS.)

But it is interesting that you can multiply stacks of matrices that way, is it not? I haven't seen it mentioned elsewhere.

Chuck

On Tue, Apr 29, 2008 at 12:50 PM, Keith Goodman kwgoodman@gmail.com wrote:

On Tue, Apr 29, 2008 at 11:46 AM, Keith Goodman kwgoodman@gmail.com wrote:

...

On Tue, Apr 29, 2008 at 11:21 AM, Alan G Isaac aisaac@american.edu wrote:

...

On Tue, 29 Apr 2008, Keith Goodman apparently wrote:

...

I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.

1. Could you give an example of the circumstances of this

use?

In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:

array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix

list index idx = [0] x[idx,:] --> matrix?

scalar index idx = 0 x[idx,:] --> not matrix

Here's another use case:

for i in xrange(x.shape[1]):
    xi = x[:,i]
    idx = M.where(M.isfinite(xi).A)[0]
    xi = xi[idx,:]
    # more code

Also the functions I have written work on matrices. If x[i,:] did not return a matrix, let's say it returned some other vector type object that you index with a scalar, then my functions will have to test whether the input is a full matrix or a vector since the indexing is different.

And here's a use case that currently doesn't work but would be nice (for me) if it did:

x[idx,i] = M.rand(2,1)

where x is a matrix, say 4x3, idx is an array with 2 elements, and i is a scalar

Keith Goodman wrote:

It would seem odd to me if x[i,:] returned different types of objects based on the type of i:

I think this is where matrices should act like arrays:

...

...

a

array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])

...

...

a[:,[1]]

array([[1], [4], [7]])

...

...

a[:,1]

array([1, 4, 7])

in this case, it's not really different types, but totally different concepts -- indexing with a sequence is different than a scalar, just like indexing with a slice is different than a scalar.

There are many, many places in python code where we have to make distinctions between a single object and a sequence of objects that happens to be of length one.

-Chris

On Tue, Apr 29, 2008 at 2:18 PM, Anne Archibald peridot.faceted@gmail.com wrote:

On 29/04/2008, Keith Goodman kwgoodman@gmail.com wrote:

...

In my use i is most commonly an array (i = M.where(y.A)[0] where y is a nx1 matrix), sometimes a list, and in ipython when debugging or first writing the code, a scalar. It would seem odd to me if x[i,:] returned different types of objects based on the type of i:

array index idx = M.where(y.A)[0] where y is a nx1 matrix x[dx,:] --> matrix

list index idx = [0] x[idx,:] --> matrix?

scalar index idx = 0 x[idx,:] --> not matrix

It is actually pretty unreasonable to hope that

A[0]

and

A[[1,2,3]] or A[[True,False,True]]

should return objects of the same rank.

Why it unreasonable to hope that

x[0,:]

and

x[0, [1,2,3]]

or

x[0, [True,False,True]]

where x is a matrix, continue to return matrices?

I will answer this question too, because I feel like Keith.

On Tue, Apr 29, 2008 at 02:21:34PM -0400, Alan G Isaac wrote:

On Tue, 29 Apr 2008, Keith Goodman apparently wrote:

...

I often use x[i,:] and x[:,i] where x is a matrix and i is a scalar. I hope this continues to return a matrix.

1. Could you give an example of the circumstances of this

use?

x[i, :]*A*x[:, i]

I find this cleaner, because more explicit, than x[i]*A*x[:, i]. I'd be _very very_ upset if this broke. Breaking the explicit for the implicit is a _bad_ idea, and if you go down this path you go down the path of inconsistency.

2. Would any or all of the following be just as good

a result? a. 1d matrix

Yes.

    b. row and column vectors (1d but "oriented" for 
    linear algebra purposes)

Yes.

3. If 2a. and 2b. are no good for you, a. would e.g. ``x[i:i+1,:]`` be too burdensome?

No way. This is hard to read and the best way to make long calculations impossible to follow and debug.

    b. would ``rows`` and ``columns`` attributes that
    yielded matrics be an offsetting convenience?

Hum, I don't love this because it introduces noise in the formulas, but at least it is explicit and readable.

I would really like us to go down the way of the row/vector columns. This is consistent at least.

Gaël

On Tue, Apr 29, 2008 at 03:20:02PM -0400, Alan G Isaac wrote:

On Tue, 29 Apr 2008, Gael Varoquaux apparently wrote:

...

x[i, :]*A*x[:, i] I find this cleaner, because more explicit, than x[i]*A*x[:, i].

...

...

    a. would e.g. ``x[i:i+1,:]``

...

No way.

So *if* ``x[0]`` is going to be 1d to fix the many anomalies that have surfaced, then you would still want ``x[0]==x[0,:]``, but to have those be 1d and have ``x[i, :]*A*x[:, i]`` work would require (oriented) "vector" objects, so you would want those.

Is that a fair summary of your view?

That an extremely good summary of my views. Better than I could have formulated it. Thank you.

Given this good short paragraph summing up the problem, as I see it, this is why I think the only correct option is the oriented 1D objects. The second best option is not to change anything in order not to break backward compatibility.

Gaël

2008/4/25 Alan G Isaac aisaac@american.edu:

1. This is **not** what I understood as the agreement

(and I think the current solution is bad).

Reverted in r5084.

Cheers Stéfan

On Fri, Apr 25, 2008 at 12:02 PM, Alan G Isaac aisaac@american.edu wrote:

On Fri, 25 Apr 2008, Stéfan van der Walt apparently wrote:

...

Reverted in r5084.

Thank you.

I think we have discovered that there is a basic conflict between two behaviors:

   x[0] == x[0,:]
   vs.
    x[0][0] == x[0,0]

To my recollection, everyone has agree that the second behavior is desirable as a *basic expectation* about the behavior of 2d objects. This implies that *eventually* we will have ``x[0][0] == x[0,0]``. But then eventually we MUST eventually have x[0] != x[0,:].

Choices:

1) x[0] equiv x[0:1,:] -- current behavior 2) x[0] equiv x[0,:] -- array behavior, gives x[0][0] == x[0,0]

These are incompatible. I think 1) is more convenient in the matrix domain and should be kept, while 2) should be given up. What is the reason for wanting 2)? Maybe we could overload the () operator so that x(0) gives 2)? Or vice versa.

Chuck

Alan G Isaac wrote:

Please return 1d arrays in response to scalar indexing as the provisional arrangement.

+1 (for the provisional solution)

This has clarified it a bit for me. The current situation allows one to create "row vectors" and "column vectors" by generating matrices (in various ways) that happen to be shape (1,n) or (n,1). These objects behave as we all want in arithmetic operations (*, **). So what's missing:

***** The ability to index one of these "vector" objects by a single index, and get a scalar as a result. *****

If I'm not mistaken (and I probably am) -- that is the crux of this entire conversation.

Do we want that? I would say yes, it's a pretty natural and common thing to want.

Note that with the current version of matrix (numpy 1.1.0rc1), you can index a (n,1) matrix with a scalar, and get a (1,1) matrix, but you get an exception if you try to index a (1,n) matrix with a scalar:

...

...

rv

matrix([[3, 4, 5]])

...

...

rv[1]

Traceback (most recent call last): File "<stdin>", line 1, in <module> File "/Library/Frameworks/Python.framework/Versions/2.5/lib/python2.5/site-packages/numpy/core/defmatrix.py", line 228, in __getitem__ out = N.ndarray.__getitem__(self, index) IndexError: index out of bounds

I know why that's happening, but I think it's a demonstration of why we might want some sort of "vector" object, for a cleaner API.

How do we achieve this?

One possibility is to special case (n,1) and (1,n) matrices, but I think there is consensus that that is a BAD IDEA (tm)

Which leaves us with needing a "vector" object -- which acts a lot like a 1-d array, except that it acts like a (n,1) or (1,n) matrix with * and **, and may get printed differently.

I think there are two ways to accomplish this: (1) RowVector or ColumnVector objects that are essentially 1-d ndarrays, with a few methods overridden.

(2) A Vector object that is a matrix of either (n,1) or (1,n) in shape, with a couple methods overridden.

I think the difference is an implementation detail. They should behave about the same either way.

In either case, we should probably have factory functions for these objects: RowVector() and ColumnVector()

By the way, I think this brings us closer to the goal of matrices "acting as much like arrays as possible" -- when you index into a 2-d array: A[:,i], you get a 1-d array. These "vectors" would act a lot like 1-d arrays, so that would be more similar than the current situation.

It also might be nice to have iterators over rows and columns, like:

for row in M.rows: ... and for col in M.columns: ...

We can get the rows as the default iterator, but it feels nice and symmetric to be able to iterate over columns too (though I have no idea how useful that would really be)

Stéfan van der Walt wrote:

The whole idea of matrices was that you always work with 2-dimensional arrays, even when you just extract a row. Unless you have a proper hierarchical container (as illustrated in the other patch I sent), returning a 1D array breaks at least some expectation.

I agree -- stopgap is OK (but having people do M.A[i] is a fine stopgap too), but if we really want this, we need some sort of "vector" object.

If that's what the matrix users want, though, then we should change it.

I still want to see some linear algebra examples -- that's what matrices are for. We can iterate over 2-d arrays and get 1-d arrays just fine now, if that's what you want.

Unless we agree on something, and soon, that won't happen. Your workaround would break x[0] == x[0,:], so we're just swapping one set of broken functionality for another.

Alan G Isaac wrote:

But then we MUST eventually have x[0] != x[0,:].

I don't think so -- I think a Vector object would allow both of:

M[i,j] == M[i][j] and M[i] == M[i,:]

however,

M[i] != M[:,i] -- M[i] can only mean one thing, and this doesn't hold for arrays, either.

and (if we do the iterators):

M.rows[i] == M[i,:] M.columns[i] == M[:,i]

You can get rows a bit more conveniently than columns, which is really just an accident of how C arranges memory, but why not? Again, to make this "pure" we'd disallow M[i], but practicality beats purity, after all.

how would you index into a vector (as in http://en.wikipedia.org/wiki/Vector_(spatial)) without it?

Right -- which is why a 1-d "vector" object has its uses.

Alan G Isaac wrote:

Since a current behavior must disappear eventually, we should make it disappear as soon as possible: before the release. The question is how. I see two simple ways to move forward immediately:

    1. scalar indexing of matrices returns a 1d array
    2. scalar indexing of matrices raises a TypeError

I don't agree -- what we have is what we have, not changing it will not break any existing code, so we should only make a change if it moves us closer to what we want in the future -- it does seem that most folks do want M[i] to return some sort of 1-d object, so I'm fine with (1). I'm also fine with (2), but I'm pretty sure that one was laughed out of the discussion a good while back!

-Chris

On Fri, Apr 25, 2008 at 2:57 PM, Alan G Isaac aisaac@american.edu wrote:

On Fri, 25 Apr 2008, Christopher Barker apparently wrote:

...

I think a Vector object would allow both of: M[i,j] == M[i][j] and M[i] == M[i,:]

The problem is that it would be a crime to give up the natural production of submatrices. The NATURAL RULE is: to get a submatrix, use nonscalar indices. We should *not* give up that x[0,:] is a sub*matrix* whose first element is x[0,0] and equivalently x[0][0]. *This* is why we must have x[0]!=x[0,:] if we want, as we do, that x[0][0]==x[0,0].

I would give it up and let x(0) be a submatrix, the normal indexing is then exactly like an array. True, the result of x(0) can't be used as an lvalue (only getitem equivalence), but otherwise it should work fine.

Chuck

OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?

Those who have spoke up have generally favored letting x[0] return a 1d array, if I count correctly. And I think that is the right "provisional" behavior. (As well as ultimate.) But many have not spoken up.

So let's at least signal that this is an area of change. Right?

Thank you, Alan Isaac

On Sat, Apr 26, 2008 at 10:27 AM, Charles R Harris charlesr.harris@gmail.com wrote:

On Fri, Apr 25, 2008 at 6:38 PM, Alan G Isaac aisaac@american.edu wrote:

...

OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?

Those who have spoke up have generally favored letting x[0] return a 1d array, if I count correctly. And I think that is the right "provisional" behavior. (As well as ultimate.) But many have not spoken up.

That would be the most compatible with arrays, but won't it break a lot of code?

Any change to the behavior of x[0] for matrices is going to break a lot of code. It actually seems like a good idea to me to make it an error for a while -- or maybe some kind of settable option -- to help people figure out where they need to fix their code. Better to get an error immediately.

--bb

Alan G Isaac wrote:

OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?

I think this is wise for the time being.

At this point, I'm leaning in the direction of the RowVector / ColumnVector approach (even if these are not really advertised and just used during indexing).

-Travis

On Sat, Apr 26, 2008 at 7:12 AM, Alan G Isaac aisaac@american.edu wrote:

On Fri, 25 Apr 2008, "Travis E. Oliphant" apparently wrote:

...

At this point, I'm leaning in the direction of the RowVector / ColumnVector approach (even if these are not really advertised and just used during indexing).

I believe that this conflicts with submatrix extraction. Details follow.

Suppose ``x`` is a matrix, and I say ``v=A[0]``. What is ``v``?

1. If ``v`` is a 1d array, its behavior is known.

E.g., ``v[0]`` is the first element and ``v[0,0]`` is an IndexError. If we want to keep submatrix extraction (as we should), we give up ``x[0] == x[0,:]``.

2. If ``v`` is a matrix, its behavior can be learned

by the experienced but breaks basic expectations (the x[0][0] problem).

3. If ``v`` is a "RowVector" what behavior do we get?

Suppose the idea is that this will rescue ``x[0][0]==x[0,0]`` **and** ``x[0]=x[0,:]``. But then we must give up that ``x[0,:]`` is a submatrix. Must ``v`` deviate from submatrix behavior in an important way? Yes: ``v[0][0]`` is an IndexError. (Note that the same problem arises if we just override ``__getitem__`` for the matrix class to special case row and column matrices.)

Since submatrix extraction is fundamental, I think it is a *very* bad idea to give it up. If so, then under the RowVector proposal we must give up ``x[0]==x[0,:]``. But this is just what we give up with the much simpler proposal that ``v`` be a 1d array.

I may have missed something since I've been getting these emails all out of order for some reason and it's been hard to follow. Can you clarify what you mean by submatrix extraction? It sounds like you want to be able index into an MxN array and get out a 1xN or Mx1 matrix. If that's the case, wouldn't the natural way to spell that under the RowVector/ColumnVector approach (and most others that I can think of) be:

m2 = m1[:1]

or

m2 = m1[:,:1]

Here we're extracting the first row and column respectively as a matrix. This doesn't appear to have any of the conflicts with row/vector extraction and in general meshes much better with arrays than having some magic associated with appending an extra ':', which seems ill advised to me.

Alan G Isaac wrote:

...

OK, we are not converging in time for the release. So can we at least raise a TypeError on scalar indexing of matrices, so that we remain free to choose the ultimate behavior?

On Fri, 25 Apr 2008, "Travis E. Oliphant" apparently wrote:

I think this is wise for the time being.

In response to Bill's concern, I chose a warning instead: URL:http://scipy.org/scipy/numpy/attachment/ticket/760

fwiw, Alan

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

I claim b is more important than a. IMHO, a is plain wrong: you should't be indexing x with x[0][0].

Why?? Would you say this about a 2d array? Why the difference?

The core argument has been that it is a **basic expectation** of the behavior of 2d array-like objects that you will be able to get, e.g., the first element with x[0][0]. (E.g., lists, tuples, arrays ... is there an exception?? It should not be broken!)

I teach with matrices and I thought you might too: if so, you surely have run into this expectation (which is **natural**). In fact I truly doubt anyone has not been puzzled on first encounter by this::

>>> x matrix([[1, 2], [3, 4]]) >>> x[0] matrix([[1, 2]]) >>> x[0][0] matrix([[1, 2]])

Another argument, which I agree with, has been that matrix behavior should match 2d array behavior except where there is a clear payoff from deviation.

What I do agree with, which I take to be implicit in your statement, is that if users of matrices should generally use multiple indexes. But that does not answer the question about what we should offer when x[0] is requested. (Or, more crucially, what you get when you iterate over a matrix.)

Cheers, Alan Isaac

PS If you have good arguments, please add them to URL:http://www.scipy.org/MatrixIndexing

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

For me this is wrong. list and tuples are not 2D. Numpy arrays happen to offer this feature, but you should not use it do to multiple dimension indexing.

But there is no proposal that people should index like this. The underlying issue is iteration. Currently::

>>> A matrix([[1, 2], [ 3, 4]]) >>> A = np.mat(x) >>> for row in A: ... for col in row: ... print col ... [[1 2]] [[3 4]]

So are you saying that one should not be able to iterate through to the elements of a matrix?

So many natural expectations are broken by the current behavior!

Cheers, Alan

On Sat, 26 Apr 2008, Gael Varoquaux apparently wrote:

for row in A.rows(): for col in row.cols()

Actually I am in favor of adding ``rows`` and ``cols`` attributes to allow such iteration. The only thing is, I say these should be matrices (i.e., 2d). I.e., this should provide a symmetric syntax (for rows and columns) for iteration across submatrices.

But once we have these, it is all the more reason not to object to letting iteration on the matrix yield 1d arrays. As others have observed: why break array behavior for no purpose?

Cheers, Alan Isaac

PS The possibility of adding these attributes is mentioned on the discussion page: URL:http://www.scipy.org/MatrixIndexing as are most of the issues that have been coming up again.

Gael Varoquaux wrote:

• to use a syntax similar to dictionnaries:

  for row in A.rows(): for col in row.cols()

  I actually think this is much better than the code you currently use,

• or implement row and column objects.

The problem in your code is that you do not distinguish iterating over rows and columns, and in linear algebra these are different objects. The two solutions I propose do this distinction. My favorite one is the first one (the rows and cols methods).

What would rows() and columns() return? I'd like them to return a vector object. If they return matrices, then you can't index them with a single value, if they return 1-d arrays, then you can't distinguish between rows and columns (which may be OK).

Travis E. Oliphant wrote:

I'm not personally persuaded by the iteration argument, because we can change iteration independently of mapping (__getitem__) access.

or just use:

for row in M.A: ...

How hard is that?

-Chris

On Sat, 26 Apr 2008, "Travis E. Oliphant" apparently wrote:

I'm not personally persuaded by the iteration argument, because we can change iteration independently of mapping (__getitem__) access.

Would that not impose another deviation from array behavior? I just do not see the purpose.

I support .row and .cols iterators (returning submatrices), and have included the idea for some time on the discussion page URL:http://www.scipy.org/MatrixIndexing.

Cheers, Alan

Alan G Isaac wrote:

On Thu, 24 Apr 2008, Christopher Barker apparently wrote:

...

I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.

This bothers me.

So, if X is 2 by 2, then X[0] will be a row vector. But if X is 1 by 2, then X[0] will be a scalar? Ouch! Bye bye generic code.

exactly one of the things that really bugged me with MATLAB!

Or are you still having separate classes and not just changing ``__getitem__``? OK, I see that is most likely what you meant.

yes, it would still be a vector class, so a 1 by 2 matrix (or array) would still act like a 2-d object.

The only issue here is if you need a separate "row vector" and "column vector" class -- I think it's a matter of what whether these vectors are special versions of 1-d arrays or 2-d matrices. I kind of like them as special versions of 1-d arrays, as a vector is a 1-d object, but I see that the implementation is a bit cleaner if you just make them 2-d matrices with a different scalar indexing.

-Chris

Christopher Barker wrote: [cjop]

Bruce Southey wrote:

...

Hi, I would like to use the matrix functionality but I, like others, get by without it.

What does it lack that keeps you from using it? That's the key question?

I mainly do matrix-matrix and matrix-vector multiplications like in linear models and predictions. The main reason I have already learnt how to do to what I want using Numpy via mainly Numarray. So at the moment I am not interested in doing something else.

...

+1 to Tim's and Nadav's comments. As Tim said, there should be seamless integration between concepts of vectors and matrices - after all there really no distinction between them in linear algebra.

This is all about being able to index a "vector" with a single index -- that's it, I think. Matlab handles this by special casing matrices that have a single dimension of one. It also has an easy way to spell a literal for column vector, though I suppose:

np.Matrix((1,2,3,4,5)).T

isn't so bad.

...

-2 for having rowvector and columnvector - both numpy and I should know the orientation, if not, I deserve garbage or an error.

But HOW can you know the orientation? a 1-d array has no orientation --

I fully agree and which is why there is absolutely no sense in saying row or column vector/array! But vectors always have an orientation in the context of a matrix (as you also indicate below) or 2-d array so in that situation vectors are matrices. Numpy and other options have an implicit orientation because they usually always check that an operation is feasible first so errors are returned if the dimensions are not appropriate. So, as I understand it, Numpy's array do have an implicit orientation because Numpy's 1-d arrays act like 1 by n arrays when used with m by n arrays: import numpy a=numpy.ones((4)) b=numpy.ones((4,3)) numpy.dot(a,b) # gives array([ 4., 4., 4.]) numpy.dot(b,a) # gives error numpy.dot(b.T,a) #gives array([ 4., 4., 4.])

Thus, Numpy's 1-d arrays do have an orientation otherwise all three multiplications should return an error!

which is why the current version always returns a matrix when indexing.

As you point out, that is inconsistent with 1-d arrays having NO orientation. But this is fully consistent with the view that a vector is a special case of a matrix when we commonly refer to an n by 1 or 1 by n vector (as you also indicate below).

...

0 for vector - I don't see the need for it as a unique entity but understand the simplicity of saying vector.

I don't think we can have vector without ColumnVector, though I suppose a "Vector" can be either a (n,1) or a (1,n) matrix that allows single indexing.

But didn't you say that '1-d array has no orientation'... :-) This why I am -2 on row and column vector because one should just be the transpose of the other.

By the way, maybe a column vector could be handy for doing array broadcasting, as an cleaner way to do:

x = np.arange(10) y = np.arange(20).reshape((-1,1))

z = FunctionOfXandY(x,y)

Final note:

If no one has an linear algebra examples, then we're wasting out time with even this cheap talk.

-Chris

Bruce