basis-splines.lisp 3.67 KB
Newer Older
1 2
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
3
;; Time-stamp: <2008-11-30 23:43:58EST basis-splines.lisp>
lhealy's avatar
lhealy committed
4
;; $Id$
5 6 7 8 9 10 11 12 13 14

(in-package :gsl)

;;; Initializing the B-splines solver

(defgo-s (basis-spline order number-of-breakpoints)
    allocate-basis-spline free-basis-spline nil 2)

(defmfun allocate-basis-spline (order number-of-breakpoints)
  "gsl_bspline_alloc"
Liam Healy's avatar
Liam Healy committed
15
  ((order sizet) (number-of-breakpoints sizet))
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
  :c-return :pointer
  :export nil
  :index (letm basis-spline)
  :documentation			; FDL
  "Allocate a workspace for computing B-splines. The number of
   breakpoints is given by number-of-breakpoints.  This leads to n =
   nbreak + k - 2 basis functions where k = order. Cubic B-splines are
   specified by k = 4. The size of the workspace is O(5k + nbreak)."  )

(defmfun free-basis-spline (bspline)
  "gsl_bspline_free"
  ((bspline :pointer))
  :c-return :void
  :export nil
  :index (letm basis-spline)
  :documentation			; FDL
  "Free the memory associated with the workspace of the bspline.")

;;; Constructing the knots vector

(defmfun knots (breakpoints workspace)
  "gsl_bspline_knots"
  ((breakpoints :pointer) (workspace :pointer))
  :documentation			; FDL
  "Compute the knots associated with the given breakpoints and store
   them in the workspace.")

(defmfun uniform-knots (a b workspace)
  "gsl_bspline_knots_uniform"
  ((a :double) (b :double) (workspace :pointer))
  :documentation			; FDL
  "Compute knots uniformly on the interval [a, b] and store
   them in the workspace.")

;;; Evaluation of B-splines

(defmfun evaluate-bspline (x B workspace)
  "gsl_bspline_eval"
Liam Healy's avatar
Liam Healy committed
54 55
  ((x :double) ((mpointer B) :pointer) (workspace :pointer))
  :outputs (B)
56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77
  :documentation			; FDL
  "Evaluate all B-spline basis functions at the position x and store
   them in the GSL vector B, so that the ith element of B is B_i(x).
   B must be of length n = nbreak + k - 2. This value is
   also stored in the workspace. It is far more
   efficient to compute all of the basis functions at once than to
   compute them individually, due to the nature of the defining
   recurrence relation.")

;;; Examples and unit test

#|

|#

(defun bspline-example (&optional (ncoeffs 8))
  (letm ((order 4)			; cubic
	 (ndata 200)
	 (nbreak (+ ncoeffs 2 (- order)))
	 (bw (basis-spline order nbreak))
	 (mw (fit-workspace ndata ncoeffs))
	 (rng (random-number-generator *mt19937* 0))
78 79 80 81 82 83 84
	 (B (make-array* 'double-float :dimensions ncoeffs))
	 (c (make-array* 'double-float :dimensions ncoeffs))
	 (cov (make-array* 'double-float :dimensions (list ncoeffs ncoeffs)))
	 (w (make-array* 'double-float :dimensions ndata))
	 (x (make-array* 'double-float :dimensions ndata))
	 (y (make-array* 'double-float :dimensions ndata))
	 (Xmatrix (make-array* 'double-float :dimensions (list ndata ncoeffs)))
85 86 87
	 (sigma 0.1d0))
    ;; The data to be fitted.
    (dotimes (i ndata)
88
      (let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float))
89 90
	     (yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
		    (gaussian rng sigma))))
91 92 93
	(setf (maref x i) xi
	      (maref y i) yi
	      (maref w i) (/ (expt sigma 2)))))
94 95 96 97 98
    ;; Uniform breakpoints [0, 15]
    (uniform-knots 0.0d0 15.0d0 bw)
    ;; Fit matrix
    (dotimes (i ndata)
      ;; Compute B_j
99
      (evaluate-bspline (maref x i) B bw)
100 101
      ;; Fill in row i of X
      (dotimes (j ncoeffs)
102
	(setf (maref Xmatrix i j) (maref B j))))
103 104 105 106 107 108 109 110 111 112 113 114
    ;; Do the fit
    (weighted-linear-mfit Xmatrix w y c cov mw)
    ;; Return the smoothed curve
    (loop for xi from 0.0d0 to 15.0d0 by 0.1d0
	  with yval and yerr
	  do
	  (evaluate-bspline xi B bw)
	  (multiple-value-setq (yval yerr)
	      (multi-linear-estimate B c cov))
	  collect yval into yvals 
	  collect yerr into yerrs
	  finally (return (values yvals yerrs)))))