[Numpy-discussion] matrix indexing question

[Bill Baxter]

On 3/27/07, Alan Isaac <aisaac at american.edu> wrote:
> > On 3/27/07, Alan Isaac <aisaac at american.edu> wrote:
> >> May I see a use case where the desired
> >> return when iterating through a matrix
> >> is rows as matrices?  That has never
> >> been what I wanted.
>
>
> On Tue, 27 Mar 2007, Bill Baxter wrote:
> > AllMyPoints = mat(rand(100,2)) # 100 two-d points
> > for pt in AllMyPoints:
> >     xformedPt = pt * TransformMatrix
> >     # do something to transformed point
>
>
> This seems artificial to me for several reasons,
> but here is one reason:
> AllxformedPts = AllMyPoints * TransformMatrix

Yeh, I was afraid you'd say that.  :-)
But maybe I've got a lot of points, and I don't feel like making a
copy of the whole set.
Or maybe it's not a linear transform, but instead
     xformedPt = someComplicatedNonLinearThing(pt)

I do stuff like the above quite frequently in my code, although with
arrays rather than matrices.  :-)
For instance in finite elements there's assembling the global
stiffness matrix step where for each node (point) in your mesh you set
some entries in the big matrix K.  Something like

for i,pt in enumerate(points):
    K[shape_fn_indices(i)] = stiffness_fn(pt)

That's cartoon code, but I think you get the idea.  I don't think
there's any good way to make that into one vectorized expression.  The
indices of K that get set depend on the connectivity of your mesh.

> Note that I am no longer disputing the convention,
> just trying to understand its usefulness.

--bb