Commit 741fd237 authored by Liam Healy's avatar Liam Healy

ODE solver functional with funcallables

Changed the defmobject definitions to handle a 'callback-included-cl
object.  Changed the ODE solver definitions to take separate arguments
for the function and Jacobian, as well as scalarsp.  Not completely
finished, but tests
(integrate-vanderpol 1.0d0 1.d-4 +step-rk2+ nil)
(integrate-vanderpol 1.0d0 1.d-4 +step-rkf45+ nil)
work correctly.
parent 1349c73a
;; The mobject that defines callbacks
;; Liam Healy 2009-03-14 11:20:03EDT callback-included.lisp
;; Time-stamp: <2009-03-29 09:52:12EDT callback-included.lisp>
;; Time-stamp: <2009-03-29 16:12:31EDT callback-included.lisp>
;; $Id: $
(in-package :gsl)
......@@ -59,8 +59,11 @@
;;;;****************************************************************************
;;; This is expanded in defmfun to set the dynamic variables.
(defun callback-set-dynamic (callback-object arglist)
(defun callback-set-dynamic (callback-object &optional arglist)
"Make a form to set the dynamic variable defining callbacks."
(when (listp callback-object)
(setf arglist callback-object
callback-object (caar callback-object)))
`((setf
,@(loop for symb in
(let ((class (category-for-argument arglist callback-object)))
......@@ -68,4 +71,4 @@
class
(number-of-callbacks (get-callbacks-for-class class))))
for n from 0
append `(,symb (nth ,n (funcallables ,callback-object)))))))
append `(,symb (nth ,n (funcallables ,callback-object))))))))
;; Macro for defining GSL functions.
;; Liam Healy 2008-04-16 20:49:50EDT defmfun.lisp
;; Time-stamp: <2009-03-14 21:54:00EDT defmfun.lisp>
;; Time-stamp: <2009-03-29 16:28:26EDT defmfun.lisp>
;; $Id$
(in-package :gsl)
......@@ -106,7 +106,9 @@
;; A list of variable names, and a list of callback names
(when (or callbacks callback-object)
(if callback-object
(let ((class (category-for-argument arglist callback-object)))
(let ((class (if (listp callback-object)
(second (first callback-object))
(category-for-argument arglist callback-object))))
(list (mobject-fnvnames
class
(number-of-callbacks (get-callbacks-for-class class)))
......
;; Generate a lambda that calls the user function; will be called by callback.
;; Liam Healy
;; Time-stamp: <2009-03-29 12:04:35EDT funcallable.lisp>
;; Time-stamp: <2009-03-29 22:37:19EDT funcallable.lisp>
;; $Id$
(in-package :gsl)
......@@ -42,12 +42,14 @@
(and (listp arg)
(eql (parse-callback-argspec arg 'io) direction)))))
(defun array-direction-specs (argspecs direction)
(defun vspecs-direction (argspecs direction &optional array-only)
"Find the specs for all variables, or all array variables, with the
specified direction."
(remove-if-not
(lambda (arg)
(when
(and (eql (parse-callback-argspec arg 'io) direction)
(parse-callback-argspec arg 'dimensions))
(if array-only (parse-callback-argspec arg 'dimensions) t))
arg))
argspecs))
......@@ -130,13 +132,13 @@
expects and returns scalars, and dimension-values should be a list
of number(s), (dim0) or (dim0 dim1), or NIL."
(let* ((argspecs (remove :slug (parse-callback-fnspec fnspec 'arguments-spec)))
(inargs-specs (array-direction-specs argspecs :input))
(inargs-specs (vspecs-direction argspecs :input))
(inargs-names
(make-symbol-cardinals
'input
(length (remove nil (mapcar (all-io :input nil) argspecs)))))
(outargs-specs (array-direction-specs argspecs :output))
(outargs-names (make-symbol-cardinals 'output (length outargs-specs)))
(outarrayp (vspecs-direction argspecs :output t))
(outargs-names (make-symbol-cardinals 'output (length outarrayp)))
(lambda-args
(loop for arg in argspecs
with oarg = (copy-list outargs-names)
......@@ -152,16 +154,20 @@
(let ((uf user-function)) `',uf)
user-function)))
`(lambda ,lambda-args
,(if (and scalarsp (or inargs-specs outargs-specs))
,(if (and scalarsp (or inargs-specs outarrayp))
(let ((call-form
`(funcall
,function-designator
,@(array-element-refs inargs-names inargs-specs dimension-values))))
(if outargs-specs
(if outarrayp
(callback-set-mvb outargs-names call-form fnspec dimension-values)
;; no specified output, return what the function returns
call-form))
`(funcall ,function-designator ,@(append inargs-names outargs-names)))
`(funcall ,function-designator
,@(if outarrayp
(append inargs-names outargs-names)
;; no arrays to return, just return the value
inargs-names)))
,@(case
(parse-callback-fnspec fnspec 'return-spec)
(:success-failure
......
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2009-03-29 12:52:22EDT mobject.lisp>
;; Time-stamp: <2009-03-29 23:04:45EDT mobject.lisp>
;;; GSL objects are represented in GSLL as and instance of a 'mobject.
;;; The macro demobject takes care of defining the appropriate
......@@ -106,7 +106,8 @@
class cl-initialize-args initialize-name prefix
initialize-suffix initialize-args inputs
(and (not (callback-arg-p class-slots-instance callbacks))
callbacks)))
callbacks)
superclasses))
(export '(,maker ,class))
,@(when callbacks `((record-callbacks-for-class ',class ',callbacks)))
,@(when callbacks (make-mobject-defmcallbacks callbacks class))
......@@ -142,7 +143,7 @@
(defun make-reinitialize-instance
(class cl-initialize-args initialize-name prefix
initialize-suffix initialize-args inputs
callbacks)
callbacks superclasses)
"Expand the reinitialize-instance form."
(let ((cbstruct (make-symbol "CBSTRUCT")))
`((defmfun reinitialize-instance
......@@ -167,7 +168,9 @@
,@(when inputs `(:inputs ,inputs))
,@(when callbacks
`(:before
((make-funcallables-for-object object))
(,@(when (member 'callback-included-cl superclasses)
`((setf (slot-value object 'callback) ,cbstruct)))
(make-funcallables-for-object object))
:after
((trivial-garbage:finalize
object
......@@ -226,7 +229,7 @@
(tg:finalize object (lambda () (foreign-free ,cbs)))))))
;; There is an initialization step
,@(when initializerp
(if initialize-args ; with arguments
(if (or initialize-args arglists) ; with arguments
(let ((reii
`(reinitialize-instance
object
......
;; Evolution functions for ODE integration.
;; Liam Healy, Sun Sep 30 2007 - 14:31
;; Time-stamp: <2009-02-15 08:58:47EST evolution.lisp>
;; Time-stamp: <2009-03-29 16:25:04EDT evolution.lisp>
;; $Id$
(in-package :gsl)
......@@ -22,6 +22,7 @@
((c-pointer step-size) :pointer) ((c-pointer y) :pointer))
:inputs (time step-size y)
:outputs (time step-size y)
:callback-object ((stepper ode-stepper))
:documentation ; FDL
"Advance the system (e, dydt) from time
and position y using the stepping function step.
......@@ -33,4 +34,3 @@
changed the value of step-size will be modified on output. The maximum
time max-time is guaranteed not to be exceeded by the time-step. On the
final time-step the value of time will be set to t1 exactly.")
;; Example ODE
;; Liam Healy Sat Sep 29 2007 - 17:49
;; Time-stamp: <2009-03-29 12:43:02EDT ode-example.lisp>
;; Time-stamp: <2009-03-29 22:51:51EDT ode-example.lisp>
;; $Id$
;;; van der Pol as given in Section 25.5 of the GSL manual. To
......@@ -33,8 +33,8 @@
(let ((mu 10.0d0) (initial-time 0.0d0) (iter 0))
(declare (special mu))
(with-ode-integration
((vanderpol vanderpol-jacobian)
time step max-time (dep0 dep1) 2 stepper)
(vanderpol time step max-time (dep0 dep1) 2
:jacobian vanderpol-jacobian :stepper stepper)
(setf dep0 1.0d0 dep1 0.0d0 step step-size time initial-time)
(loop
(when (or (>= time max-time) (> iter *max-iter*))
......
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2009-03-22 19:54:05EDT ode-system.lisp>
;; Time-stamp: <2009-03-29 22:54:21EDT ode-system.lisp>
;; $Id$
(in-package :gsl)
......@@ -8,15 +8,19 @@
(export '(with-ode-integration))
(defmacro with-ode-integration
((function time step-size max-time
dependent dimensions &optional (stepper '+step-rk8pd+)
((function time step-size max-time dependent dimensions
&key jacobian (scalarsp t) (stepper '+step-rk8pd+)
(absolute-error 1.0d-6) (relative-error 0.0d0))
&body body)
"Environment for integration of ordinary differential equations."
;; Note: the case jacobian=nil is not properly handled yet; it
;; should put a null pointer into the struct that is passed to GSL.
(let ((dep (make-symbol "DEP"))
(ctime (make-symbol "CTIME"))
(cstep (make-symbol "CSTEP")))
`(let ((stepperobj (make-ode-stepper ,stepper ,dimensions ',function))
`(let ((stepperobj
(make-ode-stepper
,stepper ,dimensions ',function ',jacobian ,scalarsp))
(control (make-y-control ,absolute-error ,relative-error))
(evolve (make-ode-evolution ,dimensions))
(,dep
......
;; Stepping functions for ODE systems.
;; Liam Healy, Mon Sep 24 2007 - 21:33
;; Time-stamp: <2009-03-29 12:34:49EDT stepping.lisp>
;; Time-stamp: <2009-03-29 17:19:19EDT stepping.lisp>
;; $Id$
;; /usr/include/gsl/gsl_odeiv.h
......@@ -31,10 +31,13 @@
(:output :double :cvector dim0 dim0) ; dfdy
(:output :double :cvector dim0) ; dydt
:slug))
:class-slots-instance (callback)
:initialize-suffix "reset"
:initialize-args nil
:singular (dimension))
:arglists-function
(lambda (set)
`((type dimension &optional (function nil ,set) jacobian (scalarsp t))
(:type type :dimensions (list dimension))
(:functions (list function jacobian) :scalarsp scalarsp))))
(cffi:defcstruct ode-system ; gsl_odeiv_system
"The definition of an ordinary differential equation system for GSL."
......
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