Re: matrix library?

[Damien Mattei]

that is:

(define-method (* (M1 <matrix>) (M2 <matrix>))

  {(n1 p1) <+ (dim-matrix M1)}
  {(n2 p2) <+ (dim-matrix M2)}

  (when {p1 ≠ n2} (error "matrix.* : matrix product impossible,
incompatible dimensions"))

  {v1 <+ (matrix-v M1)}
  {v2 <+ (matrix-v M2)}

  (define (res i j)
    {sum <+ 0}
    (for ({k <+ 0} {k < p1} {k <- k + 1})
     {sum <- sum + v1[i][k] * v2[k][j]}))

  {v <+ (create-vector-2d res n1 p2)}
  ;(display "v=") (display v) (newline)

  (matrix v))

and the full code:
a matrix multiplication (could add + - but no needs in my code...)
with GOOPS and Scheme+ , without using the overload features of
scheme+ but the native ones of guile scheme:

(define-module (matrix)

  #:use-module (oop goops)

  #:use-module ((guile)
        #:select (*)
        #:prefix orig:)

  ;;#:use-module ((srfi srfi-43)) ;; vector

  #:export (<matrix>
        matrix?
        matrix
        matrix-v
        create-matrix-by-function
        dim-matrix)

  )


(include "/Users/mattei/Dropbox/git/Scheme-PLUS-for-Guile/Scheme+.scm")

(include "/Users/mattei/Dropbox/git/Scheme-PLUS-for-Guile/scheme-infix.scm")

(include "/Users/mattei/Dropbox/git/Scheme-PLUS-for-Guile/assignment.scm")
(include 
"/Users/mattei/Dropbox/git/Scheme-PLUS-for-Guile/apply-square-brackets.scm")

(include "/Users/mattei/Dropbox/git/Scheme-PLUS-for-Guile/array.scm")
;; must know <- from assignment


(define-class <matrix> ()

  (v #:accessor matrix-v
     #:getter matrix-get-v
     #:setter matrix-set-v!
     #:init-keyword #:v)

  )


(define-method (matrix vinit)
  (make <matrix> #:v vinit))

(define-method (matrix? obj)
  (equal? (class-name (class-of obj)) (class-name <matrix>)))

;; (define M (create-matrix-by-function (lambda (l c) (+ l c)) 2 3))
;; $1 = #<<matrix> 1078a0a70>
(define (create-matrix-by-function fct lin col)
  (matrix (create-vector-2d fct lin col)))

;; return the line and column values of dimension
;; scheme@(guile-user)> (dim-matrix M)
;; $1 = 2
;; $2 = 3
(define (dim-matrix M)

  (when (not (matrix? M))
    (error "argument is not of type matrix"))

  {v <+ (matrix-v M)}
  {lin <+ (vector-length v)}
  {col <+ (vector-length {v[0]})}
  (values lin col))


;; scheme@(guile-user)> (define M1 (create-matrix-by-function (lambda
(l c) (+ l c)) 2 3))
;; scheme@(guile-user)> (define M2 (create-matrix-by-function (lambda
(l c) (- l c)) 3 2))
;; scheme@(guile-user)> (matrix-v M1)
;; $2 = #(#(0 1 2) #(1 2 3))
;; scheme@(guile-user)> (matrix-v M2)
;; $3 = #(#(0 -1) #(1 0) #(2 1))
;; (define M1*M2 {M1 * M2})
;; scheme@(guile-user)> (matrix-v M1*M2)
;; $2 = #(#(5 2) #(8 2))

(define-method (* (M1 <matrix>) (M2 <matrix>))

  {(n1 p1) <+ (dim-matrix M1)}
  {(n2 p2) <+ (dim-matrix M2)}

  (when {p1 ≠ n2} (error "matrix.* : matrix product impossible,
incompatible dimensions"))

  {v1 <+ (matrix-v M1)}
  {v2 <+ (matrix-v M2)}

  (define (res i j)
    {sum <+ 0}
    (for ({k <+ 0} {k < p1} {k <- k + 1})
     {sum <- sum + v1[i][k] * v2[k][j]}))

  {v <+ (create-vector-2d res n1 p2)}
  ;(display "v=") (display v) (newline)

  (matrix v))


and the test:

 (define M1 (create-matrix-by-function (lambda (l c) (+ l c)) 2 3))
;(define M2 (create-matrix-by-function (lambda (l c) (- l c)) 3 2))
(matrix-v M1)
 #(#(0 1 2) #(1 2 3))
(matrix-v M2)
 #(#(0 -1) #(1 0) #(2 1))
(define M1*M2 {M1 * M2})
 (matrix-v M1*M2)
 #(#(5 2) #(8 2))

additional function used:

;; create a vector of line and column with a function
(define (create-vector-2d fct lin col)
  {v <+ (make-vector lin)}
  ;;(display "ok") (newline)
  (for ({l <+ 0} {l < lin} {l <- l + 1})
       {v[l] <- (make-vector col)}
       (for ({c <+ 0} {c < col} {c <- c + 1})
        {v[l][c] <- (fct l c)}))
  v)

and a corrected version of :

(define-syntax for
   (lambda (stx)
     (syntax-case stx ()
       ((kwd (init test incrmt) body ...)
        (with-syntax ((BREAK (datum->syntax (syntax kwd) 'break))
                      (CONTINUE (datum->syntax (syntax kwd) 'continue)))
             (syntax
              (call/cc
               (lambda (escape)
             (let-syntax ((BREAK (identifier-syntax (escape))))
               init
               (let loop ((res 0)) ;; now we will return a result at
the end if no break but if we continue? what happens?
                 (if test
                 (begin
                   (call/cc
                    (lambda (next)
                      (set! res (let-syntax ((CONTINUE
(identifier-syntax (next))))
                          (let () ;; allow definitions
                            body ...)))))
                   incrmt
                   (loop res))
                 res)
                 ))))
              ) ;; close syntax

             )))))

note that expression

{sum <- sum + v1[i][k] * v2[k][j]}

is expanded as:

($nfx$ sum <- sum + ($bracket-apply$ ($bracket-apply$ v1 i) k) *
($bracket-apply$ ($bracket-apply$ v2 k) j))

which is not fast to evaluate,and slow down calculus on  intensive
computation, as the infix precedence procedure is evaluated each time
,even if it never change in the computation loop or never in program,
i'm working on an optimiser (written in scheme) of code before
evaluation , that modify the code (transform the infix expression with
precedence in prefix one definitively ) before it is passed to the
guile compiler.But it is another piece of code...

On Wed, Oct 4, 2023 at 10:29 PM Damien Mattei <damien.mattei@gmail.com> wrote:
>
> yes i will re-implement that from scratch because i can not find
> something easy in libraries to use.
> I already made it for Racket/Scheme+:
> (struct matrix-vect (v)) ;; matrix based on vector of vectors
>
> (define (create-matrix-vect-by-function fct lin col)
>   (matrix-vect (create-vector-2d fct lin col)))
>
>
> ;; return the line and column values of dimension
> (define (dim-matrix-vect M)
>
>   (when (not (matrix-vect? M))
>     (error "argument is not of type matrix"))
>
>   {v <+ (matrix-vect-v M)}
>   {lin <+ (vector-length v)}
>   {col <+ (vector-length {v[0]})}
>   (values lin col))
>
>
> (define (multiply-matrix-matrix M1 M2)
>
>   {(n1 p1) <+ (dim-matrix-vect M1)}
>   {(n2 p2) <+ (dim-matrix-vect M2)}
>
>   (when {p1 ≠ n2} (error "matrix-by-vectors.* : matrix product
> impossible, incompatible dimensions"))
>
>   {v1 <+ (matrix-vect-v M1)}
>   {v2 <+ (matrix-vect-v M2)}
>
>   (define (res i j)
>     (for/sum ([k (in-range p1)])
>          {v1[i][k] * v2[k][j]}))
>
>   {v <+ (create-vector-2d res n1 p2)}
>   ;(display "v=") (display v) (newline)
>
>   (matrix-vect v))
>
>
> (define (vector->matrix-column v)
>   (matrix-vect (vector-map (lambda (x) (make-vector 1 x))
>                v)))
>
> (define (matrix-column->vector Mc)
>   {v <+ (matrix-vect-v Mc)}
>   (vector-map (lambda (v2) {v2[0]})
>           v))
>
> (define (multiply-matrix-vector M v) ;; args: matrix ,vector ;  return vector
>   {Mc <+ (vector->matrix-column v)}
>   ;(display Mc)
>   (matrix-column->vector (multiply-matrix-matrix M Mc)))
>
>
>  there should be no difference for computation in Guile ,just that i
> will use GOOPS instead of struct and that it should be a guile module
> written in scheme+ syntax....
> see it when finished....
>
>
> On Wed, Oct 4, 2023 at 8:28 PM <tomas@tuxteam.de> wrote:
> >
> > On Wed, Oct 04, 2023 at 05:42:10PM +0200, Damien Mattei wrote:
> > > thank you. i will try to find a simple example in.
> > >
> > > On Wed, Oct 4, 2023 at 1:36 PM Nala Ginrut <nalaginrut@gmail.com> wrote:
> > > >
> > > > And I'd mention AIScm which bound tensorflow APIs. It needs more love.
> > > >
> > > > https://wedesoft.github.io/aiscm/
> > > >
> > > >
> > > > On Wed, Oct 4, 2023, 17:12 <tomas@tuxteam.de> wrote:
> > > >>
> > > >> On Wed, Oct 04, 2023 at 08:46:02AM +0200, Damien Mattei wrote:
> > > >> > hello,
> > > >> > does anyone know if it exists a basic matrix library for Guile?
> > > >> > just need to multiply a matrix and a vector.
> > > >>
> > > >> Perhaps this thread from guile-user [1] has something for you
> > > >>
> > > >> Cheers
> > > >>
> > > >> [1] 
> > > >> https://lists.gnu.org/archive/html/guile-user/2018-12/threads.html#00117
> >
> > That all said, if your case doesn't get much hairier,
> > it is fairly easy to navigate with Guile's arrays:
> >
> > (define (dot v1 v2)
> >   "The dot product of v1 and v2"
> >   (let ((sum 0))
> >     (array-for-each
> >       (lambda (x1 x2)
> >         (set! sum (+ sum (* x1 x2))))
> >       v1 v2)
> >     sum))
> >
> > (define (row m i)
> >   "The i-th row of m, as a vector"
> >   (make-shared-array m (lambda (j) (list i j)) (car (array-dimensions m))))
> >
> > (define (mat* m v)
> >   "Left multiply matrix m and vector v"
> >   (let* ((height (cadr (array-dimensions m)))
> >          (res (make-typed-array 'f64 0 height)))
> >     (do ((i 0 (1+ i)))
> >         ((>= i height))
> >         (array-set! res (dot (row m i) v) i))
> >     res))
> >
> >
> > (define vec #1f64(1 2 1))
> > (define mat #2f64((0.5 0.5 0) (0 0.5 0.5) (-0.5 0 0.5)))
> >
> > (mat* mat vec)
> >
> >  => (1.5 1.5 0)
> >
> > CAVEAT: only tested with this one example. Input value checking
> > left as an...
> >
> > Cheers
> > --
> > t
> >