Cubic spline coefficient computation

cubic-interp.lisp

;; First, some helper functions: 
(defun window (vec start &optional end)
  "As SUBSEQ, but return a displaced vector instead of copying."
  (sunless end (setf it (length vec)))
  (make-array (- end start)
              :displaced-to vec
              :displaced-index-offset start))

(defun windows (vec length)
  "Return a list of all overlapping subsequences length LENGTH of
VEC."
  (loop for i upto (- (length vec) length)
     collecting (window vec i (+ i length))))

(defun elt+ (sequence &rest offsets)
  (elt sequence (apply #'+ offsets)))

;; And now, the main one:

;; Stolen^H^H^H^H^H^H^H Inspiration taken from
;; https://en.wikipedia.org/wiki/Spline_(mathematics)
(defun cubic-spline-coeffs (xs ys)
  (let* ((k (- (length ys) 1))
         (len (1+ k))
         (a (make-array `(,len) :initial-contents ys))
         (b (make-array `(,k)))
         (d (make-array `(,k)))
         (mu (aprog1 (make-array `(,k))
               (setf (elt it 0) 0)))
         (h (map 'simple-vector
                 (lambda (xs) (- (elt xs 1) (elt xs 0)))
                 (windows xs 2)))
         (alpha (map 'simple-vector
                     (lambda (hs as) (- (* 3 (/ (elt hs 1)) (- (elt as 2) (elt as 1)))
                                        (* 3 (/ (elt hs 0)) (- (elt as 1) (elt as 0)))))
                     (windows h 2)
                     (windows a 3)))
         (c (aprog1 (make-array `(,len))
              (setf (elt it k) 0)))
         (l (aprog1 (make-array `(,len))
              (setf (elt it 0) 1)
              (setf (elt it k) 1)))
         (z (aprog1 (make-array `(,len))
              (setf (elt it 0) 0)
              (setf (elt it k) 0))))
    (loop for i from 1 to (- k 1)
       do
         (setf (elt l i) (- (* 2 (- (elt+ xs i 1) (elt+ xs i -1)))
                            (* (elt+ h i -1) (elt+ mu i -1))))
         (setf (elt mu i) (/ (elt h i) (elt l i)))
         (setf (elt z i) (/ (- (elt+ alpha i -1) (* (elt+ h i -1) (elt+ z i -1)))
                            (elt l i))))
    (loop for j downfrom (- k 1) to 0
       do
         (setf (elt c j) (- (elt z j) (* (elt mu j) (elt+ c j 1))))
         (setf (elt b j) (- (/ (- (elt+ a j 1) (elt a j))
                               (elt h j))
                            (/ (* (elt h j) (+ (elt+ c j 1) (* 2 (elt c j))))
                               3)))
         (setf (elt d j) (/ (- (elt+ c j 1) (elt c j)) (* 3 (elt h j)))))
    (map 'list #'list a b c d xs)))