Commit f0c7ecec authored by Liam Healy's avatar Liam Healy

Define callbacks by returning the struct pointer, without setting specials

Macro defcbstruct and all the macros built on it now simply return the
C struct needed by GSL functions instead of binding a special whose
name is then passed to the GSL functions.  Specials can still be used
of course in the standard way like defparameter, but the results can
also be bound locally with a let or passed in directly.  The name of a
callback passed to defmcallback and to defcbstruct in all the macros
that expand both is now gensymmed (since it is used only long enough
to put in a slot of the cbstruct) instead of taken to be the name of
the CL function, in order to make clear the different objects.
Unfortunately, there is appears to be a bug in CFFI using defcallback
with an uninterned symbol in CCL, so I have put a workaround in.
However, there is still another error even with this workaround.
Thus test results are:
SBCL 64: TOTAL: 1248 assertions passed, 5 failed, 0 execution errors.
CCL 64: TOTAL: 1229 assertions passed, 2 failed, 6 execution errors.
parent 42cf5a7b
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2009-01-19 16:30:06EST monte-carlo.lisp>
;; Time-stamp: <2009-01-24 12:58:12EST monte-carlo.lisp>
;; $Id$
(in-package :gsl)
......@@ -217,9 +217,9 @@
(dimensions sizet)
(parameters :pointer))
(export 'def-mc-function)
(defmacro def-mc-function (name dimensions)
`(def-single-function ,name :double :double monte-function
(export 'make-monte-carlo-function)
(defmacro make-monte-carlo-function (name dimensions)
`(make-single-function ,name :double :double monte-function
,dimensions nil nil))
;;;;****************************************************************************
......@@ -233,28 +233,29 @@
(* (/ (expt pi 3))
(/ (- 1 (* (cos x) (cos y) (cos z))))))
(def-mc-function monte-carlo-g 3)
(defparameter *monte-carlo-g*
(make-monte-carlo-function monte-carlo-g 3))
(defun random-walk-plain-example (&optional (nsamples 500000))
(let ((ws (make-monte-carlo-plain 3))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
(rng (make-random-number-generator *mt19937* 0)))
(monte-carlo-integrate-plain monte-carlo-g lower upper nsamples rng ws)))
(monte-carlo-integrate-plain *monte-carlo-g* lower upper nsamples rng ws)))
(defun random-walk-miser-example (&optional (nsamples 500000))
(let ((ws (make-monte-carlo-miser 3))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
(rng (make-random-number-generator *mt19937* 0)))
(monte-carlo-integrate-miser monte-carlo-g lower upper nsamples rng ws)))
(monte-carlo-integrate-miser *monte-carlo-g* lower upper nsamples rng ws)))
(defun random-walk-vegas-example (&optional (nsamples 500000))
(let ((ws (make-monte-carlo-vegas 3))
(lower #m(0.0d0 0.0d0 0.0d0))
(upper (make-marray 'double-float :initial-contents (list pi pi pi)))
(rng (make-random-number-generator *mt19937* 0)))
(monte-carlo-integrate-vegas monte-carlo-g lower upper nsamples rng ws)))
(monte-carlo-integrate-vegas *monte-carlo-g* lower upper nsamples rng ws)))
(save-test monte-carlo
(random-walk-plain-example)
......
;; Numerical differentiation.
;; Liam Healy Mon Nov 12 2007 - 22:07
;; Time-stamp: <2008-11-15 22:31:32EST numerical-differentiation.lisp>
;; Time-stamp: <2009-01-24 14:57:26EST numerical-differentiation.lisp>
;; $Id$
(in-package :gsl)
......@@ -70,56 +70,35 @@
;;; Examples from gsl-1.11/deriv/test.c
(defun-single deriv-f1 (x) (exp x))
(defun deriv-f1 (x) (exp x))
(defun deriv-f1-d (x) (exp x))
(defun-single deriv-f2 (x) (if (not (minusp x)) (expt x 3/2) 0.0d0))
(defun deriv-f2 (x) (if (not (minusp x)) (expt x 3/2) 0.0d0))
(defun deriv-f2-d (x) (if (not (minusp x)) (* 3/2 (sqrt x)) 0.0d0))
(defun-single deriv-f3 (x) (if (not (zerop x)) (sin (/ x)) 0.0d0))
(defun deriv-f3 (x) (if (not (zerop x)) (sin (/ x)) 0.0d0))
(defun deriv-f3-d (x) (if (not (zerop x)) (/ (- (cos (/ x))) (expt x 2)) 0.0d0))
(defun-single deriv-f4 (x) (exp (- (expt x 2))))
(defun deriv-f4 (x) (exp (- (expt x 2))))
(defun deriv-f4-d (x) (* -2 x (exp (- (expt x 2)))))
(defun-single deriv-f5 (x) (expt x 2))
(defun deriv-f5 (x) (expt x 2))
(defun deriv-f5-d (x) (* 2 x))
(defun-single deriv-f6 (x) (/ x))
(defun deriv-f6 (x) (/ x))
(defun deriv-f6-d (x) (- (expt x -2)))
(central-derivative deriv-f1 1.0d0 1.0d-4)
(forward-derivative deriv-f1 1.0d0 1.0d-4)
(backward-derivative deriv-f1 1.0d0 1.0d-4)
(central-derivative deriv-f2 0.1d0 1.0d-4)
(forward-derivative deriv-f2 0.1d0 1.0d-4)
(backward-derivative deriv-f2 0.1d0 1.0d-4)
(central-derivative deriv-f3 0.45d0 1.0d-4)
(forward-derivative deriv-f3 0.45d0 1.0d-4)
(backward-derivative deriv-f3 0.45d0 1.0d-4)
(central-derivative deriv-f4 0.5d0 1.0d-4)
(forward-derivative deriv-f4 0.5d0 1.0d-4)
(backward-derivative deriv-f4 0.5d0 1.0d-4)
(central-derivative deriv-f5 0.0d0 1.0d-4)
(forward-derivative deriv-f5 0.0d0 1.0d-4)
(backward-derivative deriv-f5 0.0d0 1.0d-4)
(central-derivative deriv-f6 10.0d0 1.0d-4)
(forward-derivative deriv-f6 10.0d0 1.0d-4)
(backward-derivative deriv-f6 10.0d0 1.0d-4)
(save-test numerical-differentiation
(central-derivative deriv-f1 1.0d0 1.0d-4)
(forward-derivative deriv-f1 1.0d0 1.0d-4)
(backward-derivative deriv-f1 1.0d0 1.0d-4)
(central-derivative deriv-f2 0.1d0 1.0d-4)
(forward-derivative deriv-f2 0.1d0 1.0d-4)
(backward-derivative deriv-f2 0.1d0 1.0d-4)
(central-derivative deriv-f3 0.45d0 1.0d-4)
(forward-derivative deriv-f3 0.45d0 1.0d-4)
(backward-derivative deriv-f3 0.45d0 1.0d-4)
(central-derivative deriv-f4 0.5d0 1.0d-4)
(forward-derivative deriv-f4 0.5d0 1.0d-4)
(backward-derivative deriv-f4 0.5d0 1.0d-4)
(central-derivative deriv-f5 0.0d0 1.0d-4)
(forward-derivative deriv-f5 0.0d0 1.0d-4)
(backward-derivative deriv-f5 0.0d0 1.0d-4)
(central-derivative deriv-f6 10.0d0 1.0d-4)
(forward-derivative deriv-f6 10.0d0 1.0d-4)
(backward-derivative deriv-f6 10.0d0 1.0d-4))
(central-derivative (make-single-function deriv-f1) 1.0d0 1.0d-4)
(forward-derivative (make-single-function deriv-f1) 1.0d0 1.0d-4)
(backward-derivative (make-single-function deriv-f1) 1.0d0 1.0d-4)
(central-derivative (make-single-function deriv-f2) 0.1d0 1.0d-4)
(forward-derivative (make-single-function deriv-f2) 0.1d0 1.0d-4)
(backward-derivative (make-single-function deriv-f2) 0.1d0 1.0d-4)
(central-derivative (make-single-function deriv-f3) 0.45d0 1.0d-4)
(forward-derivative (make-single-function deriv-f3) 0.45d0 1.0d-4)
(backward-derivative (make-single-function deriv-f3) 0.45d0 1.0d-4)
(central-derivative (make-single-function deriv-f4) 0.5d0 1.0d-4)
(forward-derivative (make-single-function deriv-f4) 0.5d0 1.0d-4)
(backward-derivative (make-single-function deriv-f4) 0.5d0 1.0d-4)
(central-derivative (make-single-function deriv-f5) 0.0d0 1.0d-4)
(forward-derivative (make-single-function deriv-f5) 0.0d0 1.0d-4)
(backward-derivative (make-single-function deriv-f5) 0.0d0 1.0d-4)
(central-derivative (make-single-function deriv-f6) 10.0d0 1.0d-4)
(forward-derivative (make-single-function deriv-f6) 10.0d0 1.0d-4)
(backward-derivative(make-single-function deriv-f6) 10.0d0 1.0d-4))
;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2008-12-26 18:54:27EST numerical-integration.lisp>
;; Time-stamp: <2009-01-24 14:48:39EST numerical-integration.lisp>
;; $Id$
;;; To do: QAWS, QAWO, QAWF, more tests
......@@ -206,18 +206,12 @@
;;;; Examples and unit test
;;;;****************************************************************************
(defun-single one-sine (x) (sin x))
#|
;;; Parameters may be defined through the lexical environment:
(let ((mult 2.0d0))
(defun-single two-sine (x) (sin (* mult x))))
|#
(defun one-sine (x) (sin x))
(defparameter *one-sine* (make-single-function one-sine))
(save-test numerical-integration
(integration-qng one-sine 0.0d0 pi)
(integration-qng *one-sine* 0.0d0 pi)
(let ((ws (make-integration-workspace 20)))
(integration-QAG one-sine 0.0d0 pi :gauss15 20 ws))
(integration-QAG *one-sine* 0.0d0 pi :gauss15 20 ws))
(let ((ws (make-integration-workspace 20)))
(integration-QAG one-sine 0.0d0 pi :gauss15 50 ws)))
(integration-QAG *one-sine* 0.0d0 pi :gauss15 50 ws)))
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
;; Time-stamp: <2008-12-26 18:21:11EST chebyshev.lisp>
;; Time-stamp: <2009-01-24 12:56:27EST chebyshev.lisp>
;; $Id$
(in-package :gsl)
......@@ -75,11 +75,13 @@
;;; From Chap. 28.5, except I have set steps = 100 instead of 10000
;;; to keep things sane.
(defun-single chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
(defun chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
(defun chebyshev-table-example ()
(let ((steps 100))
(let ((cheb (make-chebyshev 40 chebyshev-step 0.0d0 1.0d0)))
(let ((cheb
(make-chebyshev
40 (make-single-function chebyshev-step) 0.0d0 1.0d0)))
(dotimes (i steps)
(let ((x (coerce (/ i steps) 'double-float)))
(format t "~&~a ~a ~a ~a"
......@@ -90,7 +92,8 @@
(defun chebyshev-point-example (x)
(check-type x double-float)
(let ((cheb (make-chebyshev 40 chebyshev-step 0.0d0 1.0d0))
(let ((cheb (make-chebyshev
40 (make-single-function chebyshev-step) 0.0d0 1.0d0))
(deriv (make-chebyshev 40))
(integ (make-chebyshev 40)))
(derivative-chebyshev deriv cheb)
......
......@@ -112,6 +112,7 @@ combination for which the following are supported:</p>
version 0.10.0 or newer; callbacks and <code>foreign-funcall</code> must be supported
</li>
<li><a href="http://www.cliki.net/trivial-garbage">trivial-garbage</a>
<li><a href="http://common-lisp.net/project/cl-utilities/">cl-utilities</a>
</li>
<li><a href="http://www.cliki.net/asdf">ASDF</a></li></ul>
......@@ -150,10 +151,11 @@ Then in Lisp, load the system:
<ul>
<li><code>gsll get_git git://repo.or.cz/gsll.git</code></li>
<li><code>trivial-garbage get_darcs http://common-lisp.net/~loliveira/darcs/trivial-garbage</code></li>
<li><code>cl-utilities get_cvs_clnet</code></li>
</ul>
<p>and add
<ul>
<li><code>gsll cffi trivial-garbage</code></li>
<li><code>gsll cffi trivial-garbage cl-utilities</code></li>
</ul>
<p>
to <code>dependencies</code> and execute
......@@ -163,7 +165,7 @@ from within the clbuild directory:
</ul>
<h3>With Debian or Ubuntu</h3>
<ul>
<li><code>sudo aptitude install libgsl0-dev cl-cffi</code></li>
<li><code>sudo aptitude install libgsl0-dev cl-cffi cl-utilities</code></li>
<li><code>git clone git://repo.or.cz/gsll.git</code></li>
<li><code>darcs get http://common-lisp.net/~loliveira/darcs/trivial-garbage</code></li>
<li><code>clc-register-user-package trivial-garbage/trivial-garbage.asd</code></li>
......@@ -358,7 +360,7 @@ and arrays used internally or for function return.
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
<small>
Time-stamp: <2009-01-20 22:50:29EST index.html>
Time-stamp: <2009-01-24 18:40:10EST index.html>
</small>
<!-- hhmts end -->
</div>
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-01-16 11:13:17EST gsll.asd>
;; Time-stamp: <2009-01-24 16:56:13EST gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -9,7 +9,7 @@
:version "0"
:author "Liam M. Healy"
:licence "LLGPL v3, FDL"
:depends-on (cffi trivial-garbage)
:depends-on (cffi trivial-garbage cl-utilities)
:components
((:module init
:components
......
;; Foreign callback functions.
;; Liam Healy
;; Time-stamp: <2009-01-22 22:36:00EST callback.lisp>
;; Time-stamp: <2009-01-24 18:17:53EST callback.lisp>
;; $Id$
(in-package :gsl)
;;; Callback functions are functions which are passed as data; to Lisp
;;; that means they are just functions, but C makes a distinction.
;;; They are needed by several GSL tasks.
;;; Functions that take one double-float and return one double-float
;;; and are defined using the gsl_function structure (see gsl_math.h)
;;; are called scalar functions. They are used in
;;; numerical-integration, numerical-differentiation, chebyshev and
;;; definitions to aid in creating and using them are provided here.
;;; The idea behind the definitions is that the boundary between the CL
;;; and C functions is a narrow as possible; #'defun-single is
;;; macro that allows one to define a function in Lisp and make
;;; use it in these GSL tasks.
;;; They are needed by several GSL tasks. Functions that take one
;;; double-float and return one double-float and are defined using the
;;; gsl_function structure (see gsl_math.h) are called scalar
;;; functions. They are used in numerical-integration,
;;; numerical-differentiation, chebyshev and definitions to aid in
;;; creating and using them are provided here. The idea behind the
;;; definitions is that they absorb us much of the definition overhead
;;; as possible so that the user just defines a CL function and calls
;;; a macro to generate something to pass to the GSL functions.
;;; Other GSL tasks make use of callback functions with different
;;; characteristics. Since they are specific to each of the tasks,
......@@ -33,8 +32,7 @@
;;; right specification.
(export
'(def-single-function undef-cbstruct defun-single
with-c-double with-c-doubles))
'(make-single-function undef-cbstruct with-c-double with-c-doubles))
;;;;****************************************************************************
;;;; Setting slots
......@@ -161,7 +159,7 @@
(defmacro defmcallback
(name &optional (return-type :double) (argument-types :double)
additional-argument-types marray)
additional-argument-types marray (function-name name))
"Define a callback function used by GSL; the GSL function will call
it with an additional `parameters' argument that is ignored. the
argument-types is a single type or list of types of the argument(s)
......@@ -185,7 +183,7 @@
;; CL specials to do the same job.
(declare (ignore params))
,(callback-set-mvb
`(,name
`(,function-name
,@(append
(embedded-clfunc-args atl cbargs marray)
(embedded-clfunc-args aatl cbaddl marray)))
......@@ -209,15 +207,12 @@
This struct is bound to a CL special with the specified name.
This macro can be used whenever a callback is defined and
placed in a struct that has no other functions defined."
(let ((name
;; Bind a CL special under this name to the C structure.
(if (listp functions) (first functions) functions))
(let ((name (make-symbol "CBSTRUCT"))
(fnlist
;; Make a list of (function slot-name ...) for each function.
(if (listp functions) functions (list `,functions 'function))))
`(progn
`(let ((,name (cffi:foreign-alloc ',structure)))
;; Create the C structure and bind CL variable to it.
(defparameter ,name (cffi:foreign-alloc ',structure))
;; Set all the function slots.
,@(loop for (fn slot-name) on fnlist by #'cddr collect
`(set-slot-function ,name ',structure ',slot-name ',fn))
......@@ -227,16 +222,16 @@
,@(loop for slot in additional-slots
collect
`(set-structure-slot
,name ',structure ',(first slot) ,(second slot))))))
,name ',structure ',(first slot) ,(second slot)))
,name)))
(defun undef-cbstruct (name)
(defun undef-cbstruct (object)
"Free foreign callback function. It is not necessary to do this; think
of the memory taken by an unused foreign function as much
less than that used by an unused defun."
(cffi:foreign-free (symbol-value name))
(makunbound name))
(cffi:foreign-free object))
(defmacro def-single-function
(defmacro make-single-function
(name
&optional (return-type :double) (argument-type :double)
(structure 'gsl-function)
......@@ -247,23 +242,17 @@
This struct is bound to a CL special with the specified name.
This macro can be used whenever a callback is defined and
placed in a struct that has no other functions defined."
;; I don't think structure=nil is ever used, but it is permissible
(with-unique-names (cbsymb)
`(progn
(defmcallback ,name ,return-type
(defmcallback ,cbsymb ,return-type
,(if dimensions `((,argument-type ,dimensions)) argument-type)
,(if dimensions-return `((:set ,argument-type ,dimensions-return)))
,marray)
,marray ,name)
,@(when
structure
`((defcbstruct ,name ,structure
,(if dimensions `((dimensions ,dimensions))))))))
;;; Combine a defun and def-single-function in one:
(defmacro defun-single (name arglist &body body)
"Define a function of a scalar double-float argument returning
a double-float in CL and C."
`(progn
(defun ,name ,arglist ,@body)
(def-single-function ,name)))
`((defcbstruct ,cbsymb ,structure
,(if dimensions `((dimensions ,dimensions)))))))))
;;;;****************************************************************************
;;;; Vector of doubles
......
;; Load GSL
;; Liam Healy Sat Mar 4 2006 - 18:53
;; Time-stamp: <2009-01-16 11:06:20EST init.lisp>
;; Time-stamp: <2009-01-24 18:07:19EST init.lisp>
;; $Id$
(defpackage gsll
......@@ -38,3 +38,16 @@
#+sbcl
(eval-when (:compile-toplevel :load-toplevel :execute)
(pushnew :native *features*))
;;; CFFI currently doesn't allow uninterned callback names for CCL;
;;; patch here.
#-openmcl (import 'cl-utilities:with-unique-names :gsll)
#+openmcl
(defmacro gsll::with-unique-names ((&rest bindings) &body body)
`(cl-utilities:with-unique-names
(,@bindings)
(let ,(mapcar (lambda (symb) `(,symb (intern (symbol-name ,symb))))
bindings)
,@body)))
;; Example ODE
;; Liam Healy Sat Sep 29 2007 - 17:49
;; Time-stamp: <2009-01-18 18:29:53EST ode-example.lisp>
;; Time-stamp: <2009-01-24 14:34:20EST ode-example.lisp>
;; $Id$
;;; van der Pol as given in Section 25.5 of the GSL manual. To
......@@ -24,8 +24,6 @@
(- (* -2 mu y0 y1) 1.0d0)
(* -1 mu (- (* y0 y0) 1.0d0))))
(def-ode-functions vanderpol vanderpol-jacobian 2)
(defparameter *max-iter* 2000)
(defun integrate-vanderpol (max-time &optional (step-size 1.0d-6) (print-steps t))
......@@ -38,7 +36,8 @@
(loop (when (or (>= (dcref time) max-time) (> iter *max-iter*))
(return (values (dcref time) dep0 dep1)))
(apply-evolution
evolve control stepper vanderpol
evolve control stepper
(make-ode-functions vanderpol vanderpol-jacobian 2)
time max-time step-size dependent)
(incf iter)
(when print-steps
......
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2009-01-20 22:47:16EST ode-system.lisp>
;; Time-stamp: <2009-01-24 18:23:33EST ode-system.lisp>
;; $Id$
(in-package :gsl)
......@@ -13,9 +13,9 @@
(dimension sizet)
(parameters :pointer))
(export '(def-ode-functions with-ode-integration))
(export '(make-ode-functions with-ode-integration))
(defmacro def-ode-functions (name jacobian dimension)
(defmacro make-ode-functions (name jacobian dimension)
"Setup functions for ODE integrators. The variable name is used as the name of the
The CL functions name and jacobian should be defined previously
with defuns.
......@@ -29,16 +29,17 @@
the derivatives of f with respect to t."
;; set return values
;; Possible future improvements: take/set arrays easily, allow lambdas instead of named functions.
(with-unique-names (solverfn solverdf)
`(progn
(defmcallback
,name :success-failure
(:double (:double ,dimension) (:set :double ,dimension)))
,solverfn :success-failure
(:double (:double ,dimension) (:set :double ,dimension)) nil nil ,name)
(defmcallback
,jacobian :success-failure
,solverdf :success-failure
(:double (:double ,dimension) (:set :double ,(expt dimension 2))
(:set :double ,dimension)))
(defcbstruct (,name function ,jacobian jacobian) ode-system
((dimension ,dimension)))))
(:set :double ,dimension)) nil nil ,jacobian)
(defcbstruct (,solverfn function ,solverdf jacobian) ode-system
((dimension ,dimension))))))
(defmacro with-ode-integration
((time step-size dependent dimensions &optional (stepper '*step-rk8pd*)
......
;; Generic functions for optimization
;; Liam Healy 2009-01-03 12:59:07EST generic.lisp
;; Time-stamp: <2009-01-24 09:11:55EST generic.lisp>
;; Time-stamp: <2009-01-24 17:15:48EST generic.lisp>
;; $Id: $
(in-package :gsl)
......@@ -35,10 +35,10 @@
(dimensions sizet)
(parameters :pointer))
(export 'def-mfunction)
(defmacro def-mfunction (name dimensions)
(export 'make-mfunction)
(defmacro make-mfunction (name dimensions)
"Define a function for multivariate root solving."
`(def-single-function ,name :success-failure :double gsl-mfunction
`(make-single-function ,name :success-failure :double gsl-mfunction
,dimensions))
(cffi:defcstruct gsl-mfunction-fdf
......@@ -51,3 +51,52 @@
(dimensions sizet)
(parameters :pointer))
;;;;****************************************************************************
;;;; Function definition
;;;;****************************************************************************
;;; Solvers that require the function and its derivative use the
;;; following structure.
(cffi:defcstruct gsl-function-fdf
;; See /usr/include/gsl/gsl_math.h
"The definition of a function and its derivative for root finding in GSL."
(function :pointer)
(df :pointer)
(fdf :pointer)
(parameters :pointer))
(export 'make-solver-functions)
(defmacro make-solver-functions (function df fdf &optional dimensions array)
"Setup functions for solvers.
The CL functions name and derivative should be defined previously
with defuns. If dimensions is non-nil (positive fixnum), set multiroot solver
functions. If dimensions is a number, the functions should expect
dimensions scalar (double-float) arguments and return "
(let ((struct (if dimensions 'gsl-mfunction-fdf 'gsl-function-fdf))
(argtype (if dimensions (if array :pointer `((:double ,dimensions))) :double))
(rettype (if dimensions :success-failure :double))
(vecrettype (if array '(:pointer) `((:set :double ,dimensions))))
(matrettype
(if array '(:pointer) `((:set :double ,dimensions ,dimensions)))))
(with-unique-names (solverfn solverdf solverfdf)
`(progn
(defmcallback
,solverfn ,rettype ,argtype ,(when dimensions vecrettype) ,dimensions
,function)
(defmcallback
,solverdf ,rettype ,argtype
,(when dimensions matrettype)
,dimensions ,df)
(defmcallback
,solverfdf
,(if dimensions :success-failure :void)
,argtype
,(if dimensions
(append vecrettype matrettype)
'((:set :double 1) (:set :double 1)))
,dimensions ,fdf)
(defcbstruct (,solverfn function ,solverdf df ,solverfdf fdf)
,struct
,(when dimensions `((dimensions ,dimensions))))))))
;; Multivariate minimization.
;; Liam Healy <Tue Jan 8 2008 - 21:28>
;; Time-stamp: <2009-01-22 21:35:26EST minimization-multi.lisp>
;; Time-stamp: <2009-01-24 16:58:00EST minimization-multi.lisp>
;; $Id$
(in-package :gsl)
......@@ -27,25 +27,28 @@
;;; pointer is a pointer, even though the C definition the functions
;;; they point to have different signatures.
(export 'def-minimization-functions)
(defmacro def-minimization-functions (function dimensions &optional df fdf array)
(export 'make-minimization-functions)
(defmacro make-minimization-functions (function dimensions &optional df fdf array)
"Setup functions for multivariate minimization.
The CL functions name and derivative should be defined previously
with defuns."
(let ((vdec (if array :pointer `((:double ,dimensions))))
(svdec (if array :pointer `(:set :double ,dimensions))))
(with-unique-names (solverfn solverdf solverfdf)
`(progn
(defmcallback ,function :double ,vdec nil t)
,@(when df
`((defmcallback ,df :void ,vdec (,svdec) t)
(defmcallback ,fdf :void ,vdec ((:set :double -1) ,svdec) t)))
(defmcallback ,solverfn :double ,vdec nil t ,function)
,@(when
df
`((defmcallback ,solverdf :void ,vdec (,svdec) t ,df)
(defmcallback
,solverfdf :void ,vdec ((:set :double -1) ,svdec) t ,fdf)))
,(if df
`(defcbstruct (,function function ,df df ,fdf fdf)
`(defcbstruct (,solverfn function ,solverdf df ,solverfdf fdf)
gsl-mfunction-fdf
((dimensions ,dimensions)))
`(defcbstruct (,function function)
`(defcbstruct (,solverfn function)
gsl-mfunction
((dimensions ,dimensions)))))))
((dimensions ,dimensions))))))))
;;;;****************************************************************************
;;;; Initialization
......@@ -285,11 +288,11 @@
(defparameter *parabaloid-center* #(1.0d0 2.0d0))
(defun parabaloid (gsl-vector-pointer)
(defun parabaloid-vector (gsl-vector)
"A parabaloid function of two arguments, given in GSL manual Sec. 35.4.
This version takes a vector-double-float argument."
(let ((x (maref gsl-vector-pointer 0))
(y (maref gsl-vector-pointer 1))
(let ((x (maref gsl-vector 0))
(y (maref gsl-vector 1))
(dp0 (aref *parabaloid-center* 0))
(dp1 (aref *parabaloid-center* 1)))
(+ (* 10 (expt (- x dp0) 2))
......@@ -309,19 +312,19 @@
(defun parabaloid-and-derivative (arguments-gv-pointer derivative-gv-pointer)
(prog1
(parabaloid arguments-gv-pointer)
(parabaloid-vector arguments-gv-pointer)
(parabaloid-derivative
arguments-gv-pointer derivative-gv-pointer)))
(def-minimization-functions
parabaloid 2 parabaloid-derivative parabaloid-and-derivative t)
(defun multimin-example-derivative
(&optional (method *conjugate-fletcher-reeves*) (print-steps t))
(let* ((initial #m(5.0d0 7.0d0))
(minimizer
(make-multi-dimensional-minimizer-fdf
method 2 parabaloid
method 2
(make-minimization-functions
parabaloid-vector 2
parabaloid-derivative parabaloid-and-derivative t)
initial 0.01d0 1.0d-4)))
(loop with status = T
for iter from 0 below 100
......@@ -343,8 +346,9 @@
(values (maref x 0) (maref x 1) (function-value minimizer)))))))
;;; Example without derivatives, same function but now defined using
;;; scalars.
(defun parabaloid-f (x y)
(defun parabaloid-scalar (x y)
"A parabaloid function of two arguments, given in GSL manual Sec. 35.4.
This version takes scalar arguments."
(let ((dp0 (aref *parabaloid-center* 0))
......@@ -353,15 +357,15 @@
(* 20 (expt (- y dp1) 2))
30)))
(def-minimization-functions parabaloid-f 2)
(defun multimin-example-no-derivative
(&optional (method *simplex-nelder-mead*) (print-steps t))
(let ((step-size (make-marray 'double-float :dimensions 2)))
(set-all step-size 1.0d0)
(let ((minimizer
(make-multi-dimensional-minimizer-f
method 2 parabaloid-f #m(5.0d0 7.0d0) step-size)))
method 2
(make-minimization-functions parabaloid-scalar 2)
#m(5.0d0 7.0d0) step-size)))
(loop with status = T and size
for iter from 0 below 100
while status
......
;; Univariate minimization
;; Liam Healy Tue Jan 8 2008 - 21:02
;; Time-stamp: <2009-01-22 18:59:23EST minimization-one.lisp>
;; Time-stamp: <2009-01-24 14:46:49EST minimization-one.lisp>
;; $Id$
(in-package :gsl)
......@@ -175,7 +175,7 @@
;;; This is the example given in Sec. 33.8. The results are different
;;; than given there.
(defun-single minimization-one-fn (x)
(defun minimization-one-fn (x)
(1+ (cos x)))
(defun minimization-one-example
......@@ -184,7 +184,9 @@
(let ((max-iter 100)
(minimizer
(make-one-dimensional-minimizer
minimizer-type minimization-one-fn 2.0d0 0.0d0 6.0d0)))
minimizer-type
(make-single-function minimization-one-fn)
2.0d0 0.0d0 6.0d0)))
(when print-steps
(format
t
......
;; Nonlinear least squares fitting.
;; Liam Healy, 2008-02-09 12:59:16EST nonlinear-least-squares.lisp
;; Time-stamp: <2009-01-24 11:46:27EST nonlinear-least-squares.lisp>
;; Time-stamp: <2009-01-24 16:58:11EST nonlinear-least-squares.lisp>
;; $Id$
(in-package :gsl)
......@@ -85,29 +85,30 @@
(number-of-parameters sizet)
(parameters :pointer))
(export 'def-fitting-functions)
(defmacro def-fitting-functions
(function number-of-parameters &optional df fdf (number-of-observations 0))
"Setup functions for nonlinear least squares fitting.
(export 'make-fitting-functions)