On 05/14/2012 10:36 PM, Dag Sverre Seljebotn wrote:
> On 05/14/2012 06:31 PM, mark florisson wrote:
>> On 12 May 2012 22:55, Dag Sverre Seljebotn<d.s.seljeb...@astro.uio.no>
>> wrote:
>>> On 05/11/2012 03:37 PM, mark florisson wrote:
>>>>
>>>> On 11 May 2012 12:13, Dag Sverre Seljebotn<d.s.seljeb...@astro.uio.no>
>>>> wrote:
>>>>>
>>>>> (NumPy devs: I know, I get too many ideas. But this time I *really*
>>>>> believe
>>>>> in it, I think this is going to be *huge*. And if Mark F. likes it it's
>>>>> not
>>>>> going to be without manpower; and as his mentor I'd pitch in too here and
>>>>> there.)
>>>>>
>>>>> (Mark F.: I believe this is *very* relevant to your GSoC. I certainly
>>>>> don't
>>>>> want to micro-manage your GSoC, just have your take.)
>>>>>
>>>>> Travis, thank you very much for those good words in the "NA-mask
>>>>> interactions..." thread. It put most of my concerns away. If anybody is
>>>>> leaning towards for opaqueness because of its OOP purity, I want to refer
>>>>> to
>>>>> C++ and its walled-garden of ideological purity -- it has, what, 3-4
>>>>> different OOP array libraries, neither of which is able to out-compete
>>>>> the
>>>>> other. Meanwhile the rest of the world happily cooperates using pointers,
>>>>> strides, CSR and CSC.
>>>>>
>>>>> Now, there are limits to what you can do with strides and pointers.
>>>>> Noone's
>>>>> denying the need for more. In my mind that's an API where you can do
>>>>> fetch_block and put_block of cache-sized, N-dimensional blocks on an
>>>>> array;
>>>>> but it might be something slightly different.
>>>>>
>>>>> Here's what I'm asking: DO NOT simply keep extending ndarray and the
>>>>> NumPy C
>>>>> API to deal with this issue.
>>>>>
>>>>> What we need is duck-typing/polymorphism at the C level. If you keep
>>>>> extending ndarray and the NumPy C API, what we'll have is a one-to-many
>>>>> relationship: One provider of array technology, multiple consumers (with
>>>>> hooks, I'm sure, but all implementations of the hook concept in the NumPy
>>>>> world I've seen so far are a total disaster!).
>>>>>
>>>>> What I think we need instead is something like PEP 3118 for the
>>>>> "abstract"
>>>>> array that is only available block-wise with getters and setters. On the
>>>>> Cython list we've decided that what we want for CEP 1000 (for boxing
>>>>> callbacks etc.) is to extend PyTypeObject with our own fields; we could
>>>>> create CEP 1001 to solve this issue and make any Python object an
>>>>> exporter
>>>>> of "block-getter/setter-arrays" (better name needed).
>>>>>
>>>>> What would be exported is (of course) a simple vtable:
>>>>>
>>>>> typedef struct {
>>>>> int (*get_block)(void *ctx, ssize_t *upper_left, ssize_t
>>>>> *lower_right,
>>>>> ...);
>>>>> ...
>>>>> } block_getter_setter_array_vtable;
>>>>>
>>>>> Let's please discuss the details *after* the fundamentals. But the reason
>>>>> I
>>>>> put void* there instead of PyObject* is that I hope this could be used
>>>>> beyond the Python world (say, Python<->Julia); the void* would be handed
>>>>> to
>>>>> you at the time you receive the vtable (however we handle that).
>>>>
>>>>
>>>> I suppose it would also be useful to have some way of predicting the
>>>> output format polymorphically for the caller. E.g. dense *
>>>> block_diagonal results in block diagonal, but dense + block_diagonal
>>>> results in dense, etc. It might be useful for the caller to know
>>>> whether it needs to allocate a sparse, dense or block-structured
>>>> array. Or maybe the polymorphic function could even do the allocation.
>>>> This needs to happen recursively of course, to avoid intermediate
>>>> temporaries. The compiler could easily handle that, and so could numpy
>>>> when it gets lazy evaluation.
>>>
>>>
>>> Ah. But that depends too on the computation to be performed too; a)
>>> elementwise, b) axis-wise reductions, c) linear algebra...
>>>
>>> In my oomatrix code (please don't look at it, it's shameful) I do this using
>>> multiple dispatch.
>>>
>>> I'd rather ignore this for as long as we can, only implementing "a[:] = ..."
>>> -- I can't see how decisions here would trickle down to the API that's used
>>> in the kernel, it's more like a pre-phase, and better treated orthogonally.
>>>
>>>
>>>> I think if the heavy lifting of allocating output arrays and exporting
>>>> these arrays work in numpy, then support in Cython could use that (I
>>>> can already hear certain people object to more complicated array stuff
>>>> in Cython :). Even better here would be an external project that each
>>>> our projects could use (I still think the nditer sorting functionality
>>>> of arrays should be numpy-agnostic and externally available).
>>>
>>>
>>> I agree with the separate project idea. It's trivial for NumPy to
>>> incorporate that as one of its methods for exporting arrays, and I don't
>>> think it makes sense to either build it into Cython, or outright depend on
>>> NumPy.
>>>
>>> Here's what I'd like (working title: NumBridge?).
>>>
>>> - Mission: Be the "double* + shape + strides" in a world where that is no
>>> longer enough, by providing tight, focused APIs/ABIs that are usable across
>>> C/Fortran/Python.
>>>
>>> I basically want something I can quickly acquire from a NumPy array, then
>>> pass it into my C code without dragging along all the cruft that I don't
>>> need.
>>>
>>> - Written in pure C + specs, usable without Python
>>>
>>> - PEP 3118 "done right", basically semi-standardize the internal Cython
>>> memoryview ABI and get something that's passable on stack
>>>
>>> - Get block get/put API
>>>
>>> - Iterator APIs
>>>
>>> - Utility code for exporters and clients (iteration code, axis
>>> reordering,
>>> etc.)
>>>
>>> Is the scope of that insane, or is it at least worth a shot to see how bad
>>> it is? Beyond figuring out a small subset that can be done first, and
>>> whether performance considerations must be taken or not, there's two
>>> complicating factors: Pluggable dtypes, memory management. Perhaps you could
>>> come to Oslo for a couple of days to brainstorm...
>>>
>>> Dag
>>
>> The blocks are a good idea, but I think fairly complicated for
>> complicated matrix layouts. It would be nice to have something that is
>> reasonably efficient for at least most of the array storage
>> mechanisms.
>> I'm going to do a little brain dump below, let's see if anything is useful :)
>>
>> What if we basically take the CSR format and make it a little simpler,
>> easier to handle, and better suited for other layouts. Basically, keep
>> shape/strides for dense arrays, and for blocked arrays just "extend"
>> your number of dimensions, i.e. a 2D blocked array becomes a 4D array,
>> something like this:
>>
>>>>> a = np.arange(4).repeat(4).reshape(4, 4);
>>>>> a
>> array([[0, 0, 0, 0],
>> [1, 1, 1, 1],
>> [2, 2, 2, 2],
>> [3, 3, 3, 3]])
>>
>>>>> a.shape = (2, 2, 2, 2)
>>>>> itemsize = a.dtype.itemsize
>>>>> a.strides = (8 * itemsize, 2 * itemsize, 4 * itemsize, 1 * itemsize)
>>>>> a
>> array([[[[0, 0],
>> [1, 1]],
>>
>> [[0, 0],
>> [1, 1]]],
>>
>>
>> [[[2, 2],
>> [3, 3]],
>>
>> [[2, 2],
>> [3, 3]]]])
>>
>>>>> print a.flatten()
>> [0 0 1 1 0 0 1 1 2 2 3 3 2 2 3 3]
>>
>> Now, given some buffer flag (PyBUF_Sparse or something), use basically
>> suboffsets with indirect dimensions, where only ever the last
>> dimension is a row of contiguous memory (the entire thing may be
>> contiguous, but the point is that you don't care). This row may
>>
>> - be variable sized
>> - have either a single "column start" (e.g. diagonal, band/tri-
>> diagonal, block diagonal, etc), OR
>> - a list of column indices, the length of the row (like in the CSR
>> format)
>>
>> The length of each innermost row is then determined by looking at, in order:
>> - the extent as specified in the shape list
>> - if -1, and some flag is set, the extent is determined like CSR,
>> i.e. (uintptr_t) row[i + 1] - (uintptr_t) row[i]
>> -> maybe instead of pointers indices are better, for
>> serialization, GPUs, etc
>> - otherwise, use either a function pointer or perhaps a list of extents
>>
>> All these details will obviously be abstracted, allowing for easy
>> iteration, but it can also be used by ufuncs operating on contiguous
>> rows (often, since the actual storage is contiguous and stored in a 1D
>> array, some flag could indicate contiguity as an optimization for
>> unary ufuncs and flat iteration). Tiled nditer-ation could also work
>> here without too big a hassle.
>> When you slice, you add to the suboffset and manipulate the single
>> extent or entire list of extents in that dimension, and the flag to
>> determine the length using the pointer subtraction should be cleared
>> (this should probably all happen through vtable functions).
>>
>> An exporter would also be free to use different malloced pointers,
>> allowing variable sized array support with append/pop functionality in
>> multiple dimensions (if there are no active buffer views).
>>
>> Random access in the case where a column start is provided is still
>> contant time, and done though a[i][j][k - colstart], unless you have
>> discontiguous rows, in which case you are allowed a logarithmic search
>> (if the size exceeds some threshold). I see scipy.sparse does a linear
>> search, which is pretty slow in pathological cases:
>>
>> from scipy import sparse
>> a = np.random.random((1, 4000000))
>> b = sparse.csr_matrix(a)
>> %timeit a[0, 1000]
>> %timeit b[0, 1000]
>>
>> 10000000 loops, best of 3: 190 ns per loop
>> 10 loops, best of 3: 29.3 ms per loop
>
> Heh. That is *extremely* pathological though, nobody does that in real
> code :-)
>
> Here's an idea I had yesterday: To get an ND sparse array with good
> spatial locality (for local operations etc.), you could map the elements
> to a volume-filling fractal and then use a hash-map with linear probing.
> I bet it either doesn't work or has been done already :-)
>
>>
>> Now, depending on the density and operation, the caller may have some
>> idea of how to allocate an output array. I'm not sure how to handle
>> "insertions" of new elements, but I presume through vtable put/insert
>> functions. I'm also not sure how to fit this in with linear algebra
>> functionality, other than providing conversions of the view.
>>
>> I think a disadvantage of this scheme is that you can't reorder your
>> axes anymore, and many operations that are easy in the dense case are
>> suddenly harder, e.g. this scheme does not allow you to go from a
>> csr-like format into csc. But I think what this gives is reasonable
>> generality to allow easy use in C/Fortran/Cython compiled code/numpy
>> ufunc invocation, as well as allowing efficient-ish storage for
>> various kinds of arrays.
>>
>> Does this make any sense?
>
> I'll admit I didn't go through the finer details of your idea; let's
> deal with the below first and then I can re-read your post later.
>
> What I'm thinking is that this seems interesting, but perhaps it lowers
> the versatility so much that it's not really worth to consider, for the
> GSoC at least. If the impact isn't high enough, my hunch is that one may
> as well not bother and just do strided arrays.
>
> I actually believe that the *likely* outcome of this discussion is to
> stick to the original plan and focus on expressions with strided arrays.
> But let's see.
>
> I'm not sure if what brought me to the blocked case was really
> block-sparse arrays or diagonal arrays. Let's see...:
>
> 1) Memory conservation. The array would not be NOT element-sparse, it's
> just that you don't want to hold it all in memory at once.
>
> Examples:
>
> - "a_array *= np.random.normal(size=a_array.shape)". The right hand
> side can be generated on the fly (making it return the same data for
> each block every time is non-trivial but can be done). If a takes up a
> good chunk of RAM, there might not even be enough memory for the
> right-hand-side except for block-by-block. (That saves bandwidth too.)
>
> - You want to stream from one file to another file, neither of which
> will fit in RAM. Here we really don't care about speed, just code
> reuse...it's irritating to have to manually block in such a situation.
Actually, s/file/database/g, just to avoid comments about np.memmap.
Dag
>
> - You want to play with a new fancy array format (a blocked format
> that's faster for some linear algebra, say). But then you need to call
> another C library that takes data+shape+stride+dtype. Then you need to
> make a copy, which you'd rather avoid -- so an alternative is to make
> that C library be based on the block API so that it can interfaces with
> your fancy new format (and anything else).
>
> 2) Bandwidth conservation. Numexpr, Theano, Numba and
> Cython-vectorized-expressions all already deal with this problem on one
> level.
>
> However, I also believe there's a higher level in many programs where
> blocks come into play. The organization of many scientific codes
> essentially reads "do A on 8 GB of data, then do B on 8 GB of data"; and
> it's going to be a *long* time before you have full-program analysis to
> block up that in every case; the tools above will be able to deal with
> some almost-trivial cases only.
>
> A block API could be used to write "pipeline programs". For instance,
> consider
>
> a_array[:, None] * f(x_array)
>
> and f is some rather complicated routine in Fortran that is NOT a ufunc
> -- it takes all of x_array as input and provides all the output "at
> once"; but at some point it needs to do the writing of the output, and
> if writing to an API it could do the multiplication with a_array while
> the block is in cache anyway and save a memory bus round-trip. (Provided
> the output isn't needed as scratch space.)
>
> Summing up:
>
> Vectorized expressions on contiguous (or strided) memory in Cython is
> pretty nice by itself; it would bring us up to the current state of the
> art of static compilation (in Fortran compilers).
>
> But my hunch is your sparse-block proposal doesn't add enough
> flexibility to be worth the pain. Many of the cases above can't be
> covered with it.
>
> If it's possible to identify a little nugget API that is
> forward-compatible with the above usecases (it wouldn't solve them
> perhaps, but allow them to be solved with some extra supporting code),
> then it might be worth to a) implement it in NumPy, b) support it for
> Cython vectorized expressions; and use that to support block-transposition.
>
> But I'm starting to feel that we can't really know how that little
> nugget API should really look until the above has been explored in a
> little more depth; and then we are easily talking a too large scope
> without tying it to a larger project (which can't really before the
> GSoC..). For instance, 2) suggests push/pull rather than put/get.
>
> Dag
> _______________________________________________
> NumPy-Discussion mailing list
> NumPy-Discussion@scipy.org
> http://mail.scipy.org/mailman/listinfo/numpy-discussion
_______________________________________________
NumPy-Discussion mailing list
NumPy-Discussion@scipy.org
http://mail.scipy.org/mailman/listinfo/numpy-discussion
