gsllを入れる(quicklisp)

quicklispでgsllを入れる。対象はsbclclisp
うちの環境はVirtualBox上のUbuntu 11.10(32bit)。
結論から書くとインストールはできるし使えるがユニットテストはいくらか失敗する。

gsllとは

gsllはgnuの科学計算ライブラリgslをcommon lispから使えるようにバインディングしたもの、らしい。
gsll本家は http://common-lisp.net/project/gsll/ 。インストール終わった段階で一読しておくのがいいかも。

あとgsl自体のドキュメントは:
本家(英語) http://www.gnu.org/software/gsl/manual/
日本語訳 http://www.cbrc.jp/~tominaga/translations/gsl/gsl-1.15/index.html

quicklispのインストール

quicklispというのは最近話題?のcommon lispパッケージインストーラで、ql:quickloadをrequireみたいな感覚で書いておくと自動でインストールとrequireを済ませてくれるというもの。現在ベータ公開中。
本家:http://www.quicklisp.org/

デフォルトだと~/quicklispにインストールされるのだけど、インストール先のディレクトリを実装別に分けてインストールすることにする:

cd
curl -O http://beta.quicklisp.org/quicklisp.lisp
sbcl
(load "quicklisp.lisp")
(quicklisp-quickstart:install)
^D
mv quicklisp sbcl-quicklisp
sbcl
(load "sbcl-quicklisp/setuyp.lisp")
(ql:add-to-init-file)
^D
clisp
(load "quicklisp.lisp")
(quicklisp-quickstart:install)
^D
mv quicklisp clisp-quicklisp
sbcl
(load "sbcl-quicklisp/setup.lisp")
(ql:add-to-init-file)
^D

ただし^DはCtrl+D。

gsllのインストール

先にgsl自体のインストールをしておかないといけない(はず。apt-get install gslとかする)。

sbclはなぜかbabelへの依存を検出してくれないので二度に分けてインストールすれば良い。

sbcl
(ql:quickload :babel)
(ql:quickload :gsll)

clispはファイルの手直しが必要。

vim clisp-quicklisp/dists/quicklisp/software/gsll-20111105-git/init/init.lisp

とかして42行目のread-lineをコメントアウト*1

その後

clisp
(ql:quickload :gsll)

として読み込む。エラーが何度か発生するがcontinueと打ってreturnを繰り返せば完了する。

ユニットテスト

うまくインストールできているかのテスト。

(ql:quickload :lisp-unit)
(ql:quickload :gsll)
(in-package :gsll)
(lisp-unit:run-tests)

うちの場合の結果(だーっとログが出るが最後の行のみ):
sbcl

TOTAL: 2501 assertions passed, 100 failed, 69 execution errors.

clisp

TOTAL: 3247 assertions passed, 192 failed, 5 execution errors.

clispの方がエラー少ないのか。しかしexecution errorはいやだな。ログをざっと眺めてみるとfloating point overflowやdivision by zeroだったりする。sbclの方は"There is no applicable method for the generic function"と言っておりどうもforeign-arrayへの変換がうまく行っていない?*2

quickload時にPermission deniedと怒られるとき

もしql:quickloadで"Permiossion denied"とか言われる場合はこれまでにroot権限でなにかインストールしてしまっている可能性があるのでキャッシュを削除してquicklispを入れなおす。

cd
rm -rf .cache/common-lisp
rm -rf quicklisp #ディレクトリ名を変えた場合は適宜読み替える

例を見る

gsll:examplesを引数なしで呼ぶと例のトピック一覧が見られる。そのなかから選んで引数として渡すとそれについての例がリストとして返ってくる。
lispのreplを立ち上げたら

CL-USER> (ql:quickload 'gsll)
To load "gsll":
  Load 1 ASDF system:
    gsll
; Loading "gsll"
.............
(GSLL)
CL-USER> (in-package :gsll)
#<PACKAGE "GSLL">
GSL> (examples)
(BASIS-SPLINE NONLINEAR-LEAST-SQUARES LINEAR-LEAST-SQUARES MINIMIZATION-MULTI
 ROOTS-MULTI MINIMIZATION-ONE ROOTS-ONE HANKEL SERIES-ACCELERATION CHEBYSHEV
 INTERPOLATION ODE NUMERICAL-DIFFERENTIATION MONTE-CARLO NUMERICAL-INTEGRATION
 NTUPLE HISTOGRAM MEDIAN-PERCENTILE CORRELATION COVARIANCE AUTOCORRELATION
 HIGHER-MOMENTS ABSOLUTE-DEVIATION MATRIX-STANDARD-DEVIATION-WITH-FIXED-MEAN
 VECTOR-STANDARD-DEVIATION-WITH-FIXED-MEAN MATRIX-VARIANCE-WITH-FIXED-MEAN
 VECTOR-VARIANCE-WITH-FIXED-MEAN MATRIX-STANDARD-DEVIATION-WITH-MEAN
 VECTOR-STANDARD-DEVIATION-WITH-MEAN MATRIX-STANDARD-DEVIATION
 VECTOR-STANDARD-DEVIATION MATRIX-VARIANCE-WITH-MEAN VECTOR-VARIANCE-WITH-MEAN
 MATRIX-VARIANCE VECTOR-VARIANCE MATRIX-MEAN VECTOR-MEAN SHUFFLING-SAMPLING
 LOGARITHMIC HYPERGEOMETRIC-RANDIST GEOMETRIC NEGATIVE-BINOMIAL MULTINOMIAL
 BINOMIAL BERNOULLI POISSON DISCRETE DIRICHLET GUMBEL2 GUMBEL1 WEIBULL
 SPHERICAL-VECTOR PARETO LOGISTIC BETA TDIST FDIST CHI-SQUARED LOGNORMAL FLAT
 GAMMA-RANDIST LEVY LANDAU RAYLEIGH-TAIL RAYLEIGH CAUCHY EXPONENTIAL-POWER
 LAPLACE EXPONENTIAL GAUSSIAN-BIVARIATE GAUSSIAN-TAIL GAUSSIAN
 QUASI-RANDOM-NUMBER-GENERATORS RANDOM-NUMBER-GENERATORS FFT-PULSE EIGENSYSTEMS
 HOUSEHOLDER CHOLESKY SVD QRPT QR LU MATRIX-PRODUCT-NONSQUARE RANK-1-UPDATE
 MATRIX-PRODUCT-HERMITIAN MATRIX-PRODUCT-SYMMETRIC INVERSE-MATRIX-PRODUCT
 MATRIX-PRODUCT-TRIANGULAR MATRIX-PRODUCT GIVENS SCALE AXPY BLAS-COPY BLAS-SWAP
 INDEX-MAX ABSOLUTE-SUM EUCLIDEAN-NORM CDOT DOT SORT-VECTOR-LARGEST-INDEX
 SORT-MATRIX-LARGEST SORT-VECTOR-LARGEST SORT-VECTOR-SMALLEST-INDEX
 SORT-MATRIX-SMALLEST SORT-VECTOR-SMALLEST SORT-VECTOR-INDEX SORT-MATRIX
 SORT-VECTOR ZETA TRIGONOMETRY TRANSPORT SYNCHROTRON PSI POWER MATHIEU
 LOGARITHM LEGENDRE LAMBERT LAGUERRE HYPERGEOMETRIC GEGENBAUER GAMMA
 FERMI-DIRAC EXPONENTIAL-INTEGRALS EXPONENTIAL-FUNCTIONS ERROR-FUNCTIONS
 ELLIPTIC-FUNCTIONS ELLIPTIC-INTEGRALS ELEMENTARY DILOGARITHM DEBYE DAWSON
 COUPLING COULOMB CLAUSEN BESSEL AIRY POLYNOMIAL COMBINATION PERMUTATION
 MATRIX-TRANSPOSE MATRIX-TRANSPOSE* SWAP-ROW-COLUMN SWAP-COLUMNS SWAP-ROWS
 SETF-COLUMN COLUMN SETF-ROW ROW SET-IDENTITY VECTOR-REVERSE SWAP-ELEMENTS
 SET-BASIS MATRIX-MINMAX-INDEX VECTOR-MINMAX-INDEX MATRIX-MAX-INDEX
 VECTOR-MAX-INDEX MATRIX-MIN-INDEX VECTOR-MIN-INDEX MATRIX-MINMAX VECTOR-MINMAX
 MATRIX-MIN VECTOR-MIN MATRIX-MAX VECTOR-MAX MATRIX-ADD-SCALAR
 VECTOR-ADD-SCALAR MATRIX-MULT-SCALAR VECTOR-MULT-SCALAR MATRIX-DIV VECTOR-DIV
 MATRIX-MULT VECTOR-MULT MATRIX-SUB VECTOR-SUB MATRIX-ADD VECTOR-ADD
 MATRIX-SWAP VECTOR-SWAP MATRIX-SET-ZERO VECTOR-SET-ZERO MATRIX-SET-ALL
 VECTOR-SET-ALL MATHEMATICAL)
GSL> (examples 'cdot) ;複素ベクトルの内積
((LET ((V1
        (GRID:MAKE-FOREIGN-ARRAY '(COMPLEX SINGLE-FLOAT) :INITIAL-CONTENTS
                                 '(-34.5f0 8.24f0 3.29f0 -8.93f0 34.12f0
                                   -6.15f0 49.27f0 -13.49f0 32.5f0 42.73f0
                                   -17.24f0 43.31f0 -16.12f0 -8.25f0 21.44f0
                                   -49.08f0)))
       (V2
        (GRID:MAKE-FOREIGN-ARRAY '(COMPLEX SINGLE-FLOAT) :INITIAL-CONTENTS
                                 '(32.5f0 42.73f0 -17.24f0 43.31f0 -16.12f0
                                   -8.25f0 21.44f0 -49.08f0 -39.66f0 -49.46f0
                                   19.68f0 -5.55f0 -8.82f0 25.37f0 -30.58f0
                                   31.67f0))))
   (CDOT V1 V2))
 (LET ((V1
        (GRID:MAKE-FOREIGN-ARRAY '(COMPLEX DOUBLE-FLOAT) :INITIAL-CONTENTS
                                 '(-34.5 8.24 3.29 -8.93 34.12 -6.15 49.27
                                   -13.49 32.5 42.73 -17.24 43.31 -16.12 -8.25
                                   21.44 -49.08)))
       (V2
        (GRID:MAKE-FOREIGN-ARRAY '(COMPLEX DOUBLE-FLOAT) :INITIAL-CONTENTS
                                 '(32.5 42.73 -17.24 43.31 -16.12 -8.25 21.44
                                   -49.08 -39.66 -49.46 19.68 -5.55 -8.82 25.37
                                   -30.58 31.67))))
   (CDOT V1 V2)))
GSL> (LET ((V1
        (GRID:MAKE-FOREIGN-ARRAY '(COMPLEX SINGLE-FLOAT) :INITIAL-CONTENTS
                                 '(-34.5f0 8.24f0 3.29f0 -8.93f0 34.12f0
                                   -6.15f0 49.27f0 -13.49f0 32.5f0 42.73f0
                                   -17.24f0 43.31f0 -16.12f0 -8.25f0 21.44f0
                                   -49.08f0)))
       (V2
        (GRID:MAKE-FOREIGN-ARRAY '(COMPLEX SINGLE-FLOAT) :INITIAL-CONTENTS
                                 '(32.5f0 42.73f0 -17.24f0 43.31f0 -16.12f0
                                   -8.25f0 21.44f0 -49.08f0 -39.66f0 -49.46f0
                                   19.68f0 -5.55f0 -8.82f0 25.37f0 -30.58f0
                                   31.67f0))))
   (CDOT V1 V2))
#C(-6252.624f0 0.0f0)

なんかよくわからないけどエラーになる例なんかも混じってたりする、(examples 'MATRIX-PRODUCT-NONSQUARE)の後ろ2つとか(行列の掛け算なのだが左の行列の列数と右の行列の行数が違う)。
(mapcar #'(lambda (x) (list x '=> (cl:ignore-errors (eval x)))) (examples 'MATRIX-PRODUCT-NONSQUARE))とかすると見やすいかも

ドキュメントを見る

gslの関数はgsllでいうどの関数かというのはgsl-lookup。だいたいgslのドキュメントでやりたいことを調べてこれを使うという手順を踏むことになる。

GSL> (gsl-lookup "gsl_sf_elljac_e")
JACOBIAN-ELLIPTIC-FUNCTIONS
T

関数のドキュメント

GSL> (documentation #'cdot 'function)
"The complex conjugate scalar product x^H y for the vectors."

引数リストはslimeとかの支援環境じゃないとdescribe使ってみるしかないのかな?

GSL> (describe #'cdot)
#<STANDARD-GENERIC-FUNCTION CDOT (2)>
  [generic-function]

Lambda-list: (X Y)
Derived type: (FUNCTION (T T) *)
Documentation:
  The complex conjugate scalar product x^H y for the vectors.
Method-combination: STANDARD
Methods:
  (CDOT (GRID:VECTOR-COMPLEX-SINGLE-FLOAT
         GRID:VECTOR-COMPLEX-SINGLE-FLOAT))
  (CDOT (GRID:VECTOR-COMPLEX-DOUBLE-FLOAT
         GRID:VECTOR-COMPLEX-DOUBLE-FLOAT))
Source file: /build/buildd/sbcl-1.0.50.0/src/cold/warm.lisp
; No value

行列

例として、というわけでもないけど。

invert-matrixが定義されていると本家ページにはあるがなぜか関数がバインドされていないようだ。antik:invert-matrixが使える。

ところで行列を扱おうとするとantikのマニュアル*3のgridの部分を参照することになるが、(grid:make-foreign-matrixで作った行列でないと扱えないぽいので)とりあえずおまじないとして

(setf grid:*default-grid-type* 'grid:foreign-array)

としておけば#mリードマクロで書ける*4

GSL> (matrix-product #m((1 2) (3 4)) #m((2 3) (4 5)))
#m((10.000000000000000d0 13.000000000000000d0)
(22.000000000000000d0 29.000000000000000d0))

#mの中では^を改行のように使うことができて微妙に楽に書ける:

GSL> (matrix-product #m(1 2 ^ 3 4) #m(2 3 ^ 4 5))
#m((10.000000000000000d0 13.000000000000000d0)
(22.000000000000000d0 29.000000000000000d0))

#mはベクトルも表せるし、matrix-productは2つ目の引数がベクトルなら行列とベクトルの掛け算をやってくれる。

GSL> (matrix-product #m(1 1 ^ 0 1) #m( 1 1))
#m(2.000000000000000d0 1.000000000000000d0)

詳細は(documentation 'matrix-product 'function)と(describe 'matrix-product)を参照。

加減算はelt+、elt-、また要素ごとの乗除算はelt*、elt/だがこいつらは1つめの引数を破壊する(というか結果を1つ目に渡した行列に書きこむ)ので注意。行列じゃなくベクトル、ヒストグラムに対しても使える。

あと行列の掛け算はfloat・double(とその複素数)型に対してしか定義されてないので微妙に注意。そもそもgslでもそうだからね。

(matrix-product #m((1 1) (0 1)) #m((1 1) (0 1)))

はうまくいくが#mを#15m(32ビット符号付き整数)にすると未定義エラーを吐く。1,2,3,4(doubleとfloatの実数と複素数)は大丈夫*5