Commit 2309928b authored by liam's avatar liam

Added basis splines (untested), put documentation into a separate

directory and update it.

git-svn-id: svn+ssh://pop/opt/space/mathematics/gsl/trunk@3297 a3d8a0fb-c1db-0310-ace7-a616afeb9e30
parent cc753512
GSLL is distributed under the terms of the GNU General
Public License (GPL).
GSLL is distributed under the terms of the GNU Lesser General
Public License (LGPL), with the preamble defined by Franz, Inc.
to address specific issues about Lisp, except see below.
Most of the :documentation strings used as arguments to defun-gsl, and
Most of the :documentation strings used as arguments to defmfun, and
some comments in the Common Lisp source code are derived from the GSL
documentation. The following applies to these:
documentation. These are usually preceeded by a comment "; FDL"
in the source code. The following applies to these:
Permission is granted to copy, distribute and/or modify this
document under the terms of the GNU Free Documentation License,
......@@ -11,8 +14,7 @@ Version 1.2 or any later version published by the Free Software
Foundation; with the Invariant Sections being "GNU General Public
License" and "Free Software Needs Free Documentation", the
Front-Cover text being "A GNU Manual", and with the Back-Cover
Text being (a) (see below). A copy of the license is included in
the section entitled "GNU Free Documentation License."
Text being (a) (see below).
(a) The Back-Cover Text is: "You have freedom to copy and modify
this GNU Manual, like GNU software."
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2008-02-18 15:48:50EST basis-splines.lisp>
;; $Id: $
(in-package :gsl)
;;; Initializing the B-splines solver
(defgo-s (basis-spline order number-of-breakpoints)
allocate-basis-spline free-basis-spline nil 2)
(defmfun allocate-basis-spline (order number-of-breakpoints)
((order size) (number-of-breakpoints size))
:c-return :pointer
:export nil
:index (letm basis-spline)
:documentation ; FDL
"Allocate a workspace for computing B-splines. The number of
breakpoints is given by number-of-breakpoints. This leads to n =
nbreak + k - 2 basis functions where k = order. Cubic B-splines are
specified by k = 4. The size of the workspace is O(5k + nbreak)." )
(defmfun free-basis-spline (bspline)
((bspline :pointer))
:c-return :void
:export nil
:index (letm basis-spline)
:documentation ; FDL
"Free the memory associated with the workspace of the bspline.")
;;; Constructing the knots vector
(defmfun knots (breakpoints workspace)
((breakpoints :pointer) (workspace :pointer))
:documentation ; FDL
"Compute the knots associated with the given breakpoints and store
them in the workspace.")
(defmfun uniform-knots (a b workspace)
((a :double) (b :double) (workspace :pointer))
:documentation ; FDL
"Compute knots uniformly on the interval [a, b] and store
them in the workspace.")
;;; Evaluation of B-splines
(defmfun evaluate-bspline (x B workspace)
((x :double) (B :pointer) (workspace :pointer))
:documentation ; FDL
"Evaluate all B-spline basis functions at the position x and store
them in the GSL vector B, so that the ith element of B is B_i(x).
B must be of length n = nbreak + k - 2. This value is
also stored in the workspace. It is far more
efficient to compute all of the basis functions at once than to
compute them individually, due to the nature of the defining
recurrence relation.")
;;; Examples and unit test
(defun bspline-example (&optional (ncoeffs 8))
(letm ((order 4) ; cubic
(ndata 200)
(nbreak (+ ncoeffs 2 (- order)))
(bw (basis-spline order nbreak))
(mw (fit-workspace ndata ncoeffs))
(rng (random-number-generator *mt19937* 0))
(B (vector-double ncoeffs))
(c (vector-double ncoeffs))
(cov (matrix-double ncoeffs ncoeffs))
(w (vector-double ndata))
(x (vector-double ndata))
(y (vector-double ndata))
(Xmatrix (matrix-double ndata ncoeffs))
(sigma 0.1d0))
;; The data to be fitted.
(dotimes (i ndata)
(let* ((xi (* i (/ 15 (1- ndata))))
(yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
(gaussian rng sigma))))
(setf (gsl-aref x i) xi
(gsl-aref y i) yi
(gsl-aref w i) (/ (expt sigma 2)))))
;; Uniform breakpoints [0, 15]
(uniform-knots 0.0d0 15.0d0 bw)
;; Fit matrix
(dotimes (i ndata)
;; Compute B_j
(evaluate-bspline (gsl-aref x i) B bw)
;; Fill in row i of X
(dotimes (j ncoeffs)
(setf (gsl-aref Xmatrix i j) (gsl-aref B j))))
;; Do the fit
(weighted-linear-mfit Xmatrix w y c cov mw)
;; Return the smoothed curve
(loop for xi from 0.0d0 to 15.0d0 by 0.1d0
with yval and yerr
(evaluate-bspline xi B bw)
(multiple-value-setq (yval yerr)
(multi-linear-estimate B c cov))
collect yval into yvals
collect yerr into yerrs
finally (return (values yvals yerrs)))))
<!-- -*- mode: HTML; time-stamp-line-limit: -8; -*- -->
<?xml version="1.0" encoding="iso-8859-1"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "">
<html xmlns="">
<link rel="stylesheet" href="style.css" type="text/css"/>
<title>GSLL Development</title>
<div class="body">
<div class="header">
<h1>GSLL Development</h1>
<h2>GNU Scientific Library for Lisp Development Status</h2>
<p>The goal of this project is to provide as complete an interface
to the GSL as possible in as Lisp-natural a way as possible. It is
written only in Common Lisp (no C), with access to the foreign library
through CFFI.
The following is not part of the plan for the GSLL interface:
<li>Anything that is defined solely in C header files and not in the .so
library will not be a part of GSLL, because of the lack of access
through the library. For example, Chapter 38, Physical Constants,
are all defined as C macros.</li>
<li>Anything for which there are two near-equivalent GSL functions
and for which one is already provided. For example, all of the
special functions have two forms, one which provides the answer(s)
with an estimated error as a structure, and one which just provides
the answer as a double. This makes a
difference when C programming, but GSLL always uses the form with the
error estimate and returns the error as a later multiple value so
that it can easily be ignored.</li>
<li>Anything for which there is a Common Lisp
equivalent probably won't be done, unless there is a need to have it
present for some other GSL functionality that isn't present in CL.
<div class="content">
Development of GSLL will be tracked using Trac, to be set up.
<p>Most of the interface is complete, with a few chapters missing
and some functionality missing. See the <a
href="documentation.html">documentation</a> for details.
It is possible that the library API will change but probably not
substantially; most of the remaining work involves added functions.
In addition to the chapters below, it would be nice if the following
were done<ul>
<li>Native vectors like <code>sb-sys:vector-sap</code> in SBCL.
Contributions to the development of GSLL are welcome.
This make take the form of bug reports and fixes, suggestions for new
interface design, or of porting the remaining unported functions.
The supplied macros
<code>defmfun</code> and <code>defmpar</code>
define virtually all of the interface functions and variables
respectively, which simplifies the process of porting.
The arguments and defaults provide a convenient way of defining new
functions and variables.
<p>Where functions in the GSL library
do input and output to files, the remainder of the necessary
file-related operations must also be performed in C; no CL
interface is yet provided.</p>
<p>A Trac web page will soon be set up for those interested in
development, and for reporting bugs. In the meantime, send email to
me at this address: <code>lhealy</code> -at-
<code>common-lisp</code> -dot- <code>net</code>.</p>
<p>The following is the status of the project by <a
manual</a> chapter. When marked "Done" it means that all or a
substantial portion
that is planned that GSLL will cover has been defined and <i>some</i> of those
definitions have been tested. Improvements may still be
<h3>Regression tests</h3>
<div class="content">
<p>Regression tests using <a
(provided) are given for many of the functions defined. These
also serve as examples. Note that most of the result have not been verified
with an independent source; the values given are from previous
tests. To run the tests,
and expect a two failures on amd64, none? on i386.
<address><a href="mailto:">Liam Healy</a></address>
Time-stamp: <2008-02-18 00:25:04EST index.html>
<!-- Liam Healy 2008-02-18 11:36:11EST development.html -->
<!-- Time-stamp: <2008-02-18 12:50:26EST development.html> -->
<!-- $Id: $ -->
<!-- hhmts end -->
<script type="text/javascript">
var gaJsHost = (("https:" == document.location.protocol) ? "https://ssl." : "http://www.");
document.write(unescape("%3Cscript src='" + gaJsHost + "' type='text/javascript'%3E%3C/script%3E"));
<script type="text/javascript">
var pageTracker = _gat._getTracker("UA-3669275-1");
<!-- -*- mode: HTML; time-stamp-line-limit: -8; -*- -->
<?xml version="1.0" encoding="iso-8859-1"?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.1//EN" "">
<html xmlns="">
<link rel="stylesheet" href="style.css" type="text/css"/>
<title>GSLL Documentation</title>
<div class="body">
<div class="header">
<h1>GSLL Documentation</h1>
<h3>General Advice</h3>
<div class="content">
<p>There is little separate documentation for GSLL. Instead, the
following techniques for using the API are advised:</p>
<li>Find the appropriate function(s) in the <a
<li>Use the GSLL function <code>gsl-lookup</code>
to find the equivalent GSLL function, for example
(gsl-lookup "gsl_sf_elljac_e")
to find that the Lisp function name is <code>#'jacobian-elliptic-functions</code>.
<li>Look at the documentation for that Lisp function, e.g.
<pre>(documentation #'jacobian-elliptic-functions 'function)
"The Jacobian elliptic functions sn(u|m),
cn(u|m), dn(u|m) computed by descending Landen transformations."
to get an explanation of the arguments etc.
<li>Look at the end of the Lisp source file, or in an "example" file
in the same directory, for examples. Often, the examples are placed
in a <code>make-tests</code> form (which will be commented out with
a #| |# block).</li>
It is advisable to look at the examples first for calculations that
require more complex setup (generally, the later chapters in the GSL
manual) because often there are convenience macros defined to make
the setup simpler.
<h3>GSL objects and letm </h3>
There are a number of GSL structures available that must be manually
allocated, possibly set, and freed. In
order to make this as convenient as possible, the macro
<code>letm</code> is provided. This macro acts as a
with the additional feature that if the init-form is one of GSL objects
acceleration basis-spline block-complex block-double block-fixnum
block-single chebyshev combination complex-workspace discrete-random
eigen-herm eigen-hermv eigen-symm eigen-symmv fdfsolver fit-workspace
fminimizer fsolver hankel histogram integration-workspace
interpolation levin levin-truncated matrix-complex matrix-double
matrix-fixnum matrix-single mfdfminimizer mfdfsolver mfminimizer
mfsolver monte-carlo-miser monte-carlo-plain monte-carlo-vegas
nonlinear-fdffit nonlinear-ffit permutation
quasi-random-number-generator random-number-generator spline
vector-complex vector-double vector-fixnum vector-single wavelet
the appropriate object will be allocated and bound to the variable,
optionally set, and then freed when the body of the let is exited.
For Lisp environments with arglist prompting (such as SLIME), these
are functions whose symbols are exported so that an arglist prompt
will be visible; however, the function should not be used outside a
let binding.
<div class="content">
Consistent with the <a href="development.html">development</a>
philosophy, most of the interface to the library is done. Notes
particular chapters:
<li><a href="Vectors-and-Matrices.html">Vectors and Matrices</a>
are done for integer, real, double, complex.
Subvectors and views don't work. Other element types may be
difficult to port due to lack of CL definition.
<base href="">
<li>Some <a href="Sorting.html">sorting</a> functions fail on amd64 platforms,
perhaps because of lack of unsigned integer vectors.
<li><a href="BLAS-Support.html">BLAS</a> is completed, but with only very
limited testing.
<li><a href="Fast-Fourier-Transforms.html">FFTs</a> has not been done
because GSL does not supply an example and it is not clear how it
used. Contributions welcome.
<li><a href="Simulated-Annealing.html">Simulated Annealing</a> has
been done but does not work.
<li><a href="Wavelet-Transforms.html">Wavelet Transforms</a> has been
done but the results of the example do not agree with the GSL manual.
<li><a href="Discrete-Hankel-Transforms.html">Discrete Hankel
Transforms</a> has been done but there is no example given in the GSL
<li><a href="Basis-Splines.html">Basis splines</a> have yet to be tested.
<div class="footer">
<address><a href="mailto:">Liam Healy</a></address>
<!-- Created: Feb 25 2005 -->
<!-- hhmts start -->
Time-stamp: <2008-02-18 14:56:02EST documentation.html>
<!-- hhmts end -->
<script type="text/javascript">
var gaJsHost = (("https:" == document.location.protocol) ? "https://ssl." : "http://www.");
document.write(unescape("%3Cscript src='" + gaJsHost + "' type='text/javascript'%3E%3C/script%3E"));
<script type="text/javascript">
var pageTracker = _gat._getTracker("UA-3669275-1");
This diff is collapsed.
;; Definition of GSLL system
;; Liam Healy
;; Time-stamp: <2008-02-18 11:27:00EST gsll.asd>
;; Time-stamp: <2008-02-18 16:59:30EST gsll.asd>
;; $Id: $
(asdf:defsystem "gsll"
......@@ -179,4 +179,5 @@
(:file "roots-multi" :depends-on (roots-one))
(:file "minimization-multi")
(:file "linear-least-squares")
(:file "nonlinear-least-squares")))))
(:file "nonlinear-least-squares")))
(:file "basis-splines" :depends-on (init data))))
;; Definition of GSL objects and ways to use them.
;; Liam Healy, Sun Dec 3 2006 - 10:21
;; Time-stamp: <2008-02-16 10:35:20EST gsl-objects.lisp>
;; Time-stamp: <2008-02-18 11:46:20EST gsl-objects.lisp>
;; $Id: $
(in-package :gsl)
......@@ -14,6 +14,7 @@
(defvar *not-for-users* nil)
(export 'letm)
(defmacro letm (bindings &body body)
"A macro to bind variables to GSL objects, possibly setting
initial values, and freeing them when done. The binding is as
Markdown is supported
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment