diff --git a/general/functions.lisp b/general/functions.lisp index 8c0c259ce1147f8de969d9f48a8c1b8171ee842d..973ca634fab9f6bf6526c53db8eae897729b62dd 100644 --- a/general/functions.lisp +++ b/general/functions.lisp @@ -1,6 +1,6 @@ ;; Foreign callback functions. ;; Liam Healy -;; Time-stamp: <2008-01-14 22:56:45 liam functions.lisp> +;; Time-stamp: <2008-01-20 18:12:42EST functions.lisp> ;; $Id: $ (in-package :gsl) @@ -61,42 +61,95 @@ (parameters :pointer)) ;;;;**************************************************************************** -;;;; A function of a scalar double +;;;; Macros for defining a callback and placing in a structure ;;;;**************************************************************************** -;;; Used by numerical-integration, numerical-differentiation, chebyshev, ntuple. +;;; 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) +;;; (defmcallback myfn :pointer :int (:pointer)) +;;; Usage example for function and derivative +;;; (defmcallback fdf :pointer :double (:pointer :pointer)) +;;; (defmcallback fdf :success-failure :int (:pointer :pointer)) + +(defmacro defmcallback + (name &optional (return-type :double) (argument-types :double) + additional-argument-types) + "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) + 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)))))))) + +(defmacro defcbstruct + (functions &optional (structure 'gsl-function) additional-slots) + "Define a callback-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." + (let ((name + ;; Bind a CL special under this name to the C structure. + (if (listp functions) (first functions) functions)) + (fnlist + ;; Make a list of (function slot-name ...) for each function. + (if (listp functions) functions (list `,functions 'function)))) + `(progn + ;; Create the C structure and bind CL variable to it. + (defparameter ,name (cffi:foreign-alloc ',structure)) + ;; Set all the function slots. + ,@(loop for (fn slot-name) on fnlist by #'cddr collect + `(set-slot-function ,name ',structure ',slot-name ',fn)) + ;; Set the parameters. + (set-parameters ,name ',structure) + ;; Set any additional slots. + ,@(loop for slot in additional-slots + collect + `(set-structure-slot + ,name ',structure ',(first slot) ,(second slot)))))) (defmacro def-scalar-function (name &optional (return-type :double) (argument-type :double) (structure 'gsl-function) additional-slots - additional-arguments) + additional-argument-types) "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." - (let ((argument (gensym "CB"))) - `(progn - (cffi:defcallback ,name - ,(if (eq return-type :success-failure) :int return-type) - ((,argument ,argument-type) (params :pointer) ,@additional-arguments) - (declare (ignore params)) - (,name ,argument ,@(mapcar #'first additional-arguments)) - ,@(when (eq return-type :success-failure) - ;; We always return success, because if there was a - ;; problem, a CL error would be signalled. - '(success))) - ,@(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)))))))) + `(progn + (defmcallback ,name ,return-type ,argument-type ,additional-argument-types) + ,@(when + structure + `((defcbstruct ,name ,structure ,additional-slots))))) (defun undef-scalar-function (name) "Free foreign callback function. It is not necessary to do this; think diff --git a/gsll.asd b/gsll.asd index 8d3b2f22f404fc4dca5b390f5a1d7be5b79fe795..a903d7a0c3d6a75f7222899b801fd9369f9b1061 100644 --- a/gsll.asd +++ b/gsll.asd @@ -1,18 +1,14 @@ -;******************************************************** -; file: gsll.asd -; description: Definition of GSLL system -; date: -; author: Liam Healy -; modified: Tue Jan 8 2008 - 22:52 -;******************************************************** -;;; $Id: $ +;; Definition of GSLL system +;; Liam Healy +;; Time-stamp: <2008-01-20 21:12:57EST gsll.asd> +;; $Id: $ (asdf:defsystem "gsll" :name "gsll" :description "GNU Scientific Library for Lisp." :version "0" :author "Liam M. Healy" - :licence "GPL v3" + :licence "GPL v3, FDL" :depends-on (cffi) :components ((:module init @@ -175,4 +171,5 @@ (:file "hankel" :depends-on (init data)) (:file "roots-one" :depends-on (init general)) (:file "minimization-one" :depends-on (init general)) - (:file "roots-multi" :depends-on (init general)))) + (:file "roots-multi" :depends-on (init general)) + (:file "minimization-multi" :depends-on (init general)))) diff --git a/index.html b/index.html index 7a92a8c8a1d3bc6df0541c0e383ee79fc975c3a4..c7de1b2deaee0862459f411976b59a2996b665fd 100644 --- a/index.html +++ b/index.html @@ -315,7 +315,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p> <tr> <td><a href="Multidimensional-Minimization.html"> Multidimensional Minimization</a></td> - <td>Not started</td> + <td>Done</td> </tr> <tr> <td><a href="Least_002dSquares-Fitting.html"> @@ -370,7 +370,7 @@ and expect a few failures on amd64 (four), none? on i386. Often, the <!-- Created: Feb 25 2005 --> <!-- hhmts start --> <small> - Time-stamp: <2008-01-19 14:11:28EST index.html> + Time-stamp: <2008-01-20 21:01:45EST index.html> </small> <!-- hhmts end --> </div> diff --git a/init/interface.lisp b/init/interface.lisp index 9442a7d41072cdd906dbb577de6f3523ce81a8da..d86a67f17909a3e852c3098554d68a75bbf794d8 100644 --- a/init/interface.lisp +++ b/init/interface.lisp @@ -1,6 +1,6 @@ ;; Macros to interface GSL functions. ;; Liam Healy -;; Time-stamp: <2008-01-19 18:58:42EST interface.lisp> +;; Time-stamp: <2008-01-19 19:54:09EST interface.lisp> ;; $Id: $ (in-package :gsl) @@ -112,8 +112,8 @@ ) (defparameter *c-to-cl-types* - '((:double . double-float) (:float . single-float) (:int . fixnum) - (:size . (fixnum 0)))) + `((:double . double-float) (:float . single-float) (:int . fixnum) + (:size . (integer 0 ,most-positive-fixnum)))) (defun cl-argument-types (cl-arguments c-arguments-types) "Create CL argument and types from the C arguments." diff --git a/minimization-multi.lisp b/minimization-multi.lisp new file mode 100644 index 0000000000000000000000000000000000000000..03eeed1de1efcf7081e8fdb54ad90acfd76d6d20 --- /dev/null +++ b/minimization-multi.lisp @@ -0,0 +1,407 @@ +;; Multivariate minimization. +;; Liam Healy <Tue Jan 8 2008 - 21:28> +;; Time-stamp: <2008-01-20 21:01:12EST minimization-multi.lisp> +;; $Id: $ + +(in-package :gsl) + +;; In the parabaloid example, I notice that the consruct +;; (min-test-gradient (mfdfminimizer-gradient minimizer) 1.0d-3) +;; is constructing a CL vector-double (in mfdfminimizer-gradient) and +;; then immediately pulling out the pointer (in min-test-gradient). It +;; is easy enough to eliminate this, but then mfdfminimizer-gradient +;; would not be useful to a CL user. + +;;;;**************************************************************************** +;;;; Function definition +;;;;**************************************************************************** + +;;; The structures gsl-mfunction and gsl-mfunction-fdf are, from the +;;; CFFI point of view, equally valid for gsl_multimin_function and +;;; gsl_multimin_function_fdf defined in +;;; /usr/include/gsl/gsl_multimin.h as they are for +;;; gsl_multiroot_function and gsl_multiroot_function_fdf defined in +;;; /usr/include/gsl/gsl_multiroots.h. As far as CFFI is concerned, a +;;; pointer is a pointer, even though the C definition the functions +;;; they point to have different signatures. + +(export 'def-minimization-functions) +(defmacro def-minimization-functions (function dimensions &optional df fdf) + "Setup functions for multivariate minimization. + The CL functions name and derivative should be defined previously + with defuns." + `(progn + (defmcallback ,function :double :pointer) + ,@(when df + `((defmcallback ,df :pointer :pointer :pointer) + (defmcallback ,fdf :pointer :pointer (:pointer :pointer)))) + ,(if df + `(defcbstruct (,function function ,df df ,fdf fdf) + gsl-mfunction-fdf + ((dimensions ,dimensions))) + `(defcbstruct (,function function) + gsl-mfunction + ((dimensions ,dimensions)))))) + + +;;;;**************************************************************************** +;;;; Initialization +;;;;**************************************************************************** + +(defun-gsl allocate-mfminimizer (type dimension) + "gsl_multimin_fminimizer_alloc" + ((type :pointer) (dimension :size)) + :c-return :pointer + :documentation ; FDL + "Allocate an instance of a minimizer of the given for an + function of the given dimensions.") + +(defun-gsl allocate-mfdfminimizer (type dimension) + "gsl_multimin_fdfminimizer_alloc" + ((type :pointer) (dimension :size)) + :c-return :pointer + :documentation ; FDL + "Allocate an instance of a derivative-based minimizer of the given for an + function of the given dimensions.") + +(defun-gsl set-mfminimizer (minimizer function initial step-size) + "gsl_multimin_fminimizer_set" + ((minimizer :pointer) (function :pointer) + ((pointer initial) :pointer) ((pointer step-size) :pointer)) + :documentation ; FDL + "Initialize the minimizer to minimize the function + starting from the initial point. + The size of the initial trial steps is given in vector + step-size. The precise meaning of this parameter depends on the + method used.") + +(defun-gsl set-mfdfminimizer + (minimizer function-derivative initial step-size tolerance) + "gsl_multimin_fdfminimizer_set" + ((minimizer :pointer) (function-derivative :pointer) + ((pointer initial) :pointer) (step-size :double) + (tolerance :double)) + :documentation ; FDL + "Initialize the minimizer to minimize the function + starting from the initial point. The size of the + first trial step is given by step-size. The accuracy of the line + minimization is specified by tolernace. The precise meaning of this + parameter depends on the method used. Typically the line minimization + is considered successful if the gradient of the function g is + orthogonal to the current search direction p to a relative + accuracy of tolerance, where dot(p,g) < tol |p| |g|.") + +(defun-gsl free-mfminimizer (minimizer) + "gsl_multimin_fminimizer_free" + ((minimizer :pointer)) + :c-return :void + :documentation ; FDL + "Free all the memory associated with the minimizer.") + +(defun-gsl free-mfdfminimizer (minimizer) + "gsl_multimin_fdfminimizer_free" + ((minimizer :pointer)) + :c-return :void + :documentation ; FDL + "Free all the memory associated with the minimizer.") + +(defun-gsl mfminimizer-name (minimizer) + "gsl_multimin_fminimizer_name" + ((minimizer :pointer)) + :c-return :string + :documentation + "The name of the minimizer.") + +(defun-gsl mfdfminimizer-name (minimizer) + "gsl_multimin_fdfminimizer_name" + ((minimizer :pointer)) + :c-return :string + :documentation + "The name of the minimizer.") + +(export '(with-mfdfminimizer with-mfminimizer)) +(defmacro with-mfminimizer + ((minimizer minimizer-type function initial step-size dimensions) + &body body) + "Create and initialize an fdfminimizer for multi-dimensional problems, + and clean up afterwards." + `(let ((,minimizer (allocate-mfminimizer ,minimizer-type ,dimensions))) + (unwind-protect + (progn + (set-mfminimizer ,minimizer ,function ,initial ,step-size) + ,@body) + (free-mfminimizer ,minimizer)))) + +(defmacro with-mfdfminimizer + ((minimizer minimizer-type fdf initial step-size tolerance dimensions) + &body body) + "Create and initialize an fdfminimizer for multi-dimensional problems, + and clean up afterwards." + `(let ((,minimizer (allocate-mfdfminimizer ,minimizer-type ,dimensions))) + (unwind-protect + (progn + (set-mfdfminimizer ,minimizer ,fdf ,initial ,step-size ,tolerance) + ,@body) + (free-mfdfminimizer ,minimizer)))) + +;;;;**************************************************************************** +;;;; Iteration +;;;;**************************************************************************** + +(defun-gsl iterate-mfminimizer (minimizer) + "gsl_multimin_fminimizer_iterate" + ((minimizer :pointer)) + :documentation ; FDL + "Perform a single iteration of the minimizer. If the iteration + encounters an unexpected problem then an error code will be + returned.") + +(defun-gsl iterate-mfdfminimizer (minimizer) + "gsl_multimin_fdfminimizer_iterate" + ((minimizer :pointer)) + :documentation ; FDL + "Perform a single iteration of the minimizer. If the iteration + encounters an unexpected problem then an error code will be + returned.") + +(defun-gsl mfminimizer-x (solver) + "gsl_multimin_fminimizer_x" + ((solver :pointer)) + :c-return (canswer :pointer) + :return ((make-data-from-pointer canswer)) + :documentation ; FDL + "The current best estimate of the location of the minimum.") + +(defun-gsl mfdfminimizer-x (solver) + "gsl_multimin_fdfminimizer_x" + ((solver :pointer)) + :c-return (canswer :pointer) + :return ((make-data-from-pointer canswer)) + :documentation ; FDL + "The current best estimate of the location of the minimum.") + +(defun-gsl mfminimizer-minimum (solver) + "gsl_multimin_fminimizer_minimum" + ((solver :pointer)) + :c-return :double + :documentation ; FDL + "The current best estimate of the value of the minimum.") + +(defun-gsl mfdfminimizer-minimum (solver) + "gsl_multimin_fdfminimizer_minimum" + ((solver :pointer)) + :c-return :double + :documentation ; FDL + "The current best estimate of the value of the minimum.") + +(defun-gsl mfminimizer-size (solver) + "gsl_multimin_fminimizer_size" + ((solver :pointer)) + :c-return :double + :documentation ; FDL + "A minimizer-specific characteristic size for the minimizer.") + +(defun-gsl mfdfminimizer-gradient (solver) + "gsl_multimin_fdfminimizer_gradient" + ((solver :pointer)) + :c-return (canswer :pointer) + :return ((make-data-from-pointer canswer)) + :documentation ; FDL + "The current best estimate of the gradient for the minimizer.") + +(defun-gsl mfdfminimizer-restart (solver) + "gsl_multimin_fdfminimizer_restart" + ((solver :pointer)) + :documentation ; FDL + "Reset the minimizer to use the current point as a + new starting point.") + +;;;;**************************************************************************** +;;;; Stopping criteria +;;;;**************************************************************************** + +(defun-gsl min-test-gradient (gradient absolute-error) + "gsl_multimin_test_gradient" + (((pointer gradient) :pointer) (absolute-error :double)) + :c-return :success-continue + :documentation ; FDL + "Test the norm of the gradient against the + absolute tolerance absolute-error. The gradient of a multidimensional + function goes to zero at a minimum. The test returns T + if |g| < epsabs is achieved, and NIL otherwise. A suitable choice of + absolute-error can be made from the desired accuracy in the function for + small variations in x. The relationship between these quantities + \delta f = g \delta x.") + +(defun-gsl min-test-size (size absolute-error) + "gsl_multimin_test_size" + ((size :double) (absolute-error :double)) + :c-return :success-continue + :documentation ; FDL + "Test the minimizer specific characteristic size (if applicable to + the used minimizer) against absolute tolerance @var{epsabs}. The test + returns T if the size is smaller than tolerance, and NIL otherwise.") + +;;;;**************************************************************************** +;;;; Algorithms +;;;;**************************************************************************** + +(defvariable *conjugate-fletcher-reeves* + "gsl_multimin_fdfminimizer_conjugate_fr" + ;; FDL + "The Fletcher-Reeves conjugate gradient algorithm. The conjugate + gradient algorithm proceeds as a succession of line minimizations. The + sequence of search directions is used to build up an approximation to the + curvature of the function in the neighborhood of the minimum. + + An initial search direction p is chosen using the gradient, and + line minimization is carried out in that direction. The accuracy of + the line minimization is specified by the parameter tol. The minimum + along this line occurs when the function gradient g and the search + direction p are orthogonal. The line minimization terminates when + dot(p,g) < tol |p| |g|. The search direction is updated using the + Fletcher-Reeves formula p' = g' - \beta g where \beta=-|g'|^2/|g|^2, + and the line minimization is then repeated for the new search + direction.") + +(defvariable *conjugate-polak-ribiere* + "gsl_multimin_fdfminimizer_conjugate_pr" + ;; FDL + "The Polak-Ribiere conjugate gradient algorithm. It is similar to + the Fletcher-Reeves method, differing only in the choice of the + coefficient \beta. Both methods work well when the evaluation point is + close enough to the minimum of the objective function that it is well + approximated by a quadratic hypersurface.") + +(defvariable *vector-bfgs* + "gsl_multimin_fdfminimizer_vector_bfgs" + ;; FDL + "The vector Broyden-Fletcher-Goldfarb-Shanno (BFGS) conjugate + gradient algorithm. It is a quasi-Newton method which builds up an + approximation to the second derivatives of the function using the + difference between successive gradient vectors. By combining the + first and second derivatives the algorithm is able to take Newton-type + steps towards the function minimum, assuming quadratic behavior in + that region.") + +(defvariable *simplex-nelder-mead* + "gsl_multimin_fminimizer_nmsimplex" + ;; FDL + "The Simplex algorithm of Nelder and Mead. It constructs + n vectors p_i from the + starting vector initial and the vector step-size as follows: + p_0 = (x_0, x_1, ... , x_n) + p_1 = (x_0 + step_size_0, x_1, ... , x_n) + p_2 = (x_0, x_1 + step_size_1, ... , x_n) + ... = ... + p_n = (x_0, x_1, ... , x_n+step_size_n) + These vectors form the n+1 vertices of a simplex in n + dimensions. On each iteration the algorithm tries to improve + the parameter vector p_i corresponding to the highest + function value by simple geometrical transformations. These + are reflection, reflection followed by expansion, contraction and multiple + contraction. Using these transformations the simplex moves through + the parameter space towards the minimum, where it contracts itself. + + After each iteration, the best vertex is returned. Note, that due to + the nature of the algorithm not every step improves the current + best parameter vector. Usually several iterations are required. + + The routine calculates the minimizer specific characteristic size as the + average distance from the geometrical center of the simplex to all its + vertices. This size can be used as a stopping criteria, as the simplex + contracts itself near the minimum. The size is returned by the function + #'mfminimizer-size.") + +;;;;**************************************************************************** +;;;; Examples +;;;;**************************************************************************** + +;;; Examples from Sec. 35.8. + +(defparameter *parabaloid-center* #(1.0d0 2.0d0)) + +(defun parabaloid (gsl-vector-pointer) + "A parabaloid function of two arguments, given in GSL manual Sec. 35.4." + (let ((x (vref gsl-vector-pointer 0)) + (y (vref gsl-vector-pointer 1)) + (dp0 (aref *parabaloid-center* 0)) + (dp1 (aref *parabaloid-center* 1))) + (+ (* 10 (expt (- x dp0) 2)) + (* 20 (expt (- y dp1) 2)) + 30))) + +(defun parabaloid-derivative + (arguments-gv-pointer derivative-gv-pointer) + (let ((x (vref arguments-gv-pointer 0)) + (y (vref arguments-gv-pointer 1)) + (dp0 (aref *parabaloid-center* 0)) + (dp1 (aref *parabaloid-center* 1))) + (setf (vref derivative-gv-pointer 0) + (* 20 (- x dp0)) + (vref derivative-gv-pointer 1) + (* 40 (- y dp1))))) + +(defun parabaloid-and-derivative + (arguments-gv-pointer fnval derivative-gv-pointer) + (with-c-double (function fnval) + (setf fnval (parabaloid arguments-gv-pointer)) + (parabaloid-derivative + arguments-gv-pointer derivative-gv-pointer)))) + +(def-minimization-functions + parabaloid 2 parabaloid-derivative parabaloid-and-derivative) + +(defun multimin-example-fletcher-reeves () + (with-data (initial vector-double 2) + (setf (data initial) #(5.0d0 7.0d0)) + (with-mfdfminimizer + (minimizer *conjugate-fletcher-reeves* parabaloid + initial 0.01d0 1.0d-4 2) + (loop with status = T + for iter from 0 below 100 + while status + do (iterate-mfdfminimizer minimizer) + (setf status + (not (min-test-gradient + (mfdfminimizer-gradient minimizer) + 1.0d-3))) + (let ((x (mfdfminimizer-x minimizer))) + (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f" + iter (gsl-aref x 0) (gsl-aref x 1) + (mfdfminimizer-minimum minimizer))) + finally (return + (let ((x (mfdfminimizer-x minimizer))) + (values (gsl-aref x 0) (gsl-aref x 1)))))))) + +;;; Because def-minimization-functions bind a symbol +;;; of the same name as the first function, and we want both to run, +;;; we'll make an alias function so we can use both. +(eval-when (:load-toplevel :execute) + (setf (fdefinition 'parabaloid-f) #'parabaloid)) + +(def-minimization-functions parabaloid-f 2) + +(defun multimin-example-nelder-mead () + (with-data (initial vector-double 2) + (with-data (step-size vector-double 2) + (setf (data initial) #(5.0d0 7.0d0)) + (set-all step-size 1.0d0) + (with-mfminimizer + (minimizer *simplex-nelder-mead* parabaloid-f initial step-size 2) + (loop with status = T and size + for iter from 0 below 100 + while status + do (iterate-mfminimizer minimizer) + (setf size + (mfminimizer-size minimizer) + status + (not (min-test-size size 1.0d-2))) + (let ((x (mfminimizer-x minimizer))) + (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f~40t~8,3f" + iter (gsl-aref x 0) (gsl-aref x 1) + (mfminimizer-minimum minimizer) + size)) + finally (return + (let ((x (mfminimizer-x minimizer))) + (values (gsl-aref x 0) (gsl-aref x 1))))))))) diff --git a/minimization-one.lisp b/minimization-one.lisp index c367ad8f241c092d9f857aee72c843f7d403bd25..7faaef55335925b10b78c9f2f829254fbbd1883e 100644 --- a/minimization-one.lisp +++ b/minimization-one.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: minimization-one.lisp -; description: Univariate minimization -; date: Tue Jan 8 2008 - 21:02 -; author: Liam Healy -; modified: Tue Jan 8 2008 - 22:50 -;******************************************************** -;;; $Id: $ +;; Univariate minimization +;; Liam Healy Tue Jan 8 2008 - 21:02 +;; Time-stamp: <2008-01-20 17:02:12EST minimization-one.lisp> +;; $Id: $ (in-package :gsl) @@ -85,7 +81,7 @@ bounding interval.") (defun-gsl fminimizer-x-minimum (minimizer) - "gsl_min_fminimizer_x_lower" + "gsl_min_fminimizer_x_minimum" ((minimizer :pointer)) :c-return :double :documentation @@ -137,7 +133,7 @@ "gsl_min_test_interval" ((lower :double) (upper :double) (absolute-error :double) (relative-error :double)) - :c-return :success-failure + :c-return :success-continue ; guess that this is s-c, not s-f :documentation "Test for the convergence of the interval [lower,upper] with absolute error and relative error specified. diff --git a/ordinary-differential-equations/ode-example.lisp b/ordinary-differential-equations/ode-example.lisp index 18419f68a00fb1210f2d3872522f21452b46dc0c..b291438b0b7fff9774c84dcb55166023bdf71ee0 100644 --- a/ordinary-differential-equations/ode-example.lisp +++ b/ordinary-differential-equations/ode-example.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: ode-example.lisp -; description: Example ODE -; date: Sat Sep 29 2007 - 17:49 -; author: Liam Healy -; modified: Sat Jan 5 2008 - 21:41 -;******************************************************** -;;; $Id: $ +;; Example ODE +;; Liam Healy Sat Sep 29 2007 - 17:49 +;; Time-stamp: <2008-01-20 17:40:48EST 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) diff --git a/ordinary-differential-equations/ode-system.lisp b/ordinary-differential-equations/ode-system.lisp index 5d11509d81af4be5e2bb74a1080226b9e89bd182..6cb4d02d35b5d8048054ac5061310314b95785fa 100644 --- a/ordinary-differential-equations/ode-system.lisp +++ b/ordinary-differential-equations/ode-system.lisp @@ -1,6 +1,6 @@ ;; ODE system setup ;; Liam Healy, Sun Apr 15 2007 - 14:19 -;; Time-stamp: <2008-01-14 22:58:19 liam ode-system.lisp> +;; Time-stamp: <2008-01-20 18:06:12EST ode-system.lisp> ;; $Id: $ (in-package :gsl) @@ -20,42 +20,15 @@ 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. + ;; 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. - (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) - 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) - 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)))) + `(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))))) (defmacro with-ode-integration ((time step-size dependent dimensions &optional (stepper '*step-rk8pd*) diff --git a/roots-multi.lisp b/roots-multi.lisp index d643587a05e0d2c57b99793f0c31ec5c340a4713..962235302095e52416522149c0499ecbd9dd1147 100644 --- a/roots-multi.lisp +++ b/roots-multi.lisp @@ -1,6 +1,6 @@ ;;; Multivariate roots. ;;; Liam Healy 2008-01-12 12:49:08 -;;; Time-stamp: <2008-01-16 19:31:03EST roots-multi.lisp> +;;; Time-stamp: <2008-01-20 17:18:39EST roots-multi.lisp> ;;; $Id: $ (in-package :gsl) @@ -27,7 +27,7 @@ "Define a function for multivariate root solving." `(def-scalar-function ,name :success-failure :pointer gsl-mfunction ((dimensions ,dimensions)) - ((returned-value gsl-vector-c)))) + (gsl-vector-c))) (cffi:defcstruct gsl-mfunction-fdf ;; See /usr/include/gsl/gsl_multiroots.h diff --git a/roots-one.lisp b/roots-one.lisp index 1958edbd9087deba00f92886c7fbcd36a4cffa56..ebdba67ae878b3def61358a966a9f8784c95a53a 100644 --- a/roots-one.lisp +++ b/roots-one.lisp @@ -1,6 +1,6 @@ ;; One-dimensional root solver. ;; Liam Healy -;; Time-stamp: <2008-01-15 22:57:11 liam roots-one.lisp> +;; Time-stamp: <2008-01-20 17:24:01EST roots-one.lisp> ;; $Id: $ (in-package :gsl) @@ -20,51 +20,25 @@ (parameters :pointer)) (export 'def-solver-functions) + (defmacro def-solver-functions (function df fdf &optional dimensions) "Setup functions for solvers. The CL functions name and derivative should be defined previously with defuns. If dimensions is non-nil, set multiroot solver functions." - (let ((arg (make-symbol "ARG")) - (params (make-symbol "PARAMS")) - (func (make-symbol "F")) - (deriv (make-symbol "DF")) - (struct (if dimensions 'gsl-mfunction-fdf 'gsl-function-fdf)) + (let ((struct (if dimensions 'gsl-mfunction-fdf 'gsl-function-fdf)) (argtype (if dimensions :pointer :double)) - (rettype (if dimensions :int :double))) + (rettype (if dimensions :success-failure :double))) `(progn - (cffi:defcallback ,function ,rettype - ((,arg ,argtype) - (,params :pointer) - ,@(when dimensions `((,func :pointer)))) - (declare (ignore ,params)) - (,function ,arg ,@(when dimensions `(,func))) - ,@(when dimensions '(success))) - (cffi:defcallback ,df ,rettype - ((,arg ,argtype) - (,params :pointer) - ,@(when dimensions `((,deriv :pointer)))) - (declare (ignore ,params)) - (,df ,arg ,@(when dimensions `(,deriv))) - ,@(when dimensions '(success))) - (cffi:defcallback ,fdf ,(if dimensions :int :pointer) - ((,arg ,argtype) - (,params :pointer) - (,func :pointer) - (,deriv :pointer)) - (declare (ignore ,params)) - (,fdf ,arg ,func ,deriv) - ;; fdf function returns sucess code in multivariate case, void - ;; pointer in univariate - ,(if dimensions 'success - '(cffi:null-pointer))) - (defparameter ,function (cffi:foreign-alloc ',struct)) - (set-slot-function ,function ',struct 'function ',function) - (set-slot-function ,function ',struct 'df ',df) - (set-slot-function ,function ',struct 'fdf ',fdf) - ,@(when dimensions - `((set-structure-slot ,function ',struct 'dimensions ,dimensions))) - (set-parameters ,function ',struct)))) + (defmcallback ,function ,rettype ,argtype + ,(when dimensions '(:pointer))) + (defmcallback ,df ,rettype ,argtype + ,(when dimensions '(:pointer))) + (defmcallback + ,fdf ,(if dimensions :success-failure :pointer) + ,argtype (:pointer :pointer)) + (defcbstruct (,function 'function ,df 'df ,fdf 'fdf) ,struct + ,(when dimensions `((dimensions ,dimensions))))))) ;;;;**************************************************************************** ;;;; Initialization