Commit 0e6ff4a1 authored by Liam Healy's avatar Liam Healy
Browse files

Merge branch 'convert'

Trivial conflicts:

	data/array-structs.lisp
	init/init.lisp
parents 8c2cd8a6 fdb7a09f
;; Basis splines.
;; Liam Healy 2008-02-18 14:43:20EST basis-splines.lisp
;; Time-stamp: <2010-01-17 10:41:03EST basis-splines.lisp>
;; Time-stamp: <2010-06-30 19:57:28EDT basis-splines.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -56,7 +56,7 @@
(defmfun evaluate
((workspace basis-spline) x
&key (B (make-marray 'double-float
&key (B (grid:make-foreign-array 'double-float
:dimensions (number-of-coefficients workspace))))
"gsl_bspline_eval"
((x :double) ((mpointer B) :pointer) ((mpointer workspace) :pointer))
......@@ -112,31 +112,31 @@
(bw (make-basis-spline order nbreak))
(mw (make-fit-workspace ndata ncoeffs))
(rng (make-random-number-generator +mt19937+ 0))
(B (make-marray 'double-float :dimensions ncoeffs))
(c (make-marray 'double-float :dimensions ncoeffs))
(cov (make-marray 'double-float :dimensions (list ncoeffs ncoeffs)))
(w (make-marray 'double-float :dimensions ndata))
(x (make-marray 'double-float :dimensions ndata))
(y (make-marray 'double-float :dimensions ndata))
(Xmatrix (make-marray 'double-float :dimensions (list ndata ncoeffs)))
(B (grid:make-foreign-array 'double-float :dimensions ncoeffs))
(c (grid:make-foreign-array 'double-float :dimensions ncoeffs))
(cov (grid:make-foreign-array 'double-float :dimensions (list ncoeffs ncoeffs)))
(w (grid:make-foreign-array 'double-float :dimensions ndata))
(x (grid:make-foreign-array 'double-float :dimensions ndata))
(y (grid:make-foreign-array 'double-float :dimensions ndata))
(Xmatrix (grid:make-foreign-array 'double-float :dimensions (list ndata ncoeffs)))
(sigma 0.1d0))
;; The data to be fitted.
(dotimes (i ndata)
(let* ((xi (coerce (* i (/ 15 (1- ndata))) 'double-float))
(yi (+ (* (cos xi) (exp (* -0.1d0 xi)))
(sample rng :gaussian :sigma sigma))))
(setf (maref x i) xi
(maref y i) yi
(maref w i) (/ (expt sigma 2)))))
(setf (grid:gref x i) xi
(grid:gref y i) yi
(grid:gref 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 bw (maref x i) :b B)
(evaluate bw (grid:gref x i) :b B)
;; Fill in row i of X
(dotimes (j ncoeffs)
(setf (maref Xmatrix i j) (maref B j))))
(setf (grid:gref Xmatrix i j) (grid:gref B j))))
;; Do the fit
(linear-mfit Xmatrix y c w nil cov mw)
;; Return the smoothed curve
......
;; Monte Carlo Integration
;; Liam Healy Sat Feb 3 2007 - 17:42
;; Time-stamp: <2009-12-27 09:42:11EST monte-carlo.lisp>
;; Time-stamp: <2010-06-29 22:15:24EDT monte-carlo.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -51,7 +51,7 @@
(scalars t))
"gsl_monte_plain_integrate"
((callback :pointer)
((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
((foreign-pointer lower-limits) :pointer) ((foreign-pointer upper-limits) :pointer)
((dim0 lower-limits) sizet) (number-of-samples sizet)
((mpointer generator) :pointer)
((mpointer state) :pointer)
......@@ -112,7 +112,7 @@
(scalars t))
"gsl_monte_miser_integrate"
((callback :pointer)
((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
((foreign-pointer lower-limits) :pointer) ((foreign-pointer upper-limits) :pointer)
((dim0 lower-limits) sizet) (number-of-samples sizet)
((mpointer generator) :pointer)
((mpointer state) :pointer)
......@@ -173,7 +173,7 @@
(scalars t))
"gsl_monte_vegas_integrate"
((callback :pointer)
((c-pointer lower-limits) :pointer) ((c-pointer upper-limits) :pointer)
((foreign-pointer lower-limits) :pointer) ((foreign-pointer upper-limits) :pointer)
((dim0 lower-limits) sizet) (number-of-samples sizet)
((mpointer generator) :pointer)
((mpointer state) :pointer)
......@@ -211,7 +211,7 @@
(defparameter *mc-lower* #m(0.0d0 0.0d0 0.0d0))
(defparameter *mc-upper*
(make-marray 'double-float :initial-contents (list pi pi pi)))
(grid:make-foreign-array 'double-float :initial-contents (list pi pi pi)))
(defun random-walk-plain-example (&optional (nsamples 500000))
(monte-carlo-integrate-plain 'mcrw *mc-lower* *mc-upper* nsamples))
......
;; Numerical integration
;; Liam Healy, Wed Jul 5 2006 - 23:14
;; Time-stamp: <2009-12-27 09:42:09EST numerical-integration.lisp>
;; Time-stamp: <2010-07-07 14:21:55EDT numerical-integration.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -157,7 +157,7 @@
(limit 1000) (workspace (make-integration-workspace limit)))
"gsl_integration_qagp"
((callback :pointer)
((c-pointer points) :pointer) ((dim0 points) sizet)
((foreign-pointer points) :pointer) ((dim0 points) sizet)
(absolute-error :double) (relative-error :double) (limit sizet)
((mpointer workspace) :pointer)
(result (:pointer :double)) (abserr (:pointer :double)))
......@@ -367,6 +367,6 @@
(integration-QAGil (integration-test-myfn2 1.0d0) 1.0d0 1.0d-7 0.0d0 1000)
(integration-QAGp
'integration-test-f454
(copy (vector 0.0d0 1.0d0 (sqrt 2.0d0) 3.0d0) 'vector-double-float)
(grid:copy-to (vector 0.0d0 1.0d0 (sqrt 2.0d0) 3.0d0))
0.0d0 1.0d-3 1000)
(integration-QAWc 'integration-test-f459 -1.0d0 5.0d0 0.0d0 0.0d0 1.0d-3 1000))
;; Chebyshev Approximations
;; Liam Healy Sat Nov 17 2007 - 20:36
;; Time-stamp: <2010-04-26 09:48:32EDT chebyshev.lisp>
;; Time-stamp: <2010-06-29 21:47:34EDT chebyshev.lisp>
;;
;; Copyright 2007, 2008, 2009 Liam M. Healy
;; Copyright 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -60,7 +60,8 @@
:c-return (crtn :pointer)
:documentation
"The Chebyshev coefficient array as a CL array (foreign-friendly)."
:return ((c-array:copy-foreign-vector crtn (size chebyshev)))
:return ((grid:make-foreign-array-from-pointer
crtn (list (size chebyshev)) 'double-float t))
:gsl-version (1 12))
;;;;****************************************************************************
......
;; CFFI-Grovel definitions for unix systems.
;; Liam Healy 2009-06-06 09:36:29EDT array-structs.lisp
;; Time-stamp: <2010-05-23 11:36:42EDT array-structs.lisp>
;; Time-stamp: <2010-07-18 22:52:48EDT array-structs.lisp>
;;
;; Copyright 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -62,6 +62,6 @@
(include "gsl/gsl_combination.h")
(cstruct gsl-combination-c "gsl_combination"
(choice-of "n" :type sizet)
(range "n" :type sizet)
(size "k" :type sizet)
(data "data" :type :pointer))
;; Tests of array functions
;; Liam Healy 2008-10-20 22:41:48EDT array-tests.lisp
;; Time-stamp: <2009-12-27 09:42:06EST array-tests.lisp>
;; Time-stamp: <2010-07-07 14:25:02EDT array-tests.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -32,39 +32,33 @@
(generate-all-array-tests vector-set-all #+fsbv t #-fsbv :no-complex
(let ((v1 (array-default 3 t)))
(cl-array (set-all v1 (scalar-default)))))
(grid:copy-to (set-all v1 (scalar-default)))))
(generate-all-array-tests matrix-set-all #+fsbv t #-fsbv :no-complex
(let ((m1 (array-default '(3 3) t)))
(cl-array (set-all m1 (scalar-default)))))
(grid:copy-to (set-all m1 (scalar-default)))))
(generate-all-array-tests vector-set-zero t
(let ((v1 (array-default 3)))
(set-zero v1)
(cl-array v1)))
(grid:copy-to v1)))
(generate-all-array-tests matrix-set-zero t
(let ((m1 (array-default '(3 3))))
(set-zero m1)
(cl-array m1)))
(generate-all-array-tests vector-copy t
(cl-array (copy (array-default 3))))
(generate-all-array-tests matrix-copy t
(cl-array (copy (array-default '(3 3)))))
(grid:copy-to m1)))
(generate-all-array-tests vector-swap t
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(swap v2 v1)
(list (cl-array v1) (cl-array v2))))
(list (grid:copy-to v1) (grid:copy-to v2))))
(generate-all-array-tests matrix-swap t
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(swap m2 m1)
(list (cl-array m1) (cl-array m2))))
(list (grid:copy-to m1) (grid:copy-to m2))))
;;;;****************************************************************************
;;;; Arithmetic operations
......@@ -80,53 +74,53 @@
(generate-all-array-tests matrix-add :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt+ m1 m2))))
(grid:copy-to (elt+ m1 m2))))
(generate-all-array-tests vector-sub :no-complex
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(cl-array (elt- v1 v2))))
(grid:copy-to (elt- v1 v2))))
(generate-all-array-tests matrix-sub :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt- m1 m2))))
(grid:copy-to (elt- m1 m2))))
(generate-all-array-tests vector-mult :no-complex
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(cl-array (elt* v1 v2))))
(grid:copy-to (elt* v1 v2))))
(generate-all-array-tests matrix-mult :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt* m1 m2))))
(grid:copy-to (elt* m1 m2))))
(generate-all-array-tests vector-div :no-complex
(let ((v1 (array-default 3))
(v2 (array-default 3)))
(cl-array (elt/ v1 v2))))
(grid:copy-to (elt/ v1 v2))))
(generate-all-array-tests matrix-div :no-complex
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3))))
(cl-array (elt/ m1 m2))))
(grid:copy-to (elt/ m1 m2))))
(generate-all-array-tests vector-mult-scalar :no-complex
(let ((v1 (array-default 3)))
(cl-array (elt* v1 1.39d0))))
(grid:copy-to (elt* v1 1.39d0))))
(generate-all-array-tests matrix-mult-scalar :no-complex
(let ((m1 (array-default '(3 3))))
(cl-array (elt* m1 1.39d0))))
(grid:copy-to (elt* m1 1.39d0))))
(generate-all-array-tests vector-add-scalar :no-complex
(let ((v1 (array-default 3)))
(cl-array (elt+ v1 18.19d0))))
(grid:copy-to (elt+ v1 18.19d0))))
(generate-all-array-tests matrix-add-scalar :no-complex
(let ((m1 (array-default '(3 3))))
(cl-array (elt+ m1 18.19d0))))
(grid:copy-to (elt+ m1 18.19d0))))
;;;;****************************************************************************
;;;; Maximum and minimum elements
......@@ -189,16 +183,16 @@
(generate-all-array-tests set-basis t
(let ((v1 (array-default 8)))
(set-basis v1 2)
(cl-array v1)))
(grid:copy-to v1)))
(generate-all-array-tests swap-elements t
(let ((v1 (array-default 8)))
(swap-elements v1 2 5)
(cl-array v1)))
(grid:copy-to v1)))
(generate-all-array-tests vector-reverse t
(let ((v1 (array-default 8)))
(cl-array (vector-reverse v1))))
(grid:copy-to (vector-reverse v1))))
;;;;****************************************************************************
;;;; Matrix-only function definitions
......@@ -206,47 +200,47 @@
(generate-all-array-tests set-identity t
(let ((m1 (array-default '(3 3))))
(cl-array (set-identity m1))))
(grid:copy-to (set-identity m1))))
(generate-all-array-tests row t
(let ((m1 (array-default '(3 3)))
(row (array-default 3 t)))
(cl-array (row m1 1 row))))
(grid:copy-to (row m1 1 row))))
(generate-all-array-tests setf-row t
(let ((m1 (array-default '(3 3)))
(row (array-default 3)))
(setf (row m1 2) row)
(cl-array m1)))
(grid:copy-to m1)))
(generate-all-array-tests column t
(let ((m1 (array-default '(3 3)))
(col (array-default 3 t)))
(cl-array (column m1 1 col))))
(grid:copy-to (column m1 1 col))))
(generate-all-array-tests setf-column t
(let ((m1 (array-default '(3 3)))
(col (array-default 3)))
(setf (column m1 2) col)
(cl-array m1)))
(grid:copy-to m1)))
(generate-all-array-tests swap-rows t
(let ((m1 (array-default '(3 3))))
(cl-array (swap-rows m1 0 1))))
(grid:copy-to (swap-rows m1 0 1))))
(generate-all-array-tests swap-columns t
(let ((m1 (array-default '(3 3))))
(cl-array (swap-columns m1 1 2))))
(grid:copy-to (swap-columns m1 1 2))))
(generate-all-array-tests swap-row-column t
(let ((m1 (array-default '(3 3))))
(cl-array (swap-row-column m1 0 2))))
(grid:copy-to (swap-row-column m1 0 2))))
(generate-all-array-tests matrix-transpose* t
(let ((m1 (array-default '(3 3))))
(cl-array (matrix-transpose* m1))))
(grid:copy-to (matrix-transpose* m1))))
(generate-all-array-tests matrix-transpose t
(let ((m1 (array-default '(3 3)))
(m2 (array-default '(3 3) t)))
(cl-array (matrix-transpose m1 m2))))
(grid:copy-to (matrix-transpose m1 m2))))
;; Functions for both vectors and matrices.
;; Liam Healy 2008-04-26 20:48:44EDT both.lisp
;; Time-stamp: <2009-12-27 09:42:06EST both.lisp>
;; Time-stamp: <2010-07-13 21:46:11EDT both.lisp>
;;
;; Copyright 2008, 2009 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
......@@ -28,7 +28,7 @@
("gsl_" :category :type "_alloc_from_block")
((blockptr :pointer)
(0 sizet) ; offset
((size object) sizet) ; number of elements
((size object) sizet) ; number of elements
(1 sizet)) ; stride
:definition :generic
:c-return :pointer
......@@ -68,25 +68,6 @@
:c-return :void
:documentation "Set all elements to 0.")
(defmfun c-array:copy-to-destination
((source both) (destination both))
("gsl_" :category :type "_memcpy")
(((mpointer destination) :pointer) ((mpointer source) :pointer))
:definition :methods
:inputs (source)
:outputs (destination))
(defmfun c-array:copy-making-destination
((source both) &aux
(destination
(make-marray element-type :dimensions (dimensions source))))
("gsl_" :category :type "_memcpy")
(((mpointer destination) :pointer) ((mpointer source) :pointer))
:definition :methods
:inputs (source)
:outputs (destination)
:return (destination))
(defmfun swap ((a both) (b both))
("gsl_" :category :type "_swap")
(((mpointer a) :pointer) ((mpointer b) :pointer))
......@@ -97,33 +78,59 @@
"Exchange the elements of a and b
by copying. The two must have the same dimensions.")
;;; These functions are here as debugging tools only; they do not insure
;;; coherency between the C and CL sides and therefore should not be
;;; used in production.
;;; For matrix functions, there should be another index argument
#|
(defmfun set-value ((object vector) index value)
;;;;****************************************************************************
;;;; Array elements; used in callbacks scalarsp=T only
;;;;****************************************************************************
;;; Normal foreign array access is with grid:gref, but in order to
;;; avoid the overhead of instantiating a foreign-array object to
;;; access components, we use these functions.
(defmfun get-value ((class-name (eql vector)) mpointer &rest indices)
("gsl_" :category :type "_get")
((mpointer :pointer) ((first indices) sizet))
:definition :generic
:c-return :element-c-type
:export nil
:documentation
"Get the single element of the GSL vector. This is used
in callbacks.")
(defmfun get-value ((class-name (eql matrix)) mpointer &rest indices)
("gsl_" :category :type "_get")
((mpointer :pointer) ((first indices) sizet) ((second indices) sizet))
:definition :methods
:c-return :element-c-type
:export nil
:documentation
"Get the single element of the GSL matrix. This is used
in callbacks.")
(defmfun (setf get-value)
(value (class-name (eql vector)) mpointer &rest indices)
("gsl_" :category :type "_set")
(((mpointer object) :pointer) (index sizet) (value :element-c-type))
((value :element-c-type) (mpointer :pointer) ((first indices) sizet))
:definition :generic
:element-types #+fsbv t #-fsbv :no-complex
:inputs (object)
:outputs (object)
:c-return :void
:export nil ; debugging only
:return (value)
:export nil
:documentation
"Set single element to the value. For debugging only; do not use.")
"Set the single element of the GSL vector to the value. This is
used in callbacks.")
(defmfun get-value ((object vector) index)
("gsl_" :category :type "_get")
(((mpointer object) :pointer) (index sizet))
:definition :generic
:inputs (object)
:c-return :element-c-type
:export nil ; debugging only
(defmfun (setf get-value)
(value (class-name (eql matrix)) mpointer &rest indices)
("gsl_" :category :type "_set")
((value :element-c-type) (mpointer :pointer)
((first indices) sizet) ((second indices) sizet))
:definition :methods
:element-types #+fsbv t #-fsbv :no-complex
:c-return :void
:return (value)
:export nil
:documentation
"Set single element to the value. For debugging only; do not use.")
|#
"Set the single element of the GSL vector to the value. This is
used in callbacks.")
;;;;****************************************************************************
;;;; Elementwise arithmetic operations overwriting an array
......@@ -175,7 +182,7 @@
:documentation ; FDL
"Add the scalar complex x to all the elements of array a.")
(defmethod elt+ ((x float) (a marray))
(defmethod elt+ ((x float) (a foreign-array))
(elt+ a x))
(defmfun elt- ((a both) (b both))
......@@ -190,7 +197,7 @@
"Subtract the elements of b from the elements of a.
The two must have the same dimensions.")
(defmethod elt- ((a marray) (x float))
(defmethod elt- ((a foreign-array) (x float))
(elt+ a (- x)))
(defmfun elt* ((a vector) (b vector))
......@@ -235,7 +242,7 @@
:outputs (a)
:return (a))
(defmethod elt/ ((a marray) (x number))
(defmethod elt/ ((a foreign-array) (x number))
(elt* a (/ x)))
(defmfun elt* ((a both) (x float))
......@@ -262,7 +269,7 @@
:documentation ; FDL
"Multiply the elements of a by the scalar complex factor x.")
(defmethod elt* ((x float) (a marray))
(defmethod elt* ((x float) (a foreign-array))
(elt* a x))
;;;;****************************************************************************
......
;; Combinations
;; Liam Healy, Sun Mar 26 2006 - 11:51
;; Time-stamp: <2009-12-27 09:42:06EST combination.lisp>
;; Time-stamp: <2010-07-16 17:14:03EDT combination.lisp>
;;
;; Copyright 2006, 2007, 2008, 2009 Liam M. Healy
;; Copyright 2006, 2007, 2008, 2009, 2010 Liam M. Healy
;; Distributed under the terms of the GNU General Public License
;;
;; This program is free software: you can redistribute it and/or modify
......@@ -26,27 +26,22 @@
;;;; Combination structure and CL object
;;;;****************************************************************************
(defclass combination (mobject c-array:foreign-array)
((element-type
:initform
#+int64 '(unsigned-byte 64)
#+int32 '(unsigned-byte 32)
:reader element-type :allocation :class)
(choice-of :initarg :choice-of :reader choice-of :type (integer 0)
:documentation "Maximum possible value; n in the (n k) notation."))
(defclass combination
(#+int64 grid:vector-unsigned-byte-64 #+int32 grid:vector-unsigned-byte-32)
()
(:documentation "GSL combinations."))
(defmethod initialize-instance :after
((object combination) &key choice-of dimensions &allow-other-keys)
((object combination) &key range dimensions &allow-other-keys)
(let ((mptr (cffi:foreign-alloc 'gsl-combination-c)))
(setf (slot-value object 'mpointer)
(setf (grid:metadata-slot object 'mpointer)
mptr
(cffi:foreign-slot-value mptr 'gsl-combination-c 'data)
(c-pointer object)
(cffi:foreign-slot-value mptr 'gsl-combination-c 'choice-of)
choice-of
(foreign-pointer object)
(cffi:foreign-slot-value mptr 'gsl-combination-c 'range)
range
(cffi:foreign-slot-value mptr 'gsl-combination-c 'size)
dimensions)
(first dimensions))
(tg:finalize object (lambda () (cffi:foreign-free mptr)))))
(export 'make-combination)
......@@ -59,14 +54,24 @@
(let ((comb
(if (typep n 'combination)
(make-instance
'combination :choice-of (choice-of n) :dimensions (dimensions k))
(make-instance 'combination :choice-of n :dimensions k))))
'combination
:element-type '(unsigned-byte #+int64 64 #+int32 32)
:range (combination-range n) :dimensions (dimensions k))
(make-instance
'combination
:element-type '(unsigned-byte #+int64 64 #+int32 32)
:range n :dimensions (list k)))))
(when initialize
(if (typep n 'combination)
(copy comb n)
(error "not available yet") ; (copy comb n)
(init-first comb)))
comb))
(defmethod print-object ((object combination) stream)
(print-unreadable-object (object stream :type t)
(format stream "range ~d: " (combination-range object))
(princ (grid:contents object) stream)))
;;;;****************************************************************************
;;;; Setting values
;;;;****************************************************************************
......@@ -89,19 +94,28 @@
"Initialize the combination c to the lexicographically
last combination, i.e. (n-k,n-k+1,...,n-1).")
(defmfun c-array:copy-to-destination ((source combination) (destination combination))
(defmfun comb-copy (source destination)
"gsl_combination_memcpy"
(((mpointer destination) :pointer)
((mpointer source) :pointer))
:definition :method
:inputs (source)
:outputs (destination)
:return (destination)
:index copy
:export nil
:index grid:copy
:documentation ; FDL
"Copy the elements of the combination source into the
combination destination. The two combinations must have the same size.")
(defmethod grid:copy
((source combination) &key grid-type destination &allow-other-keys)
(if grid-type
(call-next-method)
(comb-copy
source