Commit f15ffdda authored by Liam Healy's avatar Liam Healy

More callbacks take and return scalars

More usage of callbacks using or computing multiple quantities
now do so on scalars because the macros in GSLL pull out the
individual elements of foreign arrays, and accept multiple values and
put them into the right foreign arrays.  In roots-multi,
#'multiroot-slot takes a mobject and extracts the mpointer before
calling the foreign slot value.  Several examples in
solve-minimize-fit now have optional 'print-steps argument.  These
examples are now the basis for additional tests.
parent c83788ed
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2008-12-26 19:34:42EST monte-carlo.lisp>
;; Time-stamp: <2009-01-19 16:30:06EST monte-carlo.lisp>
;; $Id$
(in-package :gsl)
......@@ -219,8 +219,8 @@
(export 'def-mc-function)
(defmacro def-mc-function (name dimensions)
`(def-single-function ,name :double :pointer monte-function
((dimensions ,dimensions))))
`(def-single-function ,name :double :double monte-function
,dimensions nil nil))
;;;; Examples and unit test
......@@ -229,10 +229,9 @@
;;; Example from Sec. 23.5
;;; This is a function that occurs in random walk studies.
(defun monte-carlo-g (arg)
(with-c-double (arg x y z)
(* (/ (expt pi 3))
(/ (- 1 (* (cos x) (cos y) (cos z)))))))
(defun monte-carlo-g (x y z)
(* (/ (expt pi 3))
(/ (- 1 (* (cos x) (cos y) (cos z))))))
(def-mc-function monte-carlo-g 3)
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-18 19:02:18EST gsll-tests.asd>
;; Time-stamp: <2009-01-19 16:40:15EST gsll-tests.asd>
;; $Id$
(asdf:defsystem "gsll-tests"
......@@ -131,6 +131,8 @@
(:file "random-number-generators")
(:file "rayleigh")
(:file "rayleigh-tail")
(:file "roots-multi")
(:file "roots-one")
(:file "row")
(:file "scale")
(:file "set-basis")
;; Foreign callback functions.
;; Liam Healy
;; Time-stamp: <2009-01-18 17:04:46EST callback.lisp>
;; Time-stamp: <2009-01-19 16:44:13EST callback.lisp>
;; $Id$
(in-package :gsl)
......@@ -64,6 +64,16 @@
;;;; Macros for defining a callback to wrap a CL function
;;; Callback functions are defined using demcallback by passing the
;;; name of the function and the argument list of types. Arrays can
;;; be handled in one of two ways. If they are declared :pointer, the
;;; CL function will be passed a C pointer, and it is responsible for
;;; reading or setting the array, with #'dcref or #'maref. If they
;;; are declared (type size) then size scalars will be passed as
;;; arguments to the CL function, and if they are declared (:set type
;;; size), the CL function return size values, to which the array
;;; elements will be set.
;;; Usage example for scalar function (e.g. numerical-integration,
;;; numerical-differentiation, chebyshev, ntuple).
;;; (defmcallback myfn :double :double)
......@@ -92,18 +102,28 @@
;;; (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)
(defun embedded-clfunc-args (types callback-args &optional marray)
"The arguments passed to the CL function call embedded in the callback."
(loop for spec in types
for (symbol type) in callback-args
for (symbol nil) in callback-args
(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))
(if (third spec)
;; matrix, marrays only
(loop for i from 0 below (second spec)
(loop for j from 0 below (third spec)
`(maref ,symbol ,i ,j ',(cffi-cl (first spec)))))
;; vector, marray or C array
(loop for ind from 0 below (second spec)
collect (if marray
`(maref ,symbol ,ind nil ',(cffi-cl (first spec)))
`(cffi:mem-aref ,symbol ',(first spec) ,ind))))
(list symbol)))))
(defun callback-set-mvb (form types callback-args)
(defun callback-set-mvb (form types callback-args &optional marray)
"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
......@@ -114,9 +134,8 @@
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)))
(let* ((setvbls (embedded-clfunc-args settype setcba marray))
(count (apply '+ (mapcar (lambda (inds) (apply '* (rest inds))) settype)))
(mvbvbls (loop repeat count collect (gensym "SETCB"))))
(if (zerop count)
......@@ -131,7 +150,7 @@
(defmacro defmcallback
(name &optional (return-type :double) (argument-types :double)
additional-argument-types marray)
"Define a callback function used by GSL; the GSL function will call
it with an additional `parameters' argument that is ignored. the
argument-types is a single type or list of types of the argument(s)
......@@ -157,9 +176,11 @@
(embedded-clfunc-args atl cbargs) (embedded-clfunc-args aatl cbaddl)))
(embedded-clfunc-args atl cbargs marray)
(embedded-clfunc-args aatl cbaddl marray)))
(append atl aatl)
(append cbargs cbaddl))
(append cbargs cbaddl)
......@@ -208,17 +229,22 @@
&optional (return-type :double) (argument-type :double)
(structure 'gsl-function)
(dimensions-return dimensions)
(marray t))
"Define a callback and optionally a related C struct used by GSL.
This struct is bound to a CL special with the specified name.
This macro can be used whenever a callback is defined and
placed in a struct that has no other functions defined."
(defmcallback ,name ,return-type ,argument-type ,additional-argument-types)
`((defcbstruct ,name ,structure ,additional-slots)))))
(defmcallback ,name ,return-type
,(if dimensions `((,argument-type ,dimensions)) argument-type)
,(if dimensions-return `((:set ,argument-type ,dimensions-return)))
`((defcbstruct ,name ,structure
,(if dimensions `((dimensions ,dimensions))))))))
;;; Combine a defun and def-single-function in one:
(defmacro defun-single (name arglist &body body)
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2009-01-11 22:13:07EST mobject.lisp>
;; Time-stamp: <2009-01-19 13:50:33EST mobject.lisp>
;; $Id$
;;; GSL objects are represented in GSLL as and instance of a 'mobject.
......@@ -56,7 +56,10 @@
:index ,maker)
(defmethod initialize-instance :after
((object ,class) &key mpointer ,@cl-alloc-args)
((object ,class) &key mpointer ,@(union cl-alloc-args cl-initialize-args))
,@(let ((not-alloc (set-difference cl-initialize-args cl-alloc-args)))
(when not-alloc)
`((declare (ignore ,@not-alloc))))
(unless mpointer
(setf mpointer
(allocate object ,@(symbol-keyword-symbol cl-alloc-args))
;;; Multivariate roots.
;;; Liam Healy 2008-01-12 12:49:08
;;; Time-stamp: <2009-01-03 16:07:02EST roots-multi.lisp>
;;; Time-stamp: <2009-01-19 16:10:31EST roots-multi.lisp>
;;; $Id$
(in-package :gsl)
;;; /usr/include/gsl/gsl_multiroots.h
;;; Currently, functions defined for root solving will be passed
;;; scalars and should return scalars as multiple values. A possible
;;; future enhancement is to optionally pass marrays and return
;;; marrays instead. This would allow directly manipulation of
;;; marrays by the user function. Notes Mon Jan 19 2009.
;;;; Function definition
......@@ -21,9 +27,8 @@
(export 'def-mfunction)
(defmacro def-mfunction (name dimensions)
"Define a function for multivariate root solving."
`(def-single-function ,name :success-failure :pointer gsl-mfunction
((dimensions ,dimensions))
`(def-single-function ,name :success-failure :double gsl-mfunction
(cffi:defcstruct gsl-mfunction-fdf
;; See /usr/include/gsl/gsl_multiroots.h
......@@ -70,7 +75,7 @@
(:type type
(if ,set (dim0 initial) function-or-dimension))
(:function function-or-dimension :initial initial)))
(:function-derivative function-or-dimension :initial initial)))
(defmfun name ((solver multi-dimensional-root-solver-f))
......@@ -188,7 +193,7 @@
(state :pointer))
(defun multiroot-slot (solver slot)
(cffi:foreign-slot-value solver 'gsl-multiroot-fsolver slot))
(cffi:foreign-slot-value (mpointer solver) 'gsl-multiroot-fsolver slot))
(defmfun multiroot-test-delta (solver absolute-error relative-error)
......@@ -359,32 +364,28 @@
(defparameter *powell-A* 1.0d4)
(defun powell (argument return)
(defun powell (arg0 arg1)
"Powell's test function."
(setf (maref return 0)
(- (* *powell-A* (maref argument 0) (maref argument 1))
(maref return 1)
(+ (exp (- (maref argument 0))) (exp (- (maref argument 1)))
(- (1+ (/ *powell-A*))))))
(- (* *powell-A* arg0 arg1) 1)
(+ (exp (- arg0)) (exp (- arg1)) (- (1+ (/ *powell-A*))))))
;;; (def-mfunction powell 2)
(def-mfunction powell 2)
;;; This is the example given in Sec. 34.8.
(defparameter *rosenbrock-a* 1.0d0)
(defparameter *rosenbrock-b* 10.0d0)
(defun rosenbrock (argument return)
(defun rosenbrock (arg0 arg1)
"Rosenbrock test function."
(setf (maref return 0)
(* *rosenbrock-a* (- 1 (maref argument 0)))
(maref return 1)
(* *rosenbrock-b* (- (maref argument 1) (expt (maref argument 0) 2)))))
(* *rosenbrock-a* (- 1 arg0))
(* *rosenbrock-b* (- arg1 (expt arg0 2)))))
(def-mfunction rosenbrock 2)
(defun roots-multi-example ()
(defun roots-multi-example (&optional (print-steps t))
"Solving Rosenbrock, the example given in Sec. 34.8 of the GSL manual."
(let ((max-iter 1000))
(let* ((vect #m(-10.0d0 -5.0d0))
......@@ -399,28 +400,34 @@
(iterate solver)
(setf fnval (cl-array (function-value solver))
argval (cl-array (solution solver)))
(format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g"
(aref argval 0)
(aref argval 1)
(aref fnval 0)
(aref fnval 1))
(when print-steps
(format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
(aref argval 0)
(aref argval 1)
(aref fnval 0)
(aref fnval 1)))
finally (return
(values (aref argval 0)
(aref argval 1)
(aref fnval 0)
(aref fnval 1)))))))
(defun rosenbrock-df (argument jacobian)
(defun rosenbrock-df (arg0 arg1)
"The partial derivatives of the Rosenbrock functions."
(setf (maref jacobian 0 0) (- *rosenbrock-a*)
(maref jacobian 0 1) 0.0d0
(maref jacobian 1 0) (* -2 *rosenbrock-b* (maref argument 0))
(maref jacobian 1 1) *rosenbrock-b*))
(defun rosenbrock-fdf (argument value jacobian)
(rosenbrock argument value)
(rosenbrock-df argument jacobian))
(declare (ignore arg1))
(values (- *rosenbrock-a*)
(* -2 *rosenbrock-b* arg0)
;;; Why is it necessary to define a function that calls the two other functions?
(defun rosenbrock-fdf (arg0 arg1)
(multiple-value-bind (v0 v1)
(rosenbrock arg0 arg1)
(multiple-value-bind (j0 j1 j2 j3)
(rosenbrock-df arg0 arg1)
(values v0 v1 j0 j1 j2 j3))))
;;; Because def-solver-functions and def-single-function bind a symbol
;;; of the same name as the first function, and we want both to run,
......@@ -430,16 +437,17 @@
(def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2)
(defun roots-multi-example-df ()
(defun roots-multi-example-df (&optional (print-steps t))
"Solving Rosenbrock with derivatives, the example given in Sec. 34.8
of the GSL manual."
(flet ((print-state (iter argval fnval)
(format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g"
(maref argval 0)
(maref argval 1)
(maref fnval 0)
(maref fnval 1))))
(when print-steps
(format t "iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g~&"
(maref argval 0)
(maref argval 1)
(maref fnval 0)
(maref fnval 1)))))
(let ((max-iter 1000))
(let* ((vect #m(-10.0d0 -5.0d0))
(solver (make-multi-dimensional-root-solver-fdf
......@@ -460,3 +468,7 @@
(maref argval 1)
(maref fnval 0)
(maref fnval 1))))))))
(save-test roots-multi
(roots-multi-example nil)
(roots-multi-example-df nil))
;; One-dimensional root solver.
;; Liam Healy
;; Time-stamp: <2009-01-03 13:08:27EST roots-one.lisp>
;; Time-stamp: <2009-01-19 16:31:35EST roots-one.lisp>
;; $Id$
(in-package :gsl)
......@@ -23,24 +23,35 @@
(export 'def-solver-functions)
(defmacro def-solver-functions (function df fdf &optional dimensions)
(defmacro def-solver-functions (function df fdf &optional dimensions array)
"Setup functions for solvers.
The CL functions name and derivative should be defined previously
with defuns. If dimensions is non-nil, set multiroot solver
with defuns. If dimensions is non-nil (positive fixnum), set multiroot solver
functions. If dimensions is a number, the functions should expect
dimensions scalar (double-float) arguments and return "
(let ((struct (if dimensions 'gsl-mfunction-fdf 'gsl-function-fdf))
(argtype (if dimensions :pointer :double))
(rettype (if dimensions :success-failure :double)))
(argtype (if dimensions (if array :pointer `((:double ,dimensions))) :double))
(rettype (if dimensions :success-failure :double))
(vecrettype (if array '(:pointer) `((:set :double ,dimensions))))
(if array '(:pointer) `((:set :double ,dimensions ,dimensions)))))
(defmcallback ,function ,rettype ,argtype
,(when dimensions '(:pointer)))
(defmcallback ,df ,rettype ,argtype
,(when dimensions '(:pointer)))
,fdf ,(if dimensions :success-failure :pointer)
,argtype (:pointer :pointer))
(defcbstruct (,function function ,df df ,fdf fdf) ,struct
,(when dimensions `((dimensions ,dimensions)))))))
,function ,rettype ,argtype ,(when dimensions vecrettype) ,dimensions)
,df ,rettype ,argtype
,(when dimensions matrettype)
,(if dimensions :success-failure :void)
,(if dimensions
(append vecrettype matrettype)
'((:set :double 1) (:set :double 1)))
(defcbstruct (,function function ,df df ,fdf fdf) ,struct
,(when dimensions `((dimensions ,dimensions)))))))
;;;; Initialization
......@@ -313,10 +324,9 @@
(+ (* (+ (* a x) b) x) c))
(defun quadratic-derivative (x)
(+ (* 2 a x) b))
(defun quadratic-and-derivative (x cy cdy)
(with-c-doubles ((cy y) (cdy dy))
(setf y (+ (* (+ (* a x) b) x) c)
dy (+ (* 2 a x) b)))))
(defun quadratic-and-derivative (x)
(values (+ (* (+ (* a x) b) x) c)
(+ (* 2 a x) b))))
(def-single-function quadratic)
......@@ -348,19 +358,23 @@
quadratic-df quadratic-derivative quadratic-and-derivative)
(defun roots-one-fdf-example ()
(defun roots-one-fdf-example (&optional (print-steps t))
"Solving a quadratic, the example given in Sec. 32.10 of the GSL manual."
(let* ((max-iter 100)
(initial 5.0d0)
(solver (make-one-dimensional-root-solver-fdf
*newton-fdfsolver* quadratic-df initial)))
(format t "~&iter ~6t ~8troot ~22terr ~34terr(est)")
(when print-steps
(format t "iter ~6t ~8troot ~22terr ~34terr(est)~&"))
(loop for iter from 0
for oldroot = initial then root
for root = (progn (iterate solver) (solution solver))
while (and (< iter max-iter)
(not (root-test-delta root oldroot 0.0d0 1.0d-5)))
(format t "~&~d~6t~10,8g ~18t~10,6g~34t~10,6g"
iter root (- root (sqrt 5.0d0)) (- root oldroot)))))
(when print-steps
(format t "~d~6t~10,8g ~18t~10,6g~34t~10,6g~&"
iter root (- root (sqrt 5.0d0)) (- root oldroot)))
finally (return root))))
(save-test roots-one (roots-one-fdf-example nil))
;; Regression test ROOTS-MULTI for GSLL, automatically generated
(in-package :gsl)
(LIST 0.999999999999998d0 0.9999999999999964d0
1.9984014443252818d-15 4.440892098500626d-15)
(LIST 1.0d0 1.0d0 0.0d0 0.0d0)
;; Regression test ROOTS-ONE for GSLL, automatically generated
(in-package :gsl)
(LIST 2.236067977499978d0)
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