;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Mod 2 linear algebra with sparse matrices
;; John Palmieri, 2003/02/06
;;
;; vectors are lists of basis elements
;; matrices are alists of vectors: car is row number, cdr is vector.
;; the function sort-vector will need to be changed to match whatever
;; basis elements look like.
;;
;; print-vector and print-matrix are probably broken right now.
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; contents:
;;  sort-vector
;;  insert-vector-in-matrix
;;  extract-vector-from-matrix
;;  row-switch
;;  add-vectors
;;  row-add
;;  row-reduce
;;  rank
;;  in-span
;;  kernel
;;  print-matrix
;;  print-vector
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; (load "front.lisp")
;; (in-package bss)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; sort-vector:
;; argument: vec
;; return: vector, sorted by some means.
(defun sort-vector (vec)
  (sort-poly vec))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; insert-vector-in-matrix:
;; arguments: vec, mat, n
;; replace row[n] of mat with vec, and return resulting matrix.
(defun insert-vector-in-matrix (vec mat n)
  "replace row[N] of MAT with VEC, and return resulting matrix"
  (let ((old-vec (assoc n mat))
	(temp mat))  ;; was (temp (copy-alist mat)))
    (if old-vec
	(setf temp (delete old-vec temp :test #'equal)))
    (acons n vec temp)))

(defun insert-vector-in-matrix (vec mat n)
  "replace row[N] of MAT with VEC, and return resulting matrix"
  (let ((old-vec (assoc n mat)))
    (if old-vec
        (progn (setf (cdr old-vec) vec)
               mat)
        (acons n vec mat))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; extract-vector-from-matrix:
;; arguments: mat, n
;; return row n of mat (a vector)
(defun extract-vector-from-matrix (mat n)
  "row[N] of MAT"
  (cdr (assoc n mat)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; row-switch:
;; arguments: mat, i, j
;; given matrix mat, switches rows i and j
;; return resulting matrix
;; 
;; plan: find row i and save it.  find row j and insert it in row i.
;; then insert row i in row j.
(defun row-switch (mat i j)
  "given sparse matrix MAT, switch rows I and J"
  (if (null mat)
      nil
    (let ((junk -1)
	  vec
	  (temp mat))    ;; was (temp (copy-alist mat))
      (unless (= i j)
	(when (setf vec (assoc i temp))
	  (rplaca vec junk))
	(when (setf vec (assoc j temp))
	  (rplaca vec i))
	(when (setf vec (assoc junk temp))
	  (rplaca vec j)))
      temp)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; add-vectors:
;; arguments: vec1, vec2
;; return: sum of vec1 and vec2, mod 2
(defun add-vectors (vec1 vec2)
  "mod 2 sum of VEC1 and VEC2"
  (sort-vector (set-exclusive-or vec1 vec2 :test #'equal)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; row-add:
;; arguments: mat, int i, int j
;; given matrix mat, adds row i to j, mod 2
;; return resulting matrix
(defun row-add (mat i j)
  "given matrix MAT, adds row I to J (mod 2) and return resulting matrix"
  (if (null mat)
      nil
    (let ((vec (add-vectors (extract-vector-from-matrix mat i)
			    (extract-vector-from-matrix mat j))))
      (if (null vec)
	  (remove (assoc j mat) mat :test #'equal)
	(insert-vector-in-matrix vec mat j)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; row-reduce:
;; arguments: mat, &optional mat2
;; return:
;; row-echelon form of mat
;; also do operations on mat2, if non-null
;;
;; strategy: go through column by column looking for a 1.
;; when you find it, put it in row 1.  then add it to other rows with 1s
;; in the same column to clear that column.  now look for the next 1
;; and repeat.
(defun row-reduce (mat &optional mat2)
  "row reduce MAT.  If optional arg MAT2 is present, do the same row
operations on that matrix (as if it were augmenting MAT)."
  (let ((row 0)	   ;; up to this row so far
	(column 0) ;; up to this column so far
	(answer mat)	;; progress so far
	next-row ;; row of next nonzero entry
	coord
	(done nil))
    (loop until done do
	  (setf next-row 0
		coord nil
		column
		(car (sort-vector 
		      (loop for matrix-row in answer
			    when (< row (car matrix-row))
			    when (consp (cdr matrix-row))
			    ;; find smallest leading term after current row
			    collect (car (cdr matrix-row))))))
	  (setf done (null column))
	  (unless done
	    (loop for matrix-row in answer
		  until coord
		  ;; only look in rows after the ones we've dealt with:
		  when (< row (car matrix-row))
		  when (consp (cdr matrix-row)) do
		  (setf coord (equal (cadr matrix-row) column)
			next-row (car matrix-row)))
	    (when (and (< 0 next-row) coord)
	      (setf row (1+ row)
		    answer (row-switch answer row next-row)
		    mat2 (row-switch mat2 row next-row))
	      (loop for matrix-row in answer
		    when (/= (car matrix-row) row)
		    when (member column (cdr matrix-row) :test #'equal) do
		    (setf answer (row-add answer row (car matrix-row))
			  mat2 (row-add mat2 row (car matrix-row))))
 	      (setf answer (sort answer #'< :key #'car)
 		    mat2 (sort mat2 #'< :key #'car)))))
    (if mat2
	(values answer mat2)
      answer)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; rank:
;; arguments: mat
;; return: 
;; rank of mat
;; 
;; strategy: assume that mat is already row-reduced, and find
;; last non-zero row.
(defun rank (mat)
  "rank of the row-reduced matrix MAT"
  (length mat))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; in-span:
;; arguments: vec, mat
;; return:
;; t if vec is in the span of the rows of mat
;; nil if vec is not in the span of the rows of mat
;;
;; strategy: row-reduce mat and compute its rank.  add vec after
;; the last row (which is indexed by the rank), row-reduce that,
;; and compute its rank.  see if the ranks are equal
(defun in-span (vec mat)
  "t if VEC is in the span of the rows of MAT, and nil otherwise."
  (if (null vec)
      t
    (if (null mat)
	nil
      (let (rank
	    (reduced (row-reduce (copy-alist mat)))
	    augmented)
	(setf rank (rank reduced)
	      augmented
	      (row-reduce 
	       (insert-vector-in-matrix vec reduced (1+ rank))))
	(<= (rank augmented) rank)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; kernel:
;; arguments: mat, basis
;; compute kernel of mat acting on space of spanned by basis,
;; return list of vectors which spans the kernel
;;
;; strategy:
;; 
;; augment mat by sticking an identity matrix on the right side.
;;   size of identity matrix = size of BASIS
;; row-reduce the result
;; find rows where left side is zero and right side is not; then 
;; the vectors on the right side give a basis for the kernel
;;
;; in our case, the entries in mat should all be larger than
;; the entries in basis, so sort and truncate accordingly.
(defun kernel (mat basis)
  "given matrix MAT acting on a space with basis BASIS, return
basis for the kernel -- a list of vectors"
  (let ((aug nil))	    ;; augmenting matrix

    (loop for x in basis
	count x into row do
	(setf aug (insert-vector-in-matrix (list x) aug row)))

    ;; augment matrix and row-reduce it
    (multiple-value-bind (reduced augmented) (row-reduce mat aug)
      ;; find rows where "left side" is zero and "right side" is not.
      (loop for vec in augmented
	  count vec into row
	  if (null (extract-vector-from-matrix reduced row))
	  collect (cdr vec)))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

(defvar matrix-printable-rows 6
  "number of rows of matrices to print")
(defvar vector-printable-cols 60
  "number of columns of vectors and matrices to print")

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; matrix-to-string:
;; arguments: mat
;; convert mat to string
(defun matrix-to-string (mat) 
  "convert MAT to string"
  (let ((last-row 0)
	(output-string ""))
    (dolist (row (sort (copy-alist mat) #'< :key #'car))
      (do ((i (1+ last-row) (1+ i)))
	  ((> i (min (1- (car row)) matrix-printable-rows)))
	(setf output-string (concatenate 'string
					 output-string
					 (vector-to-string nil)
					 "
")))
      (setf output-string (concatenate 'string
				       output-string
				       (vector-to-string (cdr row))
				       "
"))
      (setf last-row (car row))
      (when (< matrix-printable-rows last-row)
	(return)))
    (do ((i last-row (1+ i)))
	((> i matrix-printable-rows))
      (setf output-string (concatenate 'string
				       output-string
				       (vector-to-string nil)
				       "
")))
    output-string))
		
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; vector-to-string:
;; arguments: vec
;; convert vec to string
(defun vector-to-string (vec) 
  "convert VEC to string"
  (let ((last-entry 0)
	(output-string ""))
    (loop for entry in vec
	  when (<= entry vector-printable-cols) do
	  (do ((j (1+ last-entry) (1+ j)))
	      ((> j (min (1- entry) vector-printable-cols)))
	    (setf output-string (concatenate 'string
					     output-string
					     "0")))
	  (setf output-string (concatenate 'string
					   output-string
					   "1")
		last-entry entry))
    (do ((j last-entry (1+ j)))
	((>= j vector-printable-cols))
      (setf output-string (concatenate 'string
				       output-string
				       "0")))
    output-string))
	
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;; (load "front.fasl")
;; (in-package lambda-alg)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;