[SciPy-dev] Sparse matrix design and broadcasting

Ed Schofield schofield at ftw.at
Sat Jul 8 12:49:40 CDT 2006

On 07/07/2006, at 7:53 PM, Neilen Marais wrote:

> Hi
>
> I've been looking at the scipy.sparse routines a little, specifically
> the sparse_lil type. Currently fancy indexing seems to be handled by a
> lot of special case code.
>
> I wonder if there is not some way we can use numpy's build in
> broadcasting to make this easier. My thinking is something like this:
>
> 1) lil_matrix.__setitem__() figures out what the equivalent shape
> of the
> region being assigned to is. Perhaps some numpy routines exist that
> can
> help here?
>
> e.g.
>
> lil_mat[1,:] -> (1, lil_mat.shape[1])
> lil_mat[[0,1,2], [5,6,7]] -> (3,)
> lil_mat[ix_([0,1,2], [5,6,7])] -> (3,3)
>
> 2) get numpy to broadcast the value being set to that shape
>
> 3) use a single code path that knows how to handle assignation from 2D
> or 1D arrays.
>
> Does this make sense?

Interesting idea. I wouldn't be surprised if we could re-use some
code from NumPy for this, or perhaps from numarray (whose indexing
code is written in Python). I went through the numarray code for
interpreting slices and index arrays before writing the sparse fancy
indexing code, hoping to re-use whatever I could, but it made quite
heavy use of strides, and I didn't see how this could be applied to
sparse matrices. So I more or less reinvented the wheel. But some
ideas from numarray (in Lib/generic.py) could probably be lifted,
such as the way the _universalIndexing method handles both getting
and setting.

Maybe we could indeed use functionality from NumPy in interpreting
indices, but I don't know whether this is possible or how we'd do it.
I'd appreciate any help.

-- Ed