;; @module gsl.lsp ;; @description Selected functions from the GNU Scientific Library ;; @version 1.0 - initial release. Minimum newLISP version is 10.4.0. ;; @version 1.1 - added check for extended ffi enabled version. ;; @author Lutz Mueller 2012 ;;

## Module GSL For The GNU Scientific Library

;; The @link http://www.gnu.org/software/gsl/ GSL GNU Scientific Library ;; implements over a 1000 functions from different subject areas. In this ;; release of the 'gsl.lsp' module only a few linear algebra functions ;; are implemented. ;; To use this module, the main GSL library 'libgsl' and a supporting ;; library 'libgslcblas' must be installed. ;; ;; This module requires newLISP version 10.4.0 as a minimum. For 64-bit ;; newLISP on Mac OSX, Linux and other Unix, 10.4.2 is the minimum. ;; ;; Precompiled 32-bit libraries for the binary distributions of newLISP ;; versions are available in the ;; @link http://www.nuevatec.com/GSL/ http://www.nuevatec.com/GSL/ ;; directory. See the 'INSTALL.txt' file in that directory for instructions how ;; to install the library files. ;; ;; The module contains (test-gsl) to test all implemented functions. ;; @syntax (gsl:CholeskyD ) ;; @param A square matrix of m row vectors with n = m column elements each. ;; @return The matrix L . ;; The function performs a Cholesky decomposition of . The matrix A must ;; be symmetric and positive definite square. ;; ;; A = L * Lt ;; ;; Lt is the transposition of L. ;; ;; @example ;; (gsl:CholeskyD '((4 2 -2) (2 10 2) (-2 2 5))) ;; ;; gsl:Cholesky-L => ;; ;; ( ( 2 0 0 ) ;; ( 1 3 0 ) ;; (-1 1 1.732050808) ) ;; ;; @syntax (gsl:Cholesky-solve ) ;; @param A square matrix of m row vectors with n = m column elements each. ;; @params A vector of n elements. ;; @return Vector x containing a solution for Ax = b . ;; ;; @example ;; (gsl:Cholesky-solve '((4 2 -2) (2 10 2) (-2 2 5)) '(1 2 3)) ;; ;; => (0.8333333333 -0.1666666667 1) ;; ;; @syntax (gsl:QRD ) ;; @param A list of m row vector lists with n column elements each. ;; @return Vector tau of k = min(n, m) elements. ;; The function does QR decomposition of a matrix A. ;; ;; A = Q * R ;; ;; The function works for both, square and rectangular matrices. ;; The number of rows m in A must be equal to or greater than the number of columns n. ;; The orthogonal m * m matrix Q can be found in the variable 'gsl:QR-Q'. ;; The m * n matrix R can be found in the variable 'gsl:QR-R'. ;; ;; @example ;; (gsl:QRD '((12 -51 4) (6 167 -68) (-4 24 -41)) ) => (1.857142857 1.993846154 0) ;; ;; gsl:QR-Q => ;; ;; ( ( 0.8571428571 -0.3942857143 -0.3314285714 ) ;; ( 0.4285714286 0.9028571429 0.03428571429) ;; (-0.2857142857 0.1714285714 -0.9428571429 ) ) ;; ;; gsl:QR-R => ;; ;; ( (14 21 -14 ) ;; ( 0 175 -70 ) ;; ( 0 0 35 ) ) ;; @syntax (gsl:QR-solve ) ;; @param A matrix of m row vectors with n column elements each. ;; @params A vector of n elements. ;; @return Vector x containing a solution for Ax = b . ;; Matrix A is either square or overdetermined with m > n, more rows than columns. ;; When m > n, then the variable 'gsl:QR-residual' contains a vector of residuals. ;; For a square matrix A this vector contains 0's. ;; ;; @example ;; (gsl:QR-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3)) ;; ;; => (0.009387755102 -0.02432653061 -0.08832653061) ;; ;; (gsl:QR-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4)) ;; => (9.690821045e-16 0.5) ;; ;; gsl:QR-residual ;; => (2.512815377e-17 1.103917396e-16 -2.961679405e-16 1.606480472e-16) ;; @syntax (gsl:SVD ) ;; @param A matrix of m row vectors with n column elements each. ;; @return A vector of diagonal elements from the S matrix. ;; The function does a SVD (Singular Value Decomposition) of the into ;; its components U, S and V. Matrix U is orthogonal m*n. S is a diagonal ;; square matrix of n singular values. V is a n*n square orthogonal matrix. ;; The function works for both, square and rectangular matrices. ;; ;; A = U * S * Vt ;; ;; Vt is the transposition of V. ;; ;; The number of rows m in A must be equal to or greater than the number of columns n. ;; The m*n matrix U can be found in the variable 'gsl:SVD-U'. ;; The diagonal elements of matrix S can be found in the vector variable 'gsl:SVD-S'. ;; The n*n matrix V can be found in the variable 'gsl:SVD-V'. ;; ;; @example ;; (gsl:SVD '((1 2) (3 4) (5 6) (7 8))) => (14.2690955 0.6268282324) ;; ;; gsl:SVD-U => ;; ( (-0.1524832333 -0.8226474722 ) ;; (-0.3499183718 -0.4213752877 ) ;; (-0.5473535103 -0.02010310314) ;; (-0.7447886488 0.3811690814 ) ) ;; ;; gsl:SVD-S => ;; ;; (14.2690955 0.6268282324) ;; ;; gsl:SVD-V => ;; ;; ( (-0.641423028 0.7671873951) ;; (-0.7671873951 -0.641423028 ) ) ;; @syntax (gsl:SVD-solve ) ;; @param A matrix of m row vectors with n column elements each. ;; @params A vector of n elements ;; @return a vector x containing a solution for Ax = b . ;; The number of rows m in A must be equal to or greater than the number n of columns. ;; ;; @example ;; (gsl:SV-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4)) ;; ;; => (5.551115123e-16 0.5 0 0) ;; ; on Mac OSX when using 32-bit newLISP from the normal binary Mac OSX installer ; use the following to make the libraries: ; ; ./configure CFLAGS="-O2 -m32" ; make ; sudo make install ; ; install needs the root password ; (when (!= 1024 (& 1024 (sys-info -1))) (println "This module needs newLISP compiled for extended ffi.") (println "Must exit") (exit)) ; the following assumes the libararies installed in the system library path (set 'LIB (if (= ostype "Win32") "libgsl-0.dll" ; 32-bit (= ostype "OSX") "libgsl.dylib" ; 32-bit or 64-bit (= ostype "Linux") "/usr/local/lib/libgsl.so" ; 32-bit or 64-bit )) ; load libgslcblas which contans functions referenced by libgsl ; the symbol cblas_sdsdot is not needed but newLISP versions before ; 10.4.2 can not use the 'import' statement without a function name ; libgslcblas is mainly needed internally by libgsl. ; On windows the library is automatically loaded by libgsl-0.dll. (if (= ostype "OSX") (import "libgslcblas.dylib" "cblas_sdsdot") (= ostype "Linux") (import "/usr/local/lib/libgslcblas.so" "cblas_sdsdot") ) ; structs are defined but only needed for debugging, instead use "void*" (struct 'complex "double" "double") ; complex numbers (struct 'block "long" "void*") ; vectors and matrices allocation block, size data-ptr (struct 'vector "long" "long" "void*" "block" "int") ; size stride data block owner (struct 'matrix "long" "long" "long" "void*" "block" "int") ; size1 size2 tda data block owner (import LIB "gsl_set_error_handler_off" "void*") ; turn of error handler (import LIB "gsl_strerror" "char*" "int") ; get error text ; pointers are only used passed around internally ; for debugging (unpack vector vptr) can always be used (import LIB "gsl_vector_alloc" "void*" "long") (import LIB "gsl_vector_calloc" "void*" "long") (import LIB "gsl_vector_free" "void" "void*") (import LIB "gsl_vector_get" "double" "void*" "long") (import LIB "gsl_vector_set" "void" "void*" "long" "double") ; in matrix functions "void*" can be used instead of struct "matrix" ; because pointers are only used passed around internally ; for debugging (unpack matrix Xptr) can always be used (import LIB "gsl_matrix_alloc" "void*" "long" "long") (import LIB "gsl_matrix_calloc" "void*" "long" "long") (import LIB "gsl_matrix_free" "void" "void*") (import LIB "gsl_matrix_get" "double" "void*" "long" "long") (import LIB "gsl_matrix_set" "void" "void*" "long" "long" "double") (import LIB "gsl_matrix_scale" "int" "void*" "double") ; linear algebra functions (import LIB "gsl_linalg_cholesky_decomp" "int" "void*") (import LIB "gsl_linalg_cholesky_solve" "int" "void*" "void*" "void*") (import LIB "gsl_linalg_QR_decomp" "int" "void*" "void*") (import LIB "gsl_linalg_QR_unpack" "int" "void*" "void*" "void*" "void*") (import LIB "gsl_linalg_QR_solve" "int" "void*" "void*" "void*" "void*") (import LIB "gsl_linalg_QR_lssolve" "int" "void*" "void*" "void*" "void*" "void*") (import LIB "gsl_linalg_SV_decomp" "int" "void*" "void*" "void*" "void*") (import LIB "gsl_linalg_SV_solve" "int" "void*" "void*" "void*" "void*" "void*") ; turn of error handler ; instead check return values from gsl_xxx_xxx functions ; and do 'throw-error' retrieving error text with 'gsl_strerror' (gsl_set_error_handler_off) ; helper functions for translating C-arrays and vectors to lists and back (define (get-vector-from-list x) (when (array? x) (set 'x (array-list x))) (letn (m (length x) xptr (gsl_vector_calloc m)) (dotimes (i m) (gsl_vector_set xptr i (pop x))) (list xptr m)) ) (define (get-matrix-from-list X) (when (array? X) (set 'X (array-list X))) (let (m (length X) n (length (X 0)) Xptr 0) (set 'Xptr (gsl_matrix_calloc m n)) (dotimes (i m) (set 'row (pop X)) ; row (dotimes (j n) (gsl_matrix_set Xptr i j (pop row)))) ; col (list Xptr m n)) ) (define (get-list-from-vector xr) (let (xptr (xr 0) n (xr 1) result nil) (dotimes (i n) (push (gsl_vector_get xptr i) result -1)) result ) ) (define (get-list-from-matrix Xr) (let (row nil col nil result nil Xptr (Xr 0) m (Xr 1) n (Xr 2)) (dotimes (i m) (set 'row nil) (dotimes (j n) (push (gsl_matrix_get Xptr i j) row -1)) (push row result -1)) result ) ) ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; API functions ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (define (gsl:CholeskyD A) (letn ( Ar (get-matrix-from-list A) Aptr (Ar 0) m (Ar 1) n (Ar 2) result 0 ) (set 'result (gsl_linalg_cholesky_decomp Aptr)) (unless (zero? result) (throw-error (gsl_strerror result))) ; zero the triangle above diagonal (for (i 0 (- m 2)) (for (j (+ i 1) (- n 1)) (gsl_matrix_set Aptr i j 0.0))) (set 'result (get-list-from-matrix Ar)) (gsl_matrix_free Aptr) result) ; matix L of A = LLt ) (define (gsl:Cholesky-solve A b) (letn ( Ar (get-matrix-from-list A) Aptr (Ar 0) m (Ar 1) n (Ar 2) br (get-vector-from-list b) bptr (br 0) xptr (gsl_vector_calloc n) result 0 ) (set 'result (gsl_linalg_cholesky_decomp Aptr)) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'result (gsl_linalg_cholesky_solve Aptr bptr xptr)) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'result (get-list-from-vector (list xptr n))) (gsl_matrix_free Aptr) (gsl_vector_free bptr) (gsl_vector_free xptr) result) ) (define (gsl:QRD A flag) (letn ( Ar (get-matrix-from-list A) Aptr (Ar 0) m (Ar 1) n (Ar 2) nt (min m n) Tptr (gsl_vector_calloc nt) Qptr (gsl_matrix_calloc m m) Rptr (gsl_matrix_calloc m n) result 0 ) (set 'result (gsl_linalg_QR_decomp Aptr Tptr)) (unless (zero? result) (throw-error (gsl_strerror result))) ; Aptr contains QR (set 'result (gsl_linalg_QR_unpack Aptr Tptr Qptr Rptr)) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'gsl:QR-tau (get-list-from-vector (list Tptr nt))) ; Q and R have reverse signs compared with published ; examples but values are the same and verify A = QR. ; reverse signs for Q and R unless flag set true (unless flag (gsl_matrix_scale Qptr -1) (gsl_matrix_scale Rptr -1) ; avoid -0 in lower triangular (for (i 1 (- m 1)) (for (j 0 (- i 1)) (gsl_matrix_set Rptr i j 0)))) (set 'gsl:QR-Q (get-list-from-matrix (list Qptr m m))) (set 'gsl:QR-R (get-list-from-matrix (list Rptr m n))) (gsl_matrix_free Aptr) (gsl_vector_free Tptr) (gsl_matrix_free Qptr) (gsl_matrix_free Rptr) gsl:QR-tau); vector tau ) (define (gsl:QR-solve A b) (letn ( Ar (get-matrix-from-list A) Aptr (Ar 0) m (Ar 1) n (Ar 2) nt (min m n) tptr (gsl_vector_calloc nt) br (get-vector-from-list b) bptr (br 0) xptr (gsl_vector_calloc n) resptr (gsl_vector_calloc m) result 0 ) (set 'result (gsl_linalg_QR_decomp Aptr tptr)) (unless (zero? result) (throw-error (gsl_strerror result))) ; Aptr contains QR (if (> m n) ; more rows than cols (set 'result (gsl_linalg_QR_lssolve Aptr tptr bptr xptr resptr)) (set 'result (gsl_linalg_QR_solve Aptr tptr bptr xptr)) ) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'result (get-list-from-vector (list xptr n))) (set 'gsl:QR-residual (get-list-from-vector (list resptr m))) (gsl_matrix_free Aptr) (gsl_vector_free tptr) (gsl_vector_free bptr) (gsl_vector_free xptr) (gsl_vector_free resptr) result) ) (define (gsl:SVD A) (letn ( Ar (get-matrix-from-list A) Aptr (Ar 0) m (Ar 1) n (Ar 2) Sptr (gsl_vector_calloc n) Vptr (gsl_matrix_calloc n n) work (gsl_vector_calloc n) result 0) (set 'result (gsl_linalg_SV_decomp Aptr Vptr Sptr work)) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'gsl:SVD-U (get-list-from-matrix Ar)) (set 'gsl:SVD-S (get-list-from-vector (list Sptr n))) (set 'gsl:SVD-V (get-list-from-matrix (list Vptr n n))) (gsl_matrix_free Aptr) (gsl_vector_free Sptr) (gsl_matrix_free Vptr) (gsl_vector_free work) gsl:SVD-S) ; vector of diagonals of S of A = USVt ) (define (gsl:SV-solve A b) (letn ( Ar (get-matrix-from-list A) Aptr (Ar 0) m (Ar 1) n (Ar 2) Sptr (gsl_vector_calloc n) Vptr (gsl_matrix_calloc n n) work (gsl_vector_calloc n) xptr (gsl_vector_calloc n) br (get-vector-from-list b) bptr (br 0) result 0) (set 'result (gsl_linalg_SV_decomp Aptr Vptr Sptr work)) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'result (gsl_linalg_SV_solve Aptr Vptr Sptr bptr xptr)) (unless (zero? result) (throw-error (gsl_strerror result))) (set 'result (get-list-from-vector (list xptr n))) (gsl_matrix_free Aptr) (gsl_vector_free Sptr) (gsl_matrix_free Vptr) (gsl_vector_free work) (gsl_vector_free bptr) (gsl_vector_free xptr) result) ) ; User helper functions (define (gsl:diagonal x flag) (when (array? x) (set 'x (array-list x))) (letn (len (length x) A (array len len (dup 0 len))) (dotimes (i len) (setf (A i i) (pop x))) (unless flag A (array-list A))) ) (context MAIN) ;;;;;;;;;;;;;;;;;;;; tests ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; (define (test-gsl) ; Cholesky decomp (println) (println "========== Cholesky example A = L * Lt") (println "(gsl:CholeskyD '((4 2 -2) (2 10 2) (-2 2 5)))" ) (set 'L (gsl:CholeskyD '((4 2 -1) (2 10 2) (-2 2 5)))) (println) (println "L -> ") (map println L) (println) (println "verify LLt ->") (map println (multiply L (transpose L))) (println) (println "(gsl:Cholesky-solve '((4 2 -2) (2 10 2) (-2 2 5)) '(1 2 3))") (set 'x (gsl:Cholesky-solve '((4 2 -2) (2 10 2) (-2 2 5)) '(1 2 3))) (println) (println "x -> " x) (println "Ax -> " (multiply '((4 2 -2) (2 10 2) (-2 2 5)) (transpose (list x)))) ; QR decomp example (println) (println "========== QR example A = Q * R") (println "(gsl:QRD '((12 -51 4) (6 167 -68) (-4 24 -41)) ) tau vector ===> " (gsl:QRD '((12 -51 4) (6 167 -68) (-4 24 -41)) )) (println) (println "Q -> ") (map println gsl:QR-Q) (println) (println "R -> ") (map println gsl:QR-R) (println) (println "verify QR -> ") (map println (multiply gsl:QR-Q gsl:QR-R)) (println) (println "(gsl:QR-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))") (set 'x (gsl:QR-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))) (println) (println "x -> " x) (println "Ax -> " (multiply '((12 -51 4) (6 167 -68) (-4 24 -41)) (transpose (list x)))) (println "residual -> " gsl:QR-residual) (println) (println "(gsl:QR-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))") (set 'x (gsl:QR-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))) (println) (println "x -> " x) (println "Ax -> " (multiply '((1 2) (3 4) (5 6) (7 8)) (transpose (list x)))) (println "residual -> " gsl:QR-residual) ; SV decomp example (println) (println "========== SVD example A = U * S * Vt") (println "(gsl:SVD '((1 2) (3 4) (5 6) (7 8))) => " (gsl:SVD '((1 2) (3 4) (5 6) (7 8)))) (println) (println "U -> ") (map println gsl:SVD-U) (println) (println "diagonal of S -> ") (println gsl:SVD-S) (println) (println "V -> ") (map println gsl:SVD-V) (println) (println "verify USVt -> ") (map println (multiply (multiply gsl:SVD-U (gsl:diagonal gsl:SVD-S)) (transpose gsl:SVD-V))) (println) (println "(gsl:SV-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))") (set 'x (gsl:SV-solve '((1 2) (3 4) (5 6) (7 8)) '(1 2 3 4))) (println) (println "x -> " x) (println "Ax -> " (multiply '((1 2) (3 4) (5 6) (7 8)) (transpose (list x)))) (println) (println "(gsl:SV-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))") (set 'x (gsl:SV-solve '((12 -51 4) (6 167 -68) (-4 24 -41)) '(1 2 3))) (println) (println "x -> " x) (println "Ax -> " (multiply '((12 -51 4) (6 167 -68) (-4 24 -41)) (transpose (list x)))) true ) ;(test-gsl) (exit) ; eof ;