スプライン補間

LISP

急遽必要になったのでCommon Lispでスプライン補間の実装がないか調べたが、どうもないので以下のPDFを見ながら適当に実装してみる。

効率とかは知らぬ。途中、三重対角行列の逆行列を計算しているがO(N)で求める方法があるらしい。今回は普通にオレオレユーティリティの中のGauss-Jordan法を使う。

(ql:quickload :wiz-util) ; https://github.com/masatoi/wiz-util
(ql:quickload :clgplot) ; https://github.com/masatoi/clgplot

(defun make-spline (x-list y-list)
  (let* ((h (wiz:diff-list x-list))
         (dy (wiz:diff-list y-list))
         (v (wiz:nlet iter ((h h) (dy dy) (product nil))
              (if (null (cdr h))
                (nreverse product)
                (iter (cdr h)
                      (cdr dy)
                      (cons (* 6 (- (/ (cadr dy) (cadr h))
                                    (/ (car dy) (car h))))
                            product)))))
         (mat-size (- (length x-list) 2))
         (vvec (make-array (list mat-size 1) :element-type 'double-float))
         (mat (make-array (list mat-size mat-size) :element-type 'double-float :initial-element 0d0)))
    (loop for j from 0 to (1- mat-size)
          for hj-1 in h
          for hj in (cdr h)
          for vj in v
          do
       (setf (aref vvec j 0) vj) ; init vvec
       (cond ((= j 0)            ; init mat
              (setf (aref mat 0 0) (* 2 (+ hj-1 hj))
                    (aref mat 0 1) hj))
             ((= j (1- mat-size))
              (setf (aref mat (1- mat-size) (- mat-size 2)) hj-1
                    (aref mat (1- mat-size) (1- mat-size)) (* 2 (+ hj-1 hj))))
             (t (setf (aref mat j (1- j)) hj-1
                      (aref mat j j) (* 2 (+ hj-1 hj))
                      (aref mat j (1+ j)) hj))))
    (let* ((mat-1 (wiz:m-1 mat))
           (uvec (wiz:m* mat-1 vvec))
           (u (append (cons 0d0 (loop for i from 0 to (1- mat-size) collect (aref uvec i 0)))
                      '(0d0)))
           (b (mapcar (lambda (uj) (/ uj 2)) u))
           (a (mapcar (lambda (uj+1 bj hj) (/ (- uj+1 (* 2 bj)) (* 6 hj))) (cdr u) b h))
           (c (mapcar (lambda (hj yj+1 aj bj dj)
                        (/ (- yj+1 (* aj hj hj hj) (* bj hj hj) dj) hj))
                      h (cdr y-list) a b y-list)))
      (lambda (x)
        (let* ((pos0 (position-if (lambda (xj+1) (<= x xj+1)) x-list)) ; x-list must be sorted
               (pos (cond ((null pos0) (- (length x-list) 2))
                          ((<= pos0 0) 0)
                          (t (1- pos0))))
               (xj (nth pos x-list))
               (yj (nth pos y-list))
               (aj (nth pos a))
               (bj (nth pos b))
               (cj (nth pos c)))
          (+ (* aj (expt (- x xj) 3))
             (* bj (expt (- x xj) 2))
             (* cj (- x xj))
             yj))))))

make-spline関数でデータ点を補間した値を返す関数が生成される。適当なデータ点を与えてみる。

(defparameter spline-func
  (make-spline '(0d0 1.5d0 2.2d0 3.7d0 3.9d0 4.5d0 5.1d0 6.6d0 7.7d0) ; x-list
               '(1d0 3d0 2d0 5d0 1d0 4d0 3.1d0 1d0 5d0)))             ; y-list

結果をプロットしてみると、

(let ((x-list (wiz:seq -1d0 8d0 :by 0.1d0)))
  (clgp:plot-lists (list (mapcar spline-func x-list)
                         '(1d0 3d0 2d0 5d0 1d0 4d0 3.1d0 1d0 5d0))
                   :x-lists (list x-list
                                  '(0d0 1.5d0 2.2d0 3.7d0 3.9d0 4.5d0 5.1d0 6.6d0 7.7d0))
                   :style '(lines points)))


できているっぽい。