From 2035b727ec9f1ef066f581f2d2e767a91acef2b1 Mon Sep 17 00:00:00 2001 From: liam <liam@a3d8a0fb-c1db-0310-ace7-a616afeb9e30> Date: Mon, 28 Jan 2008 04:08:06 +0000 Subject: [PATCH] Introduction of #'with-gsl-objects macro and associated definitions, and application to many chapters. Most former with- macros have been removed in favor of this macro. All previous results have been reproduced. There is an experimental application to data (vectors, matrices, etc.) but it is clumsy and setting doesn't work; this requires more thought. git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3276 a3d8a0fb-c1db-0310-ace7-a616afeb9e30 --- README | 6 ++ chebyshev.lisp | 52 ++++++++--------- data/data.lisp | 31 ++++++++++- data/vector.lisp | 10 +++- eigensystems.lisp | 84 ++++++++++++---------------- hankel.lisp | 23 +++++--- index.html | 4 +- init/utility.lisp | 96 +++++++++++++++++++++++++++++--- interpolation/interpolation.lisp | 22 +++++--- interpolation/lookup.lisp | 6 +- linear-least-squares.lisp | 15 +---- minimization-multi.lisp | 92 ++++++++++++++---------------- minimization-one.lisp | 29 +++++----- monte-carlo.lisp | 56 ++++++------------- numerical-integration.lisp | 21 ++++--- random/discrete.lisp | 26 +++++---- roots-multi.lisp | 46 +++++++-------- roots-one.lisp | 45 ++++++--------- series-acceleration.lisp | 75 ++++++++++++++----------- wavelet.lisp | 43 +++++++------- 20 files changed, 436 insertions(+), 346 deletions(-) diff --git a/README b/README index a32017d5..fe33b407 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 0dcce83b..5aca6531 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 5a7be85e..008c1322 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 c5159706..e33c223a 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 bb2c23b5..10e303d5 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 4939ead1..19dcf926 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 cf5bba7e..84f49a6b 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 94d1d208..256d2cce 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 9147c0d9..4f392668 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 a3aadd3d..21993f8c 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 6312723d..f0971047 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 ccd9324e..7878cd07 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 83b7089f..76366bf3 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 88143676..c470824e 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 958bccfe..a92034f4 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 15e4548b..38780efc 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 782c30a5..55ccde6b 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 175bb1af..4abdf3d6 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 1074ba4b..11b7c1dc 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 72665eff..bb2db1b6 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))))) -- GitLab