Commit 26de0158 authored by Liam Healy's avatar Liam Healy

Merge in dynamically generated callbacks.

Merge branch 'refactoring-callbacks'
parents 462cecdc bcf61476
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2009-01-25 10:21:38EST monte-carlo.lisp>
;; Time-stamp: <2009-02-15 22:07:40EST monte-carlo.lisp>
;; $Id$
(in-package :gsl)
......@@ -30,13 +30,14 @@
(defmfun monte-carlo-integrate-plain
(function lower-limits upper-limits calls generator state)
"gsl_monte_plain_integrate"
((function :pointer)
((callback :pointer)
((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
((dim0 lower-limits) sizet) (calls sizet)
((mpointer generator) :pointer)
((mpointer state) :pointer)
(result :double) (abserr :double))
:inputs (lower-limits upper-limits)
:callback-struct monte-function
:documentation ; FDL
"Uses the plain Monte Carlo algorithm to integrate the
function f over the hypercubic region defined by the
......@@ -105,13 +106,14 @@
(defmfun monte-carlo-integrate-miser
(function lower-limits upper-limits calls generator state)
"gsl_monte_miser_integrate"
((function :pointer)
((callback :pointer)
((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
((dim0 lower-limits) sizet) (calls sizet)
((mpointer generator) :pointer)
((mpointer state) :pointer)
(result :double) (abserr :double))
:inputs (lower-limits upper-limits)
:callback-struct monte-function
:documentation ; FDL
"Uses the miser Monte Carlo algorithm to integrate the
function f over the hypercubic region defined by the
......@@ -190,13 +192,14 @@
(defmfun monte-carlo-integrate-vegas
(function lower-limits upper-limits calls generator state)
"gsl_monte_vegas_integrate"
((function :pointer)
((callback :pointer)
((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
((dim0 lower-limits) sizet) (calls sizet)
((mpointer generator) :pointer)
((mpointer state) :pointer)
(result :double) (abserr :double))
:inputs (lower-limits upper-limits)
:callback-struct monte-function
:documentation ; FDL
"Uses the vegas Monte Carlo algorithm to integrate the
function f over the dim-dimensional hypercubic region
......@@ -220,48 +223,46 @@
(dimensions sizet)
(parameters :pointer))
(export 'make-monte-carlo-function)
(defmacro make-monte-carlo-function (name-or-lambda &optional (dimensions 1))
"Define the function for use in the Monte Carlo functions. The function
should be specified with either the name of a CL function already defined,
or as a lambda form. In the latter case, the number of dimensions need
not be specified, they will be computed."
`(make-single-function ,name-or-lambda :double :double monte-function
,dimensions nil nil))
(def-make-callbacks monte-carlo (function dimension)
`(defmcallback ,function
:double ((:double ,dimension))
nil
nil
,function))
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
;;; Example from Sec. 23.5
;;; This is a function that occurs in random walk studies.
(defparameter *monte-carlo-cb*
(make-monte-carlo-function
(lambda (x y z)
(* (/ (expt pi 3))
(/ (- 1 (* (cos x) (cos y) (cos z))))))))
(defun mcrw (x y z)
"Example function for Monte Carlo used in random walk studies."
(* (/ (expt pi 3))
(/ (- 1 (* (cos x) (cos y) (cos z))))))
(make-callbacks monte-carlo mcrw 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-cb* lower upper nsamples rng ws)))
(monte-carlo-integrate-plain 'mcrw 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-cb* lower upper nsamples rng ws)))
(monte-carlo-integrate-miser 'mcrw 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-cb* lower upper nsamples rng ws)))
(monte-carlo-integrate-vegas 'mcrw 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: <2009-01-24 14:57:26EST numerical-differentiation.lisp>
;; Time-stamp: <2009-02-15 14:49:44EST numerical-differentiation.lisp>
;; $Id$
(in-package :gsl)
......@@ -12,7 +12,7 @@
(defmfun central-derivative (function x step)
"gsl_deriv_central"
((function :pointer) (x :double) (step :double)
((callback :pointer) (x :double) (step :double)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the numerical derivative of the function
......@@ -32,7 +32,7 @@
(defmfun forward-derivative (function x step)
"gsl_deriv_forward"
((function :pointer) (x :double) (step :double)
((callback :pointer) (x :double) (step :double)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the numerical derivative of the function
......@@ -54,7 +54,7 @@
(defmfun backward-derivative (function x step)
"gsl_deriv_backward"
((function :pointer) (x :double) (step :double)
((callback :pointer) (x :double) (step :double)
(result :double) (abserr :double))
:documentation ; FDL
"Compute the numerical derivative of the function at the point x
......@@ -83,22 +83,29 @@
(defun deriv-f6 (x) (/ x))
(defun deriv-f6-d (x) (- (expt x -2)))
(make-callbacks single-function deriv-f1)
(make-callbacks single-function deriv-f2)
(make-callbacks single-function deriv-f3)
(make-callbacks single-function deriv-f4)
(make-callbacks single-function deriv-f5)
(make-callbacks single-function deriv-f6)
(save-test numerical-differentiation
(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))
(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))
;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2009-01-25 10:16:31EST numerical-integration.lisp>
;; Time-stamp: <2009-02-15 14:43:39EST numerical-integration.lisp>
;; $Id$
;;; To do: QAWS, QAWO, QAWF, more tests
......@@ -19,7 +19,7 @@
;; what these are if they are too large, it will do a minimum number
;; of points anyway.
"gsl_integration_qng"
((function :pointer)
((callback :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double)
(result :double) (abserr :double) (neval sizet))
......@@ -50,13 +50,14 @@
:gauss41 :gauss51 :gauss61)
(defmfun integration-QAG
(function a b method limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a b method limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
;; Set absolute-error and relative-error to 1 because it apparently doesn't matter
;; what these are if they are too large, it will do a minimum number
;; of points anyway.
"gsl_integration_qag"
((function :pointer)
((callback :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double)
(limit sizet) (method integrate-method) ((mpointer workspace) :pointer)
......@@ -77,7 +78,7 @@
On each iteration the adaptive integration strategy bisects the interval
with the largest error estimate. The subintervals and their results are
stored in the memory provided by workspace. The maximum number of
subintervals is given by limit, which may not exceed the allocated
subintervals is given by 'limit, which may not exceed the allocated
size of the workspace.")
;;;;****************************************************************************
......@@ -85,10 +86,11 @@
;;;;****************************************************************************
(defmfun integration-QAGS
(function a b limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a b limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qags"
((function :pointer)
((callback :pointer)
(a :double) (b :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
......@@ -111,9 +113,11 @@
;;;;****************************************************************************
(defmfun integration-QAGP
(function points limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function points limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagp"
((function :pointer)
((callback :pointer)
((mpointer points) :pointer) ((dim0 points) sizet)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
......@@ -135,9 +139,11 @@
;;;;****************************************************************************
(defmfun integration-QAGi
(function limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagi"
((function :pointer)
((callback :pointer)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
:documentation ; FDL
......@@ -152,9 +158,11 @@
this case a lower-order rule is more efficient.")
(defmfun integration-QAGiu
(function a limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagiu"
((function :pointer) (a :double)
((callback :pointer) (a :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
:documentation ; FDL
......@@ -165,9 +173,11 @@
and then integrated using the QAGS algorithm.")
(defmfun integration-QAGil
(function b limit workspace &optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function b limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qagil"
((function :pointer) (b :double)
((callback :pointer) (b :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
:documentation ; FDL
......@@ -182,10 +192,11 @@
;;;;****************************************************************************
(defmfun integration-QAWC
(function a b c limit workspace
&optional (absolute-error 1.0d0) (relative-error 1.0d0))
(function a b c limit
&optional (absolute-error 1.0d0) (relative-error 1.0d0)
(workspace (make-integration-workspace limit)))
"gsl_integration_qawc"
((function :pointer)
((callback :pointer)
(a :double) (b :double) (c :double)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer) (result :double) (abserr :double))
......@@ -202,17 +213,14 @@
the singularity. Further away from the singularity the algorithm
uses an ordinary 15-point Gauss-Kronrod integration rule.")
;;;;****************************************************************************
;;;; Examples and unit test
;;;;****************************************************************************
(defun one-sine (x) (sin x))
(defparameter *one-sine* (make-single-function one-sine))
(make-callbacks single-function one-sine)
(save-test numerical-integration
(integration-qng *one-sine* 0.0d0 pi)
(let ((ws (make-integration-workspace 20)))
(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-qng 'one-sine 0.0d0 pi)
(integration-QAG 'one-sine 0.0d0 pi :gauss15 20)
(integration-QAG 'one-sine 0.0d0 pi :gauss21 40))
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
;; Time-stamp: <2009-01-28 19:05:30EST chebyshev.lisp>
;; Time-stamp: <2009-02-15 11:23:40EST chebyshev.lisp>
;; $Id$
(in-package :gsl)
......@@ -16,9 +16,12 @@
"Chebyshev series"
:documentation ; FDL
"Make a Chebyshev series of specified order."
:superclasses (callback-included)
:ci-class-slots (gsl-function nil (function))
:initialize-suffix "init"
:initialize-args
((function :pointer) (lower-limit :double) (upper-limit :double)))
((callback :pointer) (lower-limit :double) (upper-limit :double))
:singular (function))
;;;;****************************************************************************
;;;; Chebyshev series evaluation
......@@ -55,7 +58,7 @@
(defmfun derivative-chebyshev (derivative chebyshev)
"gsl_cheb_calc_deriv"
((derivative :pointer) ((mpointer chebyshev) :pointer))
(((mpointer derivative) :pointer) ((mpointer chebyshev) :pointer))
:documentation ; FDL
"Compute the derivative of the Chebyshev series, storing
the derivative coefficients in the previously allocated series.
......@@ -63,7 +66,7 @@
(defmfun integral-chebyshev (integral chebyshev)
"gsl_cheb_calc_integ"
((integral :pointer) ((mpointer chebyshev) :pointer))
(((mpointer integral) :pointer) ((mpointer chebyshev) :pointer))
:documentation ; FDL
"Compute the integral of the Chebyshev series, storing
the integral coefficients in the previously allocated series.
......@@ -80,11 +83,11 @@
(defun chebyshev-step (x) (if (< x 0.5d0) 0.25d0 0.75d0))
(make-callbacks single-function chebyshev-step)
(defun chebyshev-table-example ()
(let ((steps 100))
(let ((cheb
(make-chebyshev
40 (make-single-function chebyshev-step) 0.0d0 1.0d0)))
(let ((cheb (make-chebyshev 40 'chebyshev-step 0.0d0 1.0d0)))
(dotimes (i steps)
(let ((x (coerce (/ i steps) 'double-float)))
(format t "~&~a ~a ~a ~a"
......@@ -95,8 +98,7 @@
(defun chebyshev-point-example (x)
(check-type x double-float)
(let ((cheb (make-chebyshev
40 (make-single-function chebyshev-step) 0.0d0 1.0d0))
(let ((cheb (make-chebyshev 40 'chebyshev-step 0.0d0 1.0d0))
(deriv (make-chebyshev 40))
(integ (make-chebyshev 40)))
(derivative-chebyshev deriv cheb)
......
......@@ -299,29 +299,23 @@ area of memory.</p>
<h3>Passing functions</h3>
<p>
Functions that are passed to GSL functions (known as <i>callbacks</i>
in C) are defined with one of the following macros in GSLL, depending
on the application:
<ul>
<li><code>make-single-function</code>, used by numerical
integration, numerical differentiation, Chebyshev, one-dimensional
root solver without derivatives, one-dimensional minimizer;</li>
<li><code>make-monte-carlo-function</code>, used by Monte Carlo;</li>
<li><code>make-solver-functions</code>, used by the root
solvers;</li>
<li><code>make-fitting-functions</code>, used by nonlinear least squares;</li>
<li><code>make-minimization-functions</code>, used by
the multi-dimensional minimizer;</li>
<li><code>make-ode-functions</code>, used by ordinary differential
equation solvers;</li>
</ul>
in C) are defined by writing an ordinary CL function or functions and
then calling <code>make-callbacks</code> with arguments that
depend on the application. Functions or objects that require a
callback will accept the symbol representing the first or only function
defined in the <code>make-callbacks</code>. For example, after
defining an ordinary CL function <code>mcrw</code>,
<pre>
(make-callbacks monte-carlo mcrw 3)
(monte-carlo-integrate-plain 'mcrw lower upper nsamples rng ws)
</pre>
<p>
The function may be specified as the name of a CL function, or as a
lambda form. If the latter, the dimensions need not be specified in
<code>make-single-function</code>, <code>make-monte-carlo-function</code>,
or <code>make-solver-functions</code> for derivative-free solvers;
they will be computed from the arglist. Each of these macros will
return a foreign pointer that can be passed to the appropriate GSL
function that requires a function.
In some cases, the CL function can be written to accept and return
either scalars or arrays. See examples for numerical
integration, numerical differentiation, Chebyshev, root solvers,
minimizers, Monte Carlo, nonlinear least squares, and ordinary
differential equation solvers.
</p>
<h3>GSL objects</h3>
<p>
......@@ -388,7 +382,7 @@ and arrays used internally or for function return.
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
<small>
Time-stamp: <2009-01-25 21:38:05EST index.html>
Time-stamp: <2009-02-15 22:17:52EST index.html>
</small>
<!-- hhmts end -->
</div>
......
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2009-02-11 21:17:38EST gsll.asd>
;; Time-stamp: <2009-02-15 22:56:09EST gsll.asd>
;; $Id$
(asdf:defsystem "gsll"
......@@ -16,6 +16,7 @@
((:file "init")
(:file "conditions" :depends-on (init))
(:file "mobject" :depends-on (init))
(:file "callback" :depends-on (init mobject))
(:file "types" :depends-on (init))
(:file "complex-types" :depends-on (types))
(:file "element-types" :depends-on (init complex-types))
......@@ -24,8 +25,7 @@
:depends-on (init conditions element-types number-conversion))
(:file "defmfun" :depends-on (init element-types interface))
(:file "defmfun-array" :depends-on (defmfun))
(:file "defmfun-single" :depends-on (defmfun mobject))
(:file "callback" :depends-on (init))
(:file "defmfun-single" :depends-on (defmfun mobject callback))
(:file "generate-examples" :depends-on (init))))
(:module floating-point
:depends-on (init)
......
This diff is collapsed.
;; Helpers that define a single GSL function interface
;; Liam Healy 2009-01-07 22:02:20EST defmfun-single.lisp
;; Time-stamp: <2009-01-31 21:19:28EST defmfun-single.lisp>
;; Time-stamp: <2009-02-15 18:56:52EST defmfun-single.lisp>
;; $Id: $
(in-package :gsl)
......@@ -100,6 +100,9 @@
(cons
'values
(append before after
(when (callback-arg-p
(variables-used-in-c-arguments c-arguments))
'(function)) ; hardwired name to use for callback
(let ((auxstart (position '&aux arglist)))
;; &aux bindings are checked
(when auxstart
......@@ -108,7 +111,7 @@
(mapcar 'rest (subseq arglist (1+ auxstart))))))))))))
,@(when documentation (list documentation))
,(funcall body-maker name arglist gsl-name c-arguments key-args))
`(,definition
`(,definition
,@(when (and name (not (defgeneric-method-p name)))
(list name))
,@(when qualifier (list qualifier))
......@@ -164,7 +167,8 @@
allocated))
(clret (or ; better as a symbol macro
(substitute cret-name :c-return return)
(mapcan #'cl-convert-form allocated-decl)
(mapcan #'cl-convert-form
(callback-remove-arg allocated-decl 'st-symbol))
outputs
(unless (eq c-return :void)
(list cret-name)))))
......@@ -175,6 +179,7 @@
(mapcar #'wfo-declare allocated-decl)
'cffi:with-foreign-objects
`(,@before
,(callback-set-slots allocated callback-struct function)
(let ((,cret-name
(cffi:foreign-funcall
,gsl-name
......@@ -208,7 +213,8 @@
,@after
(values
,@(defmfun-return
c-return cret-name clret allocated return return-supplied-p
c-return cret-name clret allocated
return return-supplied-p
enumeration outputs)))))))))
(defun defmfun-return
......
;; Macro for defining GSL functions.
;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
;; Time-stamp: <2009-01-12 21:50:07EST defmfun.lisp>
;; Time-stamp: <2009-02-15 18:51:03EST defmfun.lisp>
;; $Id$
(in-package :gsl)
......@@ -74,12 +74,12 @@
(definition :function)
(export (not (member definition (list :method :methods))))
documentation inputs outputs before after enumeration qualifier
gsl-version switch)
gsl-version switch (callback-struct 'gsl-function))
,key-args
(declare (ignorable c-return return definition element-types
index export documentation inputs outputs
before after enumeration qualifier
gsl-version switch)
gsl-version switch callback-struct)
(special indexed-functions))
,@body))
......
;; Macros to interface GSL functions, including definitions necessary for defmfun.
;; Liam Healy
;; Time-stamp: <2009-01-25 10:30:39EST interface.lisp>
;; Time-stamp: <2009-02-15 12:45:17EST interface.lisp>
;; $Id$
(in-package :gsl)
......@@ -61,9 +61,11 @@
(defun wfo-declare (d)
`(,(st-symbol d)
,@(if (st-arrayp d)
`(',(st-eltype d) ,(st-dim d))
`(',(st-type d)))))
,@(if (st-arrayp d)
`(',(st-eltype d) ,(st-dim d))
(if (eq (st-symbol d) +callback-argument-name+)
'('gsl-function)
`(',(st-type d))))))
;;;;****************************************************************************
;;;; Checking results from GSL functions
......
This diff is collapsed.
;; Evolution functions for ODE integration.
;; Liam Healy, Sun Sep 30 2007 - 14:31
;; Time-stamp: <2009-01-26 21:25:15EST evolution.lisp>
;; Time-stamp: <2009-02-15 08:58:47EST evolution.lisp>
;; $Id$
(in-package :gsl)
......@@ -8,17 +8,16 @@
(defmobject ode-evolution "gsl_odeiv_evolve"
((dimensions sizet))
"evolution for ordinary differential equations"
:documentation
"Make an object to advance the ODE solution."
:documentation "Make an object to advance the ODE solution."
:initialize-suffix "reset"
:initialize-args nil)
(defmfun apply-evolution
(evolve control step dydt time max-time step-size y)
(evolution time y step-size control stepper max-time)
"gsl_odeiv_evolve_apply"
(((mpointer evolve) :pointer) ((mpointer control) :pointer)
((mpointer step) :pointer)
(dydt :pointer) ((c-pointer time) :pointer)
(((mpointer evolution) :pointer) ((mpointer control) :pointer)
((mpointer stepper) :pointer)
((callback-struct stepper) :pointer) ((c-pointer time) :pointer)
(max-time :double)
((c-pointer step-size) :pointer) ((c-pointer y) :pointer))
:inputs (time step-size y)
......
;; Example ODE
;; Liam Healy Sat Sep 29 2007 - 17:49
;; Time-stamp: <2009-01-26 21:45:43EST ode-example.lisp>
;; Time-stamp: <2009-02-15 09:09:55EST ode-example.lisp>
;; $Id$
;;; van der Pol as given in Section 25.5 of the GSL manual. To
......@@ -26,9 +26,7 @@
(defparameter *max-iter* 2000)
#+callback-toplevel-only
(defparameter *vanderpol-cb*
(make-ode-functions vanderpol vanderpol-jacobian 2))
(make-callbacks ode-stepper vanderpol vanderpol-jacobian 2)
(defun integrate-vanderpol
(max-time &optional (step-size 1.0d-6) (stepper *step-rk8pd*) (print-steps t))
......@@ -37,15 +35,13 @@
(let ((mu 10.0d0) (initial-time 0.0d0) (iter 0))
(declare (special mu))
(with-ode-integration
(#-callback-toplevel-only
(make-ode-functions vanderpol vanderpol-jacobian 2)
#+callback-toplevel-only *vanderpol-cb*
((vanderpol vanderpol-jacobian)
time step max-time (dep0 dep1) 2 stepper)
(setf dep0 1.0d0 dep1 0.0d0 step step-size time initial-time)
(loop
(when (or (>= time max-time) (> iter *max-iter*))
(return (values iter time dep0 dep1)))
make-next-step
(next-step)
(incf iter)
(when print-steps
(format t "~12,6f~10t~12,6f~24t~12,6f~&" time dep0 dep1))))))
......
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2009-01-25 19:01:57EST ode-system.lisp>
;; Time-stamp: <2009-02-15 09:09:55EST ode-system.lisp>
;; $Id$
(in-package :gsl)
......@@ -13,33 +13,7 @@
(dimension sizet)
(parameters :pointer))
(export '(make-ode-functions with-ode-integration))
(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.
The function and Jacobian arguments are the same: time and
dimension dependent variables as scalars.
The function returns dimension multiple values corresponding to
the derivatives of the dependent variables, f(t).
The Jacobian returns dimension^2 values corresponding to the
partial derivatives of each of the functions f with respect to each
of the dependent variables y, and dimension values corresponding to
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
,solverfn :success-failure
(:double (:double ,dimension) (:set :double ,dimension)) nil nil ,name)
(defmcallback
,solverdf :success-failure
(:double (:double ,dimension) (:set :double ,(expt dimension 2))
(:set :double ,dimension)) nil nil ,jacobian)
(defcbstruct (,solverfn function ,solverdf jacobian) ode-system
((dimension ,dimension))))))
(export '(with-ode-integration))
(defmacro with-ode-integration
((function time step-size max-time
......@@ -50,7 +24,7 @@
(let ((dep (make-symbol "DEP"))
(ctime (make-symbol "CTIME"))
(cstep (make-symbol "CSTEP")))
`(let ((stepperobj (make-ode-stepper ,stepper ,dimensions))
`(let ((stepperobj (make-ode-stepper ,stepper ,dimensions ',function))
(control (make-y-control ,absolute-error ,relative-error))
(evolve (make-ode-evolution ,dimensions))
(,dep
......@@ -62,10 +36,8 @@
(,step-size (maref ,cstep 0))
,@(loop for symb in dependent
for i from 0
collect `(,symb (maref ,dep ,i)))