Skip to content
Snippets Groups Projects
Commit c83788ed authored by Liam Healy's avatar Liam Healy
Browse files

Redefine defmcallback to provide and accept scalars from CL functions

New definition for defmcallback will now convert C array pointers
passed to the function by GSL into individual scalar arguments of the
CL function, and will take multiple value returns and assign as
elements of the return C arrays.  Applied to def-ode-functions and the
van der Pol oscillator example.  Redefined #'maref-function-picker as
a macro from a function, eliminating the need for an eval-when form.
parent 1b22bb29
No related branches found
No related tags found
No related merge requests found
;; Get/set array or elements: cl-array, maref
;; Liam Healy 2008-08-27 22:43:10EDT maref.lisp
;; Time-stamp: <2009-01-12 09:43:16EST maref.lisp>
;; Time-stamp: <2009-01-17 17:15:19EST maref.lisp>
;; $Id: $
(in-package :gsl)
......@@ -89,8 +89,8 @@
;;; callback functions in solve-minimize-fit and in #'cl-array method
;;; for pointers.
(eval-when (:compile-toplevel)
(defun maref-function-picker (type-symbol category ffrestargs &optional value-symbol)
(defmacro maref-function-picker
(type-symbol category ffrestargs &optional value-symbol)
"Generate sexp to select on the various gsl_{vector,matrix}*_{get,set} functions."
(cons 'cond
(mapcar (lambda (tp)
......@@ -103,18 +103,18 @@
,@(if value-symbol
(list (cl-cffi tp) value-symbol)
(list (cl-cffi tp))))))
*array-element-types*))))
*array-element-types*)))
(defmethod maref
((pointer #.+foreign-pointer-class+) index &optional index2
(type 'double-float))
(if index2
#.(maref-function-picker
'type 'matrix
'(:pointer pointer sizet index sizet index2))
#.(maref-function-picker
'type 'vector
'(:pointer pointer sizet index))))
(maref-function-picker
type matrix
(:pointer pointer sizet index sizet index2))
(maref-function-picker
type vector
(:pointer pointer sizet index))))
;;; Index the GSL function names to maref
#.(cons 'progn
......@@ -138,12 +138,12 @@
(value (pointer #.+foreign-pointer-class+) index &optional index2
(type 'double-float))
(if index2
#.(maref-function-picker
'type 'matrix
'(:pointer pointer sizet index sizet index2) 'value)
#.(maref-function-picker
'type 'vector
'(:pointer pointer sizet index) 'value)))
(maref-function-picker
type matrix
(:pointer pointer sizet index sizet index2) value)
(maref-function-picker
type vector
(:pointer pointer sizet index) value)))
;;; Index the GSL function names to (setf maref)
#.(cons 'progn
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-04 11:26:56EST gsll-tests.asd>
;; Time-stamp: <2009-01-18 19:02:18EST gsll-tests.asd>
;; $Id$
(asdf:defsystem "gsll-tests"
......@@ -120,6 +120,7 @@
(:file "negative-binomial")
(:file "numerical-differentiation")
(:file "numerical-integration")
(:file "ode")
(:file "pareto")
(:file "permutation")
(:file "poisson")
......
;; Foreign callback functions.
;; Liam Healy
;; Time-stamp: <2008-02-16 10:52:36EST callback.lisp>
;; Time-stamp: <2009-01-18 17:04:46EST callback.lisp>
;; $Id$
(in-package :gsl)
......@@ -61,17 +61,73 @@
(parameters :pointer))
;;;;****************************************************************************
;;;; Macros for defining a callback and placing in a structure
;;;; Macros for defining a callback to wrap a CL function
;;;;****************************************************************************
;;; Usage example for scalar function (e.g. numerical-integration,
;;; numerical-differentiation, chebyshev, ntuple).
;;; (defmcallback myfn :double :double)
;;; Usage example for gsl-vector function (e.g. roots-multi)
;;; Usage example for vector function (e.g. roots-multi)
;;; (defmcallback myfn :pointer :int (:pointer))
;;; Usage example for function and derivative
;;; (defmcallback fdf :pointer :double (:pointer :pointer))
;;; (defmcallback fdf :success-failure :int (:pointer :pointer))
;;; Usage example for def-ode-functions
;;; (defmcallback vanderpol :success-failure (:double (:double 2) (:set :double 2)))
;;; or
;;; (defmcallback vanderpol :success-failure (:double :pointer :pointer))
;;; to read and set within the CL function with #'dcref.
;;; (callback-args '(:double (:double 2) (:set :double 2)))
;;; ((#:ARG1193 :DOUBLE) (#:ARG1194 :POINTER) (#:ARG1195 :POINTER))
(defun callback-args (types)
"The arguments passed by GSL to the callback function."
(mapcar (lambda (type)
(let ((symbol (gensym "ARG")))
(list symbol
(if (listp type) ; like (:double 3)
:pointer ; C array
type))))
(if (listp types) types (list types))))
;;; (embedded-clfunc-args '(:double (:double 2) (:set :double 2)) (callback-args '(:double (:double 2) (:set :double 2))))
;;; (#:ARG1244 (MEM-AREF #:ARG1245 ':DOUBLE 0) (MEM-AREF #:ARG1245 ':DOUBLE 1))
(defun embedded-clfunc-args (types callback-args)
"The arguments passed to the CL function call embedded in the callback."
(loop for spec in types
for (symbol type) in callback-args
append
(unless (and (listp spec) (eq (first spec) :set))
(if (listp spec)
(loop for ind from 0 below (second spec)
collect `(cffi:mem-aref ,symbol ',(first spec) ,ind))
(list symbol)))))
(defun callback-set-mvb (form types callback-args)
"Create the multiple-value-bind form in the callback to set the return C arrays."
(multiple-value-bind (settype setcba)
(loop for cba in callback-args
for type in types
for setting = (and (listp type) (eq (first type) :set))
when setting
collect cba into setcba
when setting
collect type into settype
finally (return (values (mapcar 'rest settype) setcba)))
(print settype)
(let* ((setvbls (embedded-clfunc-args settype setcba))
(count (apply '+ (mapcar 'second settype)))
(mvbvbls (loop repeat count collect (gensym "SETCB"))))
(if (zerop count)
form
`(multiple-value-bind ,mvbvbls
,form
(setf ,@(loop for mvbvbl in mvbvbls
for setvbl in setvbls
append (list setvbl mvbvbl))))))))
;;; (DEFMCALLBACK VANDERPOL :SUCCESS-FAILURE (:DOUBLE (:DOUBLE 2) (:SET :DOUBLE 2)))
;;; (DEFMCALLBACK VANDERPOL :SUCCESS-FAILURE (:DOUBLE (:DOUBLE 2)))
(defmacro defmcallback
(name &optional (return-type :double) (argument-types :double)
......@@ -81,33 +137,39 @@
argument-types is a single type or list of types of the argument(s)
that appear before parameters, and the additional-argument-types
(default none) is a single type or list of types of the argument(s)
that appear after parameters. The return-type is the type that
should be returned to GSL. If :success-failure, a GSL_SUCCESS
code (0) is always returned; if :pointer, a null pointer is
returned."
(flet ((arg-type (types)
(when types
(mapcar (lambda (type) (list (gensym "ARG") type))
(if (listp types) types (list types))))))
(let ((arguments (arg-type argument-types))
(additional-arguments (arg-type additional-argument-types)))
`(cffi:defcallback ,name
,(if (eq return-type :success-failure) :int return-type)
(,@arguments (params :pointer) ,@additional-arguments)
;; Parameters as C argument are always ignored, because we have
;; CL specials to do the same job.
(declare (ignore params))
(,name ,@(mapcar #'first (append arguments additional-arguments)))
,@(case
return-type
(:success-failure
;; We always return success, because if there was a
;; problem, a CL error would be signalled.
'(success))
(:pointer
;; For unclear reasons, some GSL functions want callbacks
;; to return a void pointer which is apparently meaningless.
'((cffi:null-pointer))))))))
that appear after parameters. The argument types are C types or
a list of a C type and a length, indicating a C array of that type
for which each element will be passed as a separate argument.
The return-type is the type that should be returned to GSL.
If :success-failure, a GSL_SUCCESS code (0) is always returned;
if :pointer, a null pointer is returned."
(let* ((atl (if (listp argument-types) argument-types (list argument-types)))
(aatl (if (listp additional-argument-types) additional-argument-types
(list additional-argument-types)))
(cbargs (callback-args atl))
(cbaddl (callback-args aatl)))
`(cffi:defcallback ,name
,(if (eq return-type :success-failure) :int return-type)
(,@cbargs (params :pointer) ,@cbaddl)
;; Parameters as C argument are always ignored, because we have
;; CL specials to do the same job.
(declare (ignore params))
,(callback-set-mvb
`(,name
,@(append
(embedded-clfunc-args atl cbargs) (embedded-clfunc-args aatl cbaddl)))
(append atl aatl)
(append cbargs cbaddl))
,@(case
return-type
(:success-failure
;; We always return success, because if there was a
;; problem, a CL error would be signalled.
'(success))
(:pointer
;; For unclear reasons, some GSL functions want callbacks
;; to return a void pointer which is apparently meaningless.
'((cffi:null-pointer)))))))
(defmacro defcbstruct
(functions &optional (structure 'gsl-function) additional-slots)
......
;; Example ODE
;; Liam Healy Sat Sep 29 2007 - 17:49
;; Time-stamp: <2008-01-28 22:21:07EST ode-example.lisp>
;; Time-stamp: <2009-01-18 18:29:53EST ode-example.lisp>
;; $Id$
;;; van der Pol as given in Section 25.5 of the GSL manual. To
;;; reproduce that example, (integrate-vanderpol 100.0d0)
;;; Notice how the variable mu, which is a "parameter" in GSL's view,
;;; to be passed on to the evaluation functions, is declared special
;;; so its value is shared dynamically. Thus GSL's parameters are
;;; always empty (and not accessible through user functions in GSLL).
(in-package :gsl)
(defun vanderpol (time y dydt)
(defun vanderpol (time y0 y1)
(declare (special mu) (ignorable time))
(with-c-doubles ((y y0 y1) (dydt dydt0 dydt1))
(setf dydt0 y1
dydt1 (- (- y0) (* mu y1 (- (* y0 y0) 1))))))
(values y1 (- (- y0) (* mu y1 (- (* y0 y0) 1)))))
(defun vanderpol-jacobian (time y dfdy dfdt)
(defun vanderpol-jacobian (time y0 y1)
(declare (special mu) (ignorable time))
(with-c-doubles ((y y0 y1) (dfdy dfdy0 dfdy1 dfdy2 dfdy3) (dfdt dfdt0 dfdt1))
(setf dfdt0 0.0d0
dfdt1 0.0d0
dfdy0 0.0d0
dfdy1 1.0d0
dfdy2 (- (* -2 mu y0 y1) 1.0d0)
dfdy3 (* -1 mu (- (* y0 y0) 1.0d0)))))
(values 0.0d0 0.0d0 ; dfdt
0.0d0 1.0d0 ; dfdy in row-major order: 00, 01, 10, 11
(- (* -2 mu y0 y1) 1.0d0)
(* -1 mu (- (* y0 y0) 1.0d0))))
(def-ode-functions vanderpol vanderpol-jacobian 2)
(defparameter *max-iter* 2000)
(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6))
(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6) (print-steps t))
"Integrate the van der Pol oscillator as given in Section 25.5 of the
GSL manual. To reproduce that example, (integrate-vanderpol 100.0d0)."
(let ((mu 10.0d0) (time 0.0d0) (iter 0))
(declare (special mu))
(with-ode-integration (time step-size (dependent dep0 dep1) 2)
(setf dep0 1.0d0 dep1 0.0d0)
(loop (when (or (>= (dcref time) max-time) (> iter *max-iter*)) (return))
(apply-evolution
evolve control stepper vanderpol
time max-time step-size dependent)
(incf iter)
(format t "~&~12,6f~10t~12,6f~24t~12,6f"
(dcref time) dep0 dep1)))))
(loop (when (or (>= (dcref time) max-time) (> iter *max-iter*))
(return (values (dcref time) dep0 dep1)))
(apply-evolution
evolve control stepper vanderpol
time max-time step-size dependent)
(incf iter)
(when print-steps
(format t "~12,6f~10t~12,6f~24t~12,6f~&"
(dcref time) dep0 dep1))))))
(save-test ode (integrate-vanderpol 1.0d0 1.d-4 nil))
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2008-08-21 22:04:31EDT ode-system.lisp>
;; Time-stamp: <2009-01-18 18:23:07EST ode-system.lisp>
;; $Id$
(in-package :gsl)
......@@ -16,19 +16,29 @@
(export '(def-ode-functions with-ode-integration))
(defmacro def-ode-functions (name jacobian dimension)
"Setup functions for ODE integrators.
"Setup functions for ODE integrators. The variable name is used as the name of the
The CL functions name and jacobian should be defined previously
with defuns."
;; The function should take three arguments: time, dependent, derivatives
;; The latter two will be C arrays. To reference them, use #'with-c-doubles.
;; To make this more transparent using a normal CL function
;; would require transferring numbers back and forth between C and CL arrays,
;; which could be inefficient.
with defuns.
The function and Jacobian arguments are the same: time and
dimension dependent variables as scalars.
The function returns dimension multiple values corresponding to
the derivatives of the dependent variables, f(t).
The Jacobian returns dimension^2 values corresponding to the
partial derivatives of each of the functions f with respect to each
of the dependent variables y, and dimension values corresponding to
the derivatives of f with respect to t."
;; set return values
;; Possible future improvements: take/set arrays easily, allow lambdas instead of named functions.
`(progn
(defmcallback ,name :success-failure (:double :pointer :pointer))
(defmcallback ,jacobian :success-failure (:double :pointer :pointer :pointer))
(defcbstruct (,name function ,jacobian jacobian) ode-system
((dimension ,dimension)))))
(defmcallback
,name :success-failure
(:double (:double ,dimension) (:set :double ,dimension)))
(defmcallback
,jacobian :success-failure
(:double (:double ,dimension) (:set :double ,(expt dimension 2))
(:set :double ,dimension)))
(defcbstruct (,name function ,jacobian jacobian) ode-system
((dimension ,dimension)))))
(defmacro with-ode-integration
((time step-size dependent dimensions &optional (stepper '*step-rk8pd*)
......
;; Regression test ODE for GSLL, automatically generated
(in-package :gsl)
(LISP-UNIT:DEFINE-TEST ODE
(LISP-UNIT::ASSERT-NUMERICAL-EQUAL
(LIST 1.0d0 -1.4568622636249005d0
-11.547385179410822d0)
(MULTIPLE-VALUE-LIST
(INTEGRATE-VANDERPOL 1.0d0 1.d-4 NIL))))
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment