Matrix or array operations library

Matrix or array operations library

[Zelphir Kaltstahl]

Hello Guile Users,

Is there some library, that enables high performance matrix operations
or even n-dimensional array operations? I am thinking of something like
NumPy in the Python ecosystem. I think NumPy in turn also uses some
lower level thing to do what it does. I think OpenBLAS or MKL, depending
on the architecture. I wonder if there is any wrapper around OpenBLAS
for Guile or something similar.

I am writing a little code for matrix operations and currently I am
using Guile arrays, as they are made of vectors and have constant access
time, which is already great. My guess is, that this would be the right
choice if using pure Guile. I am writing data abstraction procedures, so
that later on I could exchange what is used to represent the data.

Maybe, if there is something like NumPy or lower level, I should use
that instead? (Would I have to learn how to use FFI first?)

Or maybe Guile's implementation is already so fast, that it would not
make that much difference to use a lower level thing?

Currently I have only a little experimental program, started today, so
no huge plan. OK, one can fantasize about stuff like Pandas data frames
in Guile, but I have no illusion, that it is a work of a few days or
even weeks. It would be nice to learn, how to use a low level thing or
maybe even Pandas, if there are any such bindings for Guile. I could
make the implementation use different representations, depending on a
parameter or something like that.

Regards,

Zelphir

Re: Matrix or array operations library

[Tk]

‐‐‐‐‐‐‐ Original Message ‐‐‐‐‐‐‐
On Monday, 24 December 2018 23:01, Zelphir Kaltstahl <zelphirkaltstahl@gmail.com> wrote:

> Hello Guile Users,
>
> Is there some library, that enables high performance matrix operations
> or even n-dimensional array operations? I am thinking of something like
> NumPy in the Python ecosystem. I think NumPy in turn also uses some
> lower level thing to do what it does. I think OpenBLAS or MKL, depending
> on the architecture. I wonder if there is any wrapper around OpenBLAS
> for Guile or something similar.
>
> I am writing a little code for matrix operations and currently I am
> using Guile arrays, as they are made of vectors and have constant access
> time, which is already great. My guess is, that this would be the right
> choice if using pure Guile. I am writing data abstraction procedures, so
> that later on I could exchange what is used to represent the data.
>
> Maybe, if there is something like NumPy or lower level, I should use
> that instead? (Would I have to learn how to use FFI first?)
>
> Or maybe Guile's implementation is already so fast, that it would not
> make that much difference to use a lower level thing?
>
> Currently I have only a little experimental program, started today, so
> no huge plan. OK, one can fantasize about stuff like Pandas data frames
> in Guile, but I have no illusion, that it is a work of a few days or
> even weeks. It would be nice to learn, how to use a low level thing or
> maybe even Pandas, if there are any such bindings for Guile. I could
> make the implementation use different representations, depending on a
> parameter or something like that.
>
> Regards,
>
> Zelphir


I took a different route. Instead of merely binding functionality in a lower level language, I use Guile to generate lean-and-mean modern Fortran code. The building blocks can be found here: https://gitlab.com/codetk/schemetran (wanted to save this for potluck, but here it goes ... ). Fortran compilers take care of efficient execution on HPC platforms. Actually, it is beyond me why anyone would bother with any other programming language when it comes to expressing maths efficiently.

I built a pseudo-spectral Navier-stokes solver that can work on MPI, shared-mem (OpenMP), and hopefully soon GPU/Xeon accelerators (OpenMP 4) atop of schemetran. I still need to see about publishing it under an Libre licence.

Re: Matrix or array operations library

[Matt Wette]

On 12/24/18 2:01 PM, Zelphir Kaltstahl wrote:

> Hello Guile Users,
>
> Is there some library, that enables high performance matrix operations
> or even n-dimensional array operations? I am thinking of something like
> NumPy in the Python ecosystem. I think NumPy in turn also uses some
> lower level thing to do what it does. I think OpenBLAS or MKL, depending
> on the architecture. I wonder if there is any wrapper around OpenBLAS
> for Guile or something similar.
>
> I am writing a little code for matrix operations and currently I am
> using Guile arrays, as they are made of vectors and have constant access
> time, which is already great. My guess is, that this would be the right
> choice if using pure Guile. I am writing data abstraction procedures, so
> that later on I could exchange what is used to represent the data.
>
> Maybe, if there is something like NumPy or lower level, I should use
> that instead? (Would I have to learn how to use FFI first?)
>
> Or maybe Guile's implementation is already so fast, that it would not
> make that much difference to use a lower level thing?
>
> Currently I have only a little experimental program, started today, so
> no huge plan. OK, one can fantasize about stuff like Pandas data frames
> in Guile, but I have no illusion, that it is a work of a few days or
> even weeks. It would be nice to learn, how to use a low level thing or
> maybe even Pandas, if there are any such bindings for Guile. I could
> make the implementation use different representations, depending on a
> parameter or something like that.
>
> Regards,
>
> Zelphir
>
IMO, Guile probably needs something like ndarray, the base for numpy,
if matrix operations are going to go anywhere.

As far as using blas or others, I have developed a package for helping
to use the FFI API for C libraries.  It may be of interest to you.

ffi-helper: https://www.nongnu.org/nyacc/ffi-help.html

Matt

ndarray [was: Matrix or array operations library]

[Matt Wette]

On 12/24/18 4:21 PM, Matt Wette wrote:

> IMO, Guile probably needs something like ndarray, the base for numpy,
> if matrix operations are going to go anywhere.

https://docs.scipy.org/doc/numpy-1.13.0/reference/internals.html

Re: Matrix or array operations library

[David Pirotte]

[-- Attachment #1: Type: text/plain, Size: 1479 bytes --]

Hello,

> Is there some library, that enables high performance matrix operations
> or even n-dimensional array operations? I am thinking of something like
> NumPy in the Python ecosystem. I think NumPy in turn also uses some
> lower level thing to do what it does. I think OpenBLAS or MKL, depending
> on the architecture. I wonder if there is any wrapper around OpenBLAS
> for Guile or something similar.
> ...

There are quite a few libs you might look at before you start to code your own:

	https://notabug.org/lloda
		...
		guile-newra
		guile-ffi-cblas
		...
		guile-ploy
		...

	AIscm
	http://wedesoft.github.io/aiscm

I also wrote a set of n-dim array math ops (including matrix ops) for Guile-CV:

	https://www.gnu.org/software/guile-cv/
		the manual is available online, math and matrix ops are described in
		the 'Processing' section

	mentioning Guile-CV because images are just a collection of n-dim arrays.
	(see 'Image Structure and Accessors' in the manual).

	in Guile-CV, all (almost all) arrays are f32vectors - see (srfi srfi-4) in
	the Guile manual - math and matrix ops are written in C (but memory
	allocation entirely kept in scheme). multi channel ops are multi-threaded

	some of these math and matrix ops also work on s32 and f64 (undocumented so
	far): it would not take much to actually fully generalize all guile-cv math
	and matrix ops for the full set of (srfi srfi-4 ) types.

Cheers,
David




	


[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 488 bytes --]

Re: Matrix or array operations library

[Matt Wette]

On 12/26/18 3:38 AM, Zelphir Kaltstahl wrote:
> On 25.12.18 18:00, guile-user-request@gnu.org wrote:
>> IMO, Guile probably needs something like ndarray, the base for numpy,
>> if matrix operations are going to go anywhere.
>>
>> As far as using blas or others, I have developed a package for helping
>> to use the FFI API for C libraries.  It may be of interest to you.
>>
>> ffi-helper: https://www.nongnu.org/nyacc/ffi-help.html
>>
>> Matt
> Hi!
>
> Thanks for the link Matt. Now I have a bunch of questions :-)
>
> I started reading the guide you linked to. I've never done any FFI stuff
> in any programming language so far, so this is completely new territory
> for me.
>
> (Q1) Can you explain what the advantage of the method you linked to is,
> in comparison to doing something like
> https://www.gnu.org/software/guile/manual/html_node/Dynamic-FFI.html ?

The advantage of the FFI Helper is time.  It takes minutes to generate
a ffi-module and seconds to compile to Guile Scheme.  It takes
significantly more time to code by hand.  While coding by hand gets you
something that is probably more palatable, using the FH can cover more
ground.

> (Q2) What is that `guild` command in the guide? (sounds like guildhall
> another Guile project, but I forgot what that does)

"guild" is the Guile compiler.  Try "guild compile foo.scm".

> (Q3) How would I set everything up, so that I can follow the tutorial at
> https://www.nongnu.org/nyacc/ffi-help.html?

download and read the README

* the tarball repository:
     https://download.savannah.gnu.org/releases/nyacc/

* the git repository:
     git://git.savannah.nongnu.org/nyacc.git

> iirc NYACC is a minimalistic compiler that was mentioned here on the
> mailing list several times in the context of bootstrapping Guile.
> (Q4) How is NYACC connected to doing FFI?

  It takes a decent C parser to read the include files.

> (Q5) Do I need NYACC installed on my system to do the stuff in the
> tutorial you linked?

Yes

> I guess, since I don't know things about FFI in general (except that I
> somehow get to magically call functions of another programming
> language), things like NumPy are out of my league for now and I should
> start with simple learning examples, to get an understanding of how FFI
> usually works and what things one touches, when using FFI for a library
> in another language and how much understanding of how the library works
> one needs to have.
>
> With regard to NumPy, especially the stuff about the "strides" (whatever
> that is) in the data buffer of a NumPy array and how what changes when
> creating a view and changing indices (For example I cannot imagine a
> single reason, why the ordering of how one writes the indices would
> matter, except for human readability?) … If I need to keep track of that
> … Let's just say I dislike working on a low level with a lot of
> pointers, addresses and segfaults ; )
>
> (Q6) Maybe I do not need to pay attention to how NumPy does things,
> except for procedure names, since NumPy handles it and I am just calling
> its procedures somehow? That's what I would be hoping for.

I'm not sure what you are asking.  If you want to work now in Guile try
some of the other packages offered.  To implement numpy would take some
work and I think the concept could go farther in Guile (e.g., some sort
of lazy eval scheme that can optimize evaluation of multi-operand
expressions, who knows).

> Maybe those are a lot of questions and maybe most of them would be
> answered, if I ever did any FFI stuff before. If anyone knows a good
> tutorial on FFI in Guile that explains this stuff, that would be great
> too. Or maybe someone knows some easy-to-ffi library to try things with?
>
> Best regards,
>
> Zelphir
>
>

Re: Matrix or array operations library

[Mike Gran]

On Thu, Dec 27, 2018 at 06:20:10AM -0800, Matt Wette wrote:
> > (Q1) Can you explain what the advantage of the method you linked to is,
> > in comparison to doing something like
> > https://www.gnu.org/software/guile/manual/html_node/Dynamic-FFI.html ?
> 
> The advantage of the FFI Helper is time.  It takes minutes to generate
> a ffi-module and seconds to compile to Guile Scheme.  It takes
> significantly more time to code by hand.  While coding by hand gets you
> something that is probably more palatable, using the FH can cover more
> ground.

As someone who has spent a lot of time coding bindings by hand, I can
suggest that one should always start with a tool like ffi-helper or
SWIG and then patch over any rough spots in scheme.  Unfortunately, some
C libraries are very uncooperative and do require C bindings, but I think
for most, ffi-helper and then a scheme library that cleans up the bindings
is the way to begin.

-Mike

Re: Matrix or array operations library

[Daniel Llorens]

>> With regard to NumPy, especially the stuff about the "strides" (whatever
>> that is) in the data buffer of a NumPy array and how what changes when
>> creating a view and changing indices (For example I cannot imagine a
>> single reason, why the ordering of how one writes the indices would
>> matter, except for human readability?) … If I need to keep track of that
>> … Let's just say I dislike working on a low level with a lot of
>> pointers, addresses and segfaults ; )

So I think there are three sides to this. 1) the array type, 2) using the array type in Scheme, 3) interfacing with C or Fortran (or whavever).

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.

People love to write their own matrix type on top of vectors or bytevectors or even lists. Sometimes there's a good reason for that but often there isn't.

2) IMO the only real flaw with the current Guile array type is that it's all in C so the array descriptors are opaque to Scheme, and the compiler isn't smart enough yet to remove the type dispatches that result from the base library being type generic [*]. This is an actual problem when you do everything in Scheme, but it's just the difference between things being slow or very slow.

AFAIK this isn't any worse than what you get in Numpy. If you operate with arrays ‘in the large’ (manipulating strides, etc.) then everything is fast. Of course Numpy gives you a ton of ready made functions that would be nice to have in Guile as well [**].

3) If you're using Guile arrays as an interface to C or Fortran, the above doesn't matter. Neither C nor Fortran (<90) have a standard way of passing multidimensional arrays to a function, so you are using pointers whether you like it or not. Every library has its own particular way of passing strides and its own quirks about what strides they accept, etc. You have to deal with that one way or another. Might as well have a nice, already available array type on the Scheme side!

[*] [**] I have started efforts to fix these two problems (both fizzled due to lack of time, although I think they contain some valid stuff — e.g. if you need to slice arrays, I prefer from and amend! in guile-ploy to the Numpy mess). I'd be happy to contribute to a common effort if anyone else is interested.

PS. Does anyone want to collaborate in wrapping HDF5? I've looked into mwette's FFI helper for like two seconds but I couldn't get it to work.

Re: Matrix or array operations library

[John Cowan]

>
> >> With regard to NumPy, especially the stuff about the "strides" (whatever
> >> that is)
>

Strides are an extremely powerful concept that allow transformations to be
done on arrays without actually moving or changing any of their data.  The
stride of a dimension is the distance in storage units (pointers, small
integers, floats, characters, whatever is in the array) between consecutive
indexes  along that dimension.

Suppose you have a 3 x 3 matrix:

0 1 2
3 4 5
6 7 8

If it is stored row-by-row (C order) then it will appear in storage as 0 1
2 3 4 5 6 7 8.  Suppose you are at location [0, 1], whose value is 1 (row
indices are written before column indices).  In order to get to [0, 2], the
next column in the same row, you have to move along just 1 element.  But to
get to  [1, 1],  the next row in the same column, you have to move along 3
elements.

So knowing the strides makes it easy to step along either rows or columns
fast.  But wait, there's more!  Suppose you create a new array object,
still 3 x 3, but with a row stride of 1 and a column stride of 3, but
sharing the same storage.  Hey presto, you now have transposed (rotated by
90 degrees) the original array, and it is now

0 3 6
1 4 7
2 5 8

In languages without strides, you'd have to copy all the elements.  Here,
no elements were copied in the making of this new array.

As another example, say you create yet another array sharing the same
storage.  This one has one dimension only and a stride of 4.  Its elements
are

0 4 8

and in fact it is the diagonal of the (original or transposed) array.
Again, no copying was done.  By playing tricks with upper and lower bounds,
strides, and the array offset (which is the storage address of point [0,
0]) you can translate a 0-based matrix to a 1-based matrix whose rows and
columns are numbered 1, 2, 3 instead of 0, 1, 2, or even one with row and
column numbers -1, 0, 1 if you want.  You can get any of the corners of the
array as 2 x 2 matrices.  Everything generalizes to any number of
dimensions, and there are many more possibilities: anything that can be
expressed as an affine transformation can be obtained in this way.

-- 
John Cowan          http://vrici.lojban.org/~cowan        cowan@ccil.org
There are books that are at once excellent and boring.  Those that at
once leap to the mind are Thoreau's Walden, Emerson's Essays, George
Eliot's Adam Bede, and Landor's Dialogues.  --Somerset Maugham

Re: Matrix or array operations library

[Matt Wette]

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

Re: Matrix or array operations library

[Daniel Llorens]

> 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

Re: Matrix or array operations library

[John Cowan]

On Fri, Dec 28, 2018 at 4:34 PM Daniel Llorens <daniel.llorens@bluewin.ch>
wrote:

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).


If array objects don't have an offset (the index in the backing store of
the [0,0, ... 0] element), you can't do arbitrary translations,
unfortunately.  If Guile doesn't have that, it should.


> 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.


Layering a Guile compatibility mode over SRFI 122 might be a Good Thing.
Arrays there allow arbitrary affine transformations of the indices, provide
lazy elementwise mapping (you can get eager mapping by composing mapping
with copying), and have fast paths for arrays of up to 4 dimensions.
Read-only and read-write arrays defined by arbitrary getter and setter
functions are also provided, and work exactly like storage-based arrays.
SRFI 122's only major limitation, which IMO is not a serious one, is that
it doesn't handle 0-dimensional arrays (with one element) or degenerate
arrays with non-positive dimensional ranges (with zero elements).  The code
is in Gambit Scheme, but translating it to portable Scheme is an easy
matter (I just haven't gotten around to it).  The main Gambit-specific
dependency is define-macro (non-hygienic) macros, but nothing very bad is
done with them.

There is a post-SRFI fork at https://github.com/gambiteer/srfi-122/ which
will fairly soon become an updated SRFI.

-- 
John Cowan          http://vrici.lojban.org/~cowan        cowan@ccil.org
Yakka foob mog.  Grug pubbawup zink wattoom gazork.  Chumble spuzz.
    --Calvin, giving Newton's First Law "in his own words"

Re: Matrix or array operations library

[Matt Wette]

On 12/28/18 12:23 PM, Daniel Llorens wrote:
>
> John already wrote a nice explanation of how array strides work. You mean specifically for Guile?
>
>
Yes, for Guile.  I understand strides and offsets.  I think the functionality is basically there
with shared-arrays.  So I take back that Guile does not have this.  I will play later.

For example, interesting check will be if one can make a shared-1D array from the diagonal elements
of a 2D square array (i.e., ((1 2 3) (4 5 6) (7 8 9)) => (1 5 9).

And the transpose I see as a utility that I guess just calls make-shared-array.

Matt

Re: Matrix or array operations library

[Daniel Llorens]

> On 29 Dec 2018, at 00:16, John Cowan <cowan@ccil.org> wrote:

> If array objects don't have an offset (the index in the backing store of the [0,0, ... 0] element), you can't do arbitrary translations, unfortunately.  If Guile doesn't have that, it should.

It does; I just forgot to mention it. make-shared-array accepts any affine function of indices. There's even a function to get at the offset directly  (shared-array-offset). All this is in the manual and has been in Guile for a long time.

> Layering a Guile compatibility mode over SRFI 122 might be a Good Thing.  Arrays there allow arbitrary affine transformations of the indices, provide lazy elementwise mapping (you can get eager mapping by composing mapping with copying), and have fast paths for arrays of up to 4 dimensions.  Read-only and read-write arrays defined by arbitrary getter and setter functions are also provided, and work exactly like storage-based arrays.  SRFI 122's only major limitation, which IMO is not a serious one, is that it doesn't handle 0-dimensional arrays (with one element) or degenerate arrays with non-positive dimensional ranges (with zero elements).  The code is in Gambit Scheme, but translating it to portable Scheme is an easy matter (I just haven't gotten around to it).  The main Gambit-specific dependency is define-macro (non-hygienic) macros, but nothing very bad is done with them.

Rank 0 arrays and empty arrays happen all the time if you work with arrays. I think it would be annoying to have to look out for special cases.

To be honest I find it strange that SRFI-122 doesn't support those you since you shouldn't need to do anything special (Guile supports them for free, other than a special case in the read syntax iirc). 

Guile only has the type and a few low level functions. My library guile-ploy provides arbitrary arity and arbitrary rank rank extension on top of Guile, plus a bunch of array functions, including full APL-style slicing. Still no lazy mapping, that would be nice to have. You probably want operator overloading as well.

One thing that is needed is a linear range generator of some kind, like a:b in Python. It's not enough to provide a procedure for this, since the slicing function must be able to recognize this linear range in order to avoid copying. For instance x[a:b], should produce another array with different strides, not a copy and not a lazy map either. In guile-ploy this is written (from x (range a b)), where (range a b) is a special type. In guile-newra this linear range is one of the basic vector types (along with vectors, bytevectors, and srfi-4 vectors). Does srfi-122 provide something like this? I suppose if you have lazy maps the linear range can be done with those and then the slicing function can be written to recognize it.

I shall try and have a look.

	Daniel


> There is a post-SRFI fork at https://github.com/gambiteer/srfi-122/ <https://github.com/gambiteer/srfi-122/> which will fairly soon become an updated SRFI.
> 
> -- 
> John Cowan          http://vrici.lojban.org/~cowan <http://vrici.lojban.org/~cowan>        cowan@ccil.org <mailto:cowan@ccil.org>
> Yakka foob mog.  Grug pubbawup zink wattoom gazork.  Chumble spuzz.
>     --Calvin, giving Newton's First Law "in his own words"
> 

Re: Matrix or array operations library

[Matt Wette]

So I tried an example using Guile to perform Gaussian elimination
(in situ, no column pivoting so unstable) using shared-arrays and
using explicit computation of offsets.  The explicit approach is
much more grungy but I'm guessing it does a lot less memory allocation.

Also, is it possible to get a pointer to the first element of an array
(i.e., it's #f64() root vector), for the purpose of using with FFI?

(use-modules ((srfi srfi-1) #:select (fold)))
(use-modules (srfi srfi-11))

;; === method 1: use share arrays

(define (axpy1 a x y)
   (let ((n (car (array-dimensions x))))
     (if (not (= n (car (array-dimensions y)))) (throw 'error))
     (let loop ((i 0))
       (when (< i n)
	(array-set! y (+ (* a (array-ref x i)) (array-ref y i)) i)
	(loop (1+ i))))))

;; modify A, b in-situ so that A is upper triangular
(define (reduce1 A b n)
   (let ((n-1 (1- n)))
     (let loop ((kl (iota n)))
       (when (pair? kl)
	(let* ((k (car kl))
	       (alpha (- (/ 1.0 (array-ref A k k))))
	       (x (make-shared-array A (lambda (j) (list k j)) n)))
	  (let loop1 ((i (1+ k)))
	    (when (< i n)
	      (let ((a (* alpha (array-ref A i k)))
		    (y (make-shared-array A (lambda (j) (list i j)) n)))
		(axpy1 a x y)
		(array-set! b (+ (* a (array-ref b k)) (array-ref b i)) i))
	      (loop1 (1+ i)))))
	(loop (cdr kl))))))
   
;; solve A*x = b for x where A is upper triangular
(define (bksolv1 A b x n)
   (let loop ((k (1- n)))
     (unless (negative? k)
       (let loop1 ((r 0.0) (j (1- n)))
	(cond
	 ((< k j)
	  (loop1 (+ r (* (array-ref A k j) (array-ref x j))) (1- j)))
	 (else
	  (array-set! x (/ (- (array-ref b k) r) (array-ref A k k)) k))))
       (loop (1- k)))))


;; === method 2: use root vector and computing offsets explicitly

;; based on root vectors
;; x, y are root vectors (from shared-array-root)
;; x0, y0 are starting indices
;; xi, yi are increments
(define (axpy2 n a x x0 xi y y0 yi)
   (let loop ((ix x0) (iy y0) (i 0))
     (when (< i n)
       (array-set! y (+ (* a (array-ref x ix)) (array-ref y iy)) iy)
       (loop (+ ix xi) (+ iy yi) (1+ i)))))

;; @deffn {Procedure} array-vec-model a vec-idx . idxs
;; Generate params for accessing a vector in an array.
;; Returned expression is
;; (values rv off inc)
;; where @code{rv} is the root-vector, @code{off} is the offset of the
;; first element and @code{inc} is the increment.
;; @end deffn
(define (array-vec-model a vec-idx . idxs)
   (let ((ic (shared-array-increments a)))
     (values (shared-array-root a)
	    (fold (lambda (bd ic ix of) (+ of (* ic (- ix (car bd)))))
		  (shared-array-offset a) (array-shape a) ic idxs)
	    (list-ref ic vec-idx))))

;; @deffn {Procedure} array-mat-model a row-idx col-idx . idxs
;; Generate params for accessing a vector in an array.
;; Returned expression is
;; (values rv off inc)
;; where @code{rv} is the root-vector, @code{off} is the offset of the
;; first element and @code{inc} is the increment.
;; @end deffn
(define (array-mat-model a row-idx col-idx . idxs)
   (let ((ic (shared-array-increments a)))
     (values (shared-array-root a)
	    (fold (lambda (bd ic ix of) (+ of (* ic (- ix (car bd)))))
		(shared-array-offset a) (array-shape a) ic idxs)
	    (list-ref ic row-idx)
	    (list-ref ic col-idx))))

;; modify A, b in-situ so that A is upper triangular
(define (reduce2 A b n)
   (let-values
       (((av a0 arinc acinc) (array-mat-model A 0 1 0 0))
        ((bv b0 binc) (array-vec-model b 0 0)))
     (let loop ((kl (iota n)) (akk a0) (bk b0))
       (when (pair? kl)
	(let* ((k (car kl))
	       (alpha (- (/ 1.0 (array-ref av akk)))))
	  (let loop1 ((i (1+ k)) (aik (+ akk arinc)) (bi (+ bk binc)))
	    (when (< i n)
	      (let ((a (* alpha (array-ref av aik))))
		(axpy2 (- n k) a av akk acinc av aik acinc)
		(array-set! bv (+ (* a (array-ref bv bk)) (array-ref bv bi)) bi))
	      (loop1 (1+ i) (+ aik arinc) (+ bi binc)))))
	(loop (cdr kl) (+ akk arinc acinc) (+ bk binc))))))

;; solve A*x = b for x where A is upper triangular
(define (bksolv2 A b x n)
   (let*-values
       (((av a00 arinc acinc) (array-mat-model A 0 1 0 0))
        ((bv b0 binc) (array-vec-model b 0 0))
        ((xv x0 xinc) (array-vec-model x 0 0))
        ((bn) (+ b0 (* (1- n) binc))))
     (let loop ((k (1- n))
	       (akn (+ a00 (* (1- n) arinc) (* (1- n) acinc)))
	       (xk (+ x0 (* (1- n) xinc))))
       (unless (negative? k)
	(let loop1 ((r 0.0) (j (1- n)) (akj akn) (bj bn))
	  (cond
	   ((< k j)
	    (loop1 (+ r (* (array-ref av akj) (array-ref x j)))
		   (1- j) (- akj acinc) (- bj binc)))
	   (else
	    (array-set! xv (/ (- (array-ref b j) r) (array-ref av akj)) xk))))
	(loop (1- k) (- akn arinc) (- xk xinc))))))


;; === example problem

(define A
   (list->typed-array
    'f64 2 '((11.0 1.2 1.3 1.4)
	    (2.1 22.0 2.3 2.4)
	    (3.1 3.2 33.0 3.4)
	    (4.1 4.2 4.3 44.0))))
(define b
   (list->typed-array
    'f64 1 '(4.0 3.0 2.0 1.0)))

(define x (make-typed-array 'f64 0.0 4))

;; expected result
(define x-soln
   (list->typed-array
    'f64 1 '(0.352856 0.103010 0.019728 -0.021913)))

(display "expect:\n") (vec-disp x-soln) (newline)

(when #f
   (reduce1 A b 4)
   ;;(display "A:\n") (mat-disp A)
   ;;(display "b:\n") (vec-disp b)
   (bksolv1 A b x 4)
   (display "method1:\n") (vec-disp x)
   )

(when #t
   (reduce2 A b 4)
   ;;(display "A:\n") (mat-disp A)
   ;;(display "b:\n") (vec-disp b)
   (bksolv2 A b x 4)
   (display "method2:\n") (vec-disp x)
   )

Re: Matrix or array operations library

[Linas Vepstas]

Me too.

With one significant difference: I have extremely sparse arrays. Like, only
one-in-a-million array entries are non-zero.  And my arrays are hug -- say
2M by 2M, for a total of 4 tera-entries, of which only one in a million are
non-zero, so in fact, my data might fit in a gigabyte or less. (instead of
32 terabytes of RAM)

I solved this by writing my own ad-hoc library. The library user needs to
declare a half-dozen methods, indicating how to fetch the matrix entries,
what the dimensions of the array are, and methods to loop over the non-zero
row and column entries.  The library does the rest.

My code is quasi-special-purpose, and, well, just right for me, but a
serious dedicated library programmer could do a much better job than I
can.  My code is here:

https://github.com/opencog/atomspace/tree/master/opencog/matrix

The low-level API that the user must declare is here:
https://github.com/opencog/atomspace/blob/master/opencog/matrix/object-api.scm

Hmm, maybe I should pull this out of the project, and make it generic....

-- Linas


On Mon, Dec 24, 2018 at 4:01 PM Zelphir Kaltstahl <
zelphirkaltstahl@gmail.com> wrote:

> Hello Guile Users,
>
> Is there some library, that enables high performance matrix operations
> or even n-dimensional array operations? I am thinking of something like
> NumPy in the Python ecosystem. I think NumPy in turn also uses some
> lower level thing to do what it does. I think OpenBLAS or MKL, depending
> on the architecture. I wonder if there is any wrapper around OpenBLAS
> for Guile or something similar.
>
> I am writing a little code for matrix operations and currently I am
> using Guile arrays, as they are made of vectors and have constant access
> time, which is already great. My guess is, that this would be the right
> choice if using pure Guile. I am writing data abstraction procedures, so
> that later on I could exchange what is used to represent the data.
>
> Maybe, if there is something like NumPy or lower level, I should use
> that instead? (Would I have to learn how to use FFI first?)
>
> Or maybe Guile's implementation is already so fast, that it would not
> make that much difference to use a lower level thing?
>
> Currently I have only a little experimental program, started today, so
> no huge plan. OK, one can fantasize about stuff like Pandas data frames
> in Guile, but I have no illusion, that it is a work of a few days or
> even weeks. It would be nice to learn, how to use a low level thing or
> maybe even Pandas, if there are any such bindings for Guile. I could
> make the implementation use different representations, depending on a
> parameter or something like that.
>
> Regards,
>
> Zelphir
>
>
>

-- 
cassette tapes - analog TV - film cameras - you