diff --git a/README b/README index a32017d5410b729cd17adc1f3b4b2348394687a7..fe33b4073c7af4662617d138b6914847cc7eb456 100644 --- a/README +++ b/README @@ -1,3 +1,9 @@ +TODO +remove cffi-array +rename defun-gsl, defvariable. +Document with-gsl-objects + + Features: - Designed so that usage of GSL is as if GSL were written in CL: diff --git a/chebyshev.lisp b/chebyshev.lisp index 0dcce83bffdd1bafbcf4b0cc7165602b59acd0ab..5aca6531e6175d8154610215da97d8dd9bbeab75 100644 --- a/chebyshev.lisp +++ b/chebyshev.lisp @@ -1,20 +1,22 @@ ;; Chebyshev Approximations ;; Liam Healy Sat Nov 17 2007 - 20:36 -;; Time-stamp: <2008-01-20 22:37:14EST chebyshev.lisp> +;; Time-stamp: <2008-01-26 18:17:32EST chebyshev.lisp> ;; $Id: $ (in-package :gsl) -;;; Needs with- macro to allocate, initialize, and free object. - ;;;;**************************************************************************** ;;;; Creation and calculation of Chebyshev series ;;;;**************************************************************************** +(set-asf chebyshev allocate-chebyshev free-chebyshev initialize-chebyshev) + (defun-gsl allocate-chebyshev (order) "gsl_cheb_alloc" ((order :size)) :c-return :pointer + :export nil + :index (with-gsl-objects chebyshev) :documentation "Allocate a Chebyshev series of specified order and return a pointer to it.") @@ -23,6 +25,8 @@ "gsl_cheb_free" ((chebyshev :pointer)) :c-return :void + :export nil + :index (with-gsl-objects chebyshev) :documentation "Free a previously allocated Chebyshev series.") @@ -30,6 +34,8 @@ "gsl_cheb_init" ((chebyshev :pointer) (function :pointer) (lower-limit :double) (upper-limit :double)) + :export nil + :index (with-gsl-objects chebyshev) :documentation "Compute the Chebyshev approximation for the function over the range (lower-limit, upper-limit) to the previously specified order. The @@ -118,34 +124,28 @@ (defun-single chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0)) (defun chebyshev-table-example () - (let ((cheb (allocate-chebyshev 40)) - (steps 100)) - (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))) + (let ((steps 100)) + (with-gsl-objects ((chebyshev cheb 40 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))))))) (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) + (with-gsl-objects + ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0) + (chebyshev deriv 40) + (chebyshev integ 40)) (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)))) + (list + (evaluate-chebyshev cheb x) + (evaluate-chebyshev deriv x) + (evaluate-chebyshev integ x)))) (lisp-unit:define-test chebyshev (lisp-unit:assert-equal diff --git a/data/data.lisp b/data/data.lisp index 5a7be85ef1013f2751d7e18f33a797fce7e816df..008c1322b8f93cf091e730d7c570124edebe9a1e 100644 --- a/data/data.lisp +++ b/data/data.lisp @@ -1,6 +1,6 @@ ;; Using GSL storage. ;; Liam Healy, Sun Mar 26 2006 - 16:32 -;; Time-stamp: <2008-01-14 21:09:15 liam data.lisp> +;; Time-stamp: <2008-01-27 18:50:53EST data.lisp> ;; $Id: $ (in-package :gsl) @@ -198,6 +198,8 @@ (export '(make-data with-data)) +(set-asf generic make-data free (setf data) 3) + (defgeneric alloc (object) (:documentation "Allocate GSL data; used internally.")) @@ -210,13 +212,35 @@ (:method :after ((object gsl-data)) (setf (pointer object) nil))) (defun make-data (type zero &rest size) + "Make the GSL data object, including the allocation of space. + The user is responsible for calling #'free to free the foreign + memory when done." + (let ((obj + (make-instance + (make-symbol-from-strings *gsl-prefix* type) + :storage-size size))) + (if zero (calloc obj) (alloc obj)) + obj)) + +(defmacro with-data ((symbol type size &optional zero) &body body) + "Allocate GSL data, bind to pointer, + and then deallocated it when done. If zero is T, zero the + contents when allocating." + `(let ((,symbol + (make-data ',type ,zero ,@(if (listp size) size (list size))))) + (unwind-protect + (progn ,@body) + (free ,symbol)))) + +#| +(defun make-data (type zero size) "Make the GSL data object, including the allocation of space. The user is responsible for calling #'free to free the foreign memory when done." (let ((obj (make-instance (make-symbol-from-strings *gsl-prefix* type) - :storage-size size))) + :storage-size (if (listp size) size (list size))))) (if zero (calloc obj) (alloc obj)) obj)) @@ -225,10 +249,11 @@ and then deallocated it when done. If zero is T, zero the contents when allocating." `(let ((,symbol - (make-data ',type ,zero ,@(if (listp size) size (list size))))) + (make-data ',type ,zero ,size))) (unwind-protect (progn ,@body) (free ,symbol)))) +|# ;;;;**************************************************************************** ;;;; Getting values into CL diff --git a/data/vector.lisp b/data/vector.lisp index c51597060f2e41ac2aef6c6dbc1988de274b8469..e33c223a3f10997503d3fe36e1ecc00333103dce 100644 --- a/data/vector.lisp +++ b/data/vector.lisp @@ -1,6 +1,6 @@ ;; Vectors ;; Liam Healy, Sun Mar 26 2006 - 11:51 -;; Time-stamp: <2008-01-15 19:01:37 liam vector.lisp> +;; Time-stamp: <2008-01-26 22:20:43EST vector.lisp> ;; $Id: $ (in-package :gsl) @@ -469,6 +469,14 @@ same length.") (setf (data vec) #(-3.21d0 1.0d0 12.8d0)) (data vec)) +(with-gsl-objects ((generic vec 'vector-double nil 3)) + (setf (data vec) #(-3.21d0 1.0d0 12.8d0)) + (data vec)) + + + + + (with-data (vec vector 3) (set-basis vec 1) (format t "~&~a ~a ~a" (gsl-aref vec 0) (gsl-aref vec 1) (gsl-aref vec 2))) diff --git a/eigensystems.lisp b/eigensystems.lisp index bb2c23b5c1332b1911aef6158d4451ae2f431098..10e303d5ef4d395bfda90b1754e534f5b492e731 100644 --- a/eigensystems.lisp +++ b/eigensystems.lisp @@ -1,68 +1,49 @@ ;; Eigenvectors and eigenvalues ;; Liam Healy, Sun May 21 2006 - 19:52 -;; Time-stamp: <2008-01-21 11:30:57EST eigensystems.lisp> +;; Time-stamp: <2008-01-27 18:53:17EST eigensystems.lisp> ;; $Id: $ (in-package :gsl) -;;; Eventually make a with-* macro for the eigen-herm functions. - -;;;;**************************************************************************** -;;;; Workspace macro -;;;;**************************************************************************** - -(export 'with-eigensystem-workspace) -(defmacro with-eigensystem-workspace - ((ws n &optional vectors hermitian) &body body) - "Create a workspace to use in finding eigenvalues and eigenvectors. - If eigenvectors are desired, set 'vectors to T." - `(let ((,ws ,(if vectors - (if hermitian - `(eigen-hermv-alloc ,n) - `(eigen-symmv-alloc ,n)) - (if hermitian - `(eigen-herm-alloc ,n) - `(eigen-symm-alloc ,n))))) - (check-null-pointer ,ws :ENOMEM "workspace allocation failed") - (unwind-protect - (progn ,@body) - ,(if vectors - (if hermitian - `(eigen-hermv-free ,n) - `(eigen-symmv-free ,ws)) - (if hermitian - `(eigen-herm-free ,n) - `(eigen-symm-free ,ws)))))) - - ;;;;**************************************************************************** ;;;; Real Symmetric Matrices ;;;;**************************************************************************** +(set-asf eigen-symm eigen-symm-alloc eigen-symm-free) + (defun-gsl eigen-symm-alloc (n) "gsl_eigen_symm_alloc" ((n :size)) + :c-return :pointer + :export nil + :index '(with-gsl-objects eigen-symm) :documentation ; FDL "Allocate a workspace for computing eigenvalues of n-by-n real symmetric matrices. The size of the workspace - is O(2n)." - :c-return :pointer) + is O(2n).") (defun-gsl eigen-symm-free (w) "gsl_eigen_symm_free" ((w :pointer)) :c-return :void + :export nil + :index '(with-gsl-objects eigen-symm) :documentation ; FDL "Free the memory associated with the workspace w.") +(set-asf eigen-symmv eigen-symmv-alloc eigen-symmv-free) + (defun-gsl eigen-symmv-alloc (n) "gsl_eigen_symmv_alloc" ((n :size)) + :c-return :pointer + :export nil + :index '(with-gsl-objects eigen-symmv) :documentation ; FDL "Allocate a workspace for computing eigenvalues and eigenvectors of n-by-n real symmetric matrices. The size of - the workspace is O(4n)." - :c-return :pointer) + the workspace is O(4n).") (defun-gsl eigen-symmv-free (w) "gsl_eigen_symmv_free" ((w :pointer)) + :index '(with-gsl-objects eigen-symmv) :c-return :void :documentation ; FDL "Free the memory associated with the workspace w.") @@ -81,7 +62,8 @@ (defun-gsl eigenvalues-eigenvectors-symmetric (A eval evec ws) "gsl_eigen_symmv" - ((A gsl-matrix-c) (eval gsl-vector-c) (evec gsl-matrix-c) (ws :pointer)) + (((pointer A) :pointer) ((pointer eval) :pointer) + ((pointer evec) :pointer) (ws :pointer)) :documentation ; FDL "The eigenvalues and eigenvectors of the real symmetric matrix A. Additional workspace of the appropriate size @@ -100,6 +82,8 @@ ;;;; Complex Hermitian Matrices ;;;;**************************************************************************** +(set-asf eigen-herm eigen-herm-alloc eigen-herm-free) + (defun-gsl eigen-herm-alloc (n) "gsl_eigen_herm_alloc" ((n :size)) :documentation ; FDL @@ -128,6 +112,8 @@ :invalidate (eval A) :return (eval)) +(set-asf eigen-hermv eigen-hermv-alloc eigen-hermv-free) + (defun-gsl eigen-hermv-alloc (n) "gsl_eigen_hermv_alloc" ((n :size)) :documentation ; FDL @@ -196,19 +182,23 @@ ;;;; Example ;;;;**************************************************************************** -#| -(with-data (mat matrix (3 3)) - (setf (data mat) - #2A((20.0d0 -10.0d0 0.0d0) - (-10.0d0 30.0d0 0.0d0) - (0.0d0 0.0d0 40.0d0))) - (with-data (evecs matrix (3 3)) - (with-data (evals vector 3) - (with-eigensystem-workspace (ws 3 t) - (eigenvalues-eigenvectors-symmetric mat evals evecs ws) - (values (data evals) (data evecs)))))) +(defun eigenvalue-eigenvectors-example () + (with-data (mat matrix-double (3 3)) + (setf (data mat) + #2A((20.0d0 -10.0d0 0.0d0) + (-10.0d0 30.0d0 0.0d0) + (0.0d0 0.0d0 40.0d0))) + (with-data (evecs matrix-double (3 3)) + (with-data (evals vector-double 3) + (with-gsl-objects ((eigen-symmv ws 3)) + (eigenvalues-eigenvectors-symmetric mat evals evecs ws) + (values (data evals) (data evecs))))))) + + + +#| #(13.819660112501051d0 36.180339887498945d0 40.0d0) #2A((0.8506508083520399d0 -0.5257311121191336d0 0.0d0) (0.5257311121191336d0 0.8506508083520399d0 0.0d0) diff --git a/hankel.lisp b/hankel.lisp index 4939ead1806a357b19c9b5ad81e6f863b498f6be..19dcf926d3397eb3bec9c14be8e5ecefa0d23e58 100644 --- a/hankel.lisp +++ b/hankel.lisp @@ -1,30 +1,35 @@ -;******************************************************** -; file: hankel.lisp -; description: Discrete Hankel Transforms. -; date: Sat Dec 8 2007 - 16:50 -; author: Liam Healy -; modified: Sat Dec 8 2007 - 18:38 -;******************************************************** -;;; $Id: $ +;; Discrete Hankel Transforms. +;; Liam Healy, Sat Dec 8 2007 - 16:50 +;; Time-stamp: <2008-01-26 18:36:00EST hankel.lisp> +;; $Id: $ (in-package :gsl) ;;; Everything compiles, but not tested -- need example. +(set-asf hankel allocate-hankel free-hankel init-hankel) + (defun-gsl allocate-hankel (size) "gsl_dht_alloc" ((size :size)) :c-return :pointer + :export nil + :index (with-gsl-objects hankel) :documentation "Allocate a Discrete Hankel transform object of size @var{size}.") (defun-gsl init-hankel (hankel nu xmax) "gsl_dht_new" ((hankel :pointer) (nu :double) (xmax :double)) + :export nil + :index (with-gsl-objects hankel) :documentation "Initialize the transform @var{t} for the given values of @var{nu} and @var{x}.") +;; This seems redundant; does it effectively combine allocate-hankel +;; and init-hankel? In this case we don't really need, and shouldn't +;; export it. (defun-gsl new-hankel (size nu xmax) "gsl_dht_new" ((size :size) (nu :double) (xmax :double)) @@ -38,6 +43,8 @@ "gsl_dht_free" ((hankel :pointer)) :c-return :void + :export nil + :index (with-gsl-objects hankel) :documentation "Free the Hankel object.") diff --git a/index.html b/index.html index cf5bba7e77010331291caac36e28cea91b681398..84f49a6bcf9ae1965b11116bf0e804514a195c9e 100644 --- a/index.html +++ b/index.html @@ -206,7 +206,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>.</p> </tr> <tr> <td><a href="Vectors-and-Matrices.html">Vectors and Matrices</a></td> - <td>Done for integer, real, double, complex</td> + <td>Done for integer, real, double, complex; subvectors don't work</td> </tr> <tr> <td><a href="Permutations.html">Permutations</a></td> @@ -374,7 +374,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-21 17:58:20EST index.html> + Time-stamp: <2008-01-26 22:03:36EST index.html> </small> <!-- hhmts end --> </div> diff --git a/init/utility.lisp b/init/utility.lisp index 94d1d208b95c726a0ce34dd5ba25c487fe57d5bf..256d2cce8c7c3205a287828431b353ebaf08f2ff 100644 --- a/init/utility.lisp +++ b/init/utility.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: utility.lisp -; description: Utility definitions -; date: Sun Dec 3 2006 - 10:21 -; author: Liam M. Healy -; modified: Sun Dec 3 2006 - 10:22 -;******************************************************** -;;; $Id: $ +;; Utility definitions +;; Liam Healy, Sun Dec 3 2006 - 10:21 +;; Time-stamp: <2008-01-27 13:28:29EST utility.lisp> +;; $Id: $ (in-package :gsl) @@ -14,3 +10,87 @@ (intern (format nil "~{~:@(~a~)~^-~}" (mapcar #'string strings)))) (defparameter *gsl-prefix* 'gsl) + +(defparameter *gsl-objects* nil) + +(defstruct gsl-object name allocate free set num-alloc-args) + +(defmacro set-asf (name allocate free &optional set (num-alloc-args 1)) + "Define an expansion of this name in #'with-gsl-objects + and use the appropriate functions for allocating, freeing, and setting." + `(eval-when (:compile-toplevel :load-toplevel :execute) + (pushnew + (make-gsl-object :name ',name :allocate ',allocate :free ',free + :set ',set :num-alloc-args ,num-alloc-args) + *gsl-objects* + :test #'equalp))) + +;; (set-asf poly-complex complex-workspace-alloc complex-workspace-free nil 1) +;; (set-asf mfminimizer allocate-mfminimizer free-mfminimizer set-mfminimizer 2) + +(defun gslo-letform (specs specfns) + "Create the let form for GSL objects." + (mapcar + (lambda (s sf) + `(,(second s) + (,(gsl-object-allocate sf) + ,@(subseq (cddr s) 0 (gsl-object-num-alloc-args sf))))) + specs specfns)) + +(defun gslo-setform (specs specfns) + "Create the set forms for GSL objects." + (loop for s in specs for sf in specfns + for addl-set-args + = (subseq (cddr s) (gsl-object-num-alloc-args sf)) + when (and (gsl-object-set sf) addl-set-args) + collect + `(,(gsl-object-set sf) + ;; presume the thing is the first arg to the set function, + ;; then add the other arguments + ,@(cons (second s) addl-set-args)))) + +(export 'with-gsl-objects) +(defmacro with-gsl-objects (specs &body body) + "Allocate, optionally set, and free any foreign object + used in GSL calculations." + (let ((specfns + (mapcar + (lambda (s) + (or (find (first s) *gsl-objects* :key #'gsl-object-name) + (error "Could not find GSL object named ~s" (first s)))) + specs))) + `(let (,@(gslo-letform specs specfns)) + (unwind-protect + (progn ; call set functions + ,@(gslo-setform specs specfns) + ,@body) + ,@(loop for s in specs for sf in specfns + collect `(,(gsl-object-free sf) ,(second s))))))) + +;;; roots-multi needs lambda list +;;; (solver solver-type function initial) +;;; allocate: solver-type (dim0 ,initial) +;;; set: solver function initial + + +#| +(vector-double vec init) +-> +(make-data vector-double nil (length init)) +(setf (data vec) init) + +;;; Example +(set-asf chebyshev allocate-chebyshev free-chebyshev initialize-chebyshev) + +(with-gsl-objects ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0)) + ...) + +(LET ((CHEB (ALLOCATE-CHEBYSHEV 40))) + (UNWIND-PROTECT + (PROGN + (INITIALIZE-CHEBYSHEV CHEB CHEBYSHEV-STEP 0.0d0 1.0d0) + ...) + (FREE-CHEBYSHEV CHEB))) + + +|# diff --git a/interpolation/interpolation.lisp b/interpolation/interpolation.lisp index 9147c0d94565f44ecb4614b38ab5a7ee0b6b55b6..4f392668e037e81ed1a9322e49ead64699f0bdc4 100644 --- a/interpolation/interpolation.lisp +++ b/interpolation/interpolation.lisp @@ -1,24 +1,24 @@ -;******************************************************** -; file: interpolation.lisp -; description: Interpolation allocation, initialization, -; and freeing. -; date: Sun Nov 4 2007 - 17:24 -; author: Liam Healy -; modified: Sun Nov 4 2007 - 19:00 -;******************************************************** -;;; $Id: $ +;; Interpolation allocation, initialization, and freeing. +;; Liam Healy, Sun Nov 4 2007 - 17:24 +;; Time-stamp: <2008-01-26 17:40:04EST interpolation.lisp> +;; $Id: $ (in-package :gsl) ;;; A spline is an interpolation that also stores the arrays xa and ya, ;;; so they need not be supplied on each call. +(set-asf + interpolation allocate-interpolation free-interpolation + initialize-interpolation 2) ;;; Interpolation (defun-gsl allocate-interpolation (type size) "gsl_interp_alloc" ((type :pointer) (size :size)) :c-return :pointer + :export nil + :index (with-gsl-objects interpolation) :documentation "Allocate an interpolation object of type for size data-points, and return the pointer to it.") @@ -38,9 +38,13 @@ "gsl_interp_free" ((interpolation :pointer)) :c-return :void + :export nil + :index (with-gsl-objects interpolation) :documentation "Frees the interpolation object @var{interp}.") +(set-asf spline allocate-spline free-spline intialize-spline 2) + ;;; Spline (defun-gsl allocate-spline (type size) "gsl_spline_alloc" diff --git a/interpolation/lookup.lisp b/interpolation/lookup.lisp index a3aadd3d1d6fc51629044ca58d7b2120e297e895..21993f8c7a1e110a2029264e89f89d0c1893c240 100644 --- a/interpolation/lookup.lisp +++ b/interpolation/lookup.lisp @@ -9,7 +9,7 @@ (in-package :gsl) -;;; Lacks with-acceleration macro. +(set-asf acceleration allocate-acceleration free-acceleration nil 0) (defun-gsl interpolation-search (x-array x low-index high-index) "gsl_interp_bsearch" @@ -24,6 +24,8 @@ "gsl_interp_accel_alloc" () :c-return :pointer + :export nil + :index (with-gsl-objects acceleration) :documentation "Allocate an accelerator object, which is a kind of iterator for interpolation lookups. It tracks the state of @@ -44,6 +46,8 @@ (defun-gsl free-acceleration (acceleration) "gsl_interp_accel_free" ((acceleration :pointer)) + :export nil + :index (with-gsl-objects acceleration) :c-return :void :documentation "Frees the accelerator object.") diff --git a/linear-least-squares.lisp b/linear-least-squares.lisp index 6312723d2f4c4a17db533d90a909b6917e154940..f0971047726a4db15b81defbb61918f42f8c9244 100644 --- a/linear-least-squares.lisp +++ b/linear-least-squares.lisp @@ -1,6 +1,6 @@ ;; Linear least squares, or linear regression ;; Liam Healy <2008-01-21 12:41:46EST linear-least-squares.lisp> -;; Time-stamp: <2008-01-21 17:55:19EST linear-least-squares.lisp> +;; Time-stamp: <2008-01-26 19:48:03EST linear-least-squares.lisp> ;; $Id: $ (in-package :gsl) @@ -125,16 +125,7 @@ ;;;; Multiparameter fitting ;;;;**************************************************************************** -(export 'with-fit-workspace) -(defmacro with-fit-workspace ((ws number-of-observations number-of-parameters) - &body body) - "Create and initialize an linear least squares workspace, - and clean up afterwards." - `(let ((,ws - (allocate-fit-workspace ,number-of-observations ,number-of-parameters))) - (unwind-protect - (progn ,@body) - (free-fit-workspace ,ws)))) +(set-asf fit-workspace allocate-fit-workspace free-fit-workspace nil 2) (defun-gsl allocate-fit-workspace (number-of-observations number-of-parameters) "gsl_multifit_linear_alloc" @@ -319,7 +310,7 @@ (gsl-aref X i 2) (expt (first row) 2) (gsl-aref y i) (second row) (gsl-aref w i) (/ (expt (third row) 2)))) - (with-fit-workspace (ws n 3) + (with-gsl-objects ((fit-workspace ws n 3)) (setf chisq (weighted-linear-mfit X w y c cov ws))) (format t "~&Best fit: Y = ~10,8f + ~10,8f X + ~10,8f X^2" diff --git a/minimization-multi.lisp b/minimization-multi.lisp index ccd9324e5d3b12516b06f1768fccf6352b517953..7878cd07326d19448ff914238a4a3ab4b5f102b6 100644 --- a/minimization-multi.lisp +++ b/minimization-multi.lisp @@ -1,6 +1,6 @@ ;; Multivariate minimization. ;; Liam Healy <Tue Jan 8 2008 - 21:28> -;; Time-stamp: <2008-01-20 23:17:16EST minimization-multi.lisp> +;; Time-stamp: <2008-01-27 22:55:36EST minimization-multi.lisp> ;; $Id: $ (in-package :gsl) @@ -48,10 +48,16 @@ ;;;; Initialization ;;;;**************************************************************************** +(set-asf mfminimizer allocate-mfminimizer free-mfminimizer set-mfminimizer 2) +(set-asf + mfdfminimizer allocate-mfdfminimizer free-mfdfminimizer set-mfdfminimizer 2) + (defun-gsl allocate-mfminimizer (type dimension) "gsl_multimin_fminimizer_alloc" ((type :pointer) (dimension :size)) :c-return :pointer + :export nil + :index (with-gsl-objects mfminimizer) :documentation ; FDL "Allocate an instance of a minimizer of the given for an function of the given dimensions.") @@ -60,6 +66,8 @@ "gsl_multimin_fdfminimizer_alloc" ((type :pointer) (dimension :size)) :c-return :pointer + :export nil + :index (with-gsl-objects mfdfminimizer) :documentation ; FDL "Allocate an instance of a derivative-based minimizer of the given for an function of the given dimensions.") @@ -68,6 +76,8 @@ "gsl_multimin_fminimizer_set" ((minimizer :pointer) (function :pointer) ((pointer initial) :pointer) ((pointer step-size) :pointer)) + :export nil + :index (with-gsl-objects mfminimizer) :documentation ; FDL "Initialize the minimizer to minimize the function starting from the initial point. @@ -81,6 +91,8 @@ ((minimizer :pointer) (function-derivative :pointer) ((pointer initial) :pointer) (step-size :double) (tolerance :double)) + :export nil + :index (with-gsl-objects mfdfminimizer) :documentation ; FDL "Initialize the minimizer to minimize the function starting from the initial point. The size of the @@ -95,6 +107,8 @@ "gsl_multimin_fminimizer_free" ((minimizer :pointer)) :c-return :void + :export nil + :index (with-gsl-objects mfminimizer) :documentation ; FDL "Free all the memory associated with the minimizer.") @@ -102,6 +116,8 @@ "gsl_multimin_fdfminimizer_free" ((minimizer :pointer)) :c-return :void + :export nil + :index (with-gsl-objects mfdfminimizer) :documentation ; FDL "Free all the memory associated with the minimizer.") @@ -119,31 +135,6 @@ :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 ;;;;**************************************************************************** @@ -164,54 +155,54 @@ encounters an unexpected problem then an error code will be returned.") -(defun-gsl mfminimizer-x (solver) +(defun-gsl mfminimizer-x (minimizer) "gsl_multimin_fminimizer_x" - ((solver :pointer)) + ((minimizer :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) +(defun-gsl mfdfminimizer-x (minimizer) "gsl_multimin_fdfminimizer_x" - ((solver :pointer)) + ((minimizer :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) +(defun-gsl mfminimizer-minimum (minimizer) "gsl_multimin_fminimizer_minimum" - ((solver :pointer)) + ((minimizer :pointer)) :c-return :double :documentation ; FDL "The current best estimate of the value of the minimum.") -(defun-gsl mfdfminimizer-minimum (solver) +(defun-gsl mfdfminimizer-minimum (minimizer) "gsl_multimin_fdfminimizer_minimum" - ((solver :pointer)) + ((minimizer :pointer)) :c-return :double :documentation ; FDL "The current best estimate of the value of the minimum.") -(defun-gsl mfminimizer-size (solver) +(defun-gsl mfminimizer-size (minimizer) "gsl_multimin_fminimizer_size" - ((solver :pointer)) + ((minimizer :pointer)) :c-return :double :documentation ; FDL "A minimizer-specific characteristic size for the minimizer.") -(defun-gsl mfdfminimizer-gradient (solver) +(defun-gsl mfdfminimizer-gradient (minimizer) "gsl_multimin_fdfminimizer_gradient" - ((solver :pointer)) + ((minimizer :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) +(defun-gsl mfdfminimizer-restart (minimizer) "gsl_multimin_fdfminimizer_restart" - ((solver :pointer)) + ((minimizer :pointer)) :documentation ; FDL "Reset the minimizer to use the current point as a new starting point.") @@ -344,9 +335,10 @@ (defun parabaloid-and-derivative (arguments-gv-pointer fnval derivative-gv-pointer) - (setf fnval (parabaloid arguments-gv-pointer)) - (parabaloid-derivative - arguments-gv-pointer derivative-gv-pointer)) + (with-c-double (fnval fv) + (setf fv (parabaloid arguments-gv-pointer)) + (parabaloid-derivative + arguments-gv-pointer derivative-gv-pointer))) (def-minimization-functions parabaloid 2 parabaloid-derivative parabaloid-and-derivative) @@ -354,13 +346,14 @@ (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) + (with-gsl-objects + ((mfdfminimizer minimizer *conjugate-fletcher-reeves* 2 parabaloid + initial 0.01d0 1.0d-4)) (loop with status = T for iter from 0 below 100 while status - do (iterate-mfdfminimizer minimizer) + do + (iterate-mfdfminimizer minimizer) (setf status (not (min-test-gradient (mfdfminimizer-gradient minimizer) @@ -386,8 +379,9 @@ (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) + (with-gsl-objects + ((mfminimizer + minimizer *simplex-nelder-mead* 2 parabaloid-f initial step-size)) (loop with status = T and size for iter from 0 below 100 while status diff --git a/minimization-one.lisp b/minimization-one.lisp index 83b7089fd34d1c79dae188f4dbda5d328bb537e2..76366bf3efe483da76aff683169f3cf60605ab81 100644 --- a/minimization-one.lisp +++ b/minimization-one.lisp @@ -1,6 +1,6 @@ ;; Univariate minimization ;; Liam Healy Tue Jan 8 2008 - 21:02 -;; Time-stamp: <2008-01-20 22:40:18EST minimization-one.lisp> +;; Time-stamp: <2008-01-26 19:33:45EST minimization-one.lisp> ;; $Id: $ (in-package :gsl) @@ -9,10 +9,14 @@ ;;;; Initialization ;;;;**************************************************************************** +(set-asf fminimizer allocate-fminimizer free-fminimizer set-fminimizer) + (defun-gsl allocate-fminimizer (type) "gsl_min_fminimizer_alloc" ((type :pointer)) :c-return :pointer + :export nil + :index (with-gsl-objects fminimizer) :documentation "Allocate an instance of a minimizer of the given type.") @@ -20,11 +24,14 @@ "gsl_min_fminimizer_set" ((minimizer :pointer) (function :pointer) (minimum :double) (lower :double) (upper :double)) + :export nil + :index (with-gsl-objects fminimizer) :documentation "Set, or reset, an existing minimizer to use the function and the initial search interval [lower, upper], with a guess for the location of the minimum.") +;;; Use this in with-gsl-objects macro in any way? (defun-gsl set-fminimizer-with-values (minimizer function x-minimum x-lower x-upper f-minimum f-lower f-upper) @@ -42,6 +49,8 @@ "gsl_min_fminimizer_free" ((minimizer :pointer)) :c-return :void + :export nil + :index (with-gsl-objects fminimizer) :documentation "Free all the memory associated with the minimizer.") @@ -52,18 +61,6 @@ :documentation "The name of the minimizer.") -(export '(with-fminimizer)) -(defmacro with-fminimizer - ((minimizer minimizer-type function minimum lower upper) &body body) - "Create and initialize an fminimizer for one-dimensional problems, - and clean up afterwards." - `(let ((,minimizer (allocate-fminimizer ,minimizer-type))) - (unwind-protect - (progn - (set-fminimizer ,minimizer ,function ,minimum ,lower ,upper) - ,@body) - (free-fminimizer ,minimizer)))) - ;;;;**************************************************************************** ;;;; Iteration ;;;;**************************************************************************** @@ -199,8 +196,9 @@ (defun minimization-one-example () "Solving a minimum, the example given in Sec. 33.8 of the GSL manual." (let ((max-iter 100)) - (with-fminimizer - (minimizer *brent-fminimizer* minimization-one-fn 2.0d0 0.0d0 6.0d0) + (with-gsl-objects + ((fminimizer + minimizer *brent-fminimizer* minimization-one-fn 2.0d0 0.0d0 6.0d0)) (format t "~&iter ~6t [lower ~24tupper] ~36tmin ~44tmin err ~54tupper-lower") (loop for iter from 0 for min = (fminimizer-x-minimum minimizer) @@ -215,3 +213,4 @@ iter lower upper min (- min pi) (- upper lower)))))) + diff --git a/monte-carlo.lisp b/monte-carlo.lisp index 8814367683df186b5f731e4d56eac9d2b8bb8547..c470824ebc742243d2d9534190945d2cdfcb44a3 100644 --- a/monte-carlo.lisp +++ b/monte-carlo.lisp @@ -1,6 +1,6 @@ ;; Monte Carlo Integration ;; Liam Healy Sat Feb 3 2007 - 17:42 -;; Time-stamp: <2008-01-20 22:39:10EST monte-carlo.lisp> +;; Time-stamp: <2008-01-26 17:20:19EST monte-carlo.lisp> ;; $Id: $ (in-package :gsl) @@ -13,12 +13,14 @@ (dim :size) (x :pointer)) +(set-asf monte-carlo-plain monte-carlo-plain-alloc monte-carlo-plain-free) + (defun-gsl monte-carlo-plain-alloc (dim) "gsl_monte_plain_alloc" ((dim :size)) :c-return :pointer :export nil - :index with-monte-carlo-plain + :index (with-gsl-objects monte-carlo-plain) :documentation "Allocate and initialize a workspace for Monte Carlo integration in @var{dim} dimensions.") @@ -36,18 +38,10 @@ ((integrator-state :pointer)) :c-return :void :export nil - :index with-monte-carlo-plain + :index (with-gsl-objects monte-carlo-plain) :documentation "Frees the memory associated with the integrator.") -(export 'with-monte-carlo-plain) -(defmacro with-monte-carlo-plain ((workspace size) &body body) - `(let ((,workspace (monte-carlo-plain-alloc ,size))) - (unwind-protect - (progn - ,@body) - (monte-carlo-plain-free ,workspace)))) - (defun-gsl monte-carlo-integrate-plain (function lower-limits upper-limits calls generator state) "gsl_monte_plain_integrate" @@ -102,12 +96,14 @@ (hits-l :pointer) (hits-r :pointer)) +(set-asf monte-carlo-miser monte-carlo-miser-alloc monte-carlo-miser-free) + (defun-gsl monte-carlo-miser-alloc (dim) "gsl_monte_miser_alloc" ((dim :size)) :c-return :pointer :export nil - :index with-monte-carlo-miser + :index (with-gsl-objects monte-carlo-miser) :documentation "Allocate and initialize a workspace for Monte Carlo integration in @var{dim} dimensions. The workspace is used to maintain @@ -126,21 +122,10 @@ ((integrator-state :pointer)) :c-return :void :export nil - :index with-monte-carlo-miser + :index (with-gsl-objects monte-carlo-miser) :documentation "Frees the memory associated with the integrator.") -(export 'with-monte-carlo-miser) -(defmacro with-monte-carlo-miser ((workspace size) &body body) - "Allocate and use the workspace for the MISER algorithm. - Parameters can be set by changing values in the structure - pointed to by the workspace variable." - `(let ((,workspace (monte-carlo-miser-alloc ,size))) - (unwind-protect - (progn - ,@body) - (monte-carlo-miser-free ,workspace)))) - (export 'miser-parameter) (defmacro miser-parameter (workspace parameter) "Get or set with setf the parameter value for the MISER Monte Carlo @@ -213,12 +198,14 @@ (calls-per-box :uint) (ostream :pointer)) +(set-asf monte-carlo-vegas monte-carlo-vegas-alloc monte-carlo-vegas-free) + (defun-gsl monte-carlo-vegas-alloc (dim) "gsl_monte_vegas_alloc" ((dim :size)) :c-return :pointer :export nil - :index with-monte-carlo-vegas + :index (with-gsl-objects monte-carlo-vegas) :documentation "Allocate and initialize a workspace for Monte Carlo integration in @var{dim} dimensions. The workspace is used to maintain @@ -237,21 +224,10 @@ ((integrator-state :pointer)) :c-return :void :export nil - :index with-monte-carlo-vegas + :index (with-gsl-objects monte-carlo-vegas) :documentation "Frees the memory associated with the integrator.") -(export 'with-monte-carlo-vegas) -(defmacro with-monte-carlo-vegas ((workspace size) &body body) - "Allocate and use the workspace for the VEGAS algorithm. - Parameters can be set by changing values in the structure - pointed to by the workspace variable." - `(let ((,workspace (monte-carlo-vegas-alloc ,size))) - (unwind-protect - (progn - ,@body) - (monte-carlo-vegas-free ,workspace)))) - (export 'vegas-parameter) (defmacro vegas-parameter (workspace parameter) "Get or set with setf the parameter value for the VEGAS Monte Carlo @@ -312,7 +288,7 @@ (def-mc-function monte-carlo-g 3) (defun random-walk-plain-example (&optional (nsamples 500000)) - (with-monte-carlo-plain (ws 3) + (with-gsl-objects ((monte-carlo-plain ws 3)) (with-data (lower vector-double 3) (with-data (upper vector-double 3) (setf (data lower) #(0.0d0 0.0d0 0.0d0) @@ -326,7 +302,7 @@ ws))))) (defun random-walk-miser-example (&optional (nsamples 500000)) - (with-monte-carlo-miser (ws 3) + (with-gsl-objects ((monte-carlo-miser ws 3)) (with-data (lower vector-double 3) (with-data (upper vector-double 3) (setf (data lower) #(0.0d0 0.0d0 0.0d0) @@ -340,7 +316,7 @@ ws))))) (defun random-walk-vegas-example (&optional (nsamples 500000)) - (with-monte-carlo-vegas (ws 3) + (with-gsl-objects ((monte-carlo-vegas ws 3)) (with-data (lower vector-double 3) (with-data (upper vector-double 3) (setf (data lower) #(0.0d0 0.0d0 0.0d0) diff --git a/numerical-integration.lisp b/numerical-integration.lisp index 958bccfe62fc88070fd51838432e4ed248a93e4c..a92034f437d3a7ab4ebab56e055967bec765f730 100644 --- a/numerical-integration.lisp +++ b/numerical-integration.lisp @@ -1,6 +1,6 @@ ;; Numerical integration ;; Liam Healy, Wed Jul 5 2006 - 23:14 -;; Time-stamp: <2008-01-21 12:34:03EST numerical-integration.lisp> +;; Time-stamp: <2008-01-26 17:11:25EST numerical-integration.lisp> ;; $Id: $ ;;; To do: QAWS, QAWO, QAWF, more tests @@ -36,25 +36,24 @@ ;;;; QAG adaptive Gauss-Kronrod integration ;;;;**************************************************************************** +(set-asf + integration-workspace integration-workspace-alloc integration-workspace-free) + (defun-gsl integration-workspace-alloc (size) "gsl_integration_workspace_alloc" ((size :size)) :c-return :pointer + :export nil + :index (with-gsl-objects integration-workspace) :documentation "Allocate a workspace sufficient to hold @var{n} double precision intervals, their integration results and error estimates.") (defun-gsl integration-workspace-free (pointer) "gsl_integration_workspace_free" ((pointer :pointer)) :c-return :void + :export nil + :index (with-gsl-objects integration-workspace) :documentation "Free the memory associated with the workspace @var{w}.") -(export 'with-integration-workspace) -(defmacro with-integration-workspace ((workspace size) &body body) - `(let ((,workspace (integration-workspace-alloc ,size))) - (unwind-protect - (progn - ,@body) - (integration-workspace-free ,workspace)))) - (cffi:defcenum integrate-method :gauss15 :gauss21 :gauss31 :gauss41 :gauss51 :gauss61) @@ -229,9 +228,9 @@ (integration-qng one-sine 0.0d0 pi)) (lisp-unit:assert-first-fp-equal "0.200000000000d+01" - (with-integration-workspace (ws 20) + (with-gsl-objects ((integration-workspace ws 20)) (integration-QAG one-sine 0.0d0 pi :gauss15 20 ws))) (lisp-unit:assert-error 'gsl-error - (with-integration-workspace (ws 20) + (with-gsl-objects ((integration-workspace ws 20)) (integration-QAG one-sine 0.0d0 pi :gauss15 50 ws)))) diff --git a/random/discrete.lisp b/random/discrete.lisp index 15e4548bb49e2ec0b7f21acdac808afaaa94b8ad..38780efc706f870814414c12ae30dd5f9bd0960a 100644 --- a/random/discrete.lisp +++ b/random/discrete.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: discrete.lisp -; description: Discrete random variables -; date: Sat Nov 11 2006 - 21:51 -; author: Liam M. Healy -; modified: Sat Nov 25 2006 - 16:46 -;******************************************************** -;;; $Id: $ +;; Discrete random variables +;; Liam Healy, Sat Nov 11 2006 - 21:51 +;; Time-stamp: <2008-01-26 17:22:14EST discrete.lisp> +;; $Id: $ (in-package :gsl) @@ -15,10 +11,14 @@ (A :pointer) (F :pointer)) +(set-asf discrete discrete-preprocess discrete-free) + (defun-gsl discrete-preprocess (probabilities) "gsl_ran_discrete_preproc" (((dim0 probabilities) :size) ((gsl-array probabilities) :pointer)) :c-return :pointer + :export nil + :index (with-gsl-objects discrete) :documentation "A pointer to a structure that contains the lookup table for the discrete random number generator. The array probabilities contains @@ -30,6 +30,8 @@ (defun-gsl discrete-free (table) "gsl_ran_discrete_free" ((table :pointer)) :c-return :void + :export nil + :index (with-gsl-objects discrete) :documentation "De-allocates the lookup table created by #'discrete-preprocess.") @@ -66,17 +68,17 @@ '(1 0 1 1 0 1 1 2 1 2 2) (with-data (probabilities vector-double 3) (setf (data probabilities) #(0.25d0 0.5d0 0.25d0)) - (with-discrete-table (probabilities table) + (with-gsl-objects ((discrete table probabilities)) (rng-set *rng-mt19937* 0) (loop for i from 0 to 10 - collect - (discrete *rng-mt19937* table)))) + collect + (discrete *rng-mt19937* table)))) (lisp-unit:assert-first-fp-equal "0.500000000000d+00" (with-data (probabilities vector-double 3) (setf (data probabilities) #(0.25d0 0.5d0 0.25d0)) - (with-discrete-table (probabilities table) + (with-gsl-objects ((discrete table probabilities)) (discrete-pdf 1 table)))))) diff --git a/roots-multi.lisp b/roots-multi.lisp index 782c30a5e604f727d1663904ba6292f511405a45..55ccde6bf76b69d3e32773813df7b0e7acf58fd0 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-20 22:40:18EST roots-multi.lisp> +;;; Time-stamp: <2008-01-27 23:04:02EST roots-multi.lisp> ;;; $Id: $ (in-package :gsl) @@ -43,10 +43,15 @@ ;;;; Initialization ;;;;**************************************************************************** +(set-asf mfsolver allocate-mfsolver free-mfsolver set-mfsolver 2) +(set-asf mfdfsolver allocate-mfdfsolver free-mfdfsolver set-mfdfsolver 2) + (defun-gsl allocate-mfsolver (type dimension) "gsl_multiroot_fsolver_alloc" ((type :pointer) (dimension :size)) :c-return :pointer + :export nil + :index (with-gsl-objects mfsolver) :documentation "Allocate an instance of a solver of the type specified for a system of the specified number of dimensions.") @@ -55,6 +60,8 @@ "gsl_multiroot_fdfsolver_alloc" ((type :pointer) (dimension :size)) :c-return :pointer + :export nil + :index (with-gsl-objects mfdfsolver) :documentation "Allocate an instance of a derivative solver of the type specified for a system of the specified number of dimensions.") @@ -62,6 +69,8 @@ (defun-gsl set-mfsolver (solver function initial) "gsl_multiroot_fsolver_set" ((solver :pointer) (function :pointer) ((pointer initial) :pointer)) + :export nil + :index (with-gsl-objects mfsolver) :documentation "Set or reset an existing solver to use the function and the initial guess gsl-vector.") @@ -70,6 +79,8 @@ "gsl_multiroot_fdfsolver_set" ((solver :pointer) (function-derivative :pointer) ((pointer initial) :pointer)) + :export nil + :index (with-gsl-objects mfdfsolver) :documentation "Set or reset an existing solver to use the function and derivative (fdf) and the initial guess.") @@ -78,6 +89,8 @@ "gsl_multiroot_fsolver_free" ((solver :pointer)) :c-return :void + :export nil + :index (with-gsl-objects mfsolver) :documentation "Free all the memory associated with the solver.") @@ -85,6 +98,8 @@ "gsl_multiroot_fdfsolver_free" ((solver :pointer)) :c-return :void + :export nil + :index (with-gsl-objects mfdfsolver) :documentation "Free all the memory associated with the solver.") @@ -102,28 +117,6 @@ :documentation "The name of the solver.") -(export '(with-mfsolver with-mfdfsolver)) -(defmacro with-mfsolver ((solver solver-type function initial) &body body) - "Create and initialize an fsolver for multi-dimensional problems, - and clean up afterwards." - `(let ((,solver (allocate-mfsolver ,solver-type (dim0 ,initial)))) - (unwind-protect - (progn - (set-mfsolver ,solver ,function ,initial) - ,@body) - (free-mfsolver ,solver)))) - -(defmacro with-mfdfsolver - ((solver solver-type f-df-fdf root-guess) &body body) - "Create and initialize an fdfsolver for one-dimensional problems, - and clean up afterwards." - `(let ((,solver (allocate-mfdfsolver ,solver-type (dim0 ,root-guess)))) - (unwind-protect - (progn - (set-mfdfsolver ,solver ,f-df-fdf ,root-guess) - ,@body) - (free-mfdfsolver ,solver)))) - ;;;;**************************************************************************** ;;;; Iteration ;;;;**************************************************************************** @@ -422,7 +415,8 @@ (let ((max-iter 1000)) (with-data (vect vector-double 2) (setf (data vect) #(-10.0d0 -5.0d0)) - (with-mfsolver (solver *hybrid-scaled* rosenbrock vect) + (with-gsl-objects + ((mfsolver solver *hybrid-scaled* (dim0 vect) rosenbrock vect)) (let ((fnval (mfsolver-f solver)) (argval (mfsolver-root solver))) (loop for iter from 0 @@ -442,6 +436,7 @@ (gsl-aref fnval 0) (gsl-aref fnval 1))))))))) + (defun rosenbrock-df (argument jacobian) "The partial derivatives of the Rosenbrock functions." (setf (mref jacobian 0 0) (- *rosenbrock-a*) @@ -474,7 +469,8 @@ (let ((max-iter 1000)) (with-data (vect vector-double 2) (setf (data vect) #(-10.0d0 -5.0d0)) - (with-mfdfsolver (solver *gnewton-mfdfsolver* rosenbrock-f vect) + (with-gsl-objects + ((mfdfsolver solver *gnewton-mfdfsolver* (dim0 vect) rosenbrock-f vect)) (let ((fnval (mfdfsolver-f solver)) (argval (mfdfsolver-root solver))) (loop for iter from 0 diff --git a/roots-one.lisp b/roots-one.lisp index 175bb1af6f94bf8e4ececda79985dbda9f1e2233..4abdf3d672e998790440168a018e08fd9aed46c5 100644 --- a/roots-one.lisp +++ b/roots-one.lisp @@ -1,6 +1,6 @@ ;; One-dimensional root solver. ;; Liam Healy -;; Time-stamp: <2008-01-20 23:14:27EST roots-one.lisp> +;; Time-stamp: <2008-01-26 19:07:22EST roots-one.lisp> ;; $Id: $ (in-package :gsl) @@ -44,10 +44,15 @@ ;;;; Initialization ;;;;**************************************************************************** +(set-asf fsolver allocate-fsolver free-fsolver set-fsolver) +(set-asf fdfsolver allocate-fdfsolver free-fdfsolver set-fdfsolver) + (defun-gsl allocate-fsolver (type) "gsl_root_fsolver_alloc" ((type :pointer)) :c-return :pointer + :export nil + :index (with-gsl-objects fsolver) :documentation "Allocate an instance of a solver of the given type.") @@ -55,6 +60,8 @@ "gsl_root_fdfsolver_alloc" ((type :pointer)) :c-return :pointer + :export nil + :index (with-gsl-objects fdfsolver) :documentation "Allocate an instance of a derivative-based solver of the given type.") @@ -62,6 +69,8 @@ "gsl_root_fsolver_set" ((solver :pointer) (function :pointer) (lower :double) (upper :double)) + :export nil + :index (with-gsl-objects fsolver) :documentation "Initialize or reinitialize an existing solver to use the function and the initial search interval @@ -70,6 +79,8 @@ (defun-gsl set-fdfsolver (solver function-derivative root-guess) "gsl_root_fdfsolver_set" ((solver :pointer) (function-derivative :pointer) (root-guess :double)) + :export nil + :index (with-gsl-objects fdfsolver) :documentation "Initialize or reinitialize an existing solver.") @@ -77,6 +88,8 @@ "gsl_root_fsolver_free" ((solver :pointer)) :c-return :void + :export nil + :index (with-gsl-objects fsolver) :documentation "Free all the memory associated with the solver.") @@ -84,6 +97,8 @@ "gsl_root_fdfsolver_free" ((solver :pointer)) :c-return :void + :export nil + :index (with-gsl-objects fdfsolver) :documentation "Free all the memory associated with the solver.") @@ -101,29 +116,6 @@ :documentation "The name of the solver.") -(export '(with-fsolver with-fdfsolver)) -(defmacro with-fsolver - ((solver solver-type function lower upper) &body body) - "Create and initialize an fsolver for one-dimensional problems, - and clean up afterwards." - `(let ((,solver (allocate-fsolver ,solver-type))) - (unwind-protect - (progn - (set-fsolver ,solver ,function ,lower ,upper) - ,@body) - (free-fsolver ,solver)))) - -(defmacro with-fdfsolver ((solver solver-type f-df-fdf root-guess) - &body body) - "Create and initialize an fdfsolver for one-dimensional problems, - and clean up afterwards." - `(let ((,solver (allocate-fdfsolver ,solver-type))) - (unwind-protect - (progn - (set-fdfsolver ,solver ,f-df-fdf ,root-guess) - ,@body) - (free-fdfsolver ,solver)))) - ;;;;**************************************************************************** ;;;; Iteration ;;;;**************************************************************************** @@ -362,7 +354,7 @@ (defun roots-one-example () "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." (let ((max-iter 50)) - (with-fsolver (solver *brent-fsolver* quadratic 0.0d0 5.0d0) + (with-gsl-objects ((fsolver solver *brent-fsolver* quadratic 0.0d0 5.0d0)) (format t "~&iter ~6t [lower ~24tupper] ~36troot ~44terr ~54terr(est)") (loop for iter from 0 for root = (fsolver-root solver) @@ -390,8 +382,7 @@ "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." (let ((max-iter 100) (initial 5.0d0)) - (with-fdfsolver - (solver *newton-fdfsolver* quadratic-df initial) + (with-gsl-objects ((fdfsolver solver *newton-fdfsolver* quadratic-df initial)) (format t "~&iter ~6t ~8troot ~22terr ~34terr(est)") (loop for iter from 0 for itres = (iterate-fdfsolver solver) diff --git a/series-acceleration.lisp b/series-acceleration.lisp index 1074ba4b164a566a0735b715f2d723c47680d99c..11b7c1dc458e380a21fc3f25e9af41b8f58e42c8 100644 --- a/series-acceleration.lisp +++ b/series-acceleration.lisp @@ -1,16 +1,10 @@ -;******************************************************** -; file: series-acceleration.lisp -; description: Series acceleration. -; date: Wed Nov 21 2007 - 18:41 -; author: Liam Healy -; modified: Wed Nov 21 2007 - 21:21 -;******************************************************** -;;; $Id: $ +;; Series acceleration. +;; Liam Healy, Wed Nov 21 2007 - 18:41 +;; Time-stamp: <2008-01-26 18:31:15EST series-acceleration.lisp> +;; $Id: $ (in-package :gsl) -;;; Needs with- macro to allocate, initialize, and free object. - ;;;;**************************************************************************** ;;;; Creation and calculation of Levin series acceleration ;;;;**************************************************************************** @@ -27,10 +21,14 @@ (dq-den :pointer) (dsum :pointer)) +(set-asf levin allocate-levin free-levin) + (defun-gsl allocate-levin (order) "gsl_sum_levin_u_alloc" ((order :size)) :c-return :pointer + :export nil + :index (with-gsl-objects levin) :documentation "Allocate a workspace for a Levin @math{u}-transform of @var{n} terms. The size of the workspace is @math{O(2n^2 + 3n)}.") @@ -38,6 +36,9 @@ (defun-gsl free-levin (levin) "gsl_sum_levin_u_free" ((levin :pointer)) + :c-return :void ; Error in GSL manual, should be void? + :export nil + :index (with-gsl-objects levin) :documentation "Free a previously allocated Levin acceleration.") @@ -60,10 +61,14 @@ ;;;; Acceleration with error estimation from truncation ;;;;**************************************************************************** +(set-asf levin-truncated allocate-levin-truncated free-levin-truncated) + (defun-gsl allocate-levin-truncated (order) "gsl_sum_levin_utrunc_alloc" ((order :size)) :c-return :pointer + :export nil + :index (with-gsl-objects levin-truncated) :documentation "Allocate a workspace for a Levin @math{u}-transform of @var{n} terms, without error estimation. The size of the workspace is @@ -72,6 +77,8 @@ (defun-gsl free-levin-truncated (levin) "gsl_sum_levin_utrunc_free" ((levin :pointer)) + :export nil + :index (with-gsl-objects levin-truncated) :documentation "Free a previously allocated Levin acceleration.") @@ -99,24 +106,30 @@ ;;; From Sec. 29.3 in the GSL manual. (defun acceleration-example () - (let* ((maxterms 20) - (levin (allocate-levin maxterms)) - (sum 0.0d0) - (zeta2 (/ (expt pi 2) 6))) - (with-data (array vector-double maxterms) - (dotimes (n maxterms) - (setf (gsl-aref array n) (coerce (/ (expt (1+ n) 2)) 'double-float)) - (incf sum (gsl-aref array n))) - (multiple-value-bind (accelerated-sum error) - (accelerate array levin) - (format t "~&term-by-term sum =~f using ~d terms" sum maxterms) - (format t "~&term-by-term sum =~f using ~d terms" - (cffi:foreign-slot-value levin 'levin 'sum-plain) - (cffi:foreign-slot-value levin 'levin 'terms-used)) - (format t "~&exact value = ~f" zeta2) - (format t "~&accelerated sum = ~f using ~d terms" - accelerated-sum - (cffi:foreign-slot-value levin 'levin 'terms-used)) - (format t "~&estimated error = ~f" error) - (format t "~&actual error = ~f" (- accelerated-sum zeta2)))) - (free-levin levin))) + (let ((maxterms 20) + (sum 0.0d0) + (zeta2 (/ (expt pi 2) 6))) + (with-gsl-objects ((levin levin maxterms)) + (with-data (array vector-double maxterms) + (dotimes (n maxterms) + (setf (gsl-aref array n) (coerce (/ (expt (1+ n) 2)) 'double-float)) + (incf sum (gsl-aref array n))) + (multiple-value-bind (accelerated-sum error) + (accelerate array levin) + (format t "~&term-by-term sum =~f using ~d terms" sum maxterms) + (format t "~&term-by-term sum =~f using ~d terms" + (cffi:foreign-slot-value levin 'levin 'sum-plain) + (cffi:foreign-slot-value levin 'levin 'terms-used)) + (format t "~&exact value = ~f" zeta2) + (format t "~&accelerated sum = ~f using ~d terms" + accelerated-sum + (cffi:foreign-slot-value levin 'levin 'terms-used)) + (format t "~&estimated error = ~f" error) + (format t "~&actual error = ~f" (- accelerated-sum zeta2))))))) + +;; term-by-term sum =1.5961632439130233 using 20 terms +;; term-by-term sum =1.5759958390005426 using 13 terms +;; exact value = 1.6449340668482264 +;; accelerated sum = 1.6449340669228176 using 13 terms +;; estimated error = 0.00000000008883604962761638 +;; actual error = 0.00000000007459122208786084 diff --git a/wavelet.lisp b/wavelet.lisp index 72665eff7296fc7124fd48fdf416e79fefa03783..bb2db1b6895f59b20168eb8785a49e31c44944b8 100644 --- a/wavelet.lisp +++ b/wavelet.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: wavelet.lisp -; description: Wavelet transforms. -; date: Mon Nov 26 2007 - 20:43 -; author: Liam Healy -; modified: Sun Dec 2 2007 - 17:10 -;******************************************************** -;;; $Id: $ +;; Wavelet transforms. +;; Liam Healy, Mon Nov 26 2007 - 20:43 +;; Time-stamp: <2008-01-27 18:53:48EST wavelet.lisp> +;; $Id: $ (in-package :gsl) @@ -15,10 +11,14 @@ ;;;; Allocation of wavelets ;;;;**************************************************************************** +(set-asf wavelet allocate-wavelet free-wavelet nil 2) + (defun-gsl allocate-wavelet (type member) "gsl_wavelet_alloc" ((type :pointer) (member :size)) :c-return :pointer + :export nil + :index (with-gsl-objects wavelet) :documentation "Allocate and initialize a wavelet object of type 'type. The parameter 'member selects the specific member of the @@ -29,6 +29,8 @@ "gsl_wavelet_free" ((wavelet :pointer)) :c-return :void + :export nil + :index (with-gsl-objects wavelet) :documentation "Free the wavelet object.") @@ -67,10 +69,14 @@ :documentation "The name of the wavelet family.") +(set-asf wavelet-workspace allocate-wavelet-workspace free-wavelet-workspace) + (defun-gsl allocate-wavelet-workspace (size) "gsl_wavelet_workspace_alloc" ((size :size)) :c-return :pointer + :export nil + :index (with-gsl-objects wavelet-workspace) :documentation "Allocate a workspace for the discrete wavelet transform. To perform a one-dimensional transform on size elements, a workspace @@ -83,6 +89,8 @@ "gsl_wavelet_workspace_free" ((workspace :pointer)) :c-return :void + :export nil + :index (with-gsl-objects wavelet-workspace) :documentation "Free the allocated workspace.") @@ -377,8 +385,9 @@ (let ((n (length cl-data))) (with-data (vector vector-double n) (setf (data vector) cl-data) - (let ((wavelet (allocate-wavelet *daubechies-wavelet* 4)) - (workspace (allocate-wavelet-workspace n))) + (with-gsl-objects + ((wavelet wavelet *daubechies-wavelet* 4) + (wavelet-workspace workspace n)) (wavelet-transform-forward wavelet vector 1 workspace) (with-data (absvector vector-double n) (dotimes (i n) @@ -392,20 +401,16 @@ (dotimes (i n) (format t "~&~a" (gsl-aref vector i))) (cl-invalidate vector) (wavelet-transform-inverse wavelet vector 1 workspace) - (prog1 (data vector) - (free-wavelet-workspace workspace) - (free-wavelet wavelet)))))) - + (data vector))))) (defun wavelet-forward-example (&optional (cl-data *wavelet-sample*)) "Simpler example, with only a Daubechies wavelet forward transformation." (let ((n (length cl-data))) (with-data (vector vector-double n) (setf (data vector) cl-data) - (let ((wavelet (allocate-wavelet *daubechies-wavelet* 4)) - (workspace (allocate-wavelet-workspace n))) + (with-gsl-objects + ((wavelet wavelet *daubechies-wavelet* 4) + (wavelet-workspace workspace n)) (wavelet-transform-forward wavelet vector 1 workspace) (cl-invalidate vector) - (prog1 (data vector) - (free-wavelet-workspace workspace) - (free-wavelet wavelet)))))) + (data vector)))))