diff --git a/chebyshev.lisp b/chebyshev.lisp index 5aca6531e6175d8154610215da97d8dd9bbeab75..c684c128414fd4f34fb334c371a40ff84525442f 100644 --- a/chebyshev.lisp +++ b/chebyshev.lisp @@ -1,6 +1,6 @@ ;; Chebyshev Approximations ;; Liam Healy Sat Nov 17 2007 - 20:36 -;; Time-stamp: <2008-01-26 18:17:32EST chebyshev.lisp> +;; Time-stamp: <2008-02-02 19:13:03EST chebyshev.lisp> ;; $Id: $ (in-package :gsl) @@ -9,14 +9,15 @@ ;;;; Creation and calculation of Chebyshev series ;;;;**************************************************************************** -(set-asf chebyshev allocate-chebyshev free-chebyshev initialize-chebyshev) +(set-asf (chebyshev order function lower-limit upper-limit) + 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) + :index (letm chebyshev) :documentation "Allocate a Chebyshev series of specified order and return a pointer to it.") @@ -26,7 +27,7 @@ ((chebyshev :pointer)) :c-return :void :export nil - :index (with-gsl-objects chebyshev) + :index (letm chebyshev) :documentation "Free a previously allocated Chebyshev series.") @@ -35,7 +36,7 @@ ((chebyshev :pointer) (function :pointer) (lower-limit :double) (upper-limit :double)) :export nil - :index (with-gsl-objects chebyshev) + :index (letm chebyshev) :documentation "Compute the Chebyshev approximation for the function over the range (lower-limit, upper-limit) to the previously specified order. The @@ -125,7 +126,7 @@ (defun chebyshev-table-example () (let ((steps 100)) - (with-gsl-objects ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0)) + (letm ((cheb (chebyshev 40 chebyshev-step 0.0d0 1.0d0))) (dotimes (i steps) (let ((x (coerce (/ i steps) 'double-float))) (format t "~&~a ~a ~a ~a" @@ -136,10 +137,9 @@ (defun chebyshev-point-example (x) (check-type x double-float) - (with-gsl-objects - ((chebyshev cheb 40 chebyshev-step 0.0d0 1.0d0) - (chebyshev deriv 40) - (chebyshev integ 40)) + (letm ((cheb (chebyshev 40 chebyshev-step 0.0d0 1.0d0)) + (deriv (chebyshev 40)) + (integ (chebyshev 40))) (derivative-chebyshev deriv cheb) (integral-chebyshev integ cheb) (list diff --git a/data/data.lisp b/data/data.lisp index 008c1322b8f93cf091e730d7c570124edebe9a1e..93bbed91c80b58d3d007f0486edf0a87e00b8a95 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-27 18:50:53EST data.lisp> +;; Time-stamp: <2008-02-02 23:15:39EST data.lisp> ;; $Id: $ (in-package :gsl) @@ -103,6 +103,18 @@ ;;(LONG-DOUBLE . "_long_double") (COMPLEX . "_complex"))) +(defmacro data-letm (type) + "Define the letm function for data types." + `(defun-letm ,type (size-or-initial &optional zero) + (let ((argsymb (gensym "ARG"))) + (list + `(make-data ',',type ,zero + (if (numberp ,argsymb) ,argsymb (length ,argsymb))) + 'free + (lambda (symb) + `(unless (numberp ,argsymb) (setf (data ,symb) ,argsymb))) + (lambda () (list argsymb size-or-initial)))))) + (defmacro defdata (c-string cl-symbol cl-base-type &optional (superclass 'gsl-data) (dimensions 1)) @@ -115,40 +127,41 @@ (flet ((gsl-name (function-name) (format nil "gsl_~a_~a" c-string function-name))) (let* ((cargs (loop for i below dimensions - collect `((nth ,i (storage-size object)) :size))) + collect `((nth ,i (storage-size object)) :size))) (object-name (make-symbol-from-strings *gsl-prefix* cl-symbol))) `(progn - (defclass ,object-name (,superclass) - ((cl-base-type :initform ',cl-base-type :reader cl-base-type - :allocation :class))) - (defun-gsl alloc ((object ,object-name)) - ,(gsl-name "alloc") ,cargs - :type :method - :c-return (cr :pointer) - :return ((assign-pointer object cr))) - (defun-gsl calloc ((object ,object-name)) - ,(gsl-name "calloc") ,cargs - :type :method - :c-return (cr :pointer) - :return ((assign-pointer object cr))) - (defun-gsl free ((object ,object-name)) - ,(gsl-name "free") (((pointer object) :pointer)) - :type :method - :c-return :void) - (defun-gsl write-binary ((object ,object-name) stream) - ,(gsl-name "fwrite") ((stream :pointer) ((pointer object) :pointer)) - :type :method) - (defun-gsl read-binary ((object ,object-name) stream) - ,(gsl-name "fread") ((stream :pointer) ((pointer object) :pointer)) - :type :method) - (defun-gsl write-formatted ((object ,object-name) stream format) - ,(gsl-name "fprintf") - ((stream :pointer) ((pointer object) :pointer) (format :string)) - :type :method) - (defun-gsl read-formatted ((object ,object-name) stream format) - ,(gsl-name "fscanf") - ((stream :pointer) ((pointer object) :pointer) (format :string)) - :type :method))))) + (data-letm ,cl-symbol) + (defclass ,object-name (,superclass) + ((cl-base-type :initform ',cl-base-type :reader cl-base-type + :allocation :class))) + (defun-gsl alloc ((object ,object-name)) + ,(gsl-name "alloc") ,cargs + :type :method + :c-return (cr :pointer) + :return ((assign-pointer object cr))) + (defun-gsl calloc ((object ,object-name)) + ,(gsl-name "calloc") ,cargs + :type :method + :c-return (cr :pointer) + :return ((assign-pointer object cr))) + (defun-gsl free ((object ,object-name)) + ,(gsl-name "free") (((pointer object) :pointer)) + :type :method + :c-return :void) + (defun-gsl write-binary ((object ,object-name) stream) + ,(gsl-name "fwrite") ((stream :pointer) ((pointer object) :pointer)) + :type :method) + (defun-gsl read-binary ((object ,object-name) stream) + ,(gsl-name "fread") ((stream :pointer) ((pointer object) :pointer)) + :type :method) + (defun-gsl write-formatted ((object ,object-name) stream format) + ,(gsl-name "fprintf") + ((stream :pointer) ((pointer object) :pointer) (format :string)) + :type :method) + (defun-gsl read-formatted ((object ,object-name) stream format) + ,(gsl-name "fscanf") + ((stream :pointer) ((pointer object) :pointer) (format :string)) + :type :method))))) (defun splice-name (base-name keyword symbol) "Make a new C name for a data function from a base name." @@ -198,8 +211,6 @@ (export '(make-data with-data)) -(set-asf generic make-data free (setf data) 3) - (defgeneric alloc (object) (:documentation "Allocate GSL data; used internally.")) diff --git a/data/vector.lisp b/data/vector.lisp index e33c223a3f10997503d3fe36e1ecc00333103dce..8ba0e5cf95efd6e2b24849dfd3065001c3818536 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-26 22:20:43EST vector.lisp> +;; Time-stamp: <2008-02-02 23:13:54EST vector.lisp> ;; $Id: $ (in-package :gsl) @@ -469,11 +469,15 @@ same length.") (setf (data vec) #(-3.21d0 1.0d0 12.8d0)) (data vec)) -(with-gsl-objects ((generic vec 'vector-double nil 3)) +(letm ((vec (vector-double 3))) (setf (data vec) #(-3.21d0 1.0d0 12.8d0)) (data vec)) +(letm ((vec (vector-double #(-3.21d0 1.0d0 12.8d0)))) + (data vec)) +(letm ((vec (vector-double '(-3.21d0 1.0d0 12.8d0)))) + (data vec)) diff --git a/eigensystems.lisp b/eigensystems.lisp index 10e303d5ef4d395bfda90b1754e534f5b492e731..cfddd543b0ae66f74687406f17e7f5b102d1ce13 100644 --- a/eigensystems.lisp +++ b/eigensystems.lisp @@ -1,6 +1,6 @@ ;; Eigenvectors and eigenvalues ;; Liam Healy, Sun May 21 2006 - 19:52 -;; Time-stamp: <2008-01-27 18:53:17EST eigensystems.lisp> +;; Time-stamp: <2008-02-02 21:08:35EST eigensystems.lisp> ;; $Id: $ (in-package :gsl) @@ -9,13 +9,13 @@ ;;;; Real Symmetric Matrices ;;;;**************************************************************************** -(set-asf eigen-symm eigen-symm-alloc eigen-symm-free) +(set-asf (eigen-symm n) 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) + :index '(letm eigen-symm) :documentation ; FDL "Allocate a workspace for computing eigenvalues of n-by-n real symmetric matrices. The size of the workspace @@ -25,17 +25,17 @@ "gsl_eigen_symm_free" ((w :pointer)) :c-return :void :export nil - :index '(with-gsl-objects eigen-symm) + :index '(letm eigen-symm) :documentation ; FDL "Free the memory associated with the workspace w.") -(set-asf eigen-symmv eigen-symmv-alloc eigen-symmv-free) +(set-asf (eigen-symmv n) 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) + :index '(letm eigen-symmv) :documentation ; FDL "Allocate a workspace for computing eigenvalues and eigenvectors of n-by-n real symmetric matrices. The size of @@ -43,7 +43,7 @@ (defun-gsl eigen-symmv-free (w) "gsl_eigen_symmv_free" ((w :pointer)) - :index '(with-gsl-objects eigen-symmv) + :index '(letm eigen-symmv) :c-return :void :documentation ; FDL "Free the memory associated with the workspace w.") @@ -82,7 +82,7 @@ ;;;; Complex Hermitian Matrices ;;;;**************************************************************************** -(set-asf eigen-herm eigen-herm-alloc eigen-herm-free) +(set-asf (eigen-herm n) eigen-herm-alloc eigen-herm-free) (defun-gsl eigen-herm-alloc (n) "gsl_eigen_herm_alloc" ((n :size)) @@ -112,7 +112,7 @@ :invalidate (eval A) :return (eval)) -(set-asf eigen-hermv eigen-hermv-alloc eigen-hermv-free) +(set-asf (eigen-hermv n) eigen-hermv-alloc eigen-hermv-free) (defun-gsl eigen-hermv-alloc (n) "gsl_eigen_hermv_alloc" ((n :size)) @@ -191,7 +191,7 @@ (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)) + (letm ((ws (eigen-symmv 3))) (eigenvalues-eigenvectors-symmetric mat evals evecs ws) (values (data evals) (data evecs))))))) diff --git a/hankel.lisp b/hankel.lisp index 19dcf926d3397eb3bec9c14be8e5ecefa0d23e58..8c4ca0eb59d976964d3e936623afa2a5a41b174c 100644 --- a/hankel.lisp +++ b/hankel.lisp @@ -1,20 +1,20 @@ ;; Discrete Hankel Transforms. ;; Liam Healy, Sat Dec 8 2007 - 16:50 -;; Time-stamp: <2008-01-26 18:36:00EST hankel.lisp> +;; Time-stamp: <2008-02-02 19:24:27EST hankel.lisp> ;; $Id: $ (in-package :gsl) ;;; Everything compiles, but not tested -- need example. -(set-asf hankel allocate-hankel free-hankel init-hankel) +(set-asf (hankel size nu xmax) 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) + :index (letm hankel) :documentation "Allocate a Discrete Hankel transform object of size @var{size}.") @@ -22,7 +22,7 @@ "gsl_dht_new" ((hankel :pointer) (nu :double) (xmax :double)) :export nil - :index (with-gsl-objects hankel) + :index (letm hankel) :documentation "Initialize the transform @var{t} for the given values of @var{nu} and @var{x}.") @@ -44,7 +44,7 @@ ((hankel :pointer)) :c-return :void :export nil - :index (with-gsl-objects hankel) + :index (letm hankel) :documentation "Free the Hankel object.") diff --git a/init/utility.lisp b/init/utility.lisp index 256d2cce8c7c3205a287828431b353ebaf08f2ff..49a79cff2963096117cd8793818733c627ab765f 100644 --- a/init/utility.lisp +++ b/init/utility.lisp @@ -1,6 +1,6 @@ ;; Utility definitions ;; Liam Healy, Sun Dec 3 2006 - 10:21 -;; Time-stamp: <2008-01-27 13:28:29EST utility.lisp> +;; Time-stamp: <2008-02-02 22:48:23EST utility.lisp> ;; $Id: $ (in-package :gsl) @@ -11,86 +11,91 @@ (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." +;;;;;;;;;;;; GSL objects + +;;; Each object has a defun named after it. This function is used +;;; only internally by #'letm, but it is exported so that +;;; arglist prompters such as slime's will show the arguments +;;; needed. + +(defparameter *letm-expanded-object* nil) + +(defvar *not-for-users* nil) + +(defmacro letm (bindings &body body) + (let* (setters ; save setting forms + freers ; save freeing forms + assigns ; save global assignements + (*not-for-users* t) + (bnd ; rewrite the bindings with allocators + (mapcar + (lambda (b) + (let ((defsymb (first b)) (defas (second b))) + (if (and (listp b) (listp defas) + (member (first defas) *letm-expanded-object*)) + ;; It is a GSL object + (destructuring-bind (allocator freer &optional setter assign) + (apply (first defas) (rest defas)) + ;; If there's a setting function, save the form + (when setter (push (funcall setter defsymb) setters)) + ;; Save the freeing form + (push `(,freer ,defsymb) freers) + ;; Save the global assignments + (when assign (push (funcall assign) assigns)) + ;; Return the allocation form to be used in let + `(,defsymb ,allocator)) + ;; Not a GSL object, just define it with let + b))) + bindings))) + (if freers + ;; There were GSL objects + `(let* (,@assigns ,@bnd) + ;; put the body in an unwind-protect + (unwind-protect + (progn + ,@setters ; setting the objects + ,@body) + ,@freers)) ; freeing the objects + ;; No GSL objects in the bindings, just make an ordinary let + `(let* (,@bindings) ,@body)))) + +;;; General definition for object creation +(defmacro defun-letm (symbol arglist &body body) + "Define a GSL object to be bound in a letm. + The body should return a list of two to four values. + The first is the allocation form, the second is + the freeing function name, the third is + form to set or initialize the object, + and the last is a function of no arguments + that gives a form to be bound in a let prior + to binding the object. The symbol will be exported + even though users should not call the function, so that + it will show up in arglist prompters." `(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))) - - -|# + (pushnew ',symbol *letm-expanded-object*) + (export ',symbol) + (defun ,symbol ,arglist + (unless *not-for-users* + (error "The form ~a should be placed in the binding of a letm." + ',(cons symbol arglist))) + ,@body))) + +;;; Specific, simple and common types of object creation +(defmacro set-asf (form allocate free &optional set (num-alloc-args 1)) + "Make an object usable in letm with the arglist is the allocate + and set arglists appended." + (let ((symb (gensym "SASF"))) + `(defun-letm ,(first form) + (,@(subseq (rest form) 0 num-alloc-args) + &optional + ,@(when set `((,(nth num-alloc-args (rest form)) + nil + settingp))) + ,@(when set (subseq (rest form) (1+ num-alloc-args)))) + (list + `(,',allocate ,,@(subseq (rest form) 0 num-alloc-args)) + ',free + ,@(when set + `((when settingp + (lambda (,symb) + `(,',set ,,symb ,,@(subseq (rest form) num-alloc-args)))))))))) diff --git a/interpolation/interpolation.lisp b/interpolation/interpolation.lisp index 4f392668e037e81ed1a9322e49ead64699f0bdc4..eeb76960f8a18f604ca6440c242b927b387582ec 100644 --- a/interpolation/interpolation.lisp +++ b/interpolation/interpolation.lisp @@ -1,6 +1,6 @@ ;; Interpolation allocation, initialization, and freeing. ;; Liam Healy, Sun Nov 4 2007 - 17:24 -;; Time-stamp: <2008-01-26 17:40:04EST interpolation.lisp> +;; Time-stamp: <2008-02-02 22:22:59EST interpolation.lisp> ;; $Id: $ (in-package :gsl) @@ -8,9 +8,13 @@ ;;; 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) +(defun-letm interpolation (type xa-or-size &optional ya) + (list + `(allocate-interpolation ,type ,(if ya `(dim0 ,ya) xa-or-size)) + 'free-spline + (when ya + (lambda (symb) + `(initialize-interpolation ,symb ,xa-or-size ,ya (dim0 ,ya)))))) ;;; Interpolation (defun-gsl allocate-interpolation (type size) @@ -18,7 +22,7 @@ ((type :pointer) (size :size)) :c-return :pointer :export nil - :index (with-gsl-objects interpolation) + :index (letm interpolation) :documentation "Allocate an interpolation object of type for size data-points, and return the pointer to it.") @@ -39,11 +43,17 @@ ((interpolation :pointer)) :c-return :void :export nil - :index (with-gsl-objects interpolation) + :index (letm interpolation) :documentation "Frees the interpolation object @var{interp}.") -(set-asf spline allocate-spline free-spline intialize-spline 2) +(defun-letm spline (type xa-or-size &optional ya) + (list + `(allocate-spline ,type ,(if ya `(dim0 ,ya) xa-or-size)) + 'free-spline + (when ya + (lambda (symb) + `(initialize-spline ,symb ,xa-or-size ,ya (dim0 ,ya)))))) ;;; Spline (defun-gsl allocate-spline (type size) diff --git a/interpolation/lookup.lisp b/interpolation/lookup.lisp index 21993f8c7a1e110a2029264e89f89d0c1893c240..c7a1362c028647d61d1b5a8a99aa605430cea6e4 100644 --- a/interpolation/lookup.lisp +++ b/interpolation/lookup.lisp @@ -9,7 +9,7 @@ (in-package :gsl) -(set-asf acceleration allocate-acceleration free-acceleration nil 0) +(set-asf (acceleration) allocate-acceleration free-acceleration nil 0) (defun-gsl interpolation-search (x-array x low-index high-index) "gsl_interp_bsearch" @@ -25,7 +25,7 @@ () :c-return :pointer :export nil - :index (with-gsl-objects acceleration) + :index (letm acceleration) :documentation "Allocate an accelerator object, which is a kind of iterator for interpolation lookups. It tracks the state of @@ -47,7 +47,7 @@ "gsl_interp_accel_free" ((acceleration :pointer)) :export nil - :index (with-gsl-objects acceleration) + :index (letm acceleration) :c-return :void :documentation "Frees the accelerator object.") diff --git a/linear-least-squares.lisp b/linear-least-squares.lisp index d6db05aff37b2c62893cf199ca6988063fad73f0..2eca08ce1adab8786564197491becaa414a00d05 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-28 22:22:22EST linear-least-squares.lisp> +;; Time-stamp: <2008-02-02 21:01:35EST linear-least-squares.lisp> ;; $Id: $ (in-package :gsl) @@ -125,12 +125,14 @@ ;;;; Multiparameter fitting ;;;;**************************************************************************** -(set-asf fit-workspace allocate-fit-workspace free-fit-workspace nil 2) +(set-asf (fit-workspace number-of-observations number-of-parameters) + allocate-fit-workspace free-fit-workspace nil 2) (defun-gsl allocate-fit-workspace (number-of-observations number-of-parameters) "gsl_multifit_linear_alloc" ((number-of-observations :size) (number-of-parameters :size)) :c-return :pointer + :index (letm fit-workspace) :documentation ; FDL "Allocate a workspace for fitting a linear model.") @@ -138,6 +140,7 @@ "gsl_multifit_linear_free" ((pointer :pointer)) :c-return :void + :index (letm fit-workspace) :documentation ; FDL "Free the memory associate with the workspace.") @@ -310,7 +313,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-gsl-objects ((fit-workspace ws n 3)) + (letm ((ws (fit-workspace 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 a50e4a29bdf2c5b7bb1ef449c941f804d1677a06..79fef62f0f22abe86c01bb768534aecf3ee38a7d 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-28 22:42:05EST minimization-multi.lisp> +;; Time-stamp: <2008-02-02 21:17:03EST minimization-multi.lisp> ;; $Id: $ (in-package :gsl) @@ -48,16 +48,18 @@ ;;;; Initialization ;;;;**************************************************************************** -(set-asf mfminimizer allocate-mfminimizer free-mfminimizer set-mfminimizer 2) -(set-asf - mfdfminimizer allocate-mfdfminimizer free-mfdfminimizer set-mfdfminimizer 2) +;;; Could have one fewer argument: dimension=(dim0 initial) +(set-asf (mfminimizer type dimension function intial step-size) + allocate-mfminimizer free-mfminimizer set-mfminimizer 2) +(set-asf (mfdfminimizer type dimension function-derivative intitial step-size tolerance) + 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) + :index (letm mfminimizer) :documentation ; FDL "Allocate an instance of a minimizer of the given for an function of the given dimensions.") @@ -67,7 +69,7 @@ ((type :pointer) (dimension :size)) :c-return :pointer :export nil - :index (with-gsl-objects mfdfminimizer) + :index (letm mfdfminimizer) :documentation ; FDL "Allocate an instance of a derivative-based minimizer of the given for an function of the given dimensions.") @@ -77,7 +79,7 @@ ((minimizer :pointer) (function :pointer) ((pointer initial) :pointer) ((pointer step-size) :pointer)) :export nil - :index (with-gsl-objects mfminimizer) + :index (letm mfminimizer) :documentation ; FDL "Initialize the minimizer to minimize the function starting from the initial point. @@ -92,7 +94,7 @@ ((pointer initial) :pointer) (step-size :double) (tolerance :double)) :export nil - :index (with-gsl-objects mfdfminimizer) + :index (letm mfdfminimizer) :documentation ; FDL "Initialize the minimizer to minimize the function starting from the initial point. The size of the @@ -108,7 +110,7 @@ ((minimizer :pointer)) :c-return :void :export nil - :index (with-gsl-objects mfminimizer) + :index (letm mfminimizer) :documentation ; FDL "Free all the memory associated with the minimizer.") @@ -117,7 +119,7 @@ ((minimizer :pointer)) :c-return :void :export nil - :index (with-gsl-objects mfdfminimizer) + :index (letm mfdfminimizer) :documentation ; FDL "Free all the memory associated with the minimizer.") @@ -345,9 +347,9 @@ (defun multimin-example-fletcher-reeves () (with-data (initial vector-double 2) (setf (data initial) #(5.0d0 7.0d0)) - (with-gsl-objects - ((mfdfminimizer minimizer *conjugate-fletcher-reeves* 2 parabaloid - initial 0.01d0 1.0d-4)) + (letm ((minimizer + (mfdfminimizer *conjugate-fletcher-reeves* 2 parabaloid + initial 0.01d0 1.0d-4))) (loop with status = T for iter from 0 below 100 while status @@ -378,9 +380,8 @@ (with-data (step-size vector-double 2) (setf (data initial) #(5.0d0 7.0d0)) (set-all step-size 1.0d0) - (with-gsl-objects - ((mfminimizer - minimizer *simplex-nelder-mead* 2 parabaloid-f initial step-size)) + (letm ((minimizer + (mfminimizer *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 76366bf3efe483da76aff683169f3cf60605ab81..bbc05a97d75cccb233ed94cdeb7a782845423d8a 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-26 19:33:45EST minimization-one.lisp> +;; Time-stamp: <2008-02-02 21:12:20EST minimization-one.lisp> ;; $Id: $ (in-package :gsl) @@ -9,14 +9,15 @@ ;;;; Initialization ;;;;**************************************************************************** -(set-asf fminimizer allocate-fminimizer free-fminimizer set-fminimizer) +(set-asf (fminimizer type function minimum lower upper) + 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) + :index (letm fminimizer) :documentation "Allocate an instance of a minimizer of the given type.") @@ -25,13 +26,13 @@ ((minimizer :pointer) (function :pointer) (minimum :double) (lower :double) (upper :double)) :export nil - :index (with-gsl-objects fminimizer) + :index (letm 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? +;;; Use this in letm macro in any way? (defun-gsl set-fminimizer-with-values (minimizer function x-minimum x-lower x-upper f-minimum f-lower f-upper) @@ -50,7 +51,7 @@ ((minimizer :pointer)) :c-return :void :export nil - :index (with-gsl-objects fminimizer) + :index (letm fminimizer) :documentation "Free all the memory associated with the minimizer.") @@ -195,22 +196,21 @@ (defun minimization-one-example () "Solving a minimum, the example given in Sec. 33.8 of the GSL manual." - (let ((max-iter 100)) - (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) - for lower = (fminimizer-x-lower minimizer) - for upper = (fminimizer-x-upper minimizer) - do (iterate-fminimizer minimizer) - while (and (< iter max-iter) - ;; abs and rel error swapped in example? - (not (min-test-interval lower upper 0.001d0 0.0d0))) - do - (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g" - iter lower upper - min (- min pi) - (- upper lower)))))) + (letm ((max-iter 100) + (minimizer + (fminimizer *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) + for lower = (fminimizer-x-lower minimizer) + for upper = (fminimizer-x-upper minimizer) + do (iterate-fminimizer minimizer) + while (and (< iter max-iter) + ;; abs and rel error swapped in example? + (not (min-test-interval lower upper 0.001d0 0.0d0))) + do + (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g" + iter lower upper + min (- min pi) + (- upper lower))))) diff --git a/monte-carlo.lisp b/monte-carlo.lisp index c470824ebc742243d2d9534190945d2cdfcb44a3..758304e3b1a30c7e74685f17048c3b10ba1fbd01 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-26 17:20:19EST monte-carlo.lisp> +;; Time-stamp: <2008-02-02 21:11:35EST monte-carlo.lisp> ;; $Id: $ (in-package :gsl) @@ -13,14 +13,14 @@ (dim :size) (x :pointer)) -(set-asf monte-carlo-plain monte-carlo-plain-alloc monte-carlo-plain-free) +(set-asf (monte-carlo-plain dim) 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-gsl-objects monte-carlo-plain) + :index (letm monte-carlo-plain) :documentation "Allocate and initialize a workspace for Monte Carlo integration in @var{dim} dimensions.") @@ -38,7 +38,7 @@ ((integrator-state :pointer)) :c-return :void :export nil - :index (with-gsl-objects monte-carlo-plain) + :index (letm monte-carlo-plain) :documentation "Frees the memory associated with the integrator.") @@ -96,14 +96,14 @@ (hits-l :pointer) (hits-r :pointer)) -(set-asf monte-carlo-miser monte-carlo-miser-alloc monte-carlo-miser-free) +(set-asf (monte-carlo-miser dim) 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-gsl-objects monte-carlo-miser) + :index (letm monte-carlo-miser) :documentation "Allocate and initialize a workspace for Monte Carlo integration in @var{dim} dimensions. The workspace is used to maintain @@ -122,7 +122,7 @@ ((integrator-state :pointer)) :c-return :void :export nil - :index (with-gsl-objects monte-carlo-miser) + :index (letm monte-carlo-miser) :documentation "Frees the memory associated with the integrator.") @@ -198,14 +198,14 @@ (calls-per-box :uint) (ostream :pointer)) -(set-asf monte-carlo-vegas monte-carlo-vegas-alloc monte-carlo-vegas-free) +(set-asf (monte-carlo-vegas dim) 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-gsl-objects monte-carlo-vegas) + :index (letm monte-carlo-vegas) :documentation "Allocate and initialize a workspace for Monte Carlo integration in @var{dim} dimensions. The workspace is used to maintain @@ -224,7 +224,7 @@ ((integrator-state :pointer)) :c-return :void :export nil - :index (with-gsl-objects monte-carlo-vegas) + :index (letm monte-carlo-vegas) :documentation "Frees the memory associated with the integrator.") @@ -288,7 +288,7 @@ (def-mc-function monte-carlo-g 3) (defun random-walk-plain-example (&optional (nsamples 500000)) - (with-gsl-objects ((monte-carlo-plain ws 3)) + (letm ((ws (monte-carlo-plain 3))) (with-data (lower vector-double 3) (with-data (upper vector-double 3) (setf (data lower) #(0.0d0 0.0d0 0.0d0) @@ -302,7 +302,7 @@ ws))))) (defun random-walk-miser-example (&optional (nsamples 500000)) - (with-gsl-objects ((monte-carlo-miser ws 3)) + (letm ((ws (monte-carlo-miser 3))) (with-data (lower vector-double 3) (with-data (upper vector-double 3) (setf (data lower) #(0.0d0 0.0d0 0.0d0) @@ -316,7 +316,7 @@ ws))))) (defun random-walk-vegas-example (&optional (nsamples 500000)) - (with-gsl-objects ((monte-carlo-vegas ws 3)) + (letm ((ws (monte-carlo-vegas 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 a92034f437d3a7ab4ebab56e055967bec765f730..217ec44fcfcad76333321898db7952e1dd100d65 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-26 17:11:25EST numerical-integration.lisp> +;; Time-stamp: <2008-02-02 21:06:00EST numerical-integration.lisp> ;; $Id: $ ;;; To do: QAWS, QAWO, QAWF, more tests @@ -36,14 +36,14 @@ ;;;; QAG adaptive Gauss-Kronrod integration ;;;;**************************************************************************** -(set-asf - integration-workspace integration-workspace-alloc integration-workspace-free) +(set-asf (integration-workspace size) + 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) + :index (letm integration-workspace) :documentation "Allocate a workspace sufficient to hold @var{n} double precision intervals, their integration results and error estimates.") @@ -51,7 +51,7 @@ "gsl_integration_workspace_free" ((pointer :pointer)) :c-return :void :export nil - :index (with-gsl-objects integration-workspace) + :index (letm integration-workspace) :documentation "Free the memory associated with the workspace @var{w}.") (cffi:defcenum integrate-method @@ -228,9 +228,9 @@ (integration-qng one-sine 0.0d0 pi)) (lisp-unit:assert-first-fp-equal "0.200000000000d+01" - (with-gsl-objects ((integration-workspace ws 20)) + (letm ((ws (integration-workspace 20))) (integration-QAG one-sine 0.0d0 pi :gauss15 20 ws))) (lisp-unit:assert-error 'gsl-error - (with-gsl-objects ((integration-workspace ws 20)) + (letm ((ws (integration-workspace ws 20))) (integration-QAG one-sine 0.0d0 pi :gauss15 50 ws)))) diff --git a/polynomial.lisp b/polynomial.lisp index 961921dfeb57c23e688bcb10acd6e0dbc73582da..c79a5c2b07fc9fb8013848a01286b7e8e586b638 100644 --- a/polynomial.lisp +++ b/polynomial.lisp @@ -130,19 +130,20 @@ ;;;; General Polynomial Equations ;;;;**************************************************************************** -(set-asf complex-workspace allocate-complex-workspace free-complex-workspace) +(set-asf (complex-workspace n) + allocate-complex-workspace free-complex-workspace) (defun-gsl allocate-complex-workspace (n) "gsl_poly_complex_workspace_alloc" ((n :size)) :c-return :pointer :export nil - :index (with-gsl-object complex-workspace)) + :index (letm complex-workspace)) (defun-gsl free-complex-workspace (ws) "gsl_poly_complex_workspace_free" ((ws :pointer)) :c-return :void :export nil - :index (with-gsl-object complex-workspace)) + :index (letm complex-workspace)) (defun polynomial-solve (coefficients) "The roots of the general polynomial @@ -157,7 +158,7 @@ (with-data (coef vector-double len) (setf (data coef) coefficients) (with-data (answer vector-double ((* 2 (1- len)))) - (with-gsl-objects ((complex-workspace ws len)) + (letm ((ws (complex-workspace len))) (values-list (polynomial-solve-ws coef ws answer))))))) (defun-gsl polynomial-solve-ws (coefficients workspace answer-pd) diff --git a/random/beta.lisp b/random/beta.lisp index 8d3defaf5060acb0f9df16e3e6c06994922a85a7..2e506106a3ce00f3213a5a4ccd221a1d507e5a04 100644 --- a/random/beta.lisp +++ b/random/beta.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: beta.lisp -; description: Beta distribution -; date: Sat Sep 30 2006 -; author: Liam M. Healy -; modified: Sun Oct 8 2006 - 17:05 -;******************************************************** -;;; $Id: $ +;; Beta distribution +;; Liam Healy, Sat Sep 30 2006 +;; Time-stamp: <2008-02-02 23:04:13EST beta.lisp> +;; $Id: $ (in-package :gsl) @@ -59,11 +55,10 @@ "0.230096194773d+00" "0.392834882565d+00" "0.514387412254d+00" "0.233783685805d+00" "0.198512886686d+00") (lisp-unit:fp-sequence - (progn - (rng-set *rng-mt19937* 0) + (letm ((rng (random-number-generator *mt19937* 0))) (loop for i from 0 to 10 collect - (beta-rd *rng-mt19937* 1.0d0 2.0d0))))) + (beta-rd rng 1.0d0 2.0d0))))) (lisp-unit:assert-first-fp-equal "0.180000000000d+01" (beta-pdf 0.1d0 1.0d0 2.0d0)) diff --git a/random/dirichlet.lisp b/random/dirichlet.lisp index c48e9215104060e78b159290898deb8868200c71..618c46f2a496b99dfbd3100e2baedad70328909b 100644 --- a/random/dirichlet.lisp +++ b/random/dirichlet.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: dirichlet.lisp -; description: Dirichlet distribution -; date: Sun Oct 29 2006 -; author: Liam M. Healy -; modified: Sat Mar 3 2007 - 22:40 -;******************************************************** -;;; $Id: $ +;; Dirichlet distribution +;; Liam Healy, Sun Oct 29 2006 +;; Time-stamp: <2008-02-02 23:00:40EST dirichlet.lisp> +;; $Id: $ (in-package :gsl) @@ -64,24 +60,18 @@ '("0.392833324564d-04" "0.468176310887d+00" "0.340750440311d+00" "0.191033965469d+00") (lisp-unit:fp-sequence - (progn - (rng-set *rng-mt19937* 0) - (with-data (alpha vector-double 4) - (with-data (theta vector-double 4) - (setf (data alpha) #(1.0d0 2.0d0 3.0d0 4.0d0)) - (dirichlet *rng-mt19937* alpha theta) - (data theta)))))) + (letm ((rng (random-number-generator *mt19937* 0)) + (alpha (vector-double #(1.0d0 2.0d0 3.0d0 4.0d0))) + (theta (vector-double 4))) + (dirichlet rng alpha theta) + (data theta)))) (lisp-unit:assert-first-fp-equal "0.288000000000d+01" - (with-data (alpha vector-double 4) - (with-data (theta vector-double 4) - (setf (data alpha) #(1.0d0 2.0d0 3.0d0 4.0d0) - (data theta) #(0.1d0 0.3d0 0.4d0 0.2d0)) - (dirichlet-pdf alpha theta)))) + (letm ((alpha (vector-double #(1.0d0 2.0d0 3.0d0 4.0d0))) + (theta (vector-double #(0.1d0 0.3d0 0.4d0 0.2d0)))) + (dirichlet-pdf alpha theta))) (lisp-unit:assert-first-fp-equal "0.105779029415d+01" - (with-data (alpha vector-double 4) - (with-data (theta vector-double 4) - (setf (data alpha) #(1.0d0 2.0d0 3.0d0 4.0d0) - (data theta) #(0.1d0 0.3d0 0.4d0 0.2d0)) - (dirichlet-log-pdf alpha theta))))) + (letm ((alpha (vector-double #(1.0d0 2.0d0 3.0d0 4.0d0))) + (theta (vector-double #(0.1d0 0.3d0 0.4d0 0.2d0)))) + (dirichlet-log-pdf alpha theta)))) diff --git a/random/discrete.lisp b/random/discrete.lisp index 38780efc706f870814414c12ae30dd5f9bd0960a..064a6ceabcc79bfcee2b9506b6d3a74a3d39330c 100644 --- a/random/discrete.lisp +++ b/random/discrete.lisp @@ -1,6 +1,6 @@ ;; Discrete random variables ;; Liam Healy, Sat Nov 11 2006 - 21:51 -;; Time-stamp: <2008-01-26 17:22:14EST discrete.lisp> +;; Time-stamp: <2008-02-02 20:59:58EST discrete.lisp> ;; $Id: $ (in-package :gsl) @@ -11,14 +11,14 @@ (A :pointer) (F :pointer)) -(set-asf discrete discrete-preprocess discrete-free) +(set-asf (discrete-random probabilities) 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) + :index (letm discrete) :documentation "A pointer to a structure that contains the lookup table for the discrete random number generator. The array probabilities contains @@ -31,7 +31,7 @@ "gsl_ran_discrete_free" ((table :pointer)) :c-return :void :export nil - :index (with-gsl-objects discrete) + :index (letm discrete) :documentation "De-allocates the lookup table created by #'discrete-preprocess.") @@ -68,8 +68,8 @@ '(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-gsl-objects ((discrete table probabilities)) - (rng-set *rng-mt19937* 0) + (letm ((table (discrete-random probabilities)) + (rng (random-number-generator *mt19937* 0))) (loop for i from 0 to 10 collect (discrete *rng-mt19937* table)))) @@ -78,7 +78,5 @@ (with-data (probabilities vector-double 3) (setf (data probabilities) #(0.25d0 0.5d0 0.25d0)) - (with-gsl-objects ((discrete table probabilities)) - (discrete-pdf - 1 - table)))))) + (letm ((table (discrete-random probabilities))) + (discrete-pdf 1 table)))))) diff --git a/random/generators.lisp b/random/generators.lisp index b2ae4182d62d430f6b06495c8907f2c4be39d33d..7b46f30b60e02631c5731f49fd5d2e2241939e24 100644 --- a/random/generators.lisp +++ b/random/generators.lisp @@ -1,13 +1,13 @@ ;; Generators of random numbers. ;; Liam Healy, Sat Jul 15 2006 - 14:43 -;; Time-stamp: <2008-01-31 22:43:53EST generators.lisp> +;; Time-stamp: <2008-02-02 23:01:02EST generators.lisp> ;; $Id: $ (in-package :gsl) ;;; Would it be useful to have a make-data like macro? -(export '(random-number-generator make-random-number-generator)) +(export '(make-random-number-generator)) (defgeneric generator (object) (:method ((object t)) @@ -38,7 +38,8 @@ :generator generator))) (if (eq generator t) (alloc instance) instance))) -(set-asf random-number-generator make-random-number-generator free rng-set 1) +(set-asf (random-number-generator type value) + make-random-number-generator free rng-set 1) ;;;;**************************************************************************** ;;;; Initialization @@ -243,14 +244,14 @@ ;;;; Examples and unit test ;;;;**************************************************************************** -(defparameter *rng-mt19937* (make-random-number-generator *mt19937*)) -(defparameter *rng-cmrg* (make-random-number-generator *cmrg*)) +;;; (defparameter *rng-mt19937* (make-random-number-generator *mt19937*)) +;;; (defparameter *rng-cmrg* (make-random-number-generator *cmrg*)) ;;; (defparameter *rng-default* (make-random-number-generator *default-type*)) (lisp-unit:define-test random-number-generators (lisp-unit:assert-equal '(999 162 282 947 231 484 957 744 540 739 759) - (with-gsl-objects ((random-number-generator rng *mt19937* 0)) + (letm ((rng (random-number-generator *mt19937* 0))) (loop for i from 0 to 10 collect (uniform-fixnum rng 1000)))) @@ -260,8 +261,5 @@ "0.554133304187d+00" "0.880695776958d+00" "0.597139396983d+00" "0.751874113340d+00" "0.931108462127d+00") (lisp-unit:fp-sequence - (with-gsl-objects ((random-number-generator rng *cmrg* 0)) + (letm ((rng (random-number-generator *cmrg* 0))) (loop for i from 0 to 10 collect (uniform rng)))))) - - - diff --git a/random/gumbel1.lisp b/random/gumbel1.lisp index a99ac855f0eb1a74d76e6dbc393793783b1d9ccb..d259d6f090d49845fc3425bf00a7bcc13995eeae 100644 --- a/random/gumbel1.lisp +++ b/random/gumbel1.lisp @@ -63,11 +63,9 @@ "0.382920819366d+01" "0.191289739318d+01" "0.117748457895d+01" "0.189323210797d+01" "0.198591186168d+01") (lisp-unit:fp-sequence - (progn - (rng-set *rng-mt19937* 0) + (letm ((rng (random-number-generator *mt19937* 0))) (loop for i from 0 to 10 - collect - (gumbel1 *rng-mt19937* 1.0d0 2.0d0))))) + collect (gumbel1 rng 1.0d0 2.0d0))))) (lisp-unit:assert-first-fp-equal "0.296257089650d+00" (gumbel1-pdf 0.1d0 1.0d0 2.0d0)) diff --git a/random/poisson.lisp b/random/poisson.lisp index fa18204dee33a274d95cf7d5f309116440535fc7..70d5bf268d666732c39489e024e7ae5c69d2f40d 100644 --- a/random/poisson.lisp +++ b/random/poisson.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: poisson.lisp -; description: Poisson distribution -; date: Sat Nov 25 2006 - 16:00 -; author: Liam M. Healy -; modified: Sat Nov 25 2006 - 17:00 -;******************************************************** -;;; $Id: $ +;; Poisson distribution +;; Liam Healy, Sat Nov 25 2006 - 16:00 +;; Time-stamp: <2008-02-02 23:02:43EST poisson.lisp> +;; $Id: $ (in-package :gsl) @@ -43,11 +39,10 @@ (lisp-unit:define-test poisson (lisp-unit:assert-equal '(15 6 9 9 5 8 11 9 11 5 10) - (progn - (rng-set *rng-mt19937* 0) + (letm ((rng (random-number-generator *mt19937* 0))) (loop for i from 0 to 10 - collect - (poisson *rng-mt19937* 10.0d0)))) + collect + (poisson rng 10.0d0)))) (lisp-unit:assert-first-fp-equal "0.112599032149d+00" (poisson-pdf 8 10.0d0)) diff --git a/random/quasi.lisp b/random/quasi.lisp index d71d94418257f688becb4b6daafca26cc2379157..91e9349ab42075b72444aa60b0f8560e05a6fffe 100644 --- a/random/quasi.lisp +++ b/random/quasi.lisp @@ -1,6 +1,6 @@ ;; Quasi-random sequences in arbitrary dimensions. ;; Liam Healy, Sun Jul 16 2006 - 15:54 -;; Time-stamp: <2008-01-31 22:52:29EST quasi.lisp> +;; Time-stamp: <2008-02-02 22:01:47EST quasi.lisp> ;; $Id: $ (in-package :gsl) @@ -11,7 +11,7 @@ (defparameter *default-quasi-random-number-generator* nil) -(set-asf quasi-random-number-generator +(set-asf (quasi-random-number-generator dimension type) make-quasi-random-number-generator free nil 2) (defun make-quasi-random-number-generator @@ -128,10 +128,8 @@ ;; This example is given in the GSL documentation (lisp-unit:assert-equal '(0.5d0 0.5d0 0.75d0 0.25d0 0.25d0 0.75d0 0.375d0 0.375d0 0.875d0 0.875d0) - (with-data (vec vector-double 2) - (with-gsl-objects ((quasi-random-number-generator gen 2 *sobol*)) - (loop repeat 5 - do (qrng-get gen vec) - append (coerce (data vec) 'list)))))) - - + (letm ((gen (quasi-random-number-generator 2 *sobol*)) + (vec (vector-double 2))) + (loop repeat 5 + do (qrng-get gen vec) + append (coerce (data vec) 'list))))) diff --git a/random/rayleigh-tail.lisp b/random/rayleigh-tail.lisp index 18f7d3c52b5db87407670929a3e346e7a75598e0..a4671180fb9161d9ef803e100dae143de6eb11ce 100644 --- a/random/rayleigh-tail.lisp +++ b/random/rayleigh-tail.lisp @@ -1,11 +1,7 @@ -;******************************************************** -; file: rayleigh-tail.lisp -; description: Rayleigh tail distribution -; date: Sat Sep 30 2006 -; author: Liam M. Healy -; modified: Sat Sep 30 2006 - 19:51 -;******************************************************** -;;; $Id: $ +;; Rayleigh tail distribution +;; Liam Healy, Sat Sep 30 2006 +;; Time-stamp: <2008-02-02 23:00:39EST rayleigh-tail.lisp> +;; $Id: $ (in-package :gsl) @@ -36,11 +32,9 @@ "0.311299291669d+01" "0.774988930120d+01" "0.111454501381d+02" "0.782519818732d+01" "0.747677468155d+01") (lisp-unit:fp-sequence - (progn - (rng-set *rng-mt19937* 0) + (letm ((rng (random-number-generator *mt19937* 0))) (loop for i from 0 to 10 - collect - (rayleigh-tail *rng-mt19937* 1.0d0 10.0d0))))) + collect (rayleigh-tail rng 1.0d0 10.0d0))))) (lisp-unit:assert-first-fp-equal "0.102243176249d+00" (rayleigh-tail-pdf 0.25d0 -2.0d0 2.0d0))) diff --git a/roots-multi.lisp b/roots-multi.lisp index 55ccde6bf76b69d3e32773813df7b0e7acf58fd0..6270607d1613ed2c409b0db431606e05fdb66855 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-27 23:04:02EST roots-multi.lisp> +;;; Time-stamp: <2008-02-02 22:22:20EST roots-multi.lisp> ;;; $Id: $ (in-package :gsl) @@ -43,15 +43,26 @@ ;;;; Initialization ;;;;**************************************************************************** -(set-asf mfsolver allocate-mfsolver free-mfsolver set-mfsolver 2) -(set-asf mfdfsolver allocate-mfdfsolver free-mfdfsolver set-mfdfsolver 2) +(defun-letm mfsolver (type function-derivative initial) + (list + `(allocate-mfsolver ,type (dim0 ,initial)) + 'free-mfsolver + (lambda (symb) + `(set-mfsolver ,symb ,function-derivative ,initial)))) + +(defun-letm mfdfsolver (type function-derivative initial) + (list + `(allocate-mfdfsolver ,type (dim0 ,initial)) + 'free-mfdfsolver + (lambda (symb) + `(set-mfdfsolver ,symb ,function-derivative ,initial)))) (defun-gsl allocate-mfsolver (type dimension) "gsl_multiroot_fsolver_alloc" ((type :pointer) (dimension :size)) :c-return :pointer :export nil - :index (with-gsl-objects mfsolver) + :index (letm mfsolver) :documentation "Allocate an instance of a solver of the type specified for a system of the specified number of dimensions.") @@ -61,7 +72,7 @@ ((type :pointer) (dimension :size)) :c-return :pointer :export nil - :index (with-gsl-objects mfdfsolver) + :index (letm mfdfsolver) :documentation "Allocate an instance of a derivative solver of the type specified for a system of the specified number of dimensions.") @@ -70,7 +81,7 @@ "gsl_multiroot_fsolver_set" ((solver :pointer) (function :pointer) ((pointer initial) :pointer)) :export nil - :index (with-gsl-objects mfsolver) + :index (letm mfsolver) :documentation "Set or reset an existing solver to use the function and the initial guess gsl-vector.") @@ -80,7 +91,7 @@ ((solver :pointer) (function-derivative :pointer) ((pointer initial) :pointer)) :export nil - :index (with-gsl-objects mfdfsolver) + :index (letm mfdfsolver) :documentation "Set or reset an existing solver to use the function and derivative (fdf) and the initial guess.") @@ -90,7 +101,7 @@ ((solver :pointer)) :c-return :void :export nil - :index (with-gsl-objects mfsolver) + :index (letm mfsolver) :documentation "Free all the memory associated with the solver.") @@ -99,7 +110,7 @@ ((solver :pointer)) :c-return :void :export nil - :index (with-gsl-objects mfdfsolver) + :index (letm mfdfsolver) :documentation "Free all the memory associated with the solver.") @@ -415,8 +426,7 @@ (let ((max-iter 1000)) (with-data (vect vector-double 2) (setf (data vect) #(-10.0d0 -5.0d0)) - (with-gsl-objects - ((mfsolver solver *hybrid-scaled* (dim0 vect) rosenbrock vect)) + (letm ((solver (mfsolver *hybrid-scaled* rosenbrock vect))) (let ((fnval (mfsolver-f solver)) (argval (mfsolver-root solver))) (loop for iter from 0 @@ -469,8 +479,8 @@ (let ((max-iter 1000)) (with-data (vect vector-double 2) (setf (data vect) #(-10.0d0 -5.0d0)) - (with-gsl-objects - ((mfdfsolver solver *gnewton-mfdfsolver* (dim0 vect) rosenbrock-f vect)) + (letm + ((solver (mfdfsolver *gnewton-mfdfsolver* 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 4abdf3d672e998790440168a018e08fd9aed46c5..8c1c1adcb1c75f892d6d66fb9e9ca34b18194590 100644 --- a/roots-one.lisp +++ b/roots-one.lisp @@ -1,6 +1,6 @@ ;; One-dimensional root solver. ;; Liam Healy -;; Time-stamp: <2008-01-26 19:07:22EST roots-one.lisp> +;; Time-stamp: <2008-02-02 19:29:21EST roots-one.lisp> ;; $Id: $ (in-package :gsl) @@ -44,15 +44,16 @@ ;;;; Initialization ;;;;**************************************************************************** -(set-asf fsolver allocate-fsolver free-fsolver set-fsolver) -(set-asf fdfsolver allocate-fdfsolver free-fdfsolver set-fdfsolver) +(set-asf (fsolver type function lower upper) allocate-fsolver free-fsolver set-fsolver) +(set-asf (fdfsolver type function-derivative root-guess) + 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) + :index (letm fsolver) :documentation "Allocate an instance of a solver of the given type.") @@ -61,7 +62,7 @@ ((type :pointer)) :c-return :pointer :export nil - :index (with-gsl-objects fdfsolver) + :index (letm fdfsolver) :documentation "Allocate an instance of a derivative-based solver of the given type.") @@ -70,7 +71,7 @@ ((solver :pointer) (function :pointer) (lower :double) (upper :double)) :export nil - :index (with-gsl-objects fsolver) + :index (letm fsolver) :documentation "Initialize or reinitialize an existing solver to use the function and the initial search interval @@ -80,7 +81,7 @@ "gsl_root_fdfsolver_set" ((solver :pointer) (function-derivative :pointer) (root-guess :double)) :export nil - :index (with-gsl-objects fdfsolver) + :index (letm fdfsolver) :documentation "Initialize or reinitialize an existing solver.") @@ -89,7 +90,7 @@ ((solver :pointer)) :c-return :void :export nil - :index (with-gsl-objects fsolver) + :index (letm fsolver) :documentation "Free all the memory associated with the solver.") @@ -98,7 +99,7 @@ ((solver :pointer)) :c-return :void :export nil - :index (with-gsl-objects fdfsolver) + :index (letm fdfsolver) :documentation "Free all the memory associated with the solver.") @@ -353,21 +354,21 @@ (defun roots-one-example () "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." - (let ((max-iter 50)) - (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) - for lower = (fsolver-lower solver) - for upper = (fsolver-upper solver) - do (iterate-fsolver solver) - while (and (< iter max-iter) - (not (root-test-interval lower upper 0.0d0 0.001d0))) - do - (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g" - iter lower upper - root (- root (sqrt 5.0d0)) - (- upper lower)))))) + (letm ((max-iter 50) + (solver (fsolver *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) + for lower = (fsolver-lower solver) + for upper = (fsolver-upper solver) + do (iterate-fsolver solver) + while (and (< iter max-iter) + (not (root-test-interval lower upper 0.0d0 0.001d0))) + do + (format t "~&~d~6t~10,6f~18t~10,6f~28t~12,9f ~44t~10,4g ~10,4g" + iter lower upper + root (- root (sqrt 5.0d0)) + (- upper lower))))) ;;; Because def-solver-functions and def-single-function bind a symbol ;;; of the same name as the first function, and we want both to run, @@ -382,7 +383,7 @@ "Solving a quadratic, the example given in Sec. 32.10 of the GSL manual." (let ((max-iter 100) (initial 5.0d0)) - (with-gsl-objects ((fdfsolver solver *newton-fdfsolver* quadratic-df initial)) + (letm ((solver (fdfsolver *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 11b7c1dc458e380a21fc3f25e9af41b8f58e42c8..2ae1a34a53024d38938f34c7e17affe6c45fc939 100644 --- a/series-acceleration.lisp +++ b/series-acceleration.lisp @@ -1,6 +1,6 @@ ;; Series acceleration. ;; Liam Healy, Wed Nov 21 2007 - 18:41 -;; Time-stamp: <2008-01-26 18:31:15EST series-acceleration.lisp> +;; Time-stamp: <2008-02-02 19:33:39EST series-acceleration.lisp> ;; $Id: $ (in-package :gsl) @@ -21,14 +21,14 @@ (dq-den :pointer) (dsum :pointer)) -(set-asf levin allocate-levin free-levin) +(set-asf (levin order) 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) + :index (letm 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,7 +38,7 @@ ((levin :pointer)) :c-return :void ; Error in GSL manual, should be void? :export nil - :index (with-gsl-objects levin) + :index (letm levin) :documentation "Free a previously allocated Levin acceleration.") @@ -61,14 +61,14 @@ ;;;; Acceleration with error estimation from truncation ;;;;**************************************************************************** -(set-asf levin-truncated allocate-levin-truncated free-levin-truncated) +(set-asf (levin-truncated order) 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) + :index (letm 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 @@ -78,7 +78,7 @@ "gsl_sum_levin_utrunc_free" ((levin :pointer)) :export nil - :index (with-gsl-objects levin-truncated) + :index (letm levin-truncated) :documentation "Free a previously allocated Levin acceleration.") @@ -109,7 +109,7 @@ (let ((maxterms 20) (sum 0.0d0) (zeta2 (/ (expt pi 2) 6))) - (with-gsl-objects ((levin levin maxterms)) + (letm ((levin (levin maxterms))) (with-data (array vector-double maxterms) (dotimes (n maxterms) (setf (gsl-aref array n) (coerce (/ (expt (1+ n) 2)) 'double-float)) diff --git a/wavelet.lisp b/wavelet.lisp index bb2db1b6895f59b20168eb8785a49e31c44944b8..d6acb359d3d27d532f883be55899b4512c63f1af 100644 --- a/wavelet.lisp +++ b/wavelet.lisp @@ -1,6 +1,6 @@ ;; Wavelet transforms. ;; Liam Healy, Mon Nov 26 2007 - 20:43 -;; Time-stamp: <2008-01-27 18:53:48EST wavelet.lisp> +;; Time-stamp: <2008-02-02 21:12:08EST wavelet.lisp> ;; $Id: $ (in-package :gsl) @@ -11,14 +11,14 @@ ;;;; Allocation of wavelets ;;;;**************************************************************************** -(set-asf wavelet allocate-wavelet free-wavelet nil 2) +(set-asf (wavelet type member) 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) + :index (letm wavelet) :documentation "Allocate and initialize a wavelet object of type 'type. The parameter 'member selects the specific member of the @@ -30,7 +30,7 @@ ((wavelet :pointer)) :c-return :void :export nil - :index (with-gsl-objects wavelet) + :index (letm wavelet) :documentation "Free the wavelet object.") @@ -69,14 +69,14 @@ :documentation "The name of the wavelet family.") -(set-asf wavelet-workspace allocate-wavelet-workspace free-wavelet-workspace) +(set-asf (wavelet-workspace size) 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) + :index (letm wavelet-workspace) :documentation "Allocate a workspace for the discrete wavelet transform. To perform a one-dimensional transform on size elements, a workspace @@ -90,7 +90,7 @@ ((workspace :pointer)) :c-return :void :export nil - :index (with-gsl-objects wavelet-workspace) + :index (letm wavelet-workspace) :documentation "Free the allocated workspace.") @@ -385,9 +385,8 @@ (let ((n (length cl-data))) (with-data (vector vector-double n) (setf (data vector) cl-data) - (with-gsl-objects - ((wavelet wavelet *daubechies-wavelet* 4) - (wavelet-workspace workspace n)) + (letm ((wavelet (wavelet *daubechies-wavelet* 4)) + (workspace (wavelet-workspace n))) (wavelet-transform-forward wavelet vector 1 workspace) (with-data (absvector vector-double n) (dotimes (i n) @@ -408,9 +407,8 @@ (let ((n (length cl-data))) (with-data (vector vector-double n) (setf (data vector) cl-data) - (with-gsl-objects - ((wavelet wavelet *daubechies-wavelet* 4) - (wavelet-workspace workspace n)) + (letm ((wavelet (wavelet *daubechies-wavelet* 4)) + (workspace (wavelet-workspace n))) (wavelet-transform-forward wavelet vector 1 workspace) (cl-invalidate vector) (data vector)))))