File: gsl.lsp

package info (click to toggle)
newlisp 10.7.1-1
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 8,460 kB
  • ctags: 4,357
  • sloc: ansic: 33,202; lisp: 7,369; java: 7,012; sh: 647; makefile: 273
file content (513 lines) | stat: -rw-r--r-- 18,405 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
;; @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.
;; @version 1.2 - changed ostype Win32 to Windows
;; @author Lutz Mueller 2012

;; <h2>Module GSL For The GNU Scientific Library</h2>
;; 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 <tt>(test-gsl)</tt> to test all implemented functions.

;; @syntax (gsl:CholeskyD <matrix-A>)
;; @param <matrix-A> A square matrix of m row vectors with n = m column elements each.
;; @return The matrix L .
;; The function performs a Cholesky decomposition of <matrix-A>. 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 <matrix-A> <vector-b>)
;; @param <matrix-A> A square matrix of m row vectors with n = m column elements each.
;; @params <vector-b> 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 <matrix-A>)
;; @param <matrix-A> A list of m row vector lists with n column elements each.
;; @return Vector <em>tau</em> 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 <matrix-A> <vector-b>)
;; @param <matrix-A> A matrix of m row vectors with n column elements each.
;; @params <vector-b> 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 <matrix-A>)
;; @param <matrix-A> 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 <matrix-A> 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 <matrix-A> <vector-b>)
;; @param <matrix-A> A matrix of m row vectors with n column elements each.
;; @params <vector-b> 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 "libgsl.so")

; 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.
(import "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 ;