From 1a7de5d214e6eec02f229f1b8f8758a65b405cbf Mon Sep 17 00:00:00 2001 From: AJ Rossini Date: Wed, 17 Oct 2007 09:11:31 +0200 Subject: [PATCH] 2 fixes from Carlos Ungil changed mode-?? to +mode-??+ in accordance with CL coding practices. reverted smoother to use C version, now that CFFI support is present. --- linalg.lsp | 160 +++++++++++++++++++++++++++++++------------------------------ 1 file changed, 82 insertions(+), 78 deletions(-) diff --git a/linalg.lsp b/linalg.lsp index 18cd7be..01c87a9 100644 --- a/linalg.lsp +++ b/linalg.lsp @@ -90,13 +90,13 @@ positive definite. Returns a list (L (max D))." (result (make-array (list n n))) (dpars (list maxoffl 0.0))) (check-real dpars) - (let ((mat (la-data-to-matrix a mode-re)) - (dp (la-data-to-vector dpars mode-re))) + (let ((mat (la-data-to-matrix a +mode-re+)) + (dp (la-data-to-vector dpars +mode-re+))) (unwind-protect (progn (chol-decomp-front mat n dp) - (la-matrix-to-data mat n n mode-re result) - (la-vector-to-data dp 2 mode-re dpars)) + (la-matrix-to-data mat n n +mode-re+ result) + (la-vector-to-data dp 2 +mode-re+ dpars)) (la-free-matrix mat n) (la-free-vector dp))) (list result (second dpars)))) @@ -121,17 +121,17 @@ is odd, -1 if even, and FLAG is T if A is numerically singular, NIL otherwise. Used bu LU-SOLVE." (check-square-matrix a) (let* ((n (array-dimension a 0)) - (mode (max mode-re (la-data-mode a))) + (mode (max +mode-re+ (la-data-mode a))) (result (list (make-array (list n n)) (make-array n) nil nil))) (let ((mat (la-data-to-matrix a mode)) - (iv (la-vector n mode-in)) - (d (la-vector 1 mode-re)) + (iv (la-vector n +mode-in+)) + (d (la-vector 1 +mode-re+)) (singular 0)) (unwind-protect (progn (setf singular (lu-decomp-front mat n iv mode d)) (la-matrix-to-data mat n n mode (first result)) - (la-vector-to-data iv n mode-in (second result)) + (la-vector-to-data iv n +mode-in+ (second result)) (setf (third result) (la-get-double d 0)) (setf (fourth result) (if (= singular 0.0) nil t))) (la-free-matrix mat n) @@ -156,9 +156,9 @@ singular." (b-mode (la-data-mode lb))) (if (/= n (length lidx)) (error "index sequence is wrong length")) (if (/= n (length lb)) (error "right hand side is wrong length")) - (let* ((mode (max mode-re a-mode b-mode)) + (let* ((mode (max +mode-re+ a-mode b-mode)) (a (la-data-to-matrix la mode)) - (indx (la-data-to-vector lidx mode-in)) + (indx (la-data-to-vector lidx +mode-in+)) (b (la-data-to-vector lb mode)) (singular 0)) (unwind-protect @@ -195,7 +195,7 @@ Returns the inverse of the the square matrix M; signals an error if M is ill conditioned or singular" (check-square-matrix a) (let ((n (num-rows a)) - (mode (max mode-re (la-data-mode a)))) + (mode (max +mode-re+ (la-data-mode a)))) (declare (fixnum n)) (let ((result (make-array (list n n) :initial-element 0))) (dotimes (i n) @@ -203,7 +203,7 @@ conditioned or singular" (setf (aref result i i) 1)) (let ((mat (la-data-to-matrix a mode)) (inv (la-data-to-matrix result mode)) - (iv (la-vector n mode-in)) + (iv (la-vector n +mode-in+)) (v (la-vector n mode)) (singular 0)) (unwind-protect @@ -231,16 +231,16 @@ if the algorithm converged, NIL otherwise." (check-matrix a) (let* ((m (num-rows a)) (n (num-cols a)) - (mode (max mode-re (la-data-mode a))) + (mode (max +mode-re+ (la-data-mode a))) (result (list (make-array (list m n)) (make-array n) (make-array (list n n)) nil))) (if (< m n) (error "number of rows less than number of columns")) - (if (= mode mode-cx) (error "complex SVD not available yet")) + (if (= mode +mode-cx+) (error "complex SVD not available yet")) (let ((mat (la-data-to-matrix a mode)) - (w (la-vector n mode-re)) - (v (la-matrix n n mode-re)) + (w (la-vector n +mode-re+)) + (v (la-matrix n n +mode-re+)) (converged 0)) (unwind-protect (progn @@ -270,7 +270,7 @@ the order in which they were used." (check-matrix a) (let* ((m (num-rows a)) (n (num-cols a)) - (mode (max mode-re (la-data-mode a))) + (mode (max +mode-re+ (la-data-mode a))) (p (if pivot 1 0)) (result (if pivot (list (make-array (list m n)) @@ -278,16 +278,16 @@ the order in which they were used." (make-array n)) (list (make-array (list m n)) (make-array (list n n)))))) (if (< m n) (error "number of rows less than number of columns")) - (if (= mode mode-cx) (error "complex QR decomposition not available yet")) + (if (= mode +mode-cx+) (error "complex QR decomposition not available yet")) (let ((mat (la-data-to-matrix a mode)) - (v (la-matrix n n mode-re)) - (jpvt (la-vector n mode-in))) + (v (la-matrix n n +mode-re+)) + (jpvt (la-vector n +mode-in+))) (unwind-protect (progn (qr-decomp-front mat m n v jpvt p) (la-matrix-to-data mat m n mode (first result)) (la-matrix-to-data v n n mode (second result)) - (if pivot (la-vector-to-data jpvt n mode-in (third result)))) + (if pivot (la-vector-to-data jpvt n +mode-in+ (third result)))) (la-free-matrix mat m) (la-free-matrix v n) (la-free-vector jpvt)) @@ -302,9 +302,9 @@ the order in which they were used." Returns an estimate of the reciprocal of the L1 condition number of an upper triangular matrix a." (check-square-matrix a) - (let ((mode (max mode-re (la-data-mode a))) + (let ((mode (max +mode-re+ (la-data-mode a))) (n (num-rows a))) - (if (= mode mode-cx) + (if (= mode +mode-cx+) (error "complex condition estimate not available yet")) (let ((mat (la-data-to-matrix a mode)) (est 0.0)) @@ -325,19 +325,19 @@ by angle ALPHA, in radians. X and Y are sequences of the same length." (check-sequence y) (if alpha (check-one-real alpha)) (let* ((n (length x)) - (mode (max mode-re (la-data-mode x) (la-data-mode y))) + (mode (max +mode-re+ (la-data-mode x) (la-data-mode y))) (use-angle (if alpha 1 0)) (angle (if alpha (float alpha 0.0) 0.0)) (result (make-array (list n n)))) (if (/= n (length y)) (error "sequences not the same length")) - (if (= mode mode-cx) (error "complex data not supported yet")) - (let ((px (la-data-to-vector x mode-re)) - (py (la-data-to-vector y mode-re)) - (rot (la-matrix n n mode-re))) + (if (= mode +mode-cx+) (error "complex data not supported yet")) + (let ((px (la-data-to-vector x +mode-re+)) + (py (la-data-to-vector y +mode-re+)) + (rot (la-matrix n n +mode-re+))) (unwind-protect (progn (make-rotation-front n rot px py use-angle angle) - (la-matrix-to-data rot n n mode-re result)) + (la-matrix-to-data rot n n +mode-re+ result)) (la-free-vector px) (la-free-vector py) (la-free-matrix rot n)) @@ -354,21 +354,21 @@ symmetric matrix A. Third element of result is NIL if algorithm converges. If the algorithm does not converge, the third element is an integer I. In this case the eigenvalues 0, ..., I are not reliable." (check-square-matrix a) - (let ((mode (max mode-re (la-data-mode a))) + (let ((mode (max +mode-re+ (la-data-mode a))) (n (num-rows a))) - (if (= mode mode-cx) (error "matrix must be real and symmetric")) + (if (= mode +mode-cx+) (error "matrix must be real and symmetric")) (let ((evals (make-array n)) (evecs (make-list (* n n))) - (pa (la-data-to-vector (compound-data-seq a) mode-re)) - (w (la-vector n mode-re)) - (z (la-vector (* n n) mode-re)) - (fv1 (la-vector n mode-re)) + (pa (la-data-to-vector (compound-data-seq a) +mode-re+)) + (w (la-vector n +mode-re+)) + (z (la-vector (* n n) +mode-re+)) + (fv1 (la-vector n +mode-re+)) (ierr 0)) (unwind-protect (progn (setf ierr (eigen-front pa n w z fv1)) - (la-vector-to-data w n mode-re evals) - (la-vector-to-data z (* n n) mode-re evecs)) + (la-vector-to-data w n +mode-re+ evals) + (la-vector-to-data z (* n n) +mode-re+ evecs)) (la-free-vector pa) (la-free-vector z) (la-free-vector w) @@ -417,17 +417,17 @@ In this case the eigenvalues 0, ..., I are not reliable." (error "sequences not the same length")) (if (and (not supplied) (< ns 2)) (error "too few points for interpolation")) - (let* ((px (la-data-to-vector ,x mode-re)) - (py (if ,is-reg (la-data-to-vector ,y mode-re))) + (let* ((px (la-data-to-vector ,x +mode-re+)) + (py (if ,is-reg (la-data-to-vector ,y +mode-re+))) (pxs (if supplied - (la-data-to-vector ,xvals mode-re) - (la-vector ns mode-re))) - (pys (la-vector ns mode-re))) + (la-data-to-vector ,xvals +mode-re+) + (la-vector ns +mode-re+))) + (pys (la-vector ns +mode-re+))) (unless supplied (la-range-to-rseq n px ns pxs)) (unwind-protect (progn ,@body - (la-vector-to-data pxs ns mode-re (first result)) - (la-vector-to-data pys ns mode-re (second result))) + (la-vector-to-data pxs ns +mode-re+ (first result)) + (la-vector-to-data pys ns +mode-re+ (second result))) (la-free-vector px) (if ,is-reg (la-free-vector py)) (la-free-vector pxs) @@ -441,7 +441,7 @@ X must be strictly increasing. XVALS can be an integer, the number of equally spaced points to use in the range of X, or it can be a sequence of points at which to interpolate." (with-smoother-data (x y xvals t) - (let ((work (la-vector (* 2 n) mode-re)) + (let ((work (la-vector (* 2 n) +mode-re+)) (error 0)) (unwind-protect (setf error (spline-front n px py ns pxs pys work)) @@ -483,8 +483,11 @@ U or B for Gaussian, triangular, uniform or bisquare. The default is B." (with-smoother-data (x y xvals t) (let ((code (kernel-type-code type)) (error 0)) - ;;(kernel-smooth-front px py n width pxs pys ns code) - (kernel-smooth-Cport px py n width pxs pys ns code) + (kernel-smooth-front px py n width pxs pys ns code) + ;; if we get the Lisp version ported from C, uncomment below and + ;; comment above. (thanks to Carlos Ungil for the initial CFFI + ;; work). + ;;(kernel-smooth-Cport px py n width pxs pys ns code) (if (/= 0 error) (error "bad kernel density data"))))) @@ -492,7 +495,8 @@ U or B for Gaussian, triangular, uniform or bisquare. The default is B." (defun kernel-smooth-Cport (px py n width ;;wts wds ;; see above for mismatch? xs ys ns ktype) "Port of kernel_smooth (Lib/kernel.c) to Lisp. -FIXME:kernel-smooth-Cport" +FIXME:kernel-smooth-Cport : This is broken. +Until this is fixed, we are using Luke's C code and CFFI as glue." (cond ((< n 1) 1.0) ((and (< n 2) (<= width 0)) 1.0) (t (let* ((xmin (min px)) @@ -575,16 +579,16 @@ x,y,w are doubles, type is an integer" (let* ((n (length s1)) (result (make-list n))) (if (/= n (length s2)) (error "sequences not the same length")) - (let ((x (la-data-to-vector s1 mode-re)) - (y (la-data-to-vector s2 mode-re)) - (ys (la-vector n mode-re)) - (rw (la-vector n mode-re)) - (res (la-vector n mode-re)) + (let ((x (la-data-to-vector s1 +mode-re+)) + (y (la-data-to-vector s2 +mode-re+)) + (ys (la-vector n +mode-re+)) + (rw (la-vector n +mode-re+)) + (res (la-vector n +mode-re+)) (error 0)) (unwind-protect (progn (setf error (base-lowess-front x y n f nsteps delta ys rw res)) - (la-vector-to-data ys n mode-re result)) + (la-vector-to-data ys n +mode-re+ result)) (la-free-vector x) (la-free-vector y) (la-free-vector ys) @@ -675,12 +679,12 @@ is true." (mode (la-data-mode x)) (isign (if inverse -1 1)) (result (if (consp x) (make-list n) (make-array n)))) - (let ((px (la-data-to-vector x mode-cx)) - (work (la-vector (+ (* 4 n) 15) mode-re))) + (let ((px (la-data-to-vector x +mode-cx+)) + (work (la-vector (+ (* 4 n) 15) +mode-re+))) (unwind-protect (progn (fft-front n px work isign) - (la-vector-to-data px n mode-cx result)) + (la-vector-to-data px n +mode-cx+ result)) (la-free-vector px) (la-free-vector work)) result))) @@ -966,7 +970,7 @@ This could probably be made more efficient." (la-put-double ptr i (get-next-element elem i))))) (defun maximize-callback (n px pfval pgrad phess pderivs) - (la-vector-to-data px n mode-re *maximize-callback-arg*) + (la-vector-to-data px n +mode-re+ *maximize-callback-arg*) (let* ((val (funcall *maximize-callback-function* *maximize-callback-arg*)) (derivs (if (consp val) (- (length val) 1) 0))) (la-put-integer pderivs 0 derivs) @@ -989,15 +993,15 @@ Computes the numerical gradient of F at X." (error "scale not the same length as x")) (let ((*maximize-callback-function* f) (*maximize-callback-arg* (make-list n))) - (let ((px (la-data-to-vector x mode-re)) - (pgrad (la-vector n mode-re)) + (let ((px (la-data-to-vector x +mode-re+)) + (pgrad (la-vector n +mode-re+)) (pscale (la-data-to-vector (if scale scale (make-list n :initial-element 1.0)) - mode-re))) + +mode-re+))) (unwind-protect (progn (numgrad-front n px pgrad h pscale) - (la-vector-to-data pgrad n mode-re result)) + (la-vector-to-data pgrad n +mode-re+ result)) (la-free-vector px) (la-free-vector pgrad) (la-free-vector pscale)))) @@ -1021,20 +1025,20 @@ Computes the numerical Hessian matrix of F at X." (let ((*maximize-callback-function* f) (*maximize-callback-arg* (make-list n))) (let ((hess-data (compound-data-seq (if all (third result) result))) - (px (la-data-to-vector x mode-re)) - (pf (la-vector 1 mode-re)) - (pgrad (la-vector n mode-re)) - (phess (la-vector (* n n) mode-re)) + (px (la-data-to-vector x +mode-re+)) + (pf (la-vector 1 +mode-re+)) + (pgrad (la-vector n +mode-re+)) + (phess (la-vector (* n n) +mode-re+)) (pscale (la-data-to-vector (if scale scale (make-list n :initial-element 1.0)) - mode-re))) + +mode-re+))) (unwind-protect (progn (numhess-front n px pf pgrad phess h pscale) (when all (setf (first result) (la-get-double pf 0)) - (la-vector-to-data pgrad n mode-re (second result))) - (la-vector-to-data phess (* n n) mode-re hess-data)) + (la-vector-to-data pgrad n +mode-re+ (second result))) + (la-vector-to-data phess (* n n) +mode-re+ hess-data)) (la-free-vector pf) (la-free-vector px) (la-free-vector pgrad) @@ -1126,18 +1130,18 @@ Computes the numerical Hessian matrix of F at X." (fvals (aref internals 5)) (ip (aref internals 8)) (dp (aref internals 9)) - (px (la-data-to-vector x mode-re)) - (pscale (la-data-to-vector scale mode-re)) - (pfvals (la-vector (length fvals) mode-re)) - (pip (la-data-to-vector ip mode-in)) - (pdp (la-data-to-vector dp mode-re))) + (px (la-data-to-vector x +mode-re+)) + (pscale (la-data-to-vector scale +mode-re+)) + (pfvals (la-vector (length fvals) +mode-re+)) + (pip (la-data-to-vector ip +mode-in+)) + (pdp (la-data-to-vector dp +mode-re+))) (unwind-protect (progn (base-minfo-maximize px pfvals pscale pip pdp v)) ;; access to C - (la-vector-to-data px n mode-re x) - (la-vector-to-data pfvals (+ 1 n (* n n)) mode-re fvals) - (la-vector-to-data pip (length ip) mode-in ip) - (la-vector-to-data pdp (length dp) mode-re dp)) + (la-vector-to-data px n +mode-re+ x) + (la-vector-to-data pfvals (+ 1 n (* n n)) +mode-re+ fvals) + (la-vector-to-data pip (length ip) +mode-in+ ip) + (la-vector-to-data pdp (length dp) +mode-re+ dp)) (get-buf))))) ;;;; -- 2.11.4.GIT