diff --git a/README b/README index 330aabcdb430dbdaa52ee98a295159c8e94580be..867105689f7bc76dba10712a03fb6d1e3e726c60 100644 --- a/README +++ b/README @@ -50,13 +50,7 @@ gsl_vector_short.h gsl_vector_uchar.h gsl_vector_char.h -I have implemented only one, double. - -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. - +I implemented only int, double, complex. It would be nice to use e.g. vector-sap for direct vectors in e.g. SBCL. diff --git a/index.html b/index.html index 4779e199dc3c97c4a3af66b11ee837d57fc57571..5eb47e674711dbe7b8c0210788e8a0872e901846 100644 --- a/index.html +++ b/index.html @@ -272,61 +272,61 @@ svn checkout svn://common-lisp.net/project/gsll/subversion/trunk</pre> </tr> <tr> <td><a href="Interpolation.html">Interpolation</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Numerical-Differentiation.html">Numerical Differentiation</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Chebyshev-Approximations.html">Chebyshev Approximations</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Series-Acceleration.html">Series Acceleration</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Wavelet-Transforms.html">Wavelet Transforms</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Discrete-Hankel-Transforms.html">Discrete Hankel Transforms</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="One-dimensional-Root_002dFinding.html">One-dimensional Root Finding</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="One-dimensional-Minimization.html">One-dimensional Minimization</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Multidimensional-Root_002dFinding.html"> Multidimensional Root Finding</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Multidimensional-Minimization.html"> Multidimensional Minimization</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Least_002dSquares-Fitting.html"> Least Squares Fitting</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Nonlinear-Least_002dSquares-Fitting.html"> Nonlinear Least Squares Fitting</a></td> - <td>Not done</td> + <td>Not started</td> </tr> <tr> <td><a href="Basis-Splines.html">Basis Splines</a></td> - <td>Not done</td> + <td>Not started</td> </tr> </tbody> </table> diff --git a/init/init.lisp b/init/init.lisp index d90dc97c04c0afbfcdc4f0203a851a57d0ce7dfd..28069bfa3c63003a8c5d314ec73e0557c4b73bdf 100644 --- a/init/init.lisp +++ b/init/init.lisp @@ -1,9 +1,9 @@ ;******************************************************** -; file: library.lisp +; file: init.lisp ; description: Load GSL ; date: Sat Mar 4 2006 - 18:53 ; author: Liam M. Healy -; modified: Fri Mar 24 2006 - 20:04 +; modified: Sat Sep 15 2007 - 18:59 ;******************************************************** (defpackage gsll diff --git a/init/number-conversion.lisp b/init/number-conversion.lisp index 4e284a8395d2c241e457049e5b12cc438fbd4bbb..d25ff8550f1a19cb7bfb09de07624c905c8a63c0 100644 --- a/init/number-conversion.lisp +++ b/init/number-conversion.lisp @@ -3,12 +3,36 @@ ; description: Conversion of numbers C->CL ; date: Sun May 28 2006 - 22:04 ; author: Liam M. Healy -; modified: Mon Jan 1 2007 - 12:10 +; modified: Sat Sep 15 2007 - 22:34 ;******************************************************** ;;; $Id: $ (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 ;;;;**************************************************************************** @@ -25,61 +49,22 @@ gsl-complex (* index (cffi:foreign-type-size 'gsl-complex))) 'gsl-complex 'dat))) - (complex (mem-aref carr :double 0) - (mem-aref carr :double 1)))) + (complex (double-to-cl carr 0) + (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) "Generate a form that calls the appropriate converter from C/GSL to CL." - (if (st-arrayp decl) ; eventually, check that it's a vector - `((map *make-sequence-type* - (function - ,(cl-convert-function (st-eltype decl))) - (cffi::foreign-array-to-lisp - ,(st-symbol decl) - ',(st-eltype decl) - (list ,(st-dim decl))))) - (case (st-type 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))))))) + (case (st-type 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)))))) diff --git a/ordinary-differential-equations.lisp b/ordinary-differential-equations.lisp new file mode 100644 index 0000000000000000000000000000000000000000..71d22e9730907a5ea3ef3160c1c3110e8fa61f91 --- /dev/null +++ b/ordinary-differential-equations.lisp @@ -0,0 +1,43 @@ +;******************************************************** +; 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))) +|# diff --git a/simulated-annealing.lisp b/simulated-annealing.lisp new file mode 100644 index 0000000000000000000000000000000000000000..4d96db51a677e1fad56bf0f4a121171e0116b7f6 --- /dev/null +++ b/simulated-annealing.lisp @@ -0,0 +1,178 @@ +;******************************************************** +; 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}) +|# +