Commit 4e83fcc4 authored by Liam Healy's avatar Liam Healy

ODE control objects

ODE control is now done with GSL objects, including a superclass
'ode-control.
parent cd6a85f8
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2009-01-25 10:35:08EST mobject.lisp>
;; Time-stamp: <2009-01-25 17:27:46EST mobject.lisp>
;; $Id$
;;; GSL objects are represented in GSLL as and instance of a 'mobject.
......@@ -22,15 +22,14 @@
(defmacro defmobject
(class prefix allocation-args description
&key documentation initialize-suffix initialize-args arglists-function
inputs gsl-version)
&key documentation initialize-suffix initialize-name initialize-args
arglists-function inputs gsl-version allocator allocate-inputs freer
(superclasses '(mobject)))
"Define the class, the allocate, initialize-instance and
reinitialize-instance methods, and the make-* function for the GSL object."
;; If prefix is a list, the first is the actual prefix, and the
;; second is the name of the allocator. I'm looking at you,
;; discrete-random. Grrr.
;; Argument 'initialize-suffix: string appended to prefix for form GSL function name
;; or a list of such a string and the c-return argument.
;; Argument 'initialize-suffix: string appended to prefix to form
;; GSL function name or a list of such a string and the c-return
;; argument.
(let* ((settingp (make-symbol "SETTINGP"))
(arglists
(when arglists-function
......@@ -39,19 +38,20 @@
(cl-alloc-args (variables-used-in-c-arguments allocation-args))
(cl-initialize-args
(variables-used-in-c-arguments initialize-args))
(realprefix (if (listp prefix) (first prefix) prefix)))
(initializerp (or initialize-name initialize-suffix))
(initargs ; arguments that are exclusively for reinitialize-instance
(remove-if (lambda (s) (member s cl-alloc-args)) cl-initialize-args)))
(if (have-at-least-gsl-version gsl-version)
`(progn
(defclass ,class (mobject)
(defclass ,class ,superclasses
()
(:documentation
,(format nil "The GSL representation of the ~a." description)))
(defmfun allocate ((object ,class) &key ,@cl-alloc-args)
,(if (listp prefix) (second prefix) (format nil "~a_alloc" prefix))
,(or allocator (format nil "~a_alloc" prefix))
,allocation-args
,@(if (and (listp prefix) (third prefix))
`(:inputs ,(third prefix)))
,@(if allocate-inputs `(:inputs ,allocate-inputs))
:definition :method
:c-return :pointer
:index ,maker)
......@@ -68,20 +68,21 @@
(tg:finalize object
(lambda ()
(cffi:foreign-funcall
,(format nil "~a_free" realprefix)
,(or freer (format nil "~a_free" prefix))
:pointer mpointer :void))))
,@(when initialize-suffix
,@(when initializerp
`((defmfun reinitialize-instance
((object ,class) &key ,@cl-initialize-args)
,(format nil "~a_~a" realprefix
(if (listp initialize-suffix)
(first initialize-suffix)
initialize-suffix))
,(or initialize-name
(format nil "~a_~a" prefix
(if (listp initialize-suffix)
(first initialize-suffix)
initialize-suffix)))
(((mpointer object) :pointer) ,@initialize-args)
:definition :method
:qualifier :after
,@(when (listp initialize-suffix)
,@(when (and initialize-suffix (listp initialize-suffix))
`(:c-return ,(second initialize-suffix)))
:return (object)
,@(when inputs `(:inputs ,inputs))
......@@ -93,11 +94,11 @@
,(if arglists
(first arglists)
`(,@cl-alloc-args
,@(when initialize-args
,@(when initargs
(append
(list '&optional
(list (first cl-initialize-args) nil settingp))
(rest cl-initialize-args)))))
(list
'&optional (list (first initargs) nil settingp))
(rest initargs)))))
,(format
nil "Create the GSL object representing a ~a (class ~a).~@[~&~a~]"
description class documentation)
......@@ -108,20 +109,26 @@
(second arglists)
(symbol-keyword-symbol cl-alloc-args)))))
;; There is an initialization step
,@(when initialize-suffix
,@(when initializerp
(if initialize-args ; with arguments
`((when ,settingp
(reinitialize-instance
object
,@(if arglists
(third arglists)
(symbol-keyword-symbol cl-initialize-args)))))
(if initargs
`((when ,settingp
(reinitialize-instance
object
,@(if arglists
(third arglists)
(symbol-keyword-symbol cl-initialize-args)))))
`((reinitialize-instance
object
,@(if arglists
(third arglists)
(symbol-keyword-symbol cl-initialize-args)))))
'((reinitialize-instance object)))) ; without arguments
object)))
`(defun ,maker (&rest args)
(declare (ignore args))
(error 'obsolete-gsl-version
:name ',class :gsl-name ,realprefix
:name ',class :gsl-name ,prefix
:gsl-version ',gsl-version)))))
(defun symbol-keyword-symbol (symbol)
......
;; Adaptive step-size control
;; Liam Healy 2008-02-17 17:30:04EST control.lisp
;; Time-stamp: <2009-01-25 10:01:56EST control.lisp>
;; Time-stamp: <2009-01-25 17:09:29EST control.lisp>
;; $Id$
(in-package :gsl)
#|
(defmobject standard-control ("gsl_odeiv_control_standard" "new")
(defclass ode-control (mobject)
()
(:documentation
"Objects used for control of ordinary differential equation integration."))
(defmobject standard-control
"gsl_odeiv_control_standard"
((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double))
"standard control for ordinary differential equations"
"")
|#
(defmfun new-standard-control (absolute-error relative-error y-scaling dydt-scaling)
"gsl_odeiv_control_standard_new"
:superclasses (ode-control)
:initialize-name "gsl_odeiv_control_init"
:initialize-args
((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double))
:c-return (ptr :pointer)
:return (ptr)
:allocator "gsl_odeiv_control_standard_new"
:freer "gsl_odeiv_control_free"
:documentation ; FDL
"The standard control object is a four parameter heuristic based on
absolute and relative errors absolute-error and relative-error, and
......@@ -46,11 +49,17 @@
uncontrolled changes in the stepsize, the overall scaling factor is
limited to the range 1/5 to 5.")
(defmfun new-y-control (absolute-error relative-error)
"gsl_odeiv_control_y_new"
(defmobject y-control
"gsl_odeiv_control_y"
((absolute-error :double) (relative-error :double))
:c-return (ptr :pointer)
:return (ptr)
"y control for ordinary differential equations"
:superclasses (ode-control)
:initialize-name "gsl_odeiv_control_init"
:initialize-args
((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double))
:allocator "gsl_odeiv_control_y_new"
:freer "gsl_odeiv_control_free"
:documentation ; FDL
"Create a new control object which will keep the local
error on each step within an absolute error of absolute-error and
......@@ -58,11 +67,17 @@
This is equivalent to the standard control object with y-scaling=1 and
dydt-scaling=0.")
(defmfun new-yp-control (absolute-error relative-error)
"gsl_odeiv_control_yp_new"
(defmobject yp-control
"gsl_odeiv_control_yp"
((absolute-error :double) (relative-error :double))
:c-return (ptr :pointer)
:return (ptr)
"yp control for ordinary differential equations"
:superclasses (ode-control)
:initialize-name "gsl_odeiv_control_init"
:initialize-args
((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double))
:allocator "gsl_odeiv_control_yp_new"
:freer "gsl_odeiv_control_free"
:documentation ; FDL
"Create a new control object which will keep the local
error on each step within an absolute error of absolute-error and
......@@ -70,14 +85,19 @@
solution y'_i(t). This is equivalent to the standard control
object with y-scaling=0 and dydt-scaling=1.")
(defmfun new-scaled-control
(absolute-error relative-error y-scaling dydt-scaling absolute-scale dimension)
"gsl_odeiv_control_scaled_new"
(defmobject scaled-control
"gsl_odeiv_control_scaled"
((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double)
(absolute-scale :pointer) (dimension sizet))
:c-return (ptr :pointer)
:return (ptr)
"scaled control for ordinary differential equations"
:superclasses (ode-control)
:initialize-name "gsl_odeiv_control_init"
:initialize-args
((absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double))
:allocator "gsl_odeiv_control_scaled_new"
:freer "gsl_odeiv_control_free"
:documentation ; FDL
"Create a new control object which uses the same algorithm
as #'new-standard-control but with an absolute error
......@@ -98,30 +118,12 @@
defining new types of control functions. For most purposes the standard
control functions described above should be sufficient.")
(defmfun initialize-control
(control absolute-error relative-error y-scaling dydt-scaling)
"gsl_odeiv_control_init"
((control :pointer) (absolute-error :double) (relative-error :double)
(y-scaling :double) (dydt-scaling :double))
:documentation ; FDL
"Initialize the control function control with the
parameters absolute-error, relative-error,
y-scaling (scaling factor for y) and dydt-scaling (scaling
factor for derivatives).")
(defmfun free-control (control)
"gsl_odeiv_control_free"
((control :pointer))
:c-return :void
:documentation ; FDL
"Free all the memory associated with the control function control.")
(cffi:defcenum step-size-adjustment
(:step-size-decreased -1) :step-size-unchanged :step-size-increased)
(defmfun adjust-stepsize (control stepper current-y y-error dydt step-size)
"gsl_odeiv_control_hadjust"
((control :pointer) (stepper :pointer) (current-y :pointer)
(((mpointer control) :pointer) (stepper :pointer) (current-y :pointer)
(y-error :pointer) (dydt :pointer) (step-size :pointer))
:c-return :enumerate
:enumeration step-size-adjustment
......@@ -138,9 +140,11 @@
step-size which satisfies the user-specified accuracy requirements for
the current point.")
(defmfun control-name (control)
;; Generalize this to all ODE control
(defmfun name ((control ode-control))
"gsl_odeiv_control_name"
((control :pointer))
(((mpointer control) :pointer))
:definition :method
:c-return :string
:documentation ; FDL
"The name of the control function.")
;; Evolution functions for ODE integration.
;; Liam Healy, Sun Sep 30 2007 - 14:31
;; Time-stamp: <2009-01-25 10:00:41EST evolution.lisp>
;; Time-stamp: <2009-01-25 17:11:41EST evolution.lisp>
;; $Id$
(in-package :gsl)
......@@ -16,7 +16,8 @@
(defmfun apply-evolution
(evolve control step dydt time max-time step-size y)
"gsl_odeiv_evolve_apply"
(((mpointer evolve) :pointer) (control :pointer) ((mpointer step) :pointer)
(((mpointer evolve) :pointer) ((mpointer control) :pointer)
((mpointer step) :pointer)
(dydt :pointer) (time :pointer) (max-time :double)
(step-size :pointer) (y :pointer))
:documentation ; FDL
......
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2009-01-24 18:23:33EST ode-system.lisp>
;; Time-stamp: <2009-01-25 17:10:32EST ode-system.lisp>
;; $Id$
(in-package :gsl)
......@@ -53,19 +53,17 @@
(let ((ctime (make-symbol "CTIME"))
(cstep (make-symbol "CSTEP")))
`(let ((stepper (make-ode-stepper ,stepper ,dimensions))
(control (new-y-control ,absolute-error ,relative-error))
(control (make-y-control ,absolute-error ,relative-error))
(evolve (make-ode-evolution ,dimensions)))
(unwind-protect
(cffi:with-foreign-objects
((,(if (listp dependent) (first dependent) dependent)
:double ,dimensions) (,ctime :double) (,cstep :double))
(setf
(dcref ,cstep) ,step-size
(dcref ,ctime) ,time)
,(if (listp dependent)
`(with-c-double ,dependent
(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body))
`(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body)))
(free-control control)))))
(cffi:with-foreign-objects
((,(if (listp dependent) (first dependent) dependent)
:double ,dimensions) (,ctime :double) (,cstep :double))
(setf
(dcref ,cstep) ,step-size
(dcref ,ctime) ,time)
,(if (listp dependent)
`(with-c-double ,dependent
(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body))
`(symbol-macrolet ((,time ,ctime) (,step-size ,cstep))
,@body))))))
;; Discrete random variables
;; Liam Healy, Sat Nov 11 2006 - 21:51
;; Time-stamp: <2009-01-25 09:54:54EST discrete.lisp>
;; Time-stamp: <2009-01-25 11:17:05EST discrete.lisp>
;; $Id$
(in-package :gsl)
;;; /usr/include/gsl/gsl_randist.h
(defmobject discrete-random
("gsl_ran_discrete" "gsl_ran_discrete_preproc" (probabilities))
(defmobject discrete-random "gsl_ran_discrete"
(((dim0 probabilities) sizet) ((c-pointer probabilities) :pointer))
"lookup table for the discrete random number generator"
:allocator "gsl_ran_discrete_preproc"
:allocate-inputs (probabilities)
:documentation ; FDL
"Make a structure that contains the lookup
table for the discrete random number generator. The array probabilities contains
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment