Commit 04f30601 authored by Liam Healy's avatar Liam Healy

Generic functions, fix interpolation

Turned many functions, particularly those having to do with
interpolation, into generic functions and methods, and thus renamed
them and changed their arglists.  Added test for basis spline; it
works in SBCL but does not in CCL.
parent 078abffa
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2009-01-25 10:11:48EST basis-splines.lisp>
;; Time-stamp: <2009-01-28 20:55:57EST basis-splines.lisp>
;; $Id$
(in-package :gsl)
......@@ -12,6 +12,7 @@
(defmobject basis-spline "gsl_bspline"
((order sizet) (number-of-breakpoints sizet))
"basis spline"
:gsl-version (1 9)
:documentation ; FDL
"Allocate a workspace for computing B-splines. The number of
breakpoints is given by number-of-breakpoints. This leads to n =
......@@ -22,6 +23,7 @@
(defmfun knots (breakpoints workspace)
"gsl_bspline_knots"
(((mpointer breakpoints) :pointer) ((mpointer workspace) :pointer))
:gsl-version (1 9)
:documentation ; FDL
"Compute the knots associated with the given breakpoints and store
them in the workspace.")
......@@ -30,16 +32,22 @@
"gsl_bspline_knots_uniform"
((a :double) (b :double) ((mpointer workspace) :pointer))
:return (workspace)
:gsl-version (1 9)
: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)
(defmfun evaluate
((workspace basis-spline) x
&key (B (make-marray 'double-float
:dimensions (number-of-coefficients workspace))))
"gsl_bspline_eval"
((x :double) ((mpointer B) :pointer) ((mpointer workspace) :pointer))
:definition :method
:outputs (B)
:gsl-version (1 9)
: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).
......@@ -53,27 +61,31 @@
;;; These are not documented but are in
;;; /usr/include/gsl/gsl_bspline.h; are they officially supported?
(defmfun bspline-ncoeffs (B)
(defmfun number-of-coefficients (bspline)
"gsl_bspline_ncoeffs"
(((mpointer B) :pointer))
:c-return sizet)
(((mpointer bspline) :pointer))
:c-return sizet
:gsl-version (1 9))
(defmfun bspline-order (B)
(defmfun order (bspline)
"gsl_bspline_order"
(((mpointer B) :pointer))
:c-return sizet)
(((mpointer bspline) :pointer))
:c-return sizet
:gsl-version (1 9))
(defmfun number-of-breakpoints (B)
(defmfun number-of-breakpoints (bspline)
"gsl_bspline_nbreak"
(((mpointer B) :pointer))
(((mpointer bspline) :pointer))
:c-return sizet
:documentation "The number of breakpoints of the basis spline B.")
:gsl-version (1 9)
:documentation "The number of breakpoints of the basis spline bspline.")
(defmfun breakpoint (i B)
(defmfun breakpoint (i bspline)
"gsl_bspline_breakpoint"
((i sizet) ((mpointer B) :pointer))
((i sizet) ((mpointer bspline) :pointer))
:c-return :double
:documentation "The ith breakpoint of the basis spline B.")
:gsl-version (1 9)
:documentation "The ith breakpoint of the basis spline bspline.")
;;; Examples and unit test
......@@ -105,7 +117,7 @@
;; Fit matrix
(dotimes (i ndata)
;; Compute B_j
(evaluate-bspline (maref x i) B bw)
(evaluate bw (maref x i) :b B)
;; Fill in row i of X
(dotimes (j ncoeffs)
(setf (maref Xmatrix i j) (maref B j))))
......@@ -115,9 +127,11 @@
(loop for xi from 0.0d0 to 15.0d0 by 0.1d0
with yval and yerr
do
(evaluate-bspline xi B bw)
(evaluate bw xi :b B)
(multiple-value-setq (yval yerr)
(multi-linear-estimate B c cov))
collect yval into yvals
collect yerr into yerrs
finally (return (values yvals yerrs)))))
(save-test basis-spline (bspline-example))
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
;; Time-stamp: <2009-01-25 10:12:45EST chebyshev.lisp>
;; Time-stamp: <2009-01-28 19:05:30EST chebyshev.lisp>
;; $Id$
(in-package :gsl)
......@@ -26,12 +26,13 @@
;;; The functions that don't return are defined, but it is recommended
;;; to use the functions that do return error (and ignore it if
;;; desired) in the form of #'evaluate-chebyshev.
;;; desired) in the form of #'evaluate.
(defmfun evaluate-chebyshev (chebyshev x &optional order)
(defmfun evaluate ((object chebyshev) x &key order)
("gsl_cheb_eval" "gsl_cheb_eval_n")
((((mpointer chebyshev) :pointer) (x :double))
(((mpointer chebyshev) :pointer) (order sizet) (x :double)))
((((mpointer object) :pointer) (x :double))
(((mpointer object) :pointer) (order sizet) (x :double)))
:definition :method
:c-return :double
:documentation ; FDL
"Evaluate the Chebyshev series at a point x. If order is supplied,
......@@ -89,8 +90,8 @@
(format t "~&~a ~a ~a ~a"
x
(chebyshev-step x)
(evaluate-chebyshev cheb x 10)
(evaluate-chebyshev cheb x)))))))
(evaluate cheb x :order 10)
(evaluate cheb x)))))))
(defun chebyshev-point-example (x)
(check-type x double-float)
......@@ -101,9 +102,9 @@
(derivative-chebyshev deriv cheb)
(integral-chebyshev integ cheb)
(list
(evaluate-chebyshev cheb x)
(evaluate-chebyshev deriv x)
(evaluate-chebyshev integ x))))
(evaluate cheb x)
(evaluate deriv x)
(evaluate integ x))))
;;; Unit test
(save-test chebyshev
......
......@@ -35,6 +35,7 @@ series-acceleration
linear-least-squares
nonlinear-least-squares
minimization-multi
basis-spline
=== GSL questions ===
Sun Nov 9 2008:
......@@ -48,3 +49,6 @@ Sun Jan 11 2009:
gsl_permutation_valid (presumably combination too) signals an error if
the permutation is invalid, rather than just returning an error code.
Is this intentional?
Wed Jan 28 2009:
What are the _e functions in interpolation for, how are the results to
be interpreted?
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-26 22:45:42EST gsll-tests.asd>
;; Time-stamp: <2009-01-28 20:59:16EST gsll-tests.asd>
;; $Id$
(asdf:defsystem "gsll-tests"
......@@ -24,6 +24,7 @@
(:file "airy")
(:file "autocorrelation")
(:file "axpy")
(:file "basis-spline")
(:file "bernoulli")
(:file "bessel")
(:file "beta")
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-24 16:56:13EST gsll.asd>
;; Time-stamp: <2009-01-28 21:05:10EST gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -196,4 +196,4 @@
(:file "minimization-multi" :depends-on (generic))
(:file "linear-least-squares")
(:file "nonlinear-least-squares" :depends-on (generic))))
(:file "basis-splines" :depends-on (init data))))
(:file "basis-splines" :depends-on (init data random))))
;; Load GSL
;; Liam Healy Sat Mar 4 2006 - 18:53
;; Time-stamp: <2009-01-26 21:39:39EST init.lisp>
;; Time-stamp: <2009-01-27 21:32:59EST init.lisp>
;; $Id$
(defpackage gsll
......@@ -44,8 +44,9 @@
(pushnew :callback-toplevel-only *features*))
;;; CFFI currently doesn't allow uninterned callback names for CCL;
;;; patch here.
;;; patch here. This patch can be removed once the fix in CFFI:
;;; http://common-lisp.net/pipermail/cffi-devel/2009-January/003046.html
;;; makes it into a released version.
#-openmcl (import 'cl-utilities:with-unique-names :gsll)
#+openmcl
......
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2009-01-25 17:27:46EST mobject.lisp>
;; $Id$
;; Time-stamp: <2009-01-28 19:19:54EST mobject.lisp>
;;; GSL objects are represented in GSLL as and instance of a 'mobject.
;;; The macro demobject takes care of defining the appropriate
......@@ -149,6 +148,10 @@
(defgeneric name (object)
(:documentation "The name given to the GSL object."))
(export 'evaluate)
(defgeneric evaluate (object point &key)
(:documentation "Evaluate the GSL object."))
;;; Pointer type
(eval-when (:compile-toplevel :load-toplevel :execute)
(defconstant +foreign-pointer-class+ (class-name (class-of (cffi:null-pointer)))
......
;; Evaluation of interpolation functions.
;; Liam Healy, Sun Nov 4 2007 - 18:40
;; Liam Healy Sun Nov 4 2007 - 18:40
;; Time-stamp: <2009-01-28 19:22:15EST evaluation.lisp>
(in-package :gsl)
(defmfun evaluate-interpolation (interpolation xa ya x acceleration)
(defmfun evaluate ((interpolation interpolation) x &key xa ya acceleration)
"gsl_interp_eval"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
:definition :method
:inputs (xa ya)
:c-return :double
:documentation ; FDL
......@@ -14,73 +16,76 @@
point x, using the interpolation object interpolation, data arrays
xa and ya and the accelerator acceleration.")
(defmfun derivative-evaluate-interpolation (interpolation xa ya x acceleration)
(export 'evaluate-derivative)
(defgeneric evaluate-derivative (object point &key)
(:documentation ; FDL
"Find the derivative of an interpolated function for a given point
x, using the interpolation object interpolation, data arrays
xa and ya and the accelerator acceleration."))
(defmfun evaluate-derivative
((interpolation interpolation) x &key xa ya acceleration)
"gsl_interp_eval_deriv"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
:definition :method
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the derivative of an interpolated function for a given point
:c-return :double)
(export 'evaluate-second-derivative)
(defgeneric evaluate-second-derivative (object point &key)
(:documentation ; FDL
"Find the second derivative of an interpolated function for a given point
x, using the interpolation object interpolation, data arrays
xa and ya and the accelerator acceleration.")
xa and ya and the accelerator acceleration."))
(defmfun second-derivative-evaluate-interpolation (interpolation xa ya x acceleration)
(defmfun evaluate-second-derivative
((interpolation interpolation) x &key xa ya acceleration)
"gsl_interp_eval_deriv2"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
((mpointer acceleration) :pointer))
:definition :method
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the second derivative of an interpolated function for a given point
x, using the interpolation object interpolation, data arrays
xa and ya and the accelerator acceleration.")
:c-return :double)
(defmfun integral-evaluate-interpolation (interpolation xa ya x low high acceleration)
(export 'evaluate-integral)
(defgeneric evaluate-integral (object lower-limit upper-limit &key)
(:documentation ; FDL
"Find the numerical integral of an interpolated function over the
range [low, high], using the interpolation object interpolation,
data arrays xa and ya and the accelerator 'acceleration."))
(defmfun evaluate-integral
((interpolation interpolation) low high &key xa ya acceleration)
"gsl_interp_eval_integ"
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer) (x :double)
(((mpointer interpolation) :pointer) (xa :pointer) (ya :pointer)
(low :double) (high :double) ((mpointer acceleration) :pointer))
:definition :method
:inputs (xa ya)
:c-return :double
:documentation ; FDL
"Find the numerical integral of an interpolated function over the
range [low, high], using the interpolation object interpolation,
data arrays xa and ya and the accelerator 'acceleration.")
:c-return :double)
;;; Spline
(defmfun evaluate-spline (spline x acceleration)
(defmfun evaluate ((spline spline) x &key acceleration)
"gsl_spline_eval"
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:c-return :double
:documentation
"Find the interpolated value of y for a given
point x, using the spline object spline, data arrays
xa and ya and the accelerator 'acceleration.")
:definition :method
:c-return :double)
(defmfun derivative-evaluate-spline (spline x acceleration)
(defmfun evaluate-derivative ((spline spline) x &key acceleration)
"gsl_spline_eval_deriv"
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:c-return :double
:documentation ; FDL
"Find the derivative of an interpolated function for a given point
x, using the spline object spline, data arrays
xa and ya and the accelerator acceleration.")
:definition :method
:c-return :double)
(defmfun second-derivative-evaluate-spline (spline x acceleration)
(defmfun evaluate-second-derivative ((spline spline) x &key acceleration)
"gsl_spline_eval_deriv2"
(((mpointer spline) :pointer) (x :double) ((mpointer acceleration) :pointer))
:c-return :double
:documentation ; FDL
"Find the second derivative of an interpolated function for a given point
x, using the spline object spline, data arrays
xa and ya and the accelerator acceleration.")
:definition :method
:c-return :double)
(defmfun integral-evaluate-spline (spline x low high acceleration)
(defmfun evaluate-integral ((spline spline) low high &key acceleration)
"gsl_spline_eval_integ"
(((mpointer spline) :pointer) (x :double) (low :double) (high :double)
(((mpointer spline) :pointer) (low :double) (high :double)
((mpointer acceleration) :pointer))
:c-return :double
:documentation ; FDL
"Find the numerical integral of an interpolated function over the
range [low, high], using the spline object spline,
data arrays xa and ya and the accelerator acceleration.")
:definition :method
:c-return :double)
;; Example spline
;; Liam Healy, Sat Nov 10 2007 - 21:18
;; Time-stamp: <2009-01-26 22:43:41EST spline-example.lisp>
;; Time-stamp: <2009-01-28 19:05:32EST spline-example.lisp>
;; $Id$
(in-package :gsl)
......@@ -22,6 +22,6 @@
collect (+ i (cos (expt i 2))))))
(spline (make-spline *cubic-spline-interpolation* xarr yarr)))
(loop for xi from (maref xarr 0) below (maref xarr 9) by step
collect (list xi (evaluate-spline spline xi acc)))))
collect (list xi (evaluate spline xi :acceleration acc)))))
(save-test interpolation (spline-example 0.1d0))
;; Interpolation types
;; Liam Healy, Sun Nov 4 2007 - 17:41
;; Time-stamp: <2008-12-26 17:23:12EST types.lisp>
;; Time-stamp: <2009-01-28 21:02:31EST types.lisp>
;; $Id$
(in-package :gsl)
......@@ -53,28 +53,29 @@
:documentation ; FDL
"The name of the interpolation type.")
(defmfun interpolation-minimum-size (interpolation)
"gsl_interp_min_size"
((interpolation :pointer))
:c-return :uint
:documentation ; FDL
"The minimum number of points required by the
(export 'minimum-size)
(defgeneric minimum-size (object)
(:documentation ; FDL
"The minimum number of points required by the
interpolation. For example, Akima spline interpolation
requires a minimum of 5 points.")
requires a minimum of 5 points."))
(defmfun minimum-size ((object interpolation))
"gsl_interp_min_size"
((object :pointer))
:definition :method
:c-return :uint)
(defmfun name ((interpolation spline))
(defmfun name ((object spline))
"gsl_spline_name"
((interpolation :pointer))
((object :pointer))
:definition :method
:c-return :string
:documentation ; FDL
"The name of the interpolation type.")
(defmfun spline-minimum-size (interpolation)
(defmfun minimum-size ((object spline))
"gsl_spline_min_size"
((interpolation :pointer))
:c-return :uint
:documentation ; FDL
"The minimum number of points required by the
interpolation. For example, Akima spline interpolation
requires a minimum of 5 points.")
((object :pointer))
:definition :method
:c-return :uint)
This diff is collapsed.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment