Re: Matrix or array operations library

[Daniel Llorens <daniel.llorens@bluewin.ch>]

> From: Matt Wette <matt.wette@gmail.com>
> Subject: Re: Matrix or array operations library
> Date: 27 December 2018 at 23:24:51 CET
> To: guile-user@gnu.org
> 
> On 12/27/18 10:43 AM, Daniel Llorens wrote:
>> 1) Guile already has a multidimensional array type with strides and all that. Has had it forever. The base library (array-map! etc) is low level and doesn't do much, and the implementation has some cruft, but the type is there, it's fairly simple stuff and and there's nothing fundamentally wrong with it. Should be safe too.
> 
> Daniel,
> 
> Can you explain how array strides works?  I didn't read it the same as in ndarray.
> 
> BTW, ndarray has been swallowed into Python 3 as updated def'n of buffer objects,
> if I read correctly.
> 
> Matt

Hi,

John already wrote a nice explanation of how array strides work. You mean specifically for Guile?

Guile's array type is a record with a reference to storage, a list of bounds, and a list of strides. When you look up an element in an array, the strides are multiplied by the indices to obtain a linear address in the storage. You can manipulate the strides to create another view into an existing array.

This should all be explained in the manual.

Guile shouldn't be doing this any differently from Fortran / APL / Python / etc. It's the only sensible way to do it (there are other ways, like using pointers to arrays of pointers, which aren't sensible for this kind of array). Guile arrays are type generic, so the array record also caches set! and ref functions appropriate for the type of the storage. Compared to Python there's the quirk that Guile supports not only an upper bound per dimension but also a lower bound. I don't remember that ndarray does that, because otherwise negative indices wouldn't work. Honestly I think that lower bounds are a mess and I'd love to remove them, but Fortran has them... When you operate on whole arrays (as with array-map!) the bounds are checked before entering any loop, so this doesn't make arrays slower. IIRC ndarray also keeps a bunch of flags to cache whether an array is in compact C or Fortran order, etc. Guile doesn't have those and I don't think they are necessary. Guile used to have a flag to cache if an array was in compact C order, but it was a source of bugs and I removed it.

(I'm afraid I don't follow Python development much, so I apologize if my understanding of ndarray is outdated — my only contact with it was when I had to call a C++ library from Python a few years ago.)

Of course all those other languages have extensive (and more or less consistent) array facilities, while base Guile only has make-shared-array (and transpose-array, for some reason). This seems to be out of some principle of minimalism, not that I agree with it especially. You can implement every higher level stride operation using make-shared-array (including transpose-array), and this is what my library guile-ploy does. It's not an ideal situation, it doesn't get you Fortran speed in Guile or anything of the sort, but the problem is not with the array type, other than the fact that the strides can only be set through make-shared-array.

For what is worth I think the array type needs to be moved to Scheme, but this should be done in a backwards compatible way. I started guile-newra (https://notabug.org/lloda/guile-newra or https://github.com/lloda/guile-newra) with that purpose. On 2.9 some functions are faster than Guile's built in ones, but most are slower, so there is still a long way to go. There's also a TON of low hanging fruit in Guile's array library implementation, like have a look at how array-copy! deals with typed arrays. But I haven't wanted to work on that much, since a lot of it would probably need to be replaced if the array type was moved to Scheme. Or maybe not and I'm just lazy.

I think I said this in my other comment, but in case that went unread, the speed of array-ref, make-shared-array, etc. doesn't matter if all you do with arrays is pass pointers and strides to C or C++ or Fortran. In that regard, Guile arrays are as good as anything else.

	Daniel