math.bnl.lisp

;; ROW MAJOR MODE MATRIX STUFF

;; matrix definition
;;        w h elems
(SETQ M '(2 2 (1 0
	       0 1)))

(SETQ MATRIX-NUMEL
      (LAMBDA (MAT)
	      (* (CAR MAT) (CAR (CDR MAT)))))

(SETQ MATRIX-PRINT
      (LAMBDA (MAT)
	      (SETQ count (MATRIX-NUMEL MAT))
	      (SETQ c 0)
	      (SETQ cc 0)
	      (SETQ tail (CAR (CDR (CDR MAT))))
	      (WHILE (NOT (= c count))
		     (SETQ c (+ c 1))
		     (IF (= cc (CAR MAT))
		     	 (PROGN
				 (SETQ cc 0)
				 (PRINC "\n"))
			 ())
		     (SETQ cc (+ cc 1))
		     (PRINC (CAR tail))
		     (PRINC "\t")
		     (SETQ tail (CDR tail)))
	      (PRINC "\n")
	      MAT))

(SETQ MATRIX-GET-IJ
      (LAMBDA (MAT I J)
	      (SETQ COUNT 0)
	      (SETQ IDX (+ (* I (CAR MAT)) J))
	      (SETQ TAIL (CAR (CDR (CDR MAT))))
	      (WHILE (NOT (= COUNT IDX))
		     (SETQ COUNT (+ COUNT 1))
		     (SETQ TAIL (CDR TAIL)))
	      (CAR TAIL)))

(SETQ APPEND (LAMBDA (A B)
		     (IF (EQ A NIL)
			 B
			 (CONS (CAR A) (APPEND (CDR A) B)))))

(SETQ MATRIX-SET-IJ
      (LAMBDA (MAT I J V)
	      (SETQ OUT ())
	      (SETQ IDX (+ (* I (CAR MAT)) J))
	      (SETQ TAIL (CAR (CDR (CDR MAT))))
	      (SETQ I 0)
	      (WHILE (NOT (EQ TAIL NIL))
		     (IF (= I IDX)
			 (SETQ OUT (APPEND OUT (CONS V NIL)))
			 (SETQ OUT (APPEND OUT (CONS (CAR TAIL) NIL))))
		     (SETQ TAIL (CDR TAIL))
		     (SETQ I (+ I 1)))
	      (CONS (CAR MAT) (CONS (CAR (CDR MAT)) (CONS OUT NIL) NIL))))



(SETQ MATRIX-TRANSPOSE
      (LAMBDA (MAT)
	      (SETQ I 0)
	      (SETQ M (CAR (CDR (CDR MAT))))
	      (SETQ OUT ())
	      (WHILE (NOT (= I (CAR MAT)))
		     (SETQ J 0)
		     (WHILE (NOT (= J (CAR (CDR MAT))))
			    (SETQ OUT (APPEND OUT (CONS (MATRIX-GET-IJ MAT J I) NIL)))
			    (SETQ J (+ J 1)))
		     (SETQ I (+ I 1)))
	      (CONS (CAR (CDR MAT)) (CONS (CAR MAT) (CONS OUT NIL) NIL))))

(SETQ MATRIX-GET-ROW
      (LAMBDA (M R) ;; get row R of matrix M
	      (SETQ OUT ())
	      (SETQ COUNT 0)
	      (WHILE (NOT (= COUNT (CAR M)))
		     (SETQ OUT (APPEND OUT (CONS (MATRIX-GET-IJ M R COUNT) NIL)))
		     (SETQ COUNT (+ COUNT 1)))
	      OUT))

(SETQ MATRIX-GET-COL
      (LAMBDA (M J) ;; get col J of matrix M
	      (SETQ OUT ())
	      (SETQ COUNT 0)
	      (WHILE (NOT (= COUNT (CAR (CDR M))))
		     (SETQ OUT (APPEND OUT (CONS (MATRIX-GET-IJ M COUNT J) NIL)))
		     (SETQ COUNT (+ COUNT 1)))
	      OUT))

(SETQ SUM
      (LAMBDA (V)
	      (SETQ S 0)
	      (SETQ TAIL V)
	      (WHILE (NOT (EQ TAIL NIL))
		     (SETQ S (+ S (CAR TAIL)))
		     (SETQ TAIL (CDR TAIL)))
      S))

(SETQ DOT-PRODUCT
      (LAMBDA (A B)
	      (SETQ TAIL-A A)
	      (SETQ TAIL-B B)
	      (SETQ OUT ())
	      (WHILE (NOT (EQ TAIL-A NIL))
		     (SETQ OUT (APPEND OUT (CONS (* (CAR TAIL-A) (CAR TAIL-B)) NIL)))
		     (SETQ TAIL-A (CDR TAIL-A))
		     (SETQ TAIL-B (CDR TAIL-B)))
	      (SUM OUT)))

;; assumes dimensions A[x,y] = B[y,z]
(SETQ MATRIX-MULTIPLY
      (LAMBDA (A B)
	      (SETQ OUT-W (CAR (CDR A)))
	      (SETQ OUT-H (CAR B))
	      (SETQ OUT ())
	      (SETQ J 0)
	      (WHILE (NOT (= J OUT-H))
		     (SETQ I 0)
		     (WHILE (NOT (= I OUT-W))
			    (SETQ OUT
				  (APPEND OUT (CONS (DOT-PRODUCT
						     (MATRIX-GET-ROW A I)
						     (MATRIX-GET-COL B J))
						    NIL)))
			    (SETQ I (+ I 1)))
		     (SETQ J (+ J 1)))
	      (CONS OUT-W (CONS OUT-H (CONS OUT NIL)))))

(SETQ MATRIX-SWAP-ROWS
      (LAMBDA (M R1 R2)
	      (SETQ NUM-ROWS (CAR (CDR M)))
	      (SETQ R 0)
	      (SETQ OUT ())
	      (WHILE (NOT (= R NUM-ROWS))
		     (IF (= R R1)
			 (SETQ OUT (APPEND OUT (MATRIX-GET-ROW M R2)))
			 (IF (= R R2)
			     (SETQ OUT (APPEND OUT (MATRIX-GET-ROW M R1)))
			     (SETQ OUT (APPEND OUT (MATRIX-GET-ROW M R)))))
		     (SETQ R (+ R 1)))
	      (CONS (CAR M) (CONS (CAR (CDR M)) (CONS OUT NIL)))))

(SETQ ARG-MAX
      (LAMBDA (L FROM)
	      (SETQ TAIL L)
	      (SETQ I-MAX 0)
	      (SETQ I-VAL 0)
	      (SETQ I 0)
	      (WHILE (NOT (= I FROM))
		     (SETQ I (+ I 1))
		     (SETQ TAIL (CDR TAIL)))
	      (WHILE (NOT (EQ TAIL NIL))
		     (IF (> (CAR TAIL) I-VAL)
			 (PROGN (SETQ I-MAX I)
				(SETQ I-VAL (CAR TAIL)))
			 NIL)
		     (SETQ I (+ I 1))
		     (SETQ TAIL (CDR TAIL)))
	      I-MAX))

(SETQ MATRIX-RREF ;; reduced row ecelon form
      (LAMBDA (M)
	      (SETQ NROW (CAR (CDR M)))
	      (SETQ NCOL (CAR M))
	      (SETQ H 0) ;; PIVOT ROW
	      (SETQ K 0) ;; PIVOT COL
	      (WHILE (NOT (= H NCOL))
		     (WHILE (NOT (= K NROW))
			    ;; BEGIN BY FINDING THE K PIVOT...
			    (SETQ I-MAX (ARG-MAX (MATRIX-GET-COL M K) H))

			    (IF (NOT (= (MATRIX-GET-IJ M I-MAX K) 0))
				(PROGN
				 (SETQ M (MATRIX-SWAP-ROWS M I-MAX H))
				 ;; FOR ALL ROWS BELOW PIVOT
				 (SETQ I (+ H 1))
				 (WHILE (NOT (= I NROW))
					(SETQ F (/ (MATRIX-GET-IJ M I K)
						   (MATRIX-GET-IJ M H K)))
					(SETQ M (MATRIX-SET-IJ M I K 0))
					(SETQ J (+ K 1))
					(WHILE (NOT (= J NCOL))
					       (SETQ M (MATRIX-SET-IJ
					       		M I J
					       		(- (MATRIX-GET-IJ M I J)
					       		   (* (MATRIX-GET-IJ M H J) F))))
					       (SETQ J (+ J 1)))
					(SETQ I (+ I 1)))
				 (SETQ H (+ H 1)))
				NIL)
			    (SETQ K (+ K 1))))
	      M))

(SETQ MATRIX-HORZ-CONCAT
      (LAMBDA (M1 M2)

	      ))

(SETQ M '(8 4 (1.0 2.0 3.0 4.0     1.0 0.0 0.0 0.0
		   4.0 3.0 2.0 1.0 0.0 1.0 0.0 0.0
		   1.0 2.0 3.0 4.0 0.0 0.0 1.0 0.0
		   4.0 3.0 2.0 1.0 0.0 0.0 0.0 1.0)))


(PRINC ">>> Here is matrix M:\n\n")
(MATRIX-PRINT M)

;; (PRINT (MATRIX-GET-IJ M 1 2))
;; (MATRIX-PRINT (MATRIX-SET-IJ M 1 2 0))

;; (PRINC "\n\n>>> Here is M * M(transposed):\n\n")
;; (SETQ M (MATRIX-MULTIPLY M (MATRIX-TRANSPOSE M)))
;; (MATRIX-PRINT M)

;;(MATRIX-PRINT (MATRIX-SWAP-ROWS M 0 1))

;;(SETQ M (MATRIX-SET-IJ M 2 0 1))


(PRINC "\n>>> Here is the row-echelon form of M:\n\n")
;; (MATRIX-PRINT M)
;; (PRINT (MATRIX-GET-IJ M 1 3))
;; (MATRIX-PRINT (MATRIX-SET-IJ M 1 2 10))
(MATRIX-PRINT (MATRIX-RREF M))