Commit 575fe965 authored by Liam Healy's avatar Liam Healy

Direct use of callbacks in functions

Simple function passing to functions which don't use mobjects is
possible with some modifications to defmfun which treat the C argument
'callback and the CL argument 'function specially, generating a
callback structure based on 'gsl-function and filling it
appropriately.  Ported numerical-differentiation and
numerical-integration.  Also simplified calls in numerical-integration
by making the workspace argument optional and at the end of the
arglist; it defaults to the right size based on 'limit.  Changed one
of the tests for numerical-integration because the limit test doesn't
make sense with the workspace automatically allocated.
parent e236d8d7
;; Numerical differentiation.
;; Liam Healy Mon Nov 12 2007 - 22:07
;; Time-stamp: <2009-01-24 14:57:26EST numerical-differentiation.lisp>
;; Time-stamp: <2009-02-15 14:49:44EST numerical-differentiation.lisp>
;; $Id$
(in-package :gsl)
......@@ -12,7 +12,7 @@
(defmfun central-derivative (function x step)
"gsl_deriv_central"
((function :pointer) (x :double) (step :double)
((callback :pointer) (x :double) (step :double)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the numerical derivative of the function
......@@ -32,7 +32,7 @@
(defmfun forward-derivative (function x step)
"gsl_deriv_forward"
((function :pointer) (x :double) (step :double)
((callback :pointer) (x :double) (step :double)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the numerical derivative of the function
......@@ -54,7 +54,7 @@
(defmfun backward-derivative (function x step)
"gsl_deriv_backward"
((function :pointer) (x :double) (step :double)
((callback :pointer) (x :double) (step :double)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the numerical derivative of the function at the point x
......@@ -83,22 +83,29 @@
(defun deriv-f6 (x) (/ x))
(defun deriv-f6-d (x) (- (expt x -2)))
(make-callbacks single-function deriv-f1)
(make-callbacks single-function deriv-f2)
(make-callbacks single-function deriv-f3)
(make-callbacks single-function deriv-f4)
(make-callbacks single-function deriv-f5)
(make-callbacks single-function deriv-f6)
(save-test numerical-differentiation
(central-derivative (make-single-function deriv-f1) 1.0d0 1.0d-4)
(forward-derivative (make-single-function deriv-f1) 1.0d0 1.0d-4)
(backward-derivative (make-single-function deriv-f1) 1.0d0 1.0d-4)
(central-derivative (make-single-function deriv-f2) 0.1d0 1.0d-4)
(forward-derivative (make-single-function deriv-f2) 0.1d0 1.0d-4)
(backward-derivative (make-single-function deriv-f2) 0.1d0 1.0d-4)
(central-derivative (make-single-function deriv-f3) 0.45d0 1.0d-4)
(forward-derivative (make-single-function deriv-f3) 0.45d0 1.0d-4)
(backward-derivative (make-single-function deriv-f3) 0.45d0 1.0d-4)
(central-derivative (make-single-function deriv-f4) 0.5d0 1.0d-4)
(forward-derivative (make-single-function deriv-f4) 0.5d0 1.0d-4)
(backward-derivative (make-single-function deriv-f4) 0.5d0 1.0d-4)
(central-derivative (make-single-function deriv-f5) 0.0d0 1.0d-4)
(forward-derivative (make-single-function deriv-f5) 0.0d0 1.0d-4)
(backward-derivative (make-single-function deriv-f5) 0.0d0 1.0d-4)
(central-derivative (make-single-function deriv-f6) 10.0d0 1.0d-4)
(forward-derivative (make-single-function deriv-f6) 10.0d0 1.0d-4)
(backward-derivative(make-single-function deriv-f6) 10.0d0 1.0d-4))
(central-derivative 'deriv-f1 1.0d0 1.0d-4)
(forward-derivative 'deriv-f1 1.0d0 1.0d-4)
(backward-derivative 'deriv-f1 1.0d0 1.0d-4)
(central-derivative 'deriv-f2 0.1d0 1.0d-4)
(forward-derivative 'deriv-f2 0.1d0 1.0d-4)
(backward-derivative 'deriv-f2 0.1d0 1.0d-4)
(central-derivative 'deriv-f3 0.45d0 1.0d-4)
(forward-derivative 'deriv-f3 0.45d0 1.0d-4)
(backward-derivative 'deriv-f3 0.45d0 1.0d-4)
(central-derivative 'deriv-f4 0.5d0 1.0d-4)
(forward-derivative 'deriv-f4 0.5d0 1.0d-4)
(backward-derivative 'deriv-f4 0.5d0 1.0d-4)
(central-derivative 'deriv-f5 0.0d0 1.0d-4)
(forward-derivative 'deriv-f5 0.0d0 1.0d-4)
(backward-derivative 'deriv-f5 0.0d0 1.0d-4)
(central-derivative 'deriv-f6 10.0d0 1.0d-4)
(forward-derivative 'deriv-f6 10.0d0 1.0d-4)
(backward-derivative 'deriv-f6 10.0d0 1.0d-4))
;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2009-01-25 10:16:31EST numerical-integration.lisp>
;; Time-stamp: <2009-02-15 14:43:39EST numerical-integration.lisp>
;; $Id$
;;; To do: QAWS, QAWO, QAWF, more tests
......@@ -19,7 +19,7 @@
;; what these are if they are too large, it will do a minimum number
;; of points anyway.
"gsl_integration_qng"
((function :pointer)
((callback :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double)
(result :double) (abserr :double) (neval sizet))
......@@ -50,13 +50,14 @@
:gauss41 :gauss51 :gauss61)
(defmfun integration-QAG
(function a b method limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a b method limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
;; what these are if they are too large, it will do a minimum number
;; of points anyway.
"gsl_integration_qag"
((function :pointer)
((callback :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double)
(limit sizet) (method integrate-method) ((mpointer workspace) :pointer)
......@@ -77,7 +78,7 @@
On each iteration the adaptive integration strategy bisects the interval
with the largest error estimate. The subintervals and their results are
stored in the memory provided by workspace. The maximum number of
subintervals is given by limit, which may not exceed the allocated
subintervals is given by 'limit, which may not exceed the allocated
size of the workspace.")
;;;;****************************************************************************
......@@ -85,10 +86,11 @@
;;;;****************************************************************************
(defmfun integration-QAGS
(function a b limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a b limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qags"
((function :pointer)
((callback :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
......@@ -111,9 +113,11 @@
;;;;****************************************************************************
(defmfun integration-QAGP
(function points limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function points limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagp"
((function :pointer)
((callback :pointer)
((mpointer points) :pointer) ((dim0 points) sizet)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
......@@ -135,9 +139,11 @@
;;;;****************************************************************************
(defmfun integration-QAGi
(function limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagi"
((function :pointer)
((callback :pointer)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
:documentation ; FDL
......@@ -152,9 +158,11 @@
this case a lower-order rule is more efficient.")
(defmfun integration-QAGiu
(function a limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagiu"
((function :pointer) (a :double)
((callback :pointer) (a :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
:documentation ; FDL
......@@ -165,9 +173,11 @@
and then integrated using the QAGS algorithm.")
(defmfun integration-QAGil
(function b limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function b limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagil"
((function :pointer) (b :double)
((callback :pointer) (b :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
:documentation ; FDL
......@@ -182,10 +192,11 @@
;;;;****************************************************************************
(defmfun integration-QAWC
(function a b c limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a b c limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qawc"
((function :pointer)
((callback :pointer)
(a :double) (b :double) (c :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
......@@ -202,17 +213,14 @@
the singularity. Further away from the singularity the algorithm
uses an ordinary 15-point Gauss-Kronrod integration rule.")
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(defun one-sine (x) (sin x))
(defparameter *one-sine* (make-single-function one-sine))
(make-callbacks single-function one-sine)
(save-test numerical-integration
(integration-qng *one-sine* 0.0d0 pi)
(let ((ws (make-integration-workspace 20)))
(integration-QAG *one-sine* 0.0d0 pi :gauss15 20 ws))
(let ((ws (make-integration-workspace 20)))
(integration-QAG *one-sine* 0.0d0 pi :gauss15 50 ws)))
(integration-qng 'one-sine 0.0d0 pi)
(integration-QAG 'one-sine 0.0d0 pi :gauss15 20)
(integration-QAG 'one-sine 0.0d0 pi :gauss21 40))
;; Helpers that define a single GSL function interface
;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
;; Time-stamp: <2009-01-31 21:19:28EST defmfun-single.lisp>
;; Time-stamp: <2009-02-15 14:17:36EST defmfun-single.lisp>
;; $Id: $
(in-package :gsl)
......@@ -100,6 +100,9 @@
(cons
'values
(append before after
(when (member +callback-argument-name+
(variables-used-in-c-arguments c-arguments))
'(function)) ; hardwired name to use for callback
(let ((auxstart (position '&aux arglist)))
;; &aux bindings are checked
(when auxstart
......@@ -108,7 +111,7 @@
(mapcar 'rest (subseq arglist (1+ auxstart))))))))))))
,@(when documentation (list documentation))
,(funcall body-maker name arglist gsl-name c-arguments key-args))
`(,definition
`(,definition
,@(when (and name (not (defgeneric-method-p name)))
(list name))
,@(when qualifier (list qualifier))
......@@ -164,7 +167,9 @@
allocated))
(clret (or ; better as a symbol macro
(substitute cret-name :c-return return)
(mapcan #'cl-convert-form allocated-decl)
(mapcan #'cl-convert-form
(remove +callback-argument-name+
allocated-decl :key 'st-symbol))
outputs
(unless (eq c-return :void)
(list cret-name)))))
......@@ -175,6 +180,12 @@
(mapcar #'wfo-declare allocated-decl)
'cffi:with-foreign-objects
`(,@before
,@(when (member +callback-argument-name+ allocated)
`((set-parameters
,+callback-argument-name+ 'gsl-function)
(set-slot-function
,+callback-argument-name+
'gsl-function 'function function)))
(let ((,cret-name
(cffi:foreign-funcall
,gsl-name
......@@ -208,7 +219,8 @@
,@after
(values
,@(defmfun-return
c-return cret-name clret allocated return return-supplied-p
c-return cret-name clret allocated
return return-supplied-p
enumeration outputs)))))))))
(defun defmfun-return
......
;; Macros to interface GSL functions, including definitions necessary for defmfun.
;; Liam Healy
;; Time-stamp: <2009-01-25 10:30:39EST interface.lisp>
;; Time-stamp: <2009-02-15 12:45:17EST interface.lisp>
;; $Id$
(in-package :gsl)
......@@ -61,9 +61,11 @@
(defun wfo-declare (d)
`(,(st-symbol d)
,@(if (st-arrayp d)
`(',(st-eltype d) ,(st-dim d))
`(',(st-type d)))))
,@(if (st-arrayp d)
`(',(st-eltype d) ,(st-dim d))
(if (eq (st-symbol d) +callback-argument-name+)
'('gsl-function)
`(',(st-type d))))))
;;;;****************************************************************************
;;;; Checking results from GSL functions
......
......@@ -6,91 +6,73 @@
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.718281828441488d0 4.123659913167868d-10)
(MULTIPLE-VALUE-LIST
(CENTRAL-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F1)
1.0d0 1.d-4)))
(CENTRAL-DERIVATIVE 'DERIV-F1 1.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.7182817825298398d0 9.540320743340577d-7)
(MULTIPLE-VALUE-LIST
(FORWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F1)
1.0d0 1.d-4)))
(FORWARD-DERIVATIVE 'DERIV-F1 1.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.7182818293827378d0 9.086969640080628d-7)
(MULTIPLE-VALUE-LIST
(BACKWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F1)
1.0d0 1.d-4)))
(BACKWARD-DERIVATIVE 'DERIV-F1 1.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.474341649024533d0 3.3922603135575853d-11)
(MULTIPLE-VALUE-LIST
(CENTRAL-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F2)
0.1d0 1.d-4)))
(CENTRAL-DERIVATIVE 'DERIV-F2 0.1d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.47434164943235285d0 9.102230432046565d-8)
(MULTIPLE-VALUE-LIST
(FORWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F2)
0.1d0 1.d-4)))
(FORWARD-DERIVATIVE 'DERIV-F2 0.1d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.4743416462049176d0 8.787851752748856d-8)
(MULTIPLE-VALUE-LIST
(BACKWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F2)
0.1d0 1.d-4)))
(BACKWARD-DERIVATIVE 'DERIV-F2 0.1d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.9941773945923913d0 6.730277498876722d-10)
(MULTIPLE-VALUE-LIST
(CENTRAL-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F3)
0.45d0 1.d-4)))
(CENTRAL-DERIVATIVE 'DERIV-F3 0.45d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.994177385772526d0 1.6952086150997503d-6)
(MULTIPLE-VALUE-LIST
(FORWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F3)
0.45d0 1.d-4)))
(FORWARD-DERIVATIVE 'DERIV-F3 0.45d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.9941775182635064d0 1.568063267060477d-6)
(MULTIPLE-VALUE-LIST
(BACKWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F3)
0.45d0 1.d-4)))
(BACKWARD-DERIVATIVE 'DERIV-F3 0.45d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -0.7788007830653751d0 1.7227267499194736d-10)
(MULTIPLE-VALUE-LIST
(CENTRAL-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F4)
0.5d0 1.d-4)))
(CENTRAL-DERIVATIVE 'DERIV-F4 0.5d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -0.7788007884386108d0 2.530006603059119d-7)
(MULTIPLE-VALUE-LIST
(FORWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F4)
0.5d0 1.d-4)))
(FORWARD-DERIVATIVE 'DERIV-F4 0.5d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -0.7788007743421791d0 2.503617930568489d-7)
(MULTIPLE-VALUE-LIST
(BACKWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F4)
0.5d0 1.d-4)))
(BACKWARD-DERIVATIVE 'DERIV-F4 0.5d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 0.0d0 6.66133814775094d-20)
(MULTIPLE-VALUE-LIST
(CENTRAL-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F5)
0.0d0 1.d-4)))
(CENTRAL-DERIVATIVE 'DERIV-F5 0.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -2.6337311035821125d-26 1.606762230663824d-11)
(MULTIPLE-VALUE-LIST
(FORWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F5)
0.0d0 1.d-4)))
(FORWARD-DERIVATIVE 'DERIV-F5 0.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.6337311035821125d-26 1.606762230663824d-11)
(MULTIPLE-VALUE-LIST
(BACKWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F5)
0.0d0 1.d-4)))
(BACKWARD-DERIVATIVE 'DERIV-F5 0.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -0.009999999999871223d0 2.6645352591887814d-12)
(MULTIPLE-VALUE-LIST
(CENTRAL-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F6)
10.0d0 1.d-4)))
(CENTRAL-DERIVATIVE 'DERIV-F6 10.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -0.010000000266238856d0 4.9950606958271945d-9)
(MULTIPLE-VALUE-LIST
(FORWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F6)
10.0d0 1.d-4)))
(FORWARD-DERIVATIVE 'DERIV-F6 10.0d0 1.d-4)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST -0.010000000076196142d0 4.641550044799095d-9)
(MULTIPLE-VALUE-LIST
(BACKWARD-DERIVATIVE (MAKE-SINGLE-FUNCTION DERIV-F6)
10.0d0 1.d-4))))
(BACKWARD-DERIVATIVE 'DERIV-F6 10.0d0 1.d-4))))
......@@ -6,19 +6,13 @@
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.0d0 2.220446049250313d-14 21)
(MULTIPLE-VALUE-LIST
(INTEGRATION-QNG *ONE-SINE* 0.0d0 PI)))
(INTEGRATION-QNG 'ONE-SINE 0.0d0 PI)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.0d0 2.220446049250313d-14)
(MULTIPLE-VALUE-LIST
(LET ((WS (MAKE-INTEGRATION-WORKSPACE 20)))
(INTEGRATION-QAG *ONE-SINE* 0.0d0 PI :GAUSS15 20
WS))))
(LISP-UNIT:ASSERT-ERROR 'INVALID-ARGUMENT
(LET ((WS
(MAKE-INTEGRATION-WORKSPACE
20)))
(INTEGRATION-QAG *ONE-SINE*
0.0d0 PI
:GAUSS15 50
WS))))
(INTEGRATION-QAG 'ONE-SINE 0.0d0 PI :GAUSS15 20)))
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 2.0d0 2.220446049250313d-14)
(MULTIPLE-VALUE-LIST
(INTEGRATION-QAG 'ONE-SINE 0.0d0 PI :GAUSS21 40))))
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