tests/matrix2.lsp

   1 (setf eps 1e-10)
   2
   3 (setf x '#2a((1 2)(3 4)(5 6)(7 8)))
   4 (setf y '#2a((9 10 11) (12 13 14)))
   5 (setf xv '#(15 16))
   6 (setf yv '#(18 19))
   7
   8
   9 (defun mm (x y)
  10   (let* ((m (array-dimension x 0))
  11          (k (array-dimension x 1))
  12          (n (array-dimension y 1))
  13          (z (make-array (list m n) :initial-element 0)))
  14     (dotimes (i m)
  15       (dotimes (j n)
  16         (dotimes (l k)
  17           (incf (aref z i j) (* (aref x i l) (aref y l j))))))
  18     z))
  19
  20 (defun mv (x y)
  21   (let* ((m (array-dimension x 0))
  22          (n (array-dimension x 1))
  23          (z (make-array m :initial-element 0)))
  24     (dotimes (i m)
  25       (dotimes (j n)
  26           (incf (aref z i) (* (aref x i j) (aref y j)))))
  27     z))
  28
  29 (defun vm (x y)
  30   (let* ((m (array-dimension y 0))
  31          (n (array-dimension y 1))
  32          (z (make-array n :initial-element 0)))
  33     (dotimes (i m)
  34       (dotimes (j n)
  35         (incf (aref z j) (* (aref x i) (aref y i j)))))
  36     z))
  37
  38 (defun ip (x y &optional (conj t))
  39   (sum (* x (if conj (conjugate y) y))))
  40
  41 (defun test-report (name val)
  42   ;(format t "~&~a: ~9t~f~%" name val)
  43   (check #'< val eps)
  44   )
  45
  46
  47 (test-report "MATMULT"
  48              (let ((cx (complex x (* 2 x)))
  49                    (cy (complex y (* 3 y)))
  50                    (cxv (complex xv (* 4 xv)))
  51                    (cyv (complex yv (* 5 yv))))
  52                (max (abs (- (matmult x y) (mm x y)))
  53                     (abs (- (matmult cx cy) (mm cx cy)))
  54                     (abs (- (matmult x yv) (mv x yv)))
  55                     (abs (- (matmult cx cyv) (mv cx cyv)))
  56                     (abs (- (matmult xv y) (vm xv y)))
  57                     (abs (- (matmult cxv cy) (vm cxv cy)))
  58                     (abs (- (matmult xv yv) (ip xv yv)))
  59                     (abs (- (matmult cxv cyv) (ip cxv cyv nil))))))
  60
  61 (test-report "CROSS-PRODUCT"
  62              (max (abs (- (cross-product x) (mm (transpose x) x)))
  63                   (let ((cx (complex x (* 2 x))))
  64                     (abs (- (cross-product cx)
  65                             (mm (transpose cx) (conjugate cx)))))
  66                   (let ((cx (complex x (* 2 x))))
  67                     (abs (- (cross-product cx nil)
  68                             (mm (transpose cx) cx))))))
  69
  70 (test-report "INNER-PRODUCT"
  71              (let ((cxv (complex xv (* 4 xv)))
  72                    (cyv (complex yv (* 5 yv))))
  73                (max (abs (- (inner-product xv yv) (ip xv yv)))
  74                     (abs (- (inner-product cxv cyv) (ip cxv cyv)))
  75                     (abs (- (inner-product cxv cyv nil) (ip cxv cyv nil))))))
  76
  77 (let* ((x (list (iseq 1 4) (^ (iseq 1 4) 2)))
  78        (y (^ (iseq 1 4) 3))
  79        (w (iseq 4 1))
  80        (m1 (regression-model x y :print nil))
  81        (mw (regression-model x y :print nil :weights w))
  82        (xm (apply #'bind-columns (repeat 1 4) x))
  83        (xmt (transpose xm))
  84        (wm (diagonal w))
  85        (b1 (matmult (inverse (matmult xmt xm)) xmt y))
  86        (bw (matmult (inverse (matmult xmt wm xm)) xmt wm y)))
  87   (test-report "REGRESS"
  88                (max (abs (- (send m1 :coef-estimates) b1))
  89                     (abs (- (send mw :coef-estimates) bw)))))