Commit 7b4aa0c7 authored by liam's avatar liam

Callback functions are now defined with def-scalar-function and

associated macros.  This places them in the gsl-function structure.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3257 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent db7e17a8
......@@ -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))))
;********************************************************
; 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))
......@@ -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))))
......@@ -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))
......@@ -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
......
......@@ -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
......
......@@ -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)))
......@@ -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))))
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