tests/matrix.lsp

   1 ; test of matrix related functions (not very thorough)
   2
   3 (setf eps 1.e-7
   4       xmat (transpose (matrix '(4 6) (iseq 1 24)))
   5       xmat2 (permute-array xmat '(1 0)))
   6
   7 (check #'= (array-dimension (matrix '(4 1) (iseq 1 4)) 0) 4)
   8 (check #'= (aref xmat2 0 1) (aref xmat 1 0))
   9 (check #'= (mapcar #'mean (column-list xmat)) (mapcar #'mean (row-list xmat2)))
  10
  11 (check #'=
  12        (bind-columns (iseq 1 3) (iseq 4 6)) (bind-rows '(1 4) '(2 5) '(3 6)))
  13 (check #'=
  14        (mapcar #'(lambda (x y) (select xmat x y)) '(0 1 2 3) '(0 1 2 3))
  15        (diagonal xmat))
  16 (check #'= xmat (transpose xmat2))
  17
  18
  19 #|
  20 ymat matrix(seq(0.1,18),nrow=6) # ray
  21 cancor(xmat,ymat)$cor!=1                        # ray begins
  22 sum(hat(longley.x,F))==dim(longley.x)[2]
  23 all(l1fit(matrix(rnorm(24),nrow=6),(iseq 1 6))$coef!=0)
  24 all(leaps(matrix(rnorm(24),nrow=6),(iseq 1 6))$Cp!=0)
  25 sum(ls.summary(ls.out_lsfit(matrix(runif(18),nrow=6),(iseq 1 6),int=F))$hat)-3 < eps
  26 all(qr.coef(ls.out$qr,(iseq 1 6))!=0)
  27 all(qr.fitted(ls.out$qr,(iseq 1 6))+qr.resid(ls.out$qr,(iseq 1 6))+qr.qty(ls.out$qr,(iseq 1 6))+qr.qy(ls.out$qr,(iseq 1 6))!=0)
  28 all(apply(scale(ymat),2,"var")-1<eps)
  29 mean(scale(ymat))<eps                           # ray ends
  30 {
  31         ymat <- matrix(c(1.01,2.02,-3.2,4.7,60,-14.3),ncol=2)
  32         ymat2 <- matrix(c(3.3,5.01,2.3,-7.4,-38,19.9),ncol=2)
  33         T
  34 }
  35 all(t(ymat)%*%ymat2-ymat%c%ymat2<eps)
  36 all(crossprod(ymat,ymat2)-ymat%c%ymat2==0)
  37 all(ymat+ymat2==ymat2+ymat)
  38 |#
  39
  40 (flet ((tchol (x)
  41          (let ((c (chol-decomp x)))
  42            (abs (- (matmult (first c) (transpose (first c))) x))))
  43        (tsvd (x)
  44          (let ((s (sv-decomp x)))
  45            (abs (- (matmult (first s)
  46                             (diagonal (second s))
  47                             (transpose (third s)))
  48                    x))))
  49        (teigen (x)
  50          (let ((e (eigen x)))
  51            (abs (- (matmult (apply #'bind-columns (second e))
  52                             (diagonal (first e))
  53                             (transpose (apply #'bind-columns (second e))))
  54                    x)))))
  55 #|
  56         tvar_function(x){(length(x)-1)*var(x)+length(x)*mean(x)^2-sum(x^2)}
  57         tvar2_function(x,y=x){n_nrow(x);(n-1)*var(x,y)-(x%c%y-n*outer(apply(x,2,mean),apply(y,2,mean)))}
  58         tcor_function(x,y=x){cor(x,y)-var(x,y)/sqrt(outer(diag(var(x)),diag(var(y))))}
  59 |#
  60   (let* ((x (matrix '(50 10) (uniform-rand 500)))
  61          (y (cross-product x))
  62          (z1 (matrix '(10 3) (normal-rand 30)))
  63          (z2 (matrix '(10 3) (uniform-rand 30))))
  64     (check #'< (tchol y) eps)
  65     (check #'< (tsvd x) eps)
  66     (check #'< (teigen y) eps)
  67     (check #'< (- (matmult (inverse y) y) (identity-matrix 10)) eps)))
  68 #|
  69 tvar(runif(100))<100*eps
  70 all(tvar2(z1)<100*eps)
  71 all(tvar2(z2)<100*eps)
  72 all(tvar2(z1,z2)<100*eps)
  73 all(tcor(x)<100*eps)
  74 all(tcor(z1)<100*eps)
  75 all(tcor(z2)<100*eps)
  76 all(tcor(z1,z2)<100*eps)
  77
  78 {
  79         m_matrix((iseq 1 12),3,4)
  80         s_sweep(m,2,apply(m,2,mean))
  81         s2_sweep(m,1,apply(m,1,prod),'*')
  82         s3_matrix(8*c(35,220,729,140,550,1458,245,880,2187,350,1210,2916),3,4)
  83         T
  84 }
  85 all(dim(s)==dim(m))
  86 all(s==rep(-(iseq 1 1),4))
  87 all(dim(s2)==dim(s3))
  88 all(s2==s3)
  89 {
  90         s1_sample((iseq 1 10))
  91         s2_sample((iseq 1 10),20,replace=T)
  92         s3_sample((iseq 1 10000),10)
  93         T
  94 }
  95 all(sort(s1)==(iseq 1 10))
  96 all(min(s2)>=1&&max(s2)<=10&&length(s2)==20)
  97 all(min(s3)>=1&&max(s3)<=10000&&length(s3)==10)
  98 |#
  99