Commit 5aafab66 authored by liam's avatar liam

Multivariate root solving with derivative now functions. Symbol macro

'success returns the GSL_SUCCESS code, required for many callback
definitions.  Introduced macros vref and mref (setfable) to get
components of GSL vector doubles and matrix doubles directly from the
pointers.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3265 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent 9f0886b5
;********************************************************
; file: data.lisp
; description: Using GSL storage.
; date: Sun Mar 26 2006 - 16:32
; author: Liam M. Healy
; modified: Sun Jan 13 2008 - 22:49
;********************************************************
;;; $Id: $
;; Using GSL storage.
;; Liam Healy, Sun Mar 26 2006 - 16:32
;; Time-stamp: <2008-01-14 21:09:15 liam data.lisp>
;; $Id: $
(in-package :gsl)
......@@ -234,14 +230,6 @@
(progn ,@body)
(free ,symbol))))
(defun make-data-from-pointer (pointer &optional (class 'gsl-vector-double) size)
"Given a C pointer to a GSL data type, make the CL object."
(make-instance
class
:pointer pointer
:storage-size
(or size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))))
;;;;****************************************************************************
;;;; Getting values into CL
;;;;****************************************************************************
......
;********************************************************
; file: matrix.lisp
; description: Matrices
; date: Sun Mar 26 2006 - 11:51
; author: Liam M. Healy
; modified: Sun Dec 3 2006 - 18:53
;********************************************************
;;; $Id: $
;; Matrices
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-01-15 19:34:56 liam matrix.lisp>
;; $Id: $
(in-package :gsl)
......@@ -51,6 +47,12 @@
(defmethod gsl-array ((object gsl-matrix))
(foreign-slot-value (pointer object) 'gsl-matrix-c 'data))
(export 'matrix-data)
(defun matrix-data (pointer)
"A pointer to the GSL array with the data contents, from the
sruct pointer."
(cffi:foreign-slot-value pointer 'gsl-matrix-c 'data))
;;;;****************************************************************************
;;;; Getting values
;;;;****************************************************************************
......@@ -63,6 +65,13 @@
:c-return :c-base-type
:documentation "The (i,j)-th element of the matrix.")
(defun-gsl mref (pointer index0 index1)
"gsl_matrix_get"
((pointer :pointer) (index0 :size) (index1 :size))
:c-return :double
:index nil
:documentation "An element of the matrix of doubles, computed from the pointer.")
(defun-gsl-mdsfc gsl-matrix-ptr ((matrix gsl-matrix) i j)
"gsl_matrix_ptr" (((pointer matrix) :pointer) (i :size) (j :size))
:c-return :pointer
......@@ -96,6 +105,13 @@
:c-return :void
:documentation "Set the (i,j)-th element of the matrix.")
(defun-gsl (setf mref) (value pointer index0 index1)
"gsl_matrix_set"
((pointer :pointer) (index0 :size) (index1 :size) (value :double))
:c-return :void
:index nil
:documentation "Set an element of the matrix of doubles, using its pointer.")
(defmethod (setf data) (array (object gsl-matrix))
(loop for i from 0
below
......
;********************************************************
; file: vector.lisp
; description: Vectors
; date: Sun Mar 26 2006 - 11:51
; author: Liam M. Healy
; modified: Sun Jan 13 2008 - 22:40
;********************************************************
;;; $Id: $
;; Vectors
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2008-01-15 19:01:37 liam vector.lisp>
;; $Id: $
(in-package :gsl)
......@@ -79,7 +75,16 @@ deallocated with the vector.
(defmethod gsl-array ((object gsl-vector))
(cffi:foreign-slot-value (pointer object) 'gsl-vector-c 'data))
(defun gsl-array-p (pointer)
(defun make-data-from-pointer (pointer &optional (class 'gsl-vector-double) size)
"Given a C pointer to a GSL data type, make the CL object."
(make-instance
class
:pointer pointer
:storage-size
(or size (cffi:foreign-slot-value pointer 'gsl-vector-c 'size))))
(export 'vector-data)
(defun vector-data (pointer)
"A pointer to the GSL array with the data contents, from the
sruct pointer."
(cffi:foreign-slot-value pointer 'gsl-vector-c 'data))
......@@ -94,6 +99,13 @@ deallocated with the vector.
:c-return :c-base-type
:documentation "The ith element of the vector.")
(defun-gsl vref (pointer index)
"gsl_vector_get"
((pointer :pointer) (index :size))
:c-return :double
:index nil
:documentation "An element of the vector of doubles, computed from the pointer.")
(defun-gsl gsl-vector-ptr (vector i)
"gsl_vector_ptr" (((pointer vector) :pointer) (i :size))
:c-return :pointer
......@@ -107,7 +119,14 @@ deallocated with the vector.
"gsl_vector_set"
(((pointer vector) :pointer) ((first indices) :size) (value :c-base-type))
:c-return :void
:documentation "Set the ith element of the vector.")
:documentation "Set an element of the vector.")
(defun-gsl (setf vref) (value pointer index)
"gsl_vector_set"
((pointer :pointer) (index :size) (value :double))
:c-return :void
:index nil
:documentation "Set an element of the vector of doubles, using its pointer.")
(defun-gsl-vdsfc set-all ((object gsl-vector) value)
"gsl_vector_set_all"
......
;********************************************************
; file: conditions.lisp
; description: GSL errors
; date: Sat Mar 4 2006 - 18:33
; author: Liam M. Healy
; modified: Sun Apr 2 2006 - 15:56
;********************************************************
;; GSL errors
;; Liam Healy Sat Mar 4 2006 - 18:33
;; Time-stamp: <2008-01-14 22:22:55 liam conditions.lisp>
;; $Id: $
(in-package :gsl)
......@@ -98,3 +95,11 @@
(cffi:foreign-funcall
"gsl_set_error_handler"
:pointer (cffi:callback gsl-error))
(defmacro gsl-errorno-sm (keyword)
`(define-symbol-macro
,(intern (symbol-name keyword) (find-package :gsl))
(cffi:foreign-enum-value 'gsl-errorno ,keyword)))
(gsl-errorno-sm :success)
;********************************************************
; file: functions.lisp
; description: Foreign callback functions.
; date: Sun Dec 9 2007 - 22:08
; author: Liam Healy
; modified: Sun Jan 13 2008 - 13:26
;********************************************************
;;; $Id: $
;; Foreign callback functions.
;; Liam Healy
;; Time-stamp: <2008-01-14 22:56:45 liam functions.lisp>
;; $Id: $
(in-package :gsl)
......@@ -90,7 +86,7 @@
,@(when (eq return-type :success-failure)
;; We always return success, because if there was a
;; problem, a CL error would be signalled.
'((cffi:foreign-enum-value 'gsl-errorno :SUCCESS))))
'(success)))
,@(when
structure
;; Assume that defcallback does not bind the variable 'name.
......
......@@ -313,7 +313,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
<tr>
<td><a href="Multidimensional-Root_Finding.html">
Multidimensional Root Finding</a></td>
<td>Started</td>
<td>Done</td>
</tr>
<tr>
<td><a href="Multidimensional-Minimization.html">
......@@ -357,7 +357,7 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
<small>
Time-stamp: <2008-01-13 16:51:23 liam index.html>
Time-stamp: <2008-01-15 22:37:31 liam index.html>
</small>
<!-- hhmts end -->
</div>
......
;********************************************************
; file: interface.lisp
; description: Macros to interface GSL functions.
; date: Mon Mar 6 2006 - 22:35
; author: Liam M. Healy
; modified: Sun Jan 6 2008 - 13:48
;********************************************************
;; Macros to interface GSL functions.
;; Liam Healy
;; Time-stamp: <2008-01-15 18:55:45 liam interface.lisp>
;; $Id: $
(in-package :gsl)
......@@ -101,7 +98,7 @@
(defun success-continue (value)
"If status indicates success, return T, othewise return NIL."
(eql value (cffi:foreign-enum-value 'gsl-errorno :SUCCESS)))
(eql value success))
;;;;****************************************************************************
;;;; Macro defun-gsl
......@@ -125,7 +122,8 @@
;;; or its derivatives.
;;; Default are allocated quantities in c-arguments, or :c-return if none.
;;; type :function or :method
;;; index Name under which this function should be cross-referenced
;;; index Name under which this function should be cross-referenced; t
;;; to use name, nil to not index.
;;; export Whether to export the symbol.
;;; null-pointer-info Return value if C function returns a null pointer.
;;; documentation
......@@ -136,7 +134,7 @@
(name arglist gsl-name c-arguments
&key (c-return :error-code)
(return nil return-supplied-p)
(type :function) index (export (not (eq type :method)))
(type :function) (index t) (export (not (eq type :method)))
null-pointer-info documentation invalidate after enumeration)
(let* ((cargs (substitute '(mode sf-mode) :mode c-arguments))
(carg-symbs
......@@ -164,69 +162,69 @@
(unless (eq c-return :void)
(list cret-name)))))
`(progn
(,(if (eq type :function) 'defun 'defmethod)
,name
,(if (member :mode c-arguments)
`(,@clargs &optional (mode :double-prec))
`(,@clargs))
,@(when documentation (list documentation))
(,@(if allocated
`(cffi:with-foreign-objects
,(mapcar #'wfo-declare allocated-decl))
'(let ()))
(let ((,cret-name
(cffi:foreign-funcall
,gsl-name
,@(mapcan
(lambda (arg)
(list (if (member (st-symbol arg) allocated)
:pointer
(st-type arg))
(st-symbol arg)))
cargs)
,cret-type)))
,@(case c-return
(:void `((declare (ignore ,cret-name))))
(:error-code ; fill in arguments
`((check-gsl-status ,cret-name ',name))))
,@(when invalidate `((cl-invalidate ,@invalidate)))
,@(when (or null-pointer-info (eq c-return :pointer))
`((check-null-pointer ,cret-name
,@(or null-pointer-info
'(:ENOMEM "No memory allocated")))))
,@after
(values
,@(case c-return
(:number-of-answers
(mapcan
(lambda (vbl seq)
`((when (> ,cret-name ,seq) ,vbl)))
clret
(loop for i below (length clret) collect i)))
(:success-failure
(if (equal clret invalidate)
;; success-failure more important than passed-in
`(success-failure ,@clret)
(remove cret-name ; don't return c-return itself
`(,@clret (success-failure ,cret-name)))))
(:success-continue
(if (equal clret invalidate)
;; success-failure more important than passed-in
`(success-continue ,@clret)
(remove cret-name ; don't return c-return itself
`(,@clret (success-continue ,cret-name)))))
(:true-false
`((not (zerop ,cret-name))))
(:enumerate
`((cffi:foreign-enum-keyword ',enumeration ,cret-name)))
(t (unless
(or
(and (eq c-return :error-code)
(not allocated))
(and (null return) return-supplied-p))
clret)))))))
(map-name ',(or index name) ,gsl-name)
,@(when export `((export ',name))))))
(,(if (eq type :function) 'defun 'defmethod)
,name
,(if (member :mode c-arguments)
`(,@clargs &optional (mode :double-prec))
`(,@clargs))
,@(when documentation (list documentation))
(,@(if allocated
`(cffi:with-foreign-objects
,(mapcar #'wfo-declare allocated-decl))
'(let ()))
(let ((,cret-name
(cffi:foreign-funcall
,gsl-name
,@(mapcan
(lambda (arg)
(list (if (member (st-symbol arg) allocated)
:pointer
(st-type arg))
(st-symbol arg)))
cargs)
,cret-type)))
,@(case c-return
(:void `((declare (ignore ,cret-name))))
(:error-code ; fill in arguments
`((check-gsl-status ,cret-name ',name))))
,@(when invalidate `((cl-invalidate ,@invalidate)))
,@(when (or null-pointer-info (eq c-return :pointer))
`((check-null-pointer ,cret-name
,@(or null-pointer-info
'(:ENOMEM "No memory allocated")))))
,@after
(values
,@(case c-return
(:number-of-answers
(mapcan
(lambda (vbl seq)
`((when (> ,cret-name ,seq) ,vbl)))
clret
(loop for i below (length clret) collect i)))
(:success-failure
(if (equal clret invalidate)
;; success-failure more important than passed-in
`(success-failure ,@clret)
(remove cret-name ; don't return c-return itself
`(,@clret (success-failure ,cret-name)))))
(:success-continue
(if (equal clret invalidate)
;; success-failure more important than passed-in
`(success-continue ,@clret)
(remove cret-name ; don't return c-return itself
`(,@clret (success-continue ,cret-name)))))
(:true-false
`((not (zerop ,cret-name))))
(:enumerate
`((cffi:foreign-enum-keyword ',enumeration ,cret-name)))
(t (unless
(or
(and (eq c-return :error-code)
(not allocated))
(and (null return) return-supplied-p))
clret)))))))
,@(when index `((map-name ',(if (eql index t) name index) ,gsl-name)))
,@(when export `((export ',name))))))
(defmacro defun-optionals
(name arglist no-optional optionals &optional documentation)
......
;********************************************************
; file: ode-system.lisp
; description: ODE system setup
; date: Sun Apr 15 2007 - 14:19
; author: Liam Healy
; modified: Sat Jan 5 2008 - 21:39
;********************************************************
;;; $Id: $
;; ODE system setup
;; Liam Healy, Sun Apr 15 2007 - 14:19
;; Time-stamp: <2008-01-14 22:58:19 liam ode-system.lisp>
;; $Id: $
(in-package :gsl)
......@@ -42,7 +38,7 @@
(,params :pointer))
(declare (ignore ,params))
(,name ,time ,dependent ,derivatives)
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS))
success)
;; The function should take four arguments: time, dependent, dfdy, dfdt
;; The last three will be arrays.
(cffi:defcallback ,jacobian :int
......@@ -53,7 +49,7 @@
(,params :pointer))
(declare (ignore ,params))
(,name ,time ,dependent ,dfdy ,dfdt)
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS))
success)
;; Assume that defcallback does not bind the variable 'name.
(defparameter ,name (cffi:foreign-alloc 'ode-system))
(set-slot-function ,name 'ode-system 'function ',name)
......
;;; Multivariate roots.
;;; Liam Healy 2008-01-12 12:49:08
;;; Time-stamp: <2008-01-13 22:51:02 liam roots-multi.lisp>
;;; Time-stamp: <2008-01-15 22:44:27 liam roots-multi.lisp>
;;; $Id: $
(in-package :gsl)
;;; I don't like using make-data-from-pointer.
;;; I don't like using make-data-from-pointer, but it's the only way
;;; to have access to the GSL functions when given a pointer.
;;; Alternatively, I could provide the GSL pointer and then the only
;;; thing the user could use is vref, or, of course
;;; make-data-from-pointer.
;;;;****************************************************************************
;;;; Function definition
......@@ -25,6 +29,16 @@
((dimensions ,dimensions))
((returned-value gsl-vector-c))))
(cffi:defcstruct gsl-mfunction-fdf
;; See /usr/include/gsl/gsl_multiroots.h
"The definition of a function and its derivatives for multiroot
finding in GSL."
(function :pointer)
(df :pointer)
(fdf :pointer)
(dimensions :size)
(parameters :pointer))
;;;;****************************************************************************
;;;; Initialization
;;;;****************************************************************************
......@@ -55,7 +69,7 @@
(defun-gsl set-mfdfsolver (solver function-derivative initial)
"gsl_multiroot_fdfsolver_set"
((solver :pointer) (function-derivative :pointer)
((gsl-array initial) :pointer))
((pointer initial) :pointer))
:documentation
"Set or reset an existing solver to use the function and derivative
(fdf) and the initial guess.")
......@@ -383,13 +397,22 @@
(defparameter *gsl-vector*
(make-instance 'gsl-vector-double :pointer nil :storage-size nil))
;;; One alternative way of writing the function, not recommended.
(defun rosenbrock (argument return)
"Rosenbrock test function."
(with-c-doubles (((gsl-array-p argument) x0 x1)
((gsl-array-p return) f0 f1))
(with-c-doubles (((vector-data argument) x0 x1)
((vector-data return) f0 f1))
(setf f0 (* *rosenbrock-a* (- 1 x0))
f1 (* *rosenbrock-b* (- x1 (expt x0 2))))))
;;; The recommended alternative
(defun rosenbrock (argument return)
"Rosenbrock test function."
(setf (vref return 0)
(* *rosenbrock-a* (- 1 (vref argument 0)))
(vref return 1)
(* *rosenbrock-b* (- (vref argument 1) (expt (vref argument 0) 2)))))
(def-mfunction rosenbrock 2)
(defun roots-multi-example ()
......@@ -416,3 +439,51 @@
(gsl-aref argval 1)
(gsl-aref fnval 0)
(gsl-aref fnval 1)))))))))
(defun rosenbrock-df (argument jacobian)
"The partial derivatives of the Rosenbrock functions."
(setf (mref jacobian 0 0) (- *rosenbrock-a*)
(mref jacobian 0 1) 0.0d0
(mref jacobian 1 0) (* -2 *rosenbrock-b* (vref argument 0))
(mref jacobian 1 1) *rosenbrock-b*))
(defun rosenbrock-fdf (argument value jacobian)
(rosenbrock argument value)
(rosenbrock-df argument jacobian))
;;; Because def-solver-functions and def-scalar-function bind a symbol
;;; of the same name as the first function, and we want both to run,
;;; we'll make an alias function so we can use both.
(eval-when (:load-toplevel :execute)
(setf (fdefinition 'rosenbrock-f) #'rosenbrock))
(def-solver-functions rosenbrock-f rosenbrock-df rosenbrock-fdf 2)
(defun roots-multi-example-df ()
"Solving Rosenbrock with derivatives, the example given in Sec. 34.8
of the GSL manual."
(flet ((print-state (iter argval fnval)
(format t "~&iter=~d~8tx0=~12,8g~24tx1=~12,8g~38tf0=~12,8g~52tf1=~12,8g"
iter
(gsl-aref argval 0)
(gsl-aref argval 1)
(gsl-aref fnval 0)
(gsl-aref fnval 1))))
(let ((max-iter 1000))
(with-data (vect vector-double 2)
(setf (data vect) #(-10.0d0 -5.0d0))
(with-mfdfsolver (solver *gnewton-mfdfsolver* rosenbrock-f vect)
(let ((fnval (mfdfsolver-f solver))
(argval (mfdfsolver-root solver)))
(loop for iter from 0
while (and (< iter max-iter)
(not (multiroot-test-residual solver 1.0d-7)))
initially (print-state iter argval fnval)
do
(iterate-mfdfsolver solver)
(print-state iter argval fnval)
finally (return
(values (gsl-aref argval 0)
(gsl-aref argval 1)
(gsl-aref fnval 0)
(gsl-aref fnval 1))))))))))
;********************************************************
; file: roots-one.lisp
; description: One-dimensional root solver.
; date: Sun Dec 9 2007 - 15:47
; author: Liam Healy
; modified: Sun Jan 13 2008 - 17:00
;********************************************************
;;; Time-stamp: <2008-01-13 17:00:49 liam roots-one.lisp>
;;; $Id: $
;; One-dimensional root solver.
;; Liam Healy
;; Time-stamp: <2008-01-15 22:57:11 liam roots-one.lisp>
;; $Id: $
(in-package :gsl)
......@@ -25,38 +20,51 @@
(parameters :pointer))
(export 'def-solver-functions)
(defmacro def-solver-functions (function df fdf)
(defmacro def-solver-functions (function df fdf &optional dimensions)
"Setup functions for solvers.
The CL functions name and derivative should be defined previously
with defuns."
with defuns. If dimensions is non-nil, set multiroot solver
functions."
(let ((arg (make-symbol "ARG"))
(params (make-symbol "PARAMS"))
(func (make-symbol "F"))
(deriv (make-symbol "DF")))
(deriv (make-symbol "DF"))
(struct (if dimensions 'gsl-mfunction-fdf 'gsl-function-fdf))
(argtype (if dimensions :pointer :double))
(rettype (if dimensions :int :double)))
`(progn
(cffi:defcallback ,function :double
((,arg :double)
(,params :pointer))
(cffi:defcallback ,function ,rettype
((,arg ,argtype)
(,params :pointer)
,@(when dimensions `((,func :pointer))))
(declare (ignore ,params))
(,function ,arg))
(cffi:defcallback ,df :double
((,arg :double)
(,params :pointer))
(,function ,arg ,@(when dimensions `(,func)))
,@(when dimensions '(success)))
(cffi:defcallback ,df ,rettype
((,arg ,argtype)
(,params :pointer)
,@(when dimensions `((,deriv :pointer))))
(declare (ignore ,params))
(,df ,arg))
(cffi:defcallback ,fdf :pointer
((,arg :double)
(,df ,arg ,@(when dimensions `(,deriv)))
,@(when dimensions '(success)))
(cffi:defcallback ,fdf ,(if dimensions :int :pointer)
((,arg ,argtype)
(,params :pointer)
(,func :pointer)
(,deriv :pointer))
(declare (ignore ,params))
(,fdf ,arg ,func ,deriv)
(cffi:null-pointer))
(defparameter ,function (cffi:foreign-alloc 'gsl-function-fdf))
(set-slot-function ,function 'gsl-function-fdf 'function ',function)
(set-slot-function ,function 'gsl-function-fdf 'df ',df)
(set-slot-function ,function 'gsl-function-fdf 'fdf ',fdf)
(set-parameters ,function 'gsl-function-fdf))))
;; fdf function returns sucess code in multivariate case, void
;; pointer in univariate
,(if dimensions 'success
'(cffi:null-pointer)))
(defparameter ,function (cffi:foreign-alloc ',struct))
(set-slot-function ,function ',struct 'function ',function)
(set-slot-function ,function ',struct 'df ',df)
(set-slot-function ,function ',struct 'fdf ',fdf)
,@(when dimensions
`((set-structure-slot ,function ',struct 'dimensions ,dimensions)))
(set-parameters ,function ',struct))))
;;;;****************************************************************************
;;;; Initialization
......@@ -395,7 +403,7 @@
root (- root (sqrt 5.0d0))
(- upper lower))))))
;;; Because def-solver-functions and def-scalar-fucntion bind a symbol
;;; Because def-solver-functions and def-scalar-function bind a symbol
;;; of the same name as the first function, and we want both to run,
;;; we'll make an alias function so we can use both.
(eval-when (:load-toplevel :execute)
......@@ -420,3 +428,4 @@
do
(format t "~&~d~6t~10,8g ~18t~10,6g~34t~10,6g"
iter root (- root (sqrt 5.0d0)) (- root oldroot))))))
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