Commit ecac342b authored by liam's avatar liam
Browse files

Improved definition of callbacks defined in CL. Arguments eliminated,

application to Monte Carlo, simulated annealing, ordinary differential
equations.  Simulated annealing still does not work.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3258 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 7b4aa0c7
......@@ -3,7 +3,7 @@
; description: Foreign callback functions.
; date: Sun Dec 9 2007 - 22:08
; author: Liam Healy
; modified: Sun Dec 30 2007 - 16:29
; modified: Sat Jan 5 2008 - 21:30
;********************************************************
;;; $Id: $
......@@ -30,10 +30,11 @@
;;; 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.
;;; #'with-c-double is provided to give named access to the elements.
(export
'(def-scalar-function undef-scalar-function defun-scalar with-c-vector))
'(def-scalar-function undef-scalar-function defun-scalar
with-c-double with-c-doubles))
;;;;****************************************************************************
;;;; Setting slots
......@@ -70,25 +71,35 @@
;;; Used by numerical-integration, numerical-differentiation, chebyshev, ntuple.
(defmacro def-scalar-function
(name argument &optional (return-type :double) (argument-type :double))
(name
&optional (return-type :double) (argument-type :double)
(structure 'gsl-function)
additional-slots)
"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)))
(let ((argument (gensym "CB")))
`(progn
(cffi:defcallback ,name ,return-type
((,argument ,argument-type) (params :pointer))
(declare (ignore params))
(,name ,argument))
,@(when
structure
;; Assume that defcallback does not bind the variable 'name.
`((defparameter ,name (cffi:foreign-alloc ',structure))
(set-slot-function ,name ',structure 'function ',name)
(set-parameters ,name ',structure)
,@(loop for slot in additional-slots
collect
`(set-structure-slot
,name ',structure ',(first slot) ,(second slot))))))))
(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)
(cffi:foreign-free (symbol-value name))
(makunbound name))
;;; Combine a defun and def-scalar-function in one:
......@@ -97,7 +108,7 @@
a double-float in CL and C."
`(progn
(defun ,name ,arglist ,@body)
(def-scalar-function ,name ,@arglist)))
(def-scalar-function ,name)))
;;;;****************************************************************************
;;;; Vector of doubles
......@@ -107,9 +118,9 @@
;;; (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.
;;; #'with-c-double is provided to make things easier.
(defmacro with-c-vector
(defmacro with-c-double
((c-vector &rest element-names) &body body)
"Provide named access to each element of a C array of doubles, for either
reading or setting."
......@@ -117,3 +128,13 @@
,(loop for i from 0 for a in element-names
collect `(,a (double-to-cl ,c-vector ,i)))
,@body))
(defmacro with-c-doubles ((&rest cvector-names) &body body)
"Provide named access to each element of a set of C arrays of doubles,
for either reading or setting."
(if (null (rest cvector-names))
`(with-c-double ,(first cvector-names)
,@body)
`(with-c-double ,(first cvector-names)
(with-c-doubles ,(rest cvector-names)
,@body))))
......@@ -3,7 +3,7 @@
; description: Monte Carlo Integration
; date: Sat Feb 3 2007 - 17:42
; author: Liam Healy
; modified: Sun Dec 30 2007 - 15:45
; modified: Sat Jan 5 2008 - 21:20
;********************************************************
;;; $Id: $
......@@ -297,17 +297,9 @@
(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)))
(defmacro def-mc-function (name dimensions)
`(def-scalar-function ,name :double :pointer monte-function
((dimensions ,dimensions))))
;;;;****************************************************************************
;;;; Examples and unit test
......@@ -317,11 +309,11 @@
;;; This is a function that occurs in random walk studies.
(defun monte-carlo-g (arg)
(with-c-vector (arg x y z)
(with-c-double (arg x y z)
(* (/ (expt pi 3))
(/ (- 1 (* (cos x) (cos y) (cos z)))))))
(def-mc-function monte-carlo-g arg 3)
(def-mc-function monte-carlo-g 3)
(defun random-walk-plain-example (&optional (nsamples 500000))
(with-monte-carlo-plain (ws 3)
......
......@@ -3,73 +3,47 @@
; description: Example ODE
; date: Sat Sep 29 2007 - 17:49
; author: Liam Healy
; modified: Thu Oct 18 2007 - 22:39
; modified: Sat Jan 5 2008 - 21:41
;********************************************************
;;; $Id: $
;;; Interfaces to GSL are still somewhat crude, as can be seen in the
;;; van der Pol example implemented below.
;;; Thu Oct 18 2007
;;; van der Pol as given in Section 25.5 of the GSL manual. To
;;; reproduce that example, (integrate-vanderpol 100.0d0)
(in-package :gsl)
(def-ode-function vanderpol-function (time y dydt)
(declare (special mu))
(setf (double-to-cl dydt 0) (double-to-cl y 1)
(double-to-cl dydt 1)
(- (- (double-to-cl y 0))
(* mu (double-to-cl y 1)
(- (* (double-to-cl y 0) (double-to-cl y 0)) 1)))))
(defun vanderpol (time y dydt)
(declare (special mu) (ignorable time))
(with-c-doubles ((y y0 y1) (dydt dydt0 dydt1))
(setf dydt0 y1
dydt1 (- (- y0) (* mu y1 (- (* y0 y0) 1))))))
(def-jacobian-function vanderpol-jacobian (time y dfdy dfdt)
(declare (special mu))
(setf (double-to-cl dfdt 0) 0.0d0
(double-to-cl dfdt 1) 0.0d0
(double-to-cl dfdy 0) 0.0d0
(double-to-cl dfdy 1) 1.0d0
(double-to-cl dfdy 2)
(- (* -2 mu (double-to-cl y 0) (double-to-cl y 1)) 1.0d0)
(double-to-cl dfdy 3)
(* -1 mu (- (* (double-to-cl y 0) (double-to-cl y 0)) 1.0d0))))
(defun vanderpol-jacobian (time y dfdy dfdt)
(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)))))
(def-ode-functions vanderpol vanderpol-jacobian 2)
(defparameter *max-iter* 2000)
(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6))
"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 ((stepper (step-allocate *step-rk8pd* 2))
(control (new-y-control 1.0d-6 0.0d0))
(evolve (allocate-evolution 2))
(mu 10.0d0)
(time 0.0d0)
(iter 0))
(let ((mu 10.0d0) (time 0.0d0) (iter 0))
(declare (special mu))
(cffi:with-foreign-objects
((vanderpol 'ode-system) (dependent :double 2)
(ctime :double) (cstep :double))
(setf
(cffi:foreign-slot-value vanderpol 'ode-system 'function)
(cffi:callback vanderpol-function)
(cffi:foreign-slot-value vanderpol 'ode-system 'jacobian)
(cffi:callback vanderpol-jacobian)
(cffi:foreign-slot-value vanderpol 'ode-system 'dimension)
2
(cffi:foreign-slot-value vanderpol 'ode-system 'parameters)
(cffi:null-pointer)
(double-to-cl dependent 0) 1.0d0
(double-to-cl dependent 1) 0.0d0
(double-to-cl ctime) time
(double-to-cl cstep) step-size)
(loop (when (or (>= time max-time) (> iter *max-iter*)) (return))
(with-ode-integration (time step-size (dependent dep0 dep1) 2)
(setf dep0 1.0d0 dep1 0.0d0)
(loop (when (or (>= (double-to-cl time) max-time) (> iter *max-iter*)) (return))
(apply-evolution
evolve control stepper vanderpol
ctime max-time cstep dependent)
(setf time (double-to-cl ctime))
time max-time step-size dependent)
(incf iter)
(format t "~&~12,6f~10t~12,6f~24t~12,6f"
time
(double-to-cl dependent 0)
(double-to-cl dependent 1))))
(free-evolution evolve)
(free-control control)
(step-free stepper)))
(double-to-cl time) dep0 dep1)))))
......@@ -3,7 +3,7 @@
; description: ODE system setup
; date: Sun Apr 15 2007 - 14:19
; author: Liam Healy
; modified: Sun Sep 30 2007 - 15:47
; modified: Sat Jan 5 2008 - 21:39
;********************************************************
;;; $Id: $
......@@ -17,54 +17,78 @@
(dimension :size)
(parameters :pointer))
(export '(def-ode-function def-jacobian-function))
(export '(def-ode-functions with-ode-integration))
(defmacro def-ode-function (name (time dependent derivatives) &body body)
"Define a function that will evaluate the right-hand sides (derivatives)
defining a set of ordinary differential equations (ODE).
The function should take as input the time (a double-float) and
dependent variables (a vector of double-floats) and fill the
derivatives vector with double-floats. It may refer to elements
of these vectors (arrays) using the macro double-to-cl.
This function may be passed to the GSL ODE integrators.
Parameters (non integration variables) may be passed by
using a lexical closure."
`(cffi:defcallback ,name :int
((,time :double)
(,dependent :pointer)
(,derivatives :pointer)
(params :pointer))
(declare (ignore params) (ignorable ,time))
,@body
;; Any errors or warnings should be signalled on the CL side;
;; if the function completes, we will always return success.
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS)))
(defmacro def-ode-functions (name jacobian dimension)
"Setup functions for ODE integrators.
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-vector.
;; 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.
(let ((time (make-symbol "TIME"))
(dependent (make-symbol "DEP"))
(derivatives (make-symbol "DERIV"))
(dfdy (make-symbol "DFDY"))
(dfdt (make-symbol "DFDT"))
(params (make-symbol "PARAMS")))
`(progn
(cffi:defcallback ,name :int
((,time :double)
(,dependent :pointer)
(,derivatives :pointer)
(,params :pointer))
(declare (ignore ,params))
(,name ,time ,dependent ,derivatives)
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS))
;; The function should take four arguments: time, dependent, dfdy, dfdt
;; The last three will be arrays.
(cffi:defcallback ,jacobian :int
((,time :double)
(,dependent :pointer)
(,dfdy :pointer)
(,dfdt :pointer)
(,params :pointer))
(declare (ignore ,params))
(,name ,time ,dependent ,dfdy ,dfdt)
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS))
;; Assume that defcallback does not bind the variable 'name.
(defparameter ,name (cffi:foreign-alloc 'ode-system))
(set-slot-function ,name 'ode-system 'function ',name)
(set-slot-function ,name 'ode-system 'jacobian ',jacobian)
(set-structure-slot ,name 'ode-system 'dimension ,dimension)
(set-parameters ,name 'ode-system))))
#|
(def-ode-function foo time y dydt
(setf (double-to-cl dydt 0) (- (double-to-cl y 1))
(double-to-cl dydt 1) (double-to-cl y 0)))
|#
(defmacro with-ode-integration
((time step-size dependent dimensions &optional (stepper '*step-rk8pd*)
(absolute-error 1.0d-6) (relative-error 0.0d0))
&body body)
"Environment for integration of ordinary differential equations.
The variables time and step-size will become C doubles in the body;
to convert back, use double-to-cl. The dependent variable may
be specified as a list being the same as the first argument to
with-c-double."
(let ((ctime (make-symbol "CTIME"))
(cstep (make-symbol "CSTEP")))
`(let ((stepper (step-allocate ,stepper ,dimensions))
(control (new-y-control ,absolute-error ,relative-error))
(evolve (allocate-evolution ,dimensions)))
(unwind-protect
(cffi:with-foreign-objects
((,(if (listp dependent) (first dependent) dependent)
:double ,dimensions) (,ctime :double) (,cstep :double))
(setf
(double-to-cl ,cstep) ,step-size
(double-to-cl ,ctime) ,time)
,(if (listp dependent)
`(with-c-double ,dependent
(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body))
`(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body)))
(free-evolution evolve)
(free-control control)
(step-free stepper)))))
(defmacro def-jacobian-function (name (time dependent dfdy dfdt) &body body)
"Define a function that will evaluate the Jacobian (partial derivative)
of the set of ordinary differential equations (ODE).
The function should take as input the time (a double-float) and
dependent variables (a vector of double-floats) and fill the
dfdy matrix and dfdt vector with double-floats.
It may refer to elements
of these vectors (arrays) using the macro double-to-cl.
This function may be passed to the GSL ODE integrators.
Parameters (non integration variables) may be passed by
using a lexical closure."
`(cffi:defcallback ,name :int
((,time :double)
(,dependent :pointer)
(,dfdy :pointer) ; This is a vector but should be CL array
(,dfdt :pointer)
(params :pointer))
(declare (ignore params) (ignorable ,time))
,@body
;; Any errors or warnings should be signalled on the CL side;
;; if the function completes, we will always return success.
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS)))
......@@ -3,22 +3,16 @@
; description: Simulated Annealing
; date: Sun Feb 11 2007 - 17:23
; author: Liam Healy
; modified: Sat Sep 15 2007 - 18:21
; modified: Sat Jan 5 2008 - 21:20
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; There are some problems. First, step-size is not successfully
;;; returned by GSL to the step function, it comes back as garbage.
;;; That is not a big problem because it's never used by GSL, it is
;;; only sent back, so we can just define and use it locally. The
;;; main problem is that it seems to run forever, even on the
;;; "trivial" example. This is apparently because n-tries is
;;; not getting to GSL either. It seems the problem is with
;;; simulated-annealing-parameters.
;;; Also, the code could use quite a bit of clean
;;; up to present a simpler interface.
;;; This does not work.
;;; Step size passed to the step function is incorrect.
;;; Print function is ignored, but probably couldn't work if it weren't.
;;; Does not converge.
(cffi:defcstruct simulated-annealing-parameters
......@@ -31,6 +25,36 @@
(mu-t :double)
(t-min :double))
(defmacro with-simulated-annealing-parameters
((name number-of-tries iterations-per-temperature
step-size &optional k t-initial mu-t t-min)
&body body)
`(cffi:with-foreign-object (,name 'simulated-annealing-parameters)
(setf
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 'n-tries)
,number-of-tries
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 'iterations-fixed-T)
,iterations-per-temperature
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 'step-size)
,step-size
;; The following parameters are for the Boltzmann distribution
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 'k)
,k
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 't-initial)
,t-initial
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 'mu-t)
,mu-t
(cffi:foreign-slot-value
,name 'simulated-annealing-parameters 't-min)
,t-min)
,@body))
(defun-gsl simulated-annealing
(generator x0-p
Ef take-step distance-function
......@@ -84,95 +108,83 @@
The simulated annealing routines require several user-specified
functions to define the configuration space and energy function.")
;; cribbed from def-gsl-function; unify?
(export 'def-sa-function)
(defmacro def-sa-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 simulated annealing functions."
(let ((argvec (gensym "MCARG")))
(defmacro def-energy-function (name)
"Define an energy or distance fuction for simulated annealing."
`(def-scalar-function ,name :double :pointer nil))
(defmacro def-step-function (name)
"Define a step fuction for simulated annealing."
(let ((generator (gensym "GEN"))
(arguments (gensym "ARGS"))
(step-size (gensym "SS")))
`(cffi:defcallback ,name :void
((,generator :pointer) (,arguments :pointer) (,step-size :double))
(,name ,generator ,arguments ,step-size))))
(defmacro def-distance-function (name)
"Define a metric distance fuction for simulated annealing."
(let ((x (gensym "X"))
(y (gensym "Y")))
`(cffi:defcallback ,name :double
(,(if (listp arg)
`(,argvec :pointer)
`(,arg :double)))
,@(if (listp arg)
`((symbol-macrolet
,(loop for i from 0 for a in arg
collect `(,a (cffi:mem-aref ,argvec :double ,i)))
,@body))
body))))
(defparameter *sa-function-calls* 0)
(defcallback e1 :double ((xp :pointer))
(incf *sa-function-calls*)
(when (> *sa-function-calls* 100000000)
(error "too much"))
(let ((x (cffi:mem-aref xp :double)))
((,x :pointer) (,y :pointer))
(,name ,x ,y))))
(defmacro def-print-function (name)
"Define a print function for simulated annealing."
`(def-scalar-function ,name :int :pointer nil))
;;;;****************************************************************************
;;;; Example
;;;;****************************************************************************
;;; Trivial example, Sec. 24.3.1
;;; This does not work.
(defun M2 (cx cy)
(with-c-double (cx x)
(with-c-double (cy y)
(abs (- x y)))))
(defun E2 (arg)
(with-c-double (arg x)
(incf *sa-function-calls*)
(when (> *sa-function-calls* 100)
(error "too much"))
(* (exp (- (expt (1- x) 2))) (sin (* 8 x)))))
(cffi:defcallback s1 :pointer
((generator :pointer) (parameters :pointer) (ss :double))
(declare (ignore ss))
(let ((step-size 10.0d0))
(let ((rand (uniform generator)))
(setf (cffi:mem-aref parameters :double)
(+ (cffi:mem-aref parameters :double)
(- (* 2 rand step-size) step-size))))
(cffi:null-pointer)))
(def-sa-function M1 (x y) (abs (- x y)))
(defcallback P1 :void ((xp :pointer))
(let ((x (cffi:mem-aref xp :double)))
(FORMAT T "~&from P1: ~a" X)))
(defun S2 (generator parameters step-size)
(with-c-double (parameters x)
;;(format t "~&~d ~d" x step-size)
(let ((step-size 10.0d0)) ; this is coming in wrong, so we fix it
(let ((rand (uniform generator)))
(setf x (+ x (- (* 2 rand step-size) step-size)))))))
(defparameter *sa-example-print* nil)
;;; Print functions are a problem because it is likely that the C
;;; stdout and the CL *standard-output* are not the same stream.
;;; Also, it seems to ignore that a function is supplied, and avoids
;;; printing anything.
(defun P2 (arg)
(with-c-double (arg x)
(when *sa-example-print*
(format T "~&from P2: ~a" x))))
(def-energy-function E2)
(def-distance-function M2)
(def-step-function S2)
(def-print-function P2)
(defun simulated-annealing-example ()
(cffi:with-foreign-object (initial :double)
(setf (cffi:mem-aref initial :double) 15.0d0)
(cffi:with-foreign-object (params 'simulated-annealing-parameters)
(setf
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'n-tries)
200
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'iterations-fixed-T)
10
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'step-size)
10.0d0
;; The following parameters are for the Boltzmann distribution
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'k)
1.0d0
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 't-initial)
0.002d0
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'mu-t)
1.005d0
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 't-min)
2.0d-6)
(simulated-annealing
*rng-mt19937* initial
'E1 'S1 'M1 'P1
(cffi:foreign-type-size :double)
params))))
#|
@deftypefun void
gsl_siman_solve
(const gsl_rng * @var{r},
void * @var{x0_p},
gsl_siman_Efunc_t @var{Ef},
gsl_siman_step_t @var{take_step},
gsl_siman_metric_t @var{distance},
gsl_siman_print_t @var{print_position},
gsl_siman_copy_t @var{copyfunc},
gsl_siman_copy_construct_t @var{copy_constructor},
gsl_siman_destroy_t @var{destructor},
size_t @var{element_size},
gsl_siman_params_t @var{params})
|#
(let ((*sa-function-calls* 0))
(rng-environment-setup)
(cffi:with-foreign-object (initial :double)
(setf (double-to-cl initial) 15.5d0)