Commit b6561c8e authored by Liam Healy's avatar Liam Healy

ODE evolution ported to callbacks-included

Mobject 'ode-evolution is now made as a subclass of 'callbacks-included-cl, a
new class that is a subclass of 'callbacks-included, but which stores the
cbstruct pointer in a CL slot rather than in a foreign struct.  This example is
correct,
(integrate-vanderpol 1.0d0 1.d-4 *step-rk4* nil)
40
1.0d0
-1.4568569264026898d0
-11.547449151779395d0
parent 59d0892f
;; Foreign callback functions.
;; Liam Healy
;; Time-stamp: <2009-02-14 11:50:31EST callback.lisp>
;; Time-stamp: <2009-02-14 19:22:53EST callback.lisp>
;; $Id$
(in-package :gsl)
......@@ -226,6 +226,7 @@
,name)))
(defun make-cbstruct (struct slots-values &rest function-slotnames)
"Make the callback structure."
(let ((cbstruct (cffi:foreign-alloc struct)))
(loop for (slot-name function) on function-slotnames by #'cddr
do (set-slot-function cbstruct struct slot-name function))
......@@ -235,6 +236,15 @@
do (set-structure-slot cbstruct struct slot-name value)))
cbstruct))
(defun make-cbstruct-object (object)
"Make the callback structure based on the mobject definition."
(apply
'make-cbstruct
(cbstruct-name object)
(when (dimension-names object)
(mapcan 'list (dimension-names object) (dimensions object)))
(mapcan 'list (callback-labels object) (functions object))))
(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
......@@ -309,11 +319,18 @@
(:documentation
"A mobject that includes a callback function or functions to GSL."))
(defclass callback-included-cl (callback-included)
((callback :initarg :callback :reader callback-struct))
(:documentation
"A mobject that includes a callback function or functions, which
the pointer to the callback structure is stored in a CL class
slot."))
(defmacro def-ci-subclass
(class-name documentation
(class-name superclasses documentation
cbstruct-name array-type callback-labels
&optional (dimension-names '(dimensions)))
`(defclass ,class-name (callback-included)
`(defclass ,class-name ,superclasses
((cbstruct-name :initform ',cbstruct-name :allocation :class)
(array-type :initform ',array-type :allocation :class)
(callback-labels :initform ',callback-labels :allocation :class)
......@@ -321,9 +338,9 @@
(:documentation ,documentation)))
(defmacro def-ci-subclass-1d
(class-name documentation
(class-name superclasses documentation
cbstruct-name array-type callback-labels)
`(defclass ,class-name (callback-included)
`(defclass ,class-name ,superclasses
((cbstruct-name :initform ',cbstruct-name :allocation :class)
(array-type :initform ',array-type :allocation :class)
(callback-labels :initform ',callback-labels :allocation :class)
......
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2009-02-14 11:52:53EST mobject.lisp>
;; Time-stamp: <2009-02-14 18:48:19EST mobject.lisp>
;;; GSL objects are represented in GSLL as and instance of a 'mobject.
;;; The macro demobject takes care of defining the appropriate
......@@ -26,14 +26,15 @@
&key documentation initialize-suffix initialize-name initialize-args
arglists-function inputs gsl-version allocator allocate-inputs freer
(superclasses '(mobject))
ci-class-slots singular)
class-slots-instance ci-class-slots singular)
"Define the class, the allocate, initialize-instance and
reinitialize-instance methods, and the make-* function for the GSL object."
;; 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"))
(callbackp (member 'callback-included superclasses))
(callbackp
(member-if (lambda (cl) (subtypep cl 'callback-included)) superclasses))
(arglists
(when arglists-function
(funcall (coerce arglists-function 'function) settingp)))
......@@ -51,6 +52,7 @@
`(,(if (member 'dimensions cl-alloc-args)
'def-ci-subclass 'def-ci-subclass-1d)
,class
,superclasses
,(format nil "The GSL representation of the ~a." description)
,@ci-class-slots)
`(defclass ,class ,superclasses
......@@ -72,12 +74,14 @@
(make-reinitialize-instance
class cl-initialize-args initialize-name prefix
initialize-suffix initialize-args inputs
callbackp (member 'dimensions cl-alloc-args)))
(and callbackp
(not (member +callback-argument-name+
class-slots-instance)))))
(export ',maker)
,(mobject-maker
maker arglists initargs class cl-alloc-args cl-initialize-args
description documentation initialize-args initializerp settingp
singular))
singular class-slots-instance))
`(progn
(export ',maker)
(defun ,maker (&rest args)
......@@ -103,10 +107,19 @@
,(or freer (format nil "~a_free" prefix))
:pointer mpointer :void)))))
(defun make-cbstruct-object (object)
"Make the callback structure based on the mobject definition."
(apply
'make-cbstruct
(cbstruct-name object)
(when (dimension-names object)
(mapcan 'list (dimension-names object) (dimensions object)))
(mapcan 'list (callback-labels object) (functions object))))
(defun make-reinitialize-instance
(class cl-initialize-args initialize-name prefix
initialize-suffix initialize-args inputs
callback dimensionsp)
callback)
"Expand the reinitialize-instance form. In the GSL arglist, the callback
structure pointer should be named the value of +callback-argument-name+."
(let ((cbstruct (make-symbol "CBSTRUCT")))
......@@ -114,17 +127,8 @@
((object ,class)
&key
,@cl-initialize-args
,@(if
callback
`(&aux
(,cbstruct
(apply
'make-cbstruct
(cbstruct-name object)
,(when
dimensionsp
`(mapcan 'list (dimension-names object) (dimensions object)))
(mapcan 'list (callback-labels object) (functions object)))))))
,@(when callback
`(&aux (,cbstruct (make-cbstruct-object object)))))
,(or initialize-name
(format nil "~a_~a" prefix
(if (listp initialize-suffix)
......@@ -152,7 +156,7 @@
(defun mobject-maker
(maker arglists initargs class cl-alloc-args cl-initialize-args
description documentation initialize-args initializerp settingp
singular)
singular class-slots-instance)
"Make the defun form that makes the mobject."
`(defun ,maker
,(if arglists
......@@ -160,6 +164,7 @@
(singularize
singular
`(,@cl-alloc-args
,@(subst 'functions +callback-argument-name+ class-slots-instance)
,@(when initargs
(append
(list
......@@ -171,26 +176,28 @@
(let ((object
(make-instance
',class
,@(symbol-keyword-symbol
(remove +callback-argument-name+ class-slots-instance))
,@(when (member +callback-argument-name+ class-slots-instance)
(symbol-keyword-symbol 'functions))
,@(if arglists
(second arglists)
(symbol-keyword-symbol cl-alloc-args singular)))))
;; There is callback slot variable
,@(when (member +callback-argument-name+ class-slots-instance)
`((setf (slot-value object '#.+callback-argument-name+)
(make-cbstruct-object object))))
;; There is an initialization step
,@(when initializerp
(if initialize-args ; with arguments
(if initargs
`((when ,settingp
(reinitialize-instance
(let ((reii
`(reinitialize-instance
object
,@(if arglists
(third arglists)
(symbol-keyword-symbol
cl-initialize-args singular)))))
`((reinitialize-instance
object
,@(if arglists
(third arglists)
(symbol-keyword-symbol
cl-initialize-args singular)))))
(if initargs `((when ,settingp ,reii)) `(,reii)))
'((reinitialize-instance object)))) ; without arguments
object)))
......
;; 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-14 19:29:21EST evolution.lisp>
;; $Id$
(in-package :gsl)
(defmobject ode-evolution "gsl_odeiv_evolve"
((dimensions sizet))
(((first dimensions) sizet))
"evolution for ordinary differential equations"
:documentation
"Make an object to advance the ODE solution."
:superclasses (callback-included-cl)
:ci-class-slots (ode-system marray (function jacobian) (dimension))
:class-slots-instance (#.+callback-argument-name+)
:initialize-suffix "reset"
:initialize-args nil)
:initialize-args nil
:singular (dimension))
(def-make-callbacks ode-evolution
(function jacobian dimension &optional (scalars t))
(if scalars
`(progn
(defmcallback ,function
:success-failure
(:double (:double ,dimension) (:set :double ,dimension))
nil nil
,function)
(defmcallback ,jacobian
:success-failure
(:double
(:double ,dimension)
(:set :double ,(expt dimension 2))
(:set :double ,dimension))
nil nil ,jacobian))
`(progn
(defmcallback ,function
:success-failure
(:double :pointer :pointer)
nil nil
,function)
(defmcallback ,jacobian
:success-failure
(:double :pointer :pointer :pointer)
nil nil ,jacobian))))
(defmfun apply-evolution
(evolve control step dydt time max-time step-size y)
(evolve control step time max-time step-size y)
"gsl_odeiv_evolve_apply"
(((mpointer evolve) :pointer) ((mpointer control) :pointer)
((mpointer step) :pointer)
(dydt :pointer) ((c-pointer time) :pointer)
((callback-struct evolve) :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-14 18:17:29EST 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-evolution vanderpol vanderpol-jacobian 2))
(defun integrate-vanderpol
(max-time &optional (step-size 1.0d-6) (stepper *step-rk8pd*) (print-steps t))
......@@ -37,9 +35,7 @@
(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
......
;; 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-14 19:31:58EST 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
......@@ -52,7 +26,7 @@
(cstep (make-symbol "CSTEP")))
`(let ((stepperobj (make-ode-stepper ,stepper ,dimensions))
(control (make-y-control ,absolute-error ,relative-error))
(evolve (make-ode-evolution ,dimensions))
(evolve (make-ode-evolution ,dimensions ',function))
(,dep
(make-marray 'double-float :dimensions ,dimensions))
(,ctime (make-marray 'double-float :dimensions 1))
......@@ -65,7 +39,5 @@
collect `(,symb (maref ,dep ,i)))
(make-next-step
(apply-evolution
evolve control stepperobj
,function ,ctime
,max-time ,cstep ,dep)))
evolve control stepperobj ,ctime ,max-time ,cstep ,dep)))
,@body))))
;; Generic functions for optimization
;; Liam Healy 2009-01-03 12:59:07EST generic.lisp
;; Time-stamp: <2009-02-08 15:15:15EST generic.lisp>
;; Time-stamp: <2009-02-14 12:25:14EST generic.lisp>
;; $Id: $
(in-package :gsl)
......@@ -58,42 +58,3 @@
(df :pointer)
(fdf :pointer)
(parameters :pointer))
;;; To be obsolete
(export 'make-solver-functions)
(defmacro make-solver-functions (function &optional df fdf dimensions array)
"Setup functions for solvers.
The CL functions can be specified as a previously-defined defun or
as a lambda. 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.
If df and fdf are non-nil, define for derivative solvers."
(if df
(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)))))))
`(make-single-function ,function :success-failure :double gsl-mfunction
,dimensions)))
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