Skip to content

Commit

Permalink
first set of changes from Brent Fulgham, a bunch more coming in the C…
Browse files Browse the repository at this point in the history
… code.
  • Loading branch information
blindglobe committed Sep 17, 2008
1 parent 981c085 commit 317cb33
Show file tree
Hide file tree
Showing 5 changed files with 78 additions and 38 deletions.
10 changes: 3 additions & 7 deletions lib/betabase.c
Expand Up @@ -27,17 +27,15 @@ double ppbeta(double p, double a, double b, int *ifault)
Static routines Static routines
*/ */


static double logbeta(p, q) static double logbeta(double p, double q)
double p, q;
{ {
return(gamma(p) + gamma(q) - gamma(p + q)); return(gamma(p) + gamma(q) - gamma(p + q));
} }


#define Min(x,y) (((x) < (y)) ? (x) : (y)) #define Min(x,y) (((x) < (y)) ? (x) : (y))
#define Max(x,y) (((x) > (y)) ? (x) : (y)) #define Max(x,y) (((x) > (y)) ? (x) : (y))


static double betai(x, pin, qin) static double betai(double x, double pin, double qin)
double x, pin, qin;
{ {
/* Translated from FORTRAN /* Translated from FORTRAN
july 1977 edition. w. fullerton, c3, los alamos scientific lab. july 1977 edition. w. fullerton, c3, los alamos scientific lab.
Expand Down Expand Up @@ -154,9 +152,7 @@ static double betai(x, pin, qin)
cause an underflow; a value of -308 or thereabouts will often be cause an underflow; a value of -308 or thereabouts will often be
*/ */


static double xinbta(p, q, beta, alpha, ifault) static double xinbta(double *p, double *q, double *beta, double *alpha, int *ifault)
double *p, *q, *beta, *alpha;
int *ifault;
{ {
/* Initialized data */ /* Initialized data */
static double sae = -30.0; /* this should be sufficient */ static double sae = -30.0; /* this should be sufficient */
Expand Down
27 changes: 27 additions & 0 deletions src/lisptests.lisp
@@ -0,0 +1,27 @@
(asdf:oos 'asdf:load-op 'lispstat)

(in-package :ls-user)

#|
(trace lisp-stat-linalg-data::la-allocate)
(trace lisp-stat-linalg-data::la-put-double)
(trace lisp-stat-linalg-data::la-data-to-matrix)
(trace lisp-stat-linalg-data::la-matrix)
(trace lisp-stat-linalg::lu-solve)
|#

(lu-solve
(lu-decomp #2A((2 3 4) (1 2 4) (2 4 5)))
#(2 3 4))
;; #(-2.333333333333333 1.3333333333333335 0.6666666666666666)


#|
(dotimes (i 3)
(declare (fixnum i))
(let ((vec (lisp-stat-linalg-data::la-get-pointer mat i)))
(format t "vec => ~A~%" vec)
(dotimes (j 3)
(format t "LA-PUT-DOUBLE ~A <- (~A ~A)~%" vec j (aref data i j))
(lisp-stat-linalg-data::la-put-double vec j (aref data i j)))))
|#
32 changes: 19 additions & 13 deletions src/numerics/ladata.lsp
Expand Up @@ -36,7 +36,13 @@


(in-package :lisp-stat-linalg-data) (in-package :lisp-stat-linalg-data)


#+openmcl
(defctype size-t :unsigned-long)
#+sbcl
(defctype size-t :unsigned-int)


;; Should we do the same with int's and long's? There is some
;; evidence that it might be useful?




;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Expand Down Expand Up @@ -95,7 +101,7 @@
(defun ptr-eq (p q) (cffi:pointer-eq p q)) (defun ptr-eq (p q) (cffi:pointer-eq p q))


(cffi:defcfun ("la_base_allocate" ccl-la-base-allocate) (cffi:defcfun ("la_base_allocate" ccl-la-base-allocate)
:pointer (n :int) (m :int)) :pointer (n size-t) (m size-t))
(defun la-base-allocate (n m) (defun la-base-allocate (n m)
(ccl-la-base-allocate n m)) (ccl-la-base-allocate n m))


Expand All @@ -105,7 +111,7 @@
(ccl-la-base-free-alloc p)) (ccl-la-base-free-alloc p))


(cffi:defcfun ("la_mode_size" ccl-la-mode-size) (cffi:defcfun ("la_mode_size" ccl-la-mode-size)
:int (x :int)) size-t (x size-t))


(defun la-mode-size (mode) (defun la-mode-size (mode)
(ccl-la-mode-size mode)) (ccl-la-mode-size mode))
Expand All @@ -114,13 +120,13 @@
;;; Callbacks for Internal Storage ;;; Callbacks for Internal Storage
;;; ;;;


(cffi:defcallback lisp-la-allocate :void ((n :long) (m :long)) (cffi:defcallback lisp-la-allocate :void ((n size-t) (m size-t))
(ccl-store-ptr (la-allocate n m))) (ccl-store-ptr (la-allocate n m)))
(cffi:defcfun ("register_la_allocate" register-la-allocate) (cffi:defcfun ("register_la_allocate" register-la-allocate)
:void (p :pointer)) :void (p :pointer))
(register-la-allocate (cffi:callback lisp-la-allocate)) (register-la-allocate (cffi:callback lisp-la-allocate))
(cffi:defcfun ("la_allocate" la) (cffi:defcfun ("la_allocate" la)
:pointer (x :int) (y :int)) :pointer (x size-t) (y size-t))


(cffi:defcallback lisp-la-free-alloc (cffi:defcallback lisp-la-free-alloc
:void ((p :pointer)) :void ((p :pointer))
Expand All @@ -138,53 +144,53 @@




(cffi:defcfun ("la_get_integer" ccl-la-get-integer) (cffi:defcfun ("la_get_integer" ccl-la-get-integer)
:int (p :pointer) (i :int)) :long (p :pointer) (i size-t)) ;; was int, not long, for first
(defun la-get-integer (p i) (defun la-get-integer (p i)
(ccl-la-get-integer p i)) (ccl-la-get-integer p i))


(cffi:defcfun ("la_get_double" ccl-la-get-double) (cffi:defcfun ("la_get_double" ccl-la-get-double)
:double (p :pointer) (i :int)) :double (p :pointer) (i size-t))
(defun la-get-double (p i) (defun la-get-double (p i)
(ccl-la-get-double p i)) (ccl-la-get-double p i))


(cffi:defcfun ("la_get_complex_real" ccl-la-get-complex-real) (cffi:defcfun ("la_get_complex_real" ccl-la-get-complex-real)
:double (p :pointer) (i :int)) :double (p :pointer) (i size-t))
(defun la-get-complex-real (p i) (defun la-get-complex-real (p i)
(ccl-la-get-complex-real p i)) (ccl-la-get-complex-real p i))


(cffi:defcfun ("la_get_complex_imag" ccl-la-get-complex-imag) (cffi:defcfun ("la_get_complex_imag" ccl-la-get-complex-imag)
:double (p :pointer) (i :int)) :double (p :pointer) (i size-t))
(defun la-get-complex-imag (p i) (defun la-get-complex-imag (p i)
(ccl-la-get-complex-imag p i)) (ccl-la-get-complex-imag p i))


(defun la-get-complex (p i) (defun la-get-complex (p i)
(complex (la-get-complex-real p i) (la-get-complex-imag p i))) (complex (la-get-complex-real p i) (la-get-complex-imag p i)))


(cffi:defcfun ("la_get_pointer" ccl-la-get-pointer) (cffi:defcfun ("la_get_pointer" ccl-la-get-pointer)
:pointer (p :pointer) (i :int)) :pointer (p :pointer) (i size-t))
(defun la-get-pointer (p i) (defun la-get-pointer (p i)
(ccl-la-get-pointer p i)) (ccl-la-get-pointer p i))




;;; CFFI glue for Storage Mutation Functions ;;; CFFI glue for Storage Mutation Functions


(cffi:defcfun ("la_put_integer" ccl-la-put-integer) (cffi:defcfun ("la_put_integer" ccl-la-put-integer)
:void (p :pointer) (i :int) (x :int)) :void (p :pointer) (i size-t) (x :long)) ;; last was :int ?
(defun la-put-integer (p i x) (defun la-put-integer (p i x)
(ccl-la-put-integer p i x)) (ccl-la-put-integer p i x))


(cffi:defcfun ("la_put_double" ccl-la-put-double) (cffi:defcfun ("la_put_double" ccl-la-put-double)
:void (p :pointer) (i :int) (x :double)) :void (p :pointer) (i size-t) (x :double))
(defun la-put-double (p i x) (defun la-put-double (p i x)
(ccl-la-put-double p i (float x 1d0))) (ccl-la-put-double p i (float x 1d0)))


(cffi:defcfun ("la_put_complex" ccl-la-put-complex) (cffi:defcfun ("la_put_complex" ccl-la-put-complex)
:void (p :pointer) (i :int) (x :double) (y :double)) :void (p :pointer) (i size-t) (x :double) (y :double))
(defun la-put-complex (p i x y) (defun la-put-complex (p i x y)
(ccl-la-put-complex p i (float x 1d0) (float y 1d0))) (ccl-la-put-complex p i (float x 1d0) (float y 1d0)))


(cffi:defcfun ("la_put_pointer" ccl-la-put-pointer) (cffi:defcfun ("la_put_pointer" ccl-la-put-pointer)
:void (p :pointer) (i :int) (q :pointer)) :void (p :pointer) (i size-t) (q :pointer))
(defun la-put-pointer (p i q) (defun la-put-pointer (p i q)
(ccl-la-put-pointer p i q)) (ccl-la-put-pointer p i q))


Expand Down
35 changes: 20 additions & 15 deletions src/numerics/linalg.lsp
Expand Up @@ -57,6 +57,11 @@


(in-package #:lisp-stat-linalg) (in-package #:lisp-stat-linalg)


#+openmcl
(defctype size-t :unsigned-long)
#+sbcl
(defctype size-t :unsigned-int)

;;; CFFI Support ;;; CFFI Support


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
Expand All @@ -70,7 +75,7 @@
;;; ;;;


(cffi:defcfun ("ccl_chol_decomp_front" ccl-chol-decomp-front) (cffi:defcfun ("ccl_chol_decomp_front" ccl-chol-decomp-front)
:int (x :pointer) (y :int) (z :pointer)) :int (x :pointer) (y size-t) (z :pointer))
(defun chol-decomp-front (x y z) (defun chol-decomp-front (x y z)
(ccl-chol-decomp-front x y z)) (ccl-chol-decomp-front x y z))


Expand All @@ -79,17 +84,17 @@
;;;; ;;;;


(cffi:defcfun ("ccl_lu_decomp_front" ccl-lu-decomp-front) (cffi:defcfun ("ccl_lu_decomp_front" ccl-lu-decomp-front)
:int (x :pointer) (y :int) (z :pointer) (u :int) (v :pointer)) :int (x :pointer) (y size-t) (z :pointer) (u :int) (v :pointer))
(defun lu-decomp-front (x y z u v) (defun lu-decomp-front (x y z u v)
(ccl-lu-decomp-front x y z u v)) (ccl-lu-decomp-front x y z u v))


(cffi:defcfun ("ccl_lu_solve_front" ccl-lu-solve-front) (cffi:defcfun ("ccl_lu_solve_front" ccl-lu-solve-front)
:int (x :pointer) (y :int) (z :pointer) (u :pointer) (v :int)) :int (x :pointer) (y size-t) (z :pointer) (u :pointer) (v :int))
(defun lu-solve-front (x y z u v) (defun lu-solve-front (x y z u v)
(ccl-lu-solve-front x y z u v)) (ccl-lu-solve-front x y z u v))


(cffi:defcfun ("ccl_lu_inverse_front" ccl-lu-inverse-front) (cffi:defcfun ("ccl_lu_inverse_front" ccl-lu-inverse-front)
:int (x :pointer) (y :int) (z :pointer) (u :pointer) (v :int) (w :pointer)) :int (x :pointer) (y size-t) (z :pointer) (u :pointer) (v :int) (w :pointer))
(defun lu-inverse-front (x y z u v w) (defun lu-inverse-front (x y z u v w)
(ccl-lu-inverse-front x y z u v w)) (ccl-lu-inverse-front x y z u v w))


Expand All @@ -98,7 +103,7 @@
;;;; ;;;;


(cffi:defcfun ("ccl_sv_decomp_front" ccl-sv-decomp-front) (cffi:defcfun ("ccl_sv_decomp_front" ccl-sv-decomp-front)
:int (x :pointer) (y :int) (z :int) (u :pointer) (v :pointer)) :int (x :pointer) (y size-t) (z size-t) (u :pointer) (v :pointer))
(defun sv-decomp-front (x y z u v) (defun sv-decomp-front (x y z u v)
(ccl-sv-decomp-front x y z u v)) (ccl-sv-decomp-front x y z u v))


Expand All @@ -107,7 +112,7 @@
;;;; ;;;;


(cffi:defcfun ("ccl_qr_decomp_front" ccl-qr-decomp-front) (cffi:defcfun ("ccl_qr_decomp_front" ccl-qr-decomp-front)
:int (x :pointer) (y :int) (z :int) (u :pointer) (v :pointer) (w :int)) :int (x :pointer) (y size-t) (z size-t) (u :pointer) (v :pointer) (w :int))
(defun qr-decomp-front (x y z u v w) (defun qr-decomp-front (x y z u v w)
(ccl-qr-decomp-front x y z u v w)) (ccl-qr-decomp-front x y z u v w))


Expand All @@ -116,7 +121,7 @@
;;; ;;;


(cffi:defcfun ("ccl_rcondest_front" ccl-rcondest-front) (cffi:defcfun ("ccl_rcondest_front" ccl-rcondest-front)
:double (x :pointer) (y :int)) :double (x :pointer) (y size-t))
(defun rcondest-front (x y) (defun rcondest-front (x y)
(ccl-rcondest-front x y)) (ccl-rcondest-front x y))


Expand Down Expand Up @@ -176,7 +181,7 @@
;;;; ;;;;


(cffi:defcfun ("ccl_make_rotation_front" ccl-make-rotation-front) (cffi:defcfun ("ccl_make_rotation_front" ccl-make-rotation-front)
:int (x :int) (y :pointer) (z :pointer) (u :pointer) (v :int) (w :double)) :int (x size-t) (y :pointer) (z :pointer) (u :pointer) (v :int) (w :double))
(defun make-rotation-front (x y z u v w) (defun make-rotation-front (x y z u v w)
(ccl-make-rotation-front x y z u v (float w 1d0))) (ccl-make-rotation-front x y z u v (float w 1d0)))


Expand All @@ -185,7 +190,7 @@
;;;; ;;;;


(cffi:defcfun ("ccl_eigen_front" ccl-eigen-front) (cffi:defcfun ("ccl_eigen_front" ccl-eigen-front)
:int (x :pointer) (y :int) (z :pointer) (u :pointer) (v :pointer)) :int (x :pointer) (y size-t) (z :pointer) (u :pointer) (v :pointer))
(defun eigen-front (x y z u v) (defun eigen-front (x y z u v)
(ccl-eigen-front x y z u v)) (ccl-eigen-front x y z u v))


Expand All @@ -194,12 +199,12 @@
;;;; ;;;;


(cffi:defcfun ("ccl_range_to_rseq" ccl-range-to-rseq) (cffi:defcfun ("ccl_range_to_rseq" ccl-range-to-rseq)
:int (x :int) (y :pointer) (z :int) (u :pointer)) :int (x size-t) (y :pointer) (z size-t) (u :pointer))
(defun la-range-to-rseq (x y z u) (defun la-range-to-rseq (x y z u)
(ccl-range-to-rseq x y z u)) (ccl-range-to-rseq x y z u))


(cffi:defcfun ("ccl_spline_front" ccl-spline-front) (cffi:defcfun ("ccl_spline_front" ccl-spline-front)
:int (x :int) (y :pointer) (z :pointer) (u :int) (v :pointer) (w :pointer) (a :pointer)) :int (x size-t) (y :pointer) (z :pointer) (u size-t) (v :pointer) (w :pointer) (a :pointer))
(defun spline-front (x y z u v w a) (defun spline-front (x y z u v w a)
(ccl-spline-front x y z u v w a)) (ccl-spline-front x y z u v w a))


Expand All @@ -208,12 +213,12 @@
;;;; ;;;;


(cffi:defcfun ("ccl_kernel_dens_front" ccl-kernel-dens-front) (cffi:defcfun ("ccl_kernel_dens_front" ccl-kernel-dens-front)
:int (x :pointer) (y :int) (z :double) (u :pointer) (v :pointer) (w :int) (a :int)) :int (x :pointer) (y size-t) (z :double) (u :pointer) (v :pointer) (w size-t) (a :int))
(defun kernel-dens-front (x y z u v w a) (defun kernel-dens-front (x y z u v w a)
(ccl-kernel-dens-front x y (float z 1d0) u v w a)) (ccl-kernel-dens-front x y (float z 1d0) u v w a))


(cffi:defcfun ("ccl_kernel_smooth_front" ccl-kernel-smooth-front) (cffi:defcfun ("ccl_kernel_smooth_front" ccl-kernel-smooth-front)
:int (x :pointer) (y :pointer) (z :int) (u :double) (v :pointer) (w :pointer) (a :int) (b :int)) :int (x :pointer) (y :pointer) (z size-t) (u :double) (v :pointer) (w :pointer) (a size-t) (b :int))
(defun kernel-smooth-front (x y z u v w a b) (defun kernel-smooth-front (x y z u v w a b)
(ccl-kernel-smooth-front x y z (float u 1d0) v w a b)) (ccl-kernel-smooth-front x y z (float u 1d0) v w a b))


Expand All @@ -222,7 +227,7 @@
;;;; ;;;;


(cffi:defcfun ("ccl_base_lowess_front" ccl-base-lowess-front) (cffi:defcfun ("ccl_base_lowess_front" ccl-base-lowess-front)
:int (x :pointer) (y :pointer) (z :int) (u :double) (v :int) (w :double) (a :pointer) (b :pointer) (c :pointer)) :int (x :pointer) (y :pointer) (z size-t) (u :double) (v size-t) (w :double) (a :pointer) (b :pointer) (c :pointer))
(defun base-lowess-front (x y z u v w a b c) (defun base-lowess-front (x y z u v w a b c)
(ccl-base-lowess-front x y z (float u 1d0) v (float w 1d0) a b c)) (ccl-base-lowess-front x y z (float u 1d0) v (float w 1d0) a b c))


Expand All @@ -231,7 +236,7 @@
;;;; ;;;;


(cffi:defcfun ("ccl_fft_front" ccl-fft-front) (cffi:defcfun ("ccl_fft_front" ccl-fft-front)
:int (x :int) (y :pointer) (z :pointer) (u :int)) :int (x size-t) (y :pointer) (z :pointer) (u :int))
(defun fft-front (x y z u) (defun fft-front (x y z u)
(ccl-fft-front x y z u)) (ccl-fft-front x y z u))


Expand Down
12 changes: 9 additions & 3 deletions src/numerics/optimize.lisp
Expand Up @@ -7,6 +7,7 @@


(defpackage :lisp-stat-optimize (defpackage :lisp-stat-optimize
(:use :common-lisp (:use :common-lisp
:cffi
:lisp-stat-ffi-int :lisp-stat-ffi-int
:lisp-stat-object-system :lisp-stat-object-system
:lisp-stat-types :lisp-stat-types
Expand Down Expand Up @@ -41,6 +42,11 @@


(in-package :lisp-stat-optimize) (in-package :lisp-stat-optimize)


#+openmcl
(defctype size-t :unsigned-long)
#+sbcl
(defctype size-t :unsigned-int)

(defvar *maximize-callback-function* nil (defvar *maximize-callback-function* nil
"Used in generic optimization to determine function name -- symbol or string?") "Used in generic optimization to determine function name -- symbol or string?")


Expand Down Expand Up @@ -68,17 +74,17 @@
(register-maximize-callback (cffi:callback ccl-maximize-callback)) (register-maximize-callback (cffi:callback ccl-maximize-callback))


(cffi:defcfun ("ccl_numgrad_front" ccl-numgrad-front) (cffi:defcfun ("ccl_numgrad_front" ccl-numgrad-front)
:int (x :int) (y :pointer) (z :pointer) (u :double) (v :pointer)) :void (x size-t) (y :pointer) (z :pointer) (u :double) (v :pointer))
(defun numgrad-front (x y z u v) (defun numgrad-front (x y z u v)
(ccl-numgrad-front x y z (float u 1d0) v)) (ccl-numgrad-front x y z (float u 1d0) v))


(cffi:defcfun ("ccl_numhess_front" ccl-numhess-front) (cffi:defcfun ("ccl_numhess_front" ccl-numhess-front)
:int (x :int) (y :pointer) (z :pointer) (u :pointer) (v :pointer) (w :double) (a :pointer)) :void (x size-t) (y :pointer) (z :pointer) (u :pointer) (v :pointer) (w :double) (a :pointer))
(defun numhess-front (x y z u v w a) (defun numhess-front (x y z u v w a)
(ccl-numhess-front x y z u v (float w 1d0) a)) (ccl-numhess-front x y z u v (float w 1d0) a))


(cffi:defcfun ("ccl_minfo_maximize" ccl-minfo-maximize) (cffi:defcfun ("ccl_minfo_maximize" ccl-minfo-maximize)
:int (x :pointer) (y :pointer) (z :pointer) (u :pointer) (v :pointer) (w :int)) :void (x :pointer) (y :pointer) (z :pointer) (u :pointer) (v :pointer) (w :int))
(defun base-minfo-maximize (x y z u v w) (defun base-minfo-maximize (x y z u v w)
(ccl-minfo-maximize x y z u v w)) (ccl-minfo-maximize x y z u v w))


Expand Down

0 comments on commit 317cb33

Please sign in to comment.