basis-splines.lisp 5.72 KB
Newer Older
1 2
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
3
;; Time-stamp: <2016-06-15 23:07:19EDT basis-splines.lisp>
4
;;
5
;; Copyright 2008, 2009, 2011, 2012, 2016 Liam M. Healy
6 7 8 9 10 11 12 13 14 15 16 17 18 19
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <http://www.gnu.org/licenses/>.
20 21 22

(in-package :gsl)

23 24
;;; /usr/include/gsl/gsl_bspline.h

25
;;; Should be subclass of interpolation?
26

27
(defmobject basis-spline "gsl_bspline"
28
  ((order :sizet) (number-of-breakpoints :sizet))
29
  "basis spline"
30
  :gsl-version (1 9)
31
  :documentation			; FDL
32 33 34 35
  "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).")
36 37 38 39

;;; Constructing the knots vector
(defmfun knots (breakpoints workspace)
  "gsl_bspline_knots"
40
  (((mpointer breakpoints) :pointer) ((mpointer workspace) :pointer))
41
  :gsl-version (1 9)
42 43 44 45 46 47
  :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"
48 49
  ((a :double) (b :double) ((mpointer workspace) :pointer))
  :return (workspace)
50
  :gsl-version (1 9)
51 52 53 54 55 56
  :documentation			; FDL
  "Compute knots uniformly on the interval [a, b] and store
   them in the workspace.")

;;; Evaluation of B-splines

57 58
(defmfun evaluate
    ((workspace basis-spline) x
59
     &key (B (grid:make-foreign-array 'double-float
60
			  :dimensions (number-of-coefficients workspace))))
61
  "gsl_bspline_eval"
62
  ((x :double) ((mpointer B) :pointer) ((mpointer workspace) :pointer))
63
  :definition :method
Liam Healy's avatar
Liam Healy committed
64
  :outputs (B)
65
  :gsl-version (1 9)
66 67 68 69 70 71 72 73 74
  :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.")

75 76 77 78
;;; Query settings of the B-spline
;;; These are not documented but are in
;;; /usr/include/gsl/gsl_bspline.h; are they officially supported?

79
(defmfun number-of-coefficients (bspline)
80
  "gsl_bspline_ncoeffs"
81
  (((mpointer bspline) :pointer))
82
  :c-return :sizet
83
  :gsl-version (1 9))
84

85
(defmfun order ((bspline basis-spline))
86
  "gsl_bspline_order"
87
  (((mpointer bspline) :pointer))
88
  :definition :method
89
  :c-return :sizet
90
  :gsl-version (1 9))
91

92
(defmfun number-of-breakpoints (bspline)
93
  "gsl_bspline_nbreak"
94
  (((mpointer bspline) :pointer))
95
  :c-return :sizet
96 97
  :gsl-version (1 9)
  :documentation "The number of breakpoints of the basis spline bspline.")
98

99
(defmfun breakpoint (i bspline)
100
  "gsl_bspline_breakpoint"
101
  ((i :sizet) ((mpointer bspline) :pointer))
102
  :c-return :double
103 104
  :gsl-version (1 9)
  :documentation "The ith breakpoint of the basis spline bspline.")
105

106 107 108 109 110 111 112
(defmfun greville-abscissa (i bspline)
  "gsl_bspline_greville_abscissa"
  ((i :sizet) ((mpointer bspline) :pointer))
  :c-return :double
  :gsl-version (1 13)
  :documentation "Returns the location of the @math{i}-th Greville abscissa for the given B-spline basis.  For the ill-defined case when @math{k=1}, the implementation chooses to return breakpoint interval midpoints.")

113 114 115
;;; Examples and unit test

(defun bspline-example (&optional (ncoeffs 8))
116
  (let* ((order 4)			; cubic
117 118
	 (ndata 200)
	 (nbreak (+ ncoeffs 2 (- order)))
119
	 (bw (make-basis-spline order nbreak))
120
	 (mw (make-fit-workspace ndata ncoeffs))
Liam Healy's avatar
Liam Healy committed
121
	 (rng (make-random-number-generator +mt19937+ 0))
122 123 124 125 126 127 128
	 (B (grid:make-foreign-array 'double-float :dimensions ncoeffs))
	 (c (grid:make-foreign-array 'double-float :dimensions ncoeffs))
	 (cov (grid:make-foreign-array 'double-float :dimensions (list ncoeffs ncoeffs)))
	 (w (grid:make-foreign-array 'double-float :dimensions ndata))
	 (x (grid:make-foreign-array 'double-float :dimensions ndata))
	 (y (grid:make-foreign-array 'double-float :dimensions ndata))
	 (Xmatrix (grid:make-foreign-array 'double-float :dimensions (list ndata ncoeffs)))
129 130 131
	 (sigma 0.1d0))
    ;; The data to be fitted.
    (dotimes (i ndata)
132
      (let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float))
133
	     (yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
134
		    (sample rng :gaussian :sigma sigma))))
135 136 137
	(setf (grid:aref x i) xi
	      (grid:aref y i) yi
	      (grid:aref w i) (/ (expt sigma 2)))))
138 139 140 141 142
    ;; Uniform breakpoints [0, 15]
    (uniform-knots 0.0d0 15.0d0 bw)
    ;; Fit matrix
    (dotimes (i ndata)
      ;; Compute B_j
143
      (evaluate bw (grid:aref x i) :b B)
144 145
      ;; Fill in row i of X
      (dotimes (j ncoeffs)
146
	(setf (grid:aref Xmatrix i j) (grid:aref B j))))
147
    ;; Do the fit
148
    (linear-mfit Xmatrix y c w nil cov mw)
149 150
    ;; Return the smoothed curve
    (loop for xi from 0.0d0 to 15.0d0 by 0.1d0
151 152
       with yval and yerr
       do
153
       (evaluate bw xi :b B)
154 155 156 157 158
       (multiple-value-setq (yval yerr)
	 (multi-linear-estimate B c cov))
       collect yval into yvals 
       collect yerr into yerrs
       finally (return (values yvals yerrs)))))
159 160

(save-test basis-spline (bspline-example))