diff --git a/chebyshev.lisp b/chebyshev.lisp index 8113f5c5b149006609104bf1d4bf4850f10cbe9a..3b8e82f8ddbcaa3535ac29b659573e4788c306e2 100644 --- a/chebyshev.lisp +++ b/chebyshev.lisp @@ -3,7 +3,7 @@ ; description: Chebyshev Approximations ; date: Sat Nov 17 2007 - 20:36 ; author: Liam Healy -; modified: Mon Dec 10 2007 - 18:53 +; modified: Sun Dec 30 2007 - 16:05 ;******************************************************** ;;; $Id: $ @@ -119,21 +119,39 @@ ;;; From Chap. 28.5, except I have set steps = 100 instead of 10000 ;;; to keep things sane. -;;; Calling a callback from CL is not possible, so we define a -;;; separate CL function to show what the answer really is. -(defun chebyshev-step (x) (if (< x 0.5) 0.25 0.75)) -(def-gsl-function chebyshev-step-c x (chebyshev-step x)) +(defun-scalar chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0)) -(defun chebyshev-example () +(defun chebyshev-table-example () (let ((cheb (allocate-chebyshev 40)) (steps 100)) - (with-integration-function (step-fn 'chebyshev-step-c) - (initialize-chebyshev cheb step-fn 0.0d0 1.0d0) - (dotimes (i steps) - (let ((x (coerce (/ i steps) 'double-float))) - (format t "~&~a ~a ~a ~a" - x - (chebyshev-step x) - (evaluate-chebyshev cheb x 10) - (evaluate-chebyshev cheb x)))) - (free-chebyshev cheb)))) + (initialize-chebyshev cheb chebyshev-step 0.0d0 1.0d0) + (dotimes (i steps) + (let ((x (coerce (/ i steps) 'double-float))) + (format t "~&~a ~a ~a ~a" + x + (chebyshev-step x) + (evaluate-chebyshev cheb x 10) + (evaluate-chebyshev cheb x)))) + (free-chebyshev cheb))) + +(defun chebyshev-point-example (x) + (check-type x double-float) + (let ((cheb (allocate-chebyshev 40)) + (deriv (allocate-chebyshev 40)) + (integ (allocate-chebyshev 40))) + (initialize-chebyshev cheb chebyshev-step 0.0d0 1.0d0) + (derivative-chebyshev deriv cheb) + (integral-chebyshev integ cheb) + (prog1 + (list + (evaluate-chebyshev cheb x) + (evaluate-chebyshev deriv x) + (evaluate-chebyshev integ x)) + (free-chebyshev cheb) + (free-chebyshev deriv) + (free-chebyshev integ)))) + +(lisp-unit:define-test chebyshev + (lisp-unit:assert-equal + '("0.715920990169d+00" "-0.150199666581d+01" "0.172397194040d+00") + (lisp-unit:fp-sequence (chebyshev-point-example 0.55d0)))) diff --git a/general/functions.lisp b/general/functions.lisp new file mode 100644 index 0000000000000000000000000000000000000000..3df5040cbf8051199e112a84e985e3a7cb7e01aa --- /dev/null +++ b/general/functions.lisp @@ -0,0 +1,119 @@ +;******************************************************** +; file: functions.lisp +; description: Foreign callback functions. +; date: Sun Dec 9 2007 - 22:08 +; author: Liam Healy +; modified: Sun Dec 30 2007 - 16:29 +;******************************************************** +;;; $Id: $ + +(in-package :gsl) + +;;; Callback functions are functions which are passed as data; to Lisp +;;; that means they are just functions, but C makes a distinction. +;;; They are needed by several GSL tasks. +;;; Functions that take one double-float and return one double-float +;;; and are defined using the gsl_function structure (see gsl_math.h) +;;; are called scalar functions. They are used in +;;; numerical-integration, numerical-differentiation, chebyshev and +;;; definitions to aid in creating and using them are provided here. +;;; The idea behind the definitions is that the boundary between the CL +;;; and C functions is a narrow as possible; #'defun-scalar is +;;; macro that allows one to define a function in Lisp and make +;;; use it in these GSL tasks. + +;;; Other GSL tasks make use of callback functions with different +;;; characteristics. Since they are specific to each of the tasks, +;;; they are defined with those tasks. A complexity encountered with +;;; using a vector of double floats not present with scalars is that +;;; there is no portable way to make a C array available directly to +;;; Lisp, so callbacks using C arrays must read them using a macro. +;;; Therefore it is necessary to define the function in a way that +;;; prevents its use in Lisp; to ameliorate this, the macro +;;; #'with-c-vector is provided to give named access to the elements. + +(export + '(def-scalar-function undef-scalar-function defun-scalar with-c-vector)) + +;;;;**************************************************************************** +;;;; Setting slots +;;;;**************************************************************************** + +(defun set-structure-slot (foreign-structure structure-name slot-name value) + (setf (cffi:foreign-slot-value foreign-structure structure-name slot-name) + value)) + +(defun set-slot-function (foreign-structure structure-name slot-name gsl-function) + (set-structure-slot + foreign-structure structure-name slot-name + (cffi:get-callback gsl-function))) + +(defun set-parameters (foreign-structure structure-name) + "Set the parameters slot to null." + (set-structure-slot foreign-structure structure-name + 'parameters (cffi:null-pointer))) + +;;;;**************************************************************************** +;;;; The GSL struct for holding callbacks +;;;;**************************************************************************** + +(cffi:defcstruct gsl-function + "Passing functions to GSL." + ;; see /usr/include/gsl/gsl_math.h + (function :pointer) + (parameters :pointer)) + +;;;;**************************************************************************** +;;;; A function of a scalar double +;;;;**************************************************************************** + +;;; Used by numerical-integration, numerical-differentiation, chebyshev, ntuple. + +(defmacro def-scalar-function + (name argument &optional (return-type :double) (argument-type :double)) + "Define the variable given by name + as a foreign gsl-function that contains the callback + of a CL function of the same name." + `(progn + (cffi:defcallback ,name ,return-type + ((,argument ,argument-type) (params :pointer)) + (declare (ignore params)) + (,name ,argument)) + ;; Assume that defcallback does not bind the variable 'name. + (defparameter ,name (cffi:foreign-alloc 'gsl-function)) + (set-slot-function ,name 'gsl-function 'function ',name) + (set-parameters ,name 'gsl-function))) + +(defun undef-scalar-function (name) + "Free foreign callback function. It is not necessary to do this; think + of the memory taken by an unused foreign function as much + less than that used by an unused defun." + (cffi:foreign-free name) + (makunbound name)) + +;;; Combine a defun and def-scalar-function in one: +(defmacro defun-scalar (name arglist &body body) + "Define a function of a scalar double-float argument returning + a double-float in CL and C." + `(progn + (defun ,name ,arglist ,@body) + (def-scalar-function ,name ,@arglist))) + +;;;;**************************************************************************** +;;;; Vector of doubles +;;;;**************************************************************************** + +;;; Unfortunately, because vectors must be copied between languages +;;; (even with vector-sap in callbacks, unless vector-sap can be +;;; setfed), there is no way to provide a function of a vector and +;;; have it work in both languages. As a consolation the macro +;;; #'with-c-vector is provided to make things easier. + +(defmacro with-c-vector + ((c-vector &rest element-names) &body body) + "Provide named access to each element of a C array of doubles, for either + reading or setting." + `(symbol-macrolet + ,(loop for i from 0 for a in element-names + collect `(,a (double-to-cl ,c-vector ,i))) + ,@body)) diff --git a/gsll.asd b/gsll.asd index 5bbf320a0eca6f2f4526513ae2cf03b9bd00035f..47673713e96439f6a2ae1e4e52f8cc8372e8e7c7 100644 --- a/gsll.asd +++ b/gsll.asd @@ -3,7 +3,7 @@ ; description: Definition of GSLL system ; date: ; author: Liam Healy -; modified: Sat Dec 8 2007 - 18:41 +; modified: Sun Dec 30 2007 - 15:46 ;******************************************************** ;;; $Id: $ @@ -20,7 +20,7 @@ ((:file "init") (:file "utility" :depends-on (init)) (:file "number-conversion" :depends-on (init)) - (:file "interface" :depends-on (init)) + (:file "interface" :depends-on (init number-conversion)) ;; http://www.cs.northwestern.edu/academics/courses/325/readings/lisp-unit.html (:file "lisp-unit") (:file "tests" :depends-on (init lisp-unit)))) @@ -28,7 +28,8 @@ :depends-on (init) :components ((:file "conditions") - (:file "mathematical"))) + (:file "mathematical") + (:file "functions"))) ;; complex numbers not necessary? Just make a struct. (:module data :depends-on (init) @@ -89,7 +90,7 @@ (:file "householder"))) (:file "eigensystems" :depends-on (init data)) ;; Skip fft for now, I'm not sure how it works in C - (:file "numerical-integration" :depends-on (init)) + (:file "numerical-integration" :depends-on (init general)) (:module random :depends-on (init) :components @@ -151,7 +152,7 @@ (:file "read-write" :depends-on (histogram)) (:file "probability-distribution" :depends-on (histogram)) (:file "ntuple"))) - (:file "monte-carlo" :depends-on (init data random numerical-integration)) + (:file "monte-carlo" :depends-on (init data random general)) (:module ordinary-differential-equations :depends-on (init) :components @@ -170,8 +171,8 @@ ;; Once common callback definitions are moved out of ;; numerical-integration to something in init, the ;; numerical-integration dependency won't be necessary. - (:file "numerical-differentiation" :depends-on (init numerical-integration)) - (:file "chebyshev" :depends-on (init numerical-integration)) + (:file "numerical-differentiation" :depends-on (init general)) + (:file "chebyshev" :depends-on (init general)) (:file "series-acceleration" :depends-on (init)) (:file "wavelet" :depends-on (init data)) (:file "hankel" :depends-on (init data)))) diff --git a/histogram/ntuple.lisp b/histogram/ntuple.lisp index 8b0a6d6c30afbf59ba4e1f0aa1a9a82e4713cd1c..337541c0839dc19e37cbb41a2f2ee54e3816d128 100644 --- a/histogram/ntuple.lisp +++ b/histogram/ntuple.lisp @@ -3,7 +3,7 @@ ; description: N-tuples ; date: Sat Feb 3 2007 - 12:53 ; author: Liam Healy -; modified: Sat Feb 3 2007 - 14:32 +; modified: Sun Dec 30 2007 - 16:33 ;******************************************************** ;;; $Id: $ @@ -81,9 +81,19 @@ the histogram, so subsequent calls can be used to accumulate further data in the same histogram.") -;;; This works for either select-function or value-function. -;;; The functions themselves will have to be defined with cffi:defcallback. -(cffi:defcstruct ntuple-function - (function :pointer) - (parameters :pointer)) +;;; Callback definitions +(export '(def-ntuple-select-function def-ntuple-value-function)) + +(defmacro def-ntuple-select-function (name arg) + "The selection function determines which ntuple rows are selected + for histogramming. The struct component function should return a + non-zero value for each ntuple row that is to be included in the + histogram. " + `(def-scalar-function ,name ,arg :int :pointer)) + +(defmacro def-ntuple-value-function (name arg) + "The value function computes scalar values for those ntuple rows + selected by the selection function which should return the value + to be added to the histogram." + `(def-scalar-function ,name ,arg :double :pointer)) diff --git a/init/interface.lisp b/init/interface.lisp index 90f15f0d86e5c32703be5eec6fe07d12dc05ff23..e2f299c15da2930f5efeb680254d35c45a815f3e 100644 --- a/init/interface.lisp +++ b/init/interface.lisp @@ -3,7 +3,7 @@ ; description: Macros to interface GSL functions. ; date: Mon Mar 6 2006 - 22:35 ; author: Liam M. Healy -; modified: Sat Nov 17 2007 - 21:35 +; modified: Sun Dec 9 2007 - 17:51 ;******************************************************** (in-package :gsl) @@ -95,7 +95,9 @@ (defun success-failure (value) "If status indicates failure, return NIL, othewise return T." - (not (eql value (cffi:foreign-enum-value 'gsl-errorno :FAILURE)))) + ;;(not (eql value (cffi:foreign-enum-value 'gsl-errorno :FAILURE))) + ;; More general, to allow :CONTINUE + (not (minusp value))) ;;;;**************************************************************************** ;;;; Macro defun-gsl diff --git a/monte-carlo.lisp b/monte-carlo.lisp index 73e5e1f54ecc3c8a73c77846d9377590c54d778f..1403c4f9f8b1edeb2c44ca4fe729f062edd712b4 100644 --- a/monte-carlo.lisp +++ b/monte-carlo.lisp @@ -3,17 +3,12 @@ ; description: Monte Carlo Integration ; date: Sat Feb 3 2007 - 17:42 ; author: Liam Healy -; modified: Sun Feb 11 2007 - 12:03 +; modified: Sun Dec 30 2007 - 15:45 ;******************************************************** ;;; $Id: $ (in-package :gsl) -(cffi:defcstruct monte-function - (function :pointer) - (dimensions :size) - (parameters :pointer)) - ;;;;**************************************************************************** ;;;; PLAIN Monte Carlo ;;;;**************************************************************************** @@ -292,6 +287,28 @@ is returned via the state struct component, @var{s->chisq}, and must be consistent with 1 for the weighted average to be reliable.") +;;;;**************************************************************************** +;;;; Callback definition +;;;;**************************************************************************** + +(cffi:defcstruct monte-function + (function :pointer) + (dimensions :size) + (parameters :pointer)) + +(export 'def-mc-function) +(defmacro def-mc-function (name arg dimensions) + `(progn + (cffi:defcallback ,name :double + ((,arg :pointer) (params :pointer)) + (declare (ignore params)) + (,name ,arg)) + ;; Assume that defcallback does not bind the variable 'name. + (defparameter ,name (cffi:foreign-alloc 'monte-function)) + (set-slot-function ,name 'monte-function 'function ',name) + (set-structure-slot ,name 'monte-function 'dimensions ,dimensions) + (set-parameters ,name 'monte-function))) + ;;;;**************************************************************************** ;;;; Examples and unit test ;;;;**************************************************************************** @@ -299,9 +316,12 @@ ;;; Example from Sec. 23.5 ;;; This is a function that occurs in random walk studies. -(def-gsl-function monte-carlo-g (x y z) - (* (/ (expt pi 3)) - (/ (- 1 (* (cos x) (cos y) (cos z)))))) +(defun monte-carlo-g (arg) + (with-c-vector (arg x y z) + (* (/ (expt pi 3)) + (/ (- 1 (* (cos x) (cos y) (cos z))))))) + +(def-mc-function monte-carlo-g arg 3) (defun random-walk-plain-example (&optional (nsamples 500000)) (with-monte-carlo-plain (ws 3) @@ -310,13 +330,12 @@ (setf (data lower) #(0.0d0 0.0d0 0.0d0) (data upper) (vector pi pi pi)) (rng-set *rng-mt19937* 0) - (with-integration-function (mcf 'monte-carlo-g 3) - (monte-carlo-integrate-plain - mcf - lower upper - nsamples - *rng-mt19937* - ws)))))) + (monte-carlo-integrate-plain + monte-carlo-g + lower upper + nsamples + *rng-mt19937* + ws))))) (defun random-walk-miser-example (&optional (nsamples 500000)) (with-monte-carlo-miser (ws 3) @@ -325,13 +344,12 @@ (setf (data lower) #(0.0d0 0.0d0 0.0d0) (data upper) (vector pi pi pi)) (rng-set *rng-mt19937* 0) - (with-integration-function (mcf 'monte-carlo-g 3) - (monte-carlo-integrate-miser - mcf - lower upper - nsamples - *rng-mt19937* - ws)))))) + (monte-carlo-integrate-miser + monte-carlo-g + lower upper + nsamples + *rng-mt19937* + ws))))) (defun random-walk-vegas-example (&optional (nsamples 500000)) (with-monte-carlo-vegas (ws 3) @@ -340,13 +358,12 @@ (setf (data lower) #(0.0d0 0.0d0 0.0d0) (data upper) (vector pi pi pi)) (rng-set *rng-mt19937* 0) - (with-integration-function (mcf 'monte-carlo-g 3) - (monte-carlo-integrate-vegas - mcf - lower upper - nsamples - *rng-mt19937* - ws)))))) + (monte-carlo-integrate-vegas + monte-carlo-g + lower upper + nsamples + *rng-mt19937* + ws))))) (lisp-unit:define-test monte-carlo (lisp-unit:assert-first-fp-equal diff --git a/numerical-differentiation.lisp b/numerical-differentiation.lisp index 555601182f2e3357fe421e29d144b24db878bca2..1d5ee971c0bf7566b1947ab9fa7cdf0199c42c97 100644 --- a/numerical-differentiation.lisp +++ b/numerical-differentiation.lisp @@ -3,7 +3,7 @@ ; description: Numerical differentiation. ; date: Mon Nov 12 2007 - 22:07 ; author: Liam Healy -; modified: Mon Nov 12 2007 - 22:33 +; modified: Sun Dec 30 2007 - 14:25 ;******************************************************** ;;; $Id: $ @@ -71,13 +71,14 @@ @var{x}. This function is equivalent to calling #'forward-derivative with a negative step-size.") -#| -;;; This is the example given in the GSL manual, Sec. 27.2. - -(def-gsl-function 3/2-power x (expt x 3/2)) +;;;; Examples and unit test -(with-integration-function (fn '3/2-power) - (central-derivative fn 2.0d0 1.0d-8)) +;;; This is the example given in the GSL manual, Sec. 27.2. +(defun-scalar 3/2-power (x) (expt x 3/2)) +;;; (3/2-power 2.0d0) -;;; Compare with (* 3/2 (sqrt 2.0d0)). -|# +(lisp-unit:define-test numerical-differentiation + (lisp-unit:assert-first-fp-equal + "0.212132031200d+01" + ;; Compare to (* 3/2 (sqrt 2.0d0)) + (central-derivative 3/2-power 2.0d0 1.d-8))) diff --git a/numerical-integration.lisp b/numerical-integration.lisp index aa86c34eeb76e6f6c9811ed1a21e0a91d82eb0de..6199d4d1473c545045cbbef05739b06255d15119 100644 --- a/numerical-integration.lisp +++ b/numerical-integration.lisp @@ -3,7 +3,7 @@ ; description: Numerical integration ; date: Wed Jul 5 2006 - 23:14 ; author: Liam M. Healy -; modified: Sat Nov 17 2007 - 21:51 +; modified: Sun Dec 30 2007 - 14:19 ;******************************************************** ;;; $Id: $ @@ -11,51 +11,6 @@ (in-package :gsl) -(cffi:defcstruct gsl-function - "Passing functions to GSL." - ;; see /usr/include/gsl/gsl_math.h - (function :pointer) - (parameters :pointer)) - -(export 'def-gsl-function) -(defmacro def-gsl-function (name arg &body body) - "Define a GSL (C) function of either one argument of type - double (if arg is a symbol), or a C array of doubles - (if arg is a list), for GSL numerical - integration functions. - Parameters (non integration variables) may be passed by - using a lexical closure." - (let ((argvec (gensym "MCARG"))) - `(cffi:defcallback ,name :double - (,(if (listp arg) - `(,argvec :pointer) - `(,arg :double)) - (params :pointer)) - (declare (ignore params)) - ,@(if (listp arg) - `((symbol-macrolet - ,(loop for i from 0 for a in arg - collect `(,a (cffi:mem-aref ,argvec :double ,i))) - ,@body)) - body)))) - -(export 'with-integration-function) -(defmacro with-integration-function - ((name function &optional number-of-arguments) &body body) - "Make a function for GSL to integrate." - ;; Is this creating and deallocating the gsl-function object? - (let ((structure (if number-of-arguments 'monte-function 'gsl-function))) - `(cffi:with-foreign-object (,name ',structure) - (setf (cffi:foreign-slot-value ,name ',structure 'function) - (cffi:get-callback ,function) - ,@(when number-of-arguments - `((cffi:foreign-slot-value ,name ',structure 'dimensions) - ,number-of-arguments)) - ;; Pass the parameters in with a closure. - (cffi:foreign-slot-value ,name ',structure 'parameters) - (cffi:null-pointer)) - ,@body))) - ;;;;**************************************************************************** ;;;; QNG non-adaptive Gauss-Kronrod integration ;;;;**************************************************************************** @@ -264,23 +219,22 @@ ;;;; Examples and unit test ;;;;**************************************************************************** -(let ((mult 2.0d0)) - (def-gsl-function two-sine x (sin (* mult x)))) +;;; defun-scalar is a combination of defun and def-scalar-function: +(defun-scalar one-sine (x) (sin x)) -(def-gsl-function one-sine x (sin x)) +;;; Parameters may be defined through the lexical environment: +(let ((mult 2.0d0)) + (defun-scalar two-sine (x) (sin (* mult x)))) (lisp-unit:define-test numerical-integration (lisp-unit:assert-first-fp-equal "0.200000000000d+01" - (with-integration-function (os 'one-sine) - (integration-qng os 0.0d0 pi))) + (integration-qng one-sine 0.0d0 pi)) (lisp-unit:assert-first-fp-equal "0.200000000000d+01" - (with-integration-function (os 'one-sine) - (with-integration-workspace (ws 20) - (integration-QAG os 0.0d0 pi :gauss15 20 ws)))) + (with-integration-workspace (ws 20) + (integration-QAG one-sine 0.0d0 pi :gauss15 20 ws))) (lisp-unit:assert-error 'gsl-error - (with-integration-function (os 'one-sine) - (with-integration-workspace (ws 20) - (integration-QAG os 0.0d0 pi :gauss15 50 ws))))) + (with-integration-workspace (ws 20) + (integration-QAG one-sine 0.0d0 pi :gauss15 50 ws))))