From 5ae086f74f44e842f1f91efcf116dbe09de464b8 Mon Sep 17 00:00:00 2001
From: tony <tony@jeeves.blindglobe.net>
Date: Sun, 15 Jun 2008 08:15:24 +0200
Subject: [PATCH] small matlisp changes, big ones staged for later.

---
 matrix-matlisp.lisp    | 57 +++++++++++++++++++++++++-------------------------
 regression-matlisp.lsp | 15 +++++++------
 2 files changed, 35 insertions(+), 37 deletions(-)

diff --git a/matrix-matlisp.lisp b/matrix-matlisp.lisp
index ffc872c..7b8ab11 100644
--- a/matrix-matlisp.lisp
+++ b/matrix-matlisp.lisp
@@ -28,16 +28,17 @@ i;;; -*- mode: lisp -*-
 ;;;
 
 (defun make-sweep-front (x y w n p mode has_w x_mean result)
-  "X is: 
-y is: 
-w is =: 
-n is : 
-p is: 
-mode is: real, complex
-has_w :
+  "Compute the sweep matrix, a generalized computation based on the regression components Y and X. 
+X is: numerical version of design/regt matrix
+y is: response
+w is: weights (might be simply 1
+n is: ?
+p is: ?
+mode is: real, complex, though complex is not yet supported.
+has_w just indicates if weights (w) are meaningful.
 x_mean: 
-result:
-."
+result: ? place to store?
+"
   (declare (fixnum n p mode has_w))
   (let ((x_data nil)
 	(result_data nil)
@@ -62,8 +63,6 @@ result:
     (setf result_data (compound-data-seq result))
   
     ;; find the mean of y
-    (setf val 0.0)
-    (setf sum_w 0.0)
     (dotimes (i n)
       (declare (fixnum i))
       (setf dyi (makedouble (aref y i)))
@@ -76,7 +75,7 @@ result:
     (if (<= sum_w 0.0) (error "non positive sum of weights"))
     (setf dy_mean (/ val sum_w))
   
-    ;; find the column means
+    ;; find the column means of X
     (dotimes (j p)
       (declare (fixnum j))
       (setf val 0.0)
@@ -201,23 +200,23 @@ result:
   "Args: (x y &optional weights)
  X is matrix, Y and WEIGHTS are sequences. Returns the sweep matrix of the
  (weighted) regression of Y on X"
-   (check-matrix x)
-   (check-sequence y)
-   (if w (check-sequence w))
-   (let ((n (num-rows x))
-      	(p (num-cols x)))
-      (if (/= n (length y)) (error "dimensions do not match"))
-    (if (and w (/= n (length w))) (error "dimensions do not match"))
-    (let ((mode (max (la-data-mode x) 
-                     (la-data-mode x) 
-                     (if w (la-data-mode w) 0)))
-	  (result (make-array (list (+ p 2) (+ p 2))))
-	  (x-mean (make-array p))
-	  (y (coerce y 'vector))
- 	  (w (if w (coerce w 'vector)))
-	  (has-w (if w 1 0)))
-      (make-sweep-front x y w n p mode has-w x-mean result)
-      result)))
+   (if (not (typep x 'real-matrix)) (error "Make Sweep Matrix: X not real-matrix"))
+   (if (not (typep y 'real-matrix)) (error "Make Sweep Matrix: Y not real-matrix"))
+   (if w (if (not (typep w 'real-matrix)) (error "Make Sweep Matrix: W not real-matrix")))
+   (let ((n (n x))
+	 (p (m x)))
+     (if (/= n (length y)) (error "dimensions do not match"))
+     (if (and w (/= n (length w))) (error "dimensions do not match"))
+     (let ((mode (max (la-data-mode x) 
+		      (la-data-mode x) 
+		      (if w (la-data-mode w) 0)))
+	   (result (make-array (list (+ p 2) (+ p 2))))
+	   (x-mean (make-array p))
+	   (y (coerce y 'vector))
+	   (w (if w (coerce w 'vector)))
+	   (has-w (if w 1 0)))
+       (make-sweep-front x y w n p mode has-w x-mean result)
+       result)))
 
  (defun sweep-in-place (a k tol)
   (check-matrix a)
diff --git a/regression-matlisp.lsp b/regression-matlisp.lsp
index b1289da..512fb2b 100644
--- a/regression-matlisp.lsp
+++ b/regression-matlisp.lsp
@@ -15,8 +15,8 @@
 ;;;
 ;;; Incorporates modifications suggested by Sandy Weisberg.
 ;;;
-
-;;; This version uses matlisp rather than homebrewed linear algebra. 
+;;; Modified to use MATLISP as the matrix package, rather than
+;;; homebrewed linear algebra.
 
 (in-package :cl-user)
 
@@ -30,6 +30,8 @@
 	:lisp-stat-descriptive-statistics)
   (:shadowing-import-from :lisp-stat-object-system
 			  slot-value call-method call-next-method)
+  (:shadowing-import-from :matlisp
+			  real)
   (:shadowing-import-from :lisp-stat-math
       expt + - * / ** mod rem abs 1+ 1- log exp sqrt sin cos tan
       asin acos atan sinh cosh tanh asinh acosh atanh float random
@@ -104,15 +106,12 @@ Example (data are in file absorbtion.lsp in the sample data directory):
   (def m (regression-model (list iron aluminum) absorbtion))
   (send m :help) (send m :plot-residuals)"
   (let ((x (cond 
-	     ((typep x real-matrix) x)   ;; matrix, or ... 
-	     ((typep x 'vector) (list x)) ;; vector, or ...
-	     ((and (consp x)  ;; what?
-		   (numberp (car x))) (list x))
-	    (t x))) ;; defaulting...
+	     ((typep x 'real-matrix) x)   ;; matrix, or ... 
+	     (t (error "Regression:  X is not a real matrix")))) ;; defaulting...
         (m (send regression-model-proto :new)))
     (format t "~%")
     (send m :doc doc)
-    (send m :x (if (typep x real-matrix) x (apply #'bind-columns x)))
+    (send m :x x)
     (send m :y y)
     (send m :intercept intercept)
     (send m :weights weights)
-- 
2.11.4.GIT