Re: Matrix or array operations library

[Daniel Llorens]

> From: Matt Wette <address@hidden>
> Subject: Re: Matrix or array operations library
> Date: 27 December 2018 at 23:24:51 CET
> To: address@hidden
> 
> 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