Skip to content
Snippets Groups Projects
Commit d5250819 authored by liam's avatar liam
Browse files

Added simulated-annealing (not working well) and

ordinary-differential-equations (started).  Change *-to-cl from
functions into macros, eliminate obsolete definitions.  Update
descriptions.


git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3226 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent c02a8b7f
No related branches found
No related tags found
No related merge requests found
...@@ -50,13 +50,7 @@ gsl_vector_short.h ...@@ -50,13 +50,7 @@ gsl_vector_short.h
gsl_vector_uchar.h gsl_vector_uchar.h
gsl_vector_char.h gsl_vector_char.h
I have implemented only one, double. I implemented only int, double, complex.
This means that the following functions have been put aside:
- #'polynomial-solve (requires vector_complex)
- The index functions in sorting that are not returning permutations,
because they need a vector_int.
It would be nice to use e.g. vector-sap for direct vectors in It would be nice to use e.g. vector-sap for direct vectors in
e.g. SBCL. e.g. SBCL.
......
...@@ -272,61 +272,61 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre> ...@@ -272,61 +272,61 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre>
</tr> </tr>
<tr> <tr>
<td><a href="Interpolation.html">Interpolation</a></td> <td><a href="Interpolation.html">Interpolation</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Numerical-Differentiation.html">Numerical Differentiation</a></td> <td><a href="Numerical-Differentiation.html">Numerical Differentiation</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Chebyshev-Approximations.html">Chebyshev Approximations</a></td> <td><a href="Chebyshev-Approximations.html">Chebyshev Approximations</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Series-Acceleration.html">Series Acceleration</a></td> <td><a href="Series-Acceleration.html">Series Acceleration</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Wavelet-Transforms.html">Wavelet Transforms</a></td> <td><a href="Wavelet-Transforms.html">Wavelet Transforms</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Discrete-Hankel-Transforms.html">Discrete Hankel Transforms</a></td> <td><a href="Discrete-Hankel-Transforms.html">Discrete Hankel Transforms</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a <td><a
href="One-dimensional-Root_002dFinding.html">One-dimensional Root Finding</a></td> href="One-dimensional-Root_002dFinding.html">One-dimensional Root Finding</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a <td><a
href="One-dimensional-Minimization.html">One-dimensional Minimization</a></td> href="One-dimensional-Minimization.html">One-dimensional Minimization</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Multidimensional-Root_002dFinding.html"> <td><a href="Multidimensional-Root_002dFinding.html">
Multidimensional Root Finding</a></td> Multidimensional Root Finding</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Multidimensional-Minimization.html"> <td><a href="Multidimensional-Minimization.html">
Multidimensional Minimization</a></td> Multidimensional Minimization</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Least_002dSquares-Fitting.html"> <td><a href="Least_002dSquares-Fitting.html">
Least Squares Fitting</a></td> Least Squares Fitting</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Nonlinear-Least_002dSquares-Fitting.html"> <td><a href="Nonlinear-Least_002dSquares-Fitting.html">
Nonlinear Least Squares Fitting</a></td> Nonlinear Least Squares Fitting</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
<tr> <tr>
<td><a href="Basis-Splines.html">Basis Splines</a></td> <td><a href="Basis-Splines.html">Basis Splines</a></td>
<td>Not done</td> <td>Not started</td>
</tr> </tr>
</tbody> </tbody>
</table> </table>
......
;******************************************************** ;********************************************************
; file: library.lisp ; file: init.lisp
; description: Load GSL ; description: Load GSL
; date: Sat Mar 4 2006 - 18:53 ; date: Sat Mar 4 2006 - 18:53
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Fri Mar 24 2006 - 20:04 ; modified: Sat Sep 15 2007 - 18:59
;******************************************************** ;********************************************************
(defpackage gsll (defpackage gsll
......
...@@ -3,12 +3,36 @@ ...@@ -3,12 +3,36 @@
; description: Conversion of numbers C->CL ; description: Conversion of numbers C->CL
; date: Sun May 28 2006 - 22:04 ; date: Sun May 28 2006 - 22:04
; author: Liam M. Healy ; author: Liam M. Healy
; modified: Mon Jan 1 2007 - 12:10 ; modified: Sat Sep 15 2007 - 22:34
;******************************************************** ;********************************************************
;;; $Id: $ ;;; $Id: $
(in-package :gsl) (in-package :gsl)
;;;;****************************************************************************
;;;; Built-in numerical types
;;;;****************************************************************************
(defmacro float-to-cl (float &optional (index 0))
`(cffi:mem-aref ,float :float ,index))
(defmacro double-to-cl (double &optional (index 0))
`(cffi:mem-aref ,double :double ,index))
(defmacro size-to-cl (size &optional (index 0))
`(cffi:mem-aref ,size :size ,index))
(defmacro int-to-cl (integer &optional (index 0))
`(cffi:mem-aref ,integer :int ,index))
(defun cl-convert-function (type)
(case type
(:float 'float-to-cl)
(:double 'double-to-cl)
(:size 'size-to-cl)
(:int 'int-to-cl)
(gsl-complex 'complex-to-cl)))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Complex numbers ;;;; Complex numbers
;;;;**************************************************************************** ;;;;****************************************************************************
...@@ -25,61 +49,22 @@ ...@@ -25,61 +49,22 @@
gsl-complex gsl-complex
(* index (cffi:foreign-type-size 'gsl-complex))) (* index (cffi:foreign-type-size 'gsl-complex)))
'gsl-complex 'dat))) 'gsl-complex 'dat)))
(complex (mem-aref carr :double 0) (complex (double-to-cl carr 0)
(mem-aref carr :double 1)))) (double-to-cl carr 1))))
;;;;**************************************************************************** ;;;;****************************************************************************
;;;; Built-in numerical types ;;;; Conversion form
;;;;**************************************************************************** ;;;;****************************************************************************
(defun float-to-cl (float &optional (index 0))
(cffi:mem-aref float :float index))
(defun double-to-cl (double &optional (index 0))
(cffi:mem-aref double :double index))
(defun size-to-cl (size &optional (index 0))
(cffi:mem-aref size :size index))
(defun int-to-cl (integer &optional (index 0))
(cffi:mem-aref integer :int index))
(defun cl-convert-function (type)
(case type
(:float 'float-to-cl)
(:double 'double-to-cl)
(:size 'size-to-cl)
(:int 'int-to-cl)
(gsl-complex 'complex-to-cl)))
(defparameter *make-sequence-type* 'list
"Whether sequences should be returned as list or vector.")
(defun items-in-sequence (element-function length)
"Make a CL sequence of the type specified by *make-sequence-type*,
computing each element with the function element-function."
(let ((ans (make-sequence *make-sequence-type* length)))
(dotimes (i length ans)
(setf (elt ans i)
(funcall element-function i)))))
(defun cl-convert-form (decl) (defun cl-convert-form (decl)
"Generate a form that calls the appropriate converter from C/GSL to CL." "Generate a form that calls the appropriate converter from C/GSL to CL."
(if (st-arrayp decl) ; eventually, check that it's a vector (case (st-type decl)
`((map *make-sequence-type* (sf-result
(function `((val ,(st-symbol decl))
,(cl-convert-function (st-eltype decl))) (err ,(st-symbol decl))))
(cffi::foreign-array-to-lisp (sf-result-e10
,(st-symbol decl) `((val ,(st-symbol decl) 'sf-result-e10)
',(st-eltype decl) (e10 ,(st-symbol decl))
(list ,(st-dim decl))))) (err ,(st-symbol decl) 'sf-result-e10)))
(case (st-type decl) (t `((,(cl-convert-function (st-type decl)) ,(st-symbol decl))))))
(sf-result
`((val ,(st-symbol decl))
(err ,(st-symbol decl))))
(sf-result-e10
`((val ,(st-symbol decl) 'sf-result-e10)
(e10 ,(st-symbol decl))
(err ,(st-symbol decl) 'sf-result-e10)))
(t `((,(cl-convert-function (st-type decl)) ,(st-symbol decl)))))))
;********************************************************
; file: ordinary-differential-equations.lisp
; description: ODE initial value problems
; date: Sun Apr 15 2007 - 14:19
; author: Liam Healy
; modified: Sat Sep 15 2007 - 19:24
;********************************************************
;;; $Id: $
(in-package :gsl)
(defun-gsl
"gsl_odeiv_system"
)
(export 'def-ode-function)
(defmacro def-ode-function (name time dependent derivatives &body body)
"Define a function that will evaluate the right-hand sides (derivatives)
defining a set of ordinary differential equations (ODE).
The function should take as input the time (a double-float) and
dependent variables (a vector of double-floats) and fill the
derivatives vector with double-floats. It may refer to elements
of these vectors (arrays) using the macro double-to-cl.
This function may be passed to the GSL ODE integrators.
Parameters (non integration variables) may be passed by
using a lexical closure."
`(cffi:defcallback ,name :int
((,time :double)
(,dependent :pointer) ; This needs to be a C array
(,derivatives :pointer) ; This needs to be a C array
(params :pointer))
(declare (ignore params) (ignorable ,time))
,@body
;; Any errors or warnings should be signalled on the CL side;
;; if the function completes, we will always return success.
(cffi:foreign-enum-value 'gsl-errorno :SUCCESS)))
#|
(def-ode-function foo time y dydt
(setf (double-to-cl dydt 0) (- (double-to-cl y 1))
(double-to-cl dydt 1) (double-to-cl y 0)))
|#
;********************************************************
; file: simulated-annealing.lisp
; description: Simulated Annealing
; date: Sun Feb 11 2007 - 17:23
; author: Liam Healy
; modified: Sat Sep 15 2007 - 18:21
;********************************************************
;;; $Id: $
(in-package :gsl)
;;; There are some problems. First, step-size is not successfully
;;; returned by GSL to the step function, it comes back as garbage.
;;; That is not a big problem because it's never used by GSL, it is
;;; only sent back, so we can just define and use it locally. The
;;; main problem is that it seems to run forever, even on the
;;; "trivial" example. This is apparently because n-tries is
;;; not getting to GSL either. It seems the problem is with
;;; simulated-annealing-parameters.
;;; Also, the code could use quite a bit of clean
;;; up to present a simpler interface.
(cffi:defcstruct simulated-annealing-parameters
(n-tries :int) ; how many points to try for each step
(iterations-fixed-T :int) ; how many iterations at each temperature?
(step-size :double) ; max step size in the random walk
;; The following parameters are for the Boltzmann distribution
(k :double)
(t-initial :double)
(mu-t :double)
(t-min :double))
(defun-gsl simulated-annealing
(generator x0-p
Ef take-step distance-function
print-position
;; copy-function copy-constructor destructor
element-size parameters)
"gsl_siman_solve"
(((generator generator) :pointer) (x0-p :pointer)
((cffi:get-callback Ef) :pointer)
((cffi:get-callback take-step) :pointer)
((cffi:get-callback distance-function) :pointer)
((cffi:get-callback print-position) :pointer)
((cffi:null-pointer) :pointer)
((cffi:null-pointer) :pointer)
((cffi:null-pointer) :pointer)
;;((cffi:get-callback copy-function) :pointer)
;;((cffi:get-callback copy-constructor) :pointer)
;;((cffi:get-callback destructor) :pointer)
(element-size :size) (parameters simulated-annealing-parameters))
:c-return :void
:documentation
"Perform a simulated annealing search through a given
space. The space is specified by providing the functions @var{Ef} and
@var{distance}. The simulated annealing steps are generated using the
random number generator @var{r} and the function @var{take_step}.
The starting configuration of the system should be given by x0-p
The routine offers two modes for updating configurations, a fixed-size
mode and a variable-size mode. In the fixed-size mode the configuration
is stored as a single block of memory of size element-size
The functions copy-function,
copy-constructor and destructor should be NIL in
fixed-size mode. In the variable-size mode the functions
copy-function, copy-constructor and destructor are used to
create, copy and destroy configurations internally. The variable
@var{element_size} should be zero in the variable-size mode.
The parameters structure (described below) controls the run by
providing the temperature schedule and other tunable parameters to the
algorithm.
On exit the best result achieved during the search is placed in
x0-p. If the annealing process has been successful this
should be a good approximation to the optimal point in the space.
If the function pointer @var{print_position} is not null, a debugging
log will be printed to @code{stdout} with the following columns:
number_of_iterations temperature x x-x0p Ef(x)
and the output of the function @var{print_position} itself. If
@var{print_position} is null then no information is printed.
The simulated annealing routines require several user-specified
functions to define the configuration space and energy function.")
;; cribbed from def-gsl-function; unify?
(export 'def-sa-function)
(defmacro def-sa-function (name arg &body body)
"Define a GSL (C) function of either one argument of type
double (if arg is a symbol), or a C array of doubles
(if arg is a list), for GSL simulated annealing functions."
(let ((argvec (gensym "MCARG")))
`(cffi:defcallback ,name :double
(,(if (listp arg)
`(,argvec :pointer)
`(,arg :double)))
,@(if (listp arg)
`((symbol-macrolet
,(loop for i from 0 for a in arg
collect `(,a (cffi:mem-aref ,argvec :double ,i)))
,@body))
body))))
(defparameter *sa-function-calls* 0)
(defcallback e1 :double ((xp :pointer))
(incf *sa-function-calls*)
(when (> *sa-function-calls* 100000000)
(error "too much"))
(let ((x (cffi:mem-aref xp :double)))
(* (exp (- (expt (1- x) 2))) (sin (* 8 x)))))
(cffi:defcallback s1 :pointer
((generator :pointer) (parameters :pointer) (ss :double))
(declare (ignore ss))
(let ((step-size 10.0d0))
(let ((rand (uniform generator)))
(setf (cffi:mem-aref parameters :double)
(+ (cffi:mem-aref parameters :double)
(- (* 2 rand step-size) step-size))))
(cffi:null-pointer)))
(def-sa-function M1 (x y) (abs (- x y)))
(defcallback P1 :void ((xp :pointer))
(let ((x (cffi:mem-aref xp :double)))
(FORMAT T "~&from P1: ~a" X)))
(defun simulated-annealing-example ()
(cffi:with-foreign-object (initial :double)
(setf (cffi:mem-aref initial :double) 15.0d0)
(cffi:with-foreign-object (params 'simulated-annealing-parameters)
(setf
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'n-tries)
200
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'iterations-fixed-T)
10
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'step-size)
10.0d0
;; The following parameters are for the Boltzmann distribution
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'k)
1.0d0
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 't-initial)
0.002d0
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 'mu-t)
1.005d0
(cffi:foreign-slot-value
params 'simulated-annealing-parameters 't-min)
2.0d-6)
(simulated-annealing
*rng-mt19937* initial
'E1 'S1 'M1 'P1
(cffi:foreign-type-size :double)
params))))
#|
@deftypefun void
gsl_siman_solve
(const gsl_rng * @var{r},
void * @var{x0_p},
gsl_siman_Efunc_t @var{Ef},
gsl_siman_step_t @var{take_step},
gsl_siman_metric_t @var{distance},
gsl_siman_print_t @var{print_position},
gsl_siman_copy_t @var{copyfunc},
gsl_siman_copy_construct_t @var{copy_constructor},
gsl_siman_destroy_t @var{destructor},
size_t @var{element_size},
gsl_siman_params_t @var{params})
|#
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment