(setf eps 1e-10)

(setf x '#2a((1 2)(3 4)(5 6)(7 8)))
(setf y '#2a((9 10 11) (12 13 14)))
(setf xv '#(15 16))
(setf yv '#(18 19))


(defun mm (x y)
  (let* ((m (array-dimension x 0))
	 (k (array-dimension x 1))
	 (n (array-dimension y 1))
	 (z (make-array (list m n) :initial-element 0)))
    (dotimes (i m)
      (dotimes (j n)
        (dotimes (l k)
          (incf (aref z i j) (* (aref x i l) (aref y l j))))))
    z))

(defun mv (x y)
  (let* ((m (array-dimension x 0))
	 (n (array-dimension x 1))
	 (z (make-array m :initial-element 0)))
    (dotimes (i m)
      (dotimes (j n)
          (incf (aref z i) (* (aref x i j) (aref y j)))))
    z))

(defun vm (x y)
  (let* ((m (array-dimension y 0))
	 (n (array-dimension y 1))
	 (z (make-array n :initial-element 0)))
    (dotimes (i m)
      (dotimes (j n)
        (incf (aref z j) (* (aref x i) (aref y i j)))))
    z))

(defun ip (x y &optional (conj t))
  (sum (* x (if conj (conjugate y) y))))

(defun test-report (name val)
  ;(format t "~&~a: ~9t~f~%" name val)
  (check #'< val eps)
  )


(test-report "MATMULT"
	     (let ((cx (complex x (* 2 x)))
		   (cy (complex y (* 3 y)))
		   (cxv (complex xv (* 4 xv)))
		   (cyv (complex yv (* 5 yv))))
	       (max (abs (- (matmult x y) (mm x y)))
		    (abs (- (matmult cx cy) (mm cx cy)))
		    (abs (- (matmult x yv) (mv x yv)))
		    (abs (- (matmult cx cyv) (mv cx cyv)))
		    (abs (- (matmult xv y) (vm xv y)))
		    (abs (- (matmult cxv cy) (vm cxv cy)))
		    (abs (- (matmult xv yv) (ip xv yv)))
		    (abs (- (matmult cxv cyv) (ip cxv cyv nil))))))

(test-report "CROSS-PRODUCT"
	     (max (abs (- (cross-product x) (mm (transpose x) x)))
		  (let ((cx (complex x (* 2 x))))
		    (abs (- (cross-product cx)
			    (mm (transpose cx) (conjugate cx)))))
		  (let ((cx (complex x (* 2 x))))
		    (abs (- (cross-product cx nil)
			    (mm (transpose cx) cx))))))

(test-report "INNER-PRODUCT"
	     (let ((cxv (complex xv (* 4 xv)))
		   (cyv (complex yv (* 5 yv))))
	       (max (abs (- (inner-product xv yv) (ip xv yv)))
		    (abs (- (inner-product cxv cyv) (ip cxv cyv)))
		    (abs (- (inner-product cxv cyv nil) (ip cxv cyv nil))))))

(let* ((x (list (iseq 1 4) (^ (iseq 1 4) 2)))
       (y (^ (iseq 1 4) 3))
       (w (iseq 4 1))
       (m1 (regression-model x y :print nil))
       (mw (regression-model x y :print nil :weights w))
       (xm (apply #'bind-columns (repeat 1 4) x))
       (xmt (transpose xm))
       (wm (diagonal w))
       (b1 (matmult (inverse (matmult xmt xm)) xmt y))
       (bw (matmult (inverse (matmult xmt wm xm)) xmt wm y)))
  (test-report "REGRESS"
	       (max (abs (- (send m1 :coef-estimates) b1))
		    (abs (- (send mw :coef-estimates) bw)))))