(use-modules (system foreign) 
	     (rnrs bytevectors))

(define libgsl (dynamic-link "libgsl"))
(define libgslcblas (dynamic-link "libgslcblas"))

;=== Misc. utility functions ========================================
; (make-guard) creates a guard closure that collects objects
; that have memory stored outside the garbage collector with their
; destructors. if the guard is called without arguments, all
; objects in this guard are destroyed.

(define (apply-list l)
  (if (null? l)
    '()
    (begin (apply (caar l) (cdar l))
	   (apply-list (cdr l)))))

(define (make-guard)
  (let ((collection '()))
    (case-lambda 
      (() (apply-list collection))
      ((destructor obj) (set! collection 
			  (cons (list destructor obj) collection))))))
;====================================================================

;=== matrices =======================================================
(define gsl-matrix-alloc
  (pointer->procedure '*
		      (dynamic-func "gsl_matrix_alloc" libgsl)
		      (list size_t size_t)))

(define gsl-matrix-free
  (pointer->procedure void
		      (dynamic-func "gsl_matrix_free" libgsl)
		      (list '*)))

(define (make-gsl-matrix n m guard)
  (let ((m (gsl-matrix-alloc n m)))
    (guard gsl-matrix-free m)
    m))

(define (gsl-matrix->array M)
  (let* ((raw (parse-c-struct M (list size_t size_t size_t '* '* int)))
	 (n   (car raw))
	 (m   (cadr raw))
	 (tda (caddr raw))
	 (ptr (cadddr raw)))
    (make-shared-array
      (pointer->bytevector ptr (* tda n) 0 'f64)
      (lambda x (list (+ (* tda (car x)) (cadr x))))
      n m)))
;====================================================================

;=== vectors ========================================================
(define gsl-vector-alloc
  (pointer->procedure '*
		      (dynamic-func "gsl_vector_alloc" libgsl)
		      (list size_t)))

(define gsl-vector-free
  (pointer->procedure void
		      (dynamic-func "gsl_vector_free" libgsl)
		      (list '*)))

(define (gsl-vector->array v)
  (let* ((raw    (parse-c-struct v (list size_t size_t '* '* int)))
	 (size   (car raw))
	 (stride (cadr raw))
	 (ptr    (caddr raw)))
    (pointer->bytevector ptr size 0 'f64)))

(define (make-gsl-vector n guard)
  (let ((v (gsl-vector-alloc n)))
    (guard gsl-vector-free v)
    v))
;====================================================================

;=== permutations ===================================================
(define permutation-alloc
  (pointer->procedure '*
		      (dynamic-func "gsl_permutation_alloc" libgsl)
		      (list size_t)))

(define permutation-free
  (pointer->procedure void
		      (dynamic-func "gsl_permutation_free" libgsl)
		      (list '*)))

(define (make-permutation n guard)
  (let ((p (permutation-alloc n)))
    (guard permutation-free p)
    p))
;====================================================================

;=== linear algebra =================================================

;=== this could be easier, need a c-pointer to a single int ===
(define (dynamic-int) 
  (make-bytevector (sizeof int)))
(define (dynamic-int->pointer di) 
  (bytevector->pointer di))
(define (dynamic-int->int di) 
  (bytevector-sint-ref di 0 (native-endianness) (sizeof int)))
;======

(define linalg-LU-decomp
  (pointer->procedure int
		      (dynamic-func "gsl_linalg_LU_decomp" libgsl)
		      (list '* '* '*)))

(define linalg-LU-solve
  (pointer->procedure int
		      (dynamic-func "gsl_linalg_LU_solve" libgsl)
		      (list '* '* '* '*)))

(define linalg-LU-det
  (pointer->procedure double
		      (dynamic-func "gsl_linalg_LU_det" libgsl)
		      (list '* int)))

(define (make-linalg-det n guard)
  (let ((perm (make-permutation n guard))
	(sgn (dynamic-int)))
    (lambda (m)
      (linalg-LU-decomp m perm (dynamic-int->pointer sgn))
      (linalg-LU-det m (dynamic-int->int sgn)))))

(define (make-linalg-solve n guard)
  (let ((perm (make-permutation n guard))
	(sgn (dynamic-int)))
    (lambda (A b x)
      (linalg-LU-decomp A perm (dynamic-int->pointer sgn))
      (linalg-LU-solve A perm b x))))

;====================================================================
(define (randomize-timer)
  (let ((time (gettimeofday)))
    (set! *random-state*
      (seed->random-state (+ (car time)	
			     (cdr time))))))
(let* ((guard (make-guard))
       (b (make-gsl-vector 5 guard))
       (b-vec (gsl-vector->array b))
       (x (make-gsl-vector 5 guard))
       (x-vec (gsl-vector->array b))
       (A (make-gsl-matrix 5 5 guard))
       (A-vec (gsl-matrix->array A))
       (det (make-linalg-det 5 guard))
       (solve (make-linalg-solve 5 guard)))

  (randomize-timer)
  (array-index-map! (array-contents A-vec) (lambda X (random:normal)))
  (array-index-map! (array-contents b-vec) (lambda X (random:normal)))

  (display A-vec) (newline)
  (display b-vec) (newline)
  (display "=================") (newline)
  (display (det A)) (newline)

  (solve A b x)
  (display x-vec) (newline)

  (guard))