;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; Lambda algebra calculator
;; John Palmieri, 2003/02/06
;;
;; computing homology of Lambda, homology of
;; Lambda/im(theta), where
;; theta: lambda_n --> lambda_{2n+1}.
;; compute the Bockstein spectral sequence associated to theta.
;; 
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;
;; contents
;;  odd-ending
;;  strip-odds-from-poly
;;  strip-odds-from-set
;;  matrix-to-set
;;  set-to-matrix
;;  boundaries
;;  boundary-matrix
;;  cycles
;;  cycle-matrix
;;  ext
;;  theta-inverse
;;  boundaries-all
;;  homological
;;  internal
;;  boundaryp
;;  boundary-preimage
;;  bss
;;  print-betti
;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; odd-ending:
;; arguments: poly
;; return:
;; subset of poly, consisting of those terms which end in odd numbers.
(defun odd-ending (poly)
  "given a polynomial POLY,  return terms of POLY which end in odd numbers.
terms which end in odd numbers."
  (loop for mono in poly
	when (oddp (car (last mono)))
	collect mono))
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; strip-odds-from-poly:
;; arguments: poly
;; return:
;; all terms in poly which do not consist entirely of odd lambdas
(defun strip-odds-from-poly (poly)
  "all terms in POLY which do not consist entirely of odd lambdas"
  (loop for mono in poly
	when (loop for n in mono
		   when (evenp n)
		   return (evenp n))
	collect mono))
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; strip-odds-from-set:
;; arguments: set
;; return:
;; set after removing all monomials which consist
;; only of odd lambdas
(defun strip-odds-from-set (set)
  "remove all monomials from all polynomials in SET which only of odd lambdas"
  (loop for poly in set
	unless (null (strip-odds-from-poly poly))
	collect (strip-odds-from-poly poly)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; matrix-to-set:
;; arguments: mat
;; convert each row of matrix to a poly
(defun matrix-to-set (mat)
  "convert each row of MAT to a polynomial by removing the row index
for each row."
  (loop for matrix-row in mat
	collect (cdr matrix-row)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; set-to-matrix:
;; arguments: set
;; convert set to a matrix by writing each poly in coords with
;; respect to basis.
(defun set-to-matrix (set)
  "convert set to matrix by adding a row index for each row."
  (loop for poly in set
	count poly into row
	collect (cons row (sort-poly poly))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; boundaries:
;; arguments: int homological, int internal, int no-odd
;; return:
;; basis for the odd-ending boundaries in bidegree (hom,int).
;; if no-odd is non-zero, remove terms which are all-odd before and
;; after computing boundaries.
;;
;; strategy: call boundary-matrix, and then convert to set of polys.
(defun boundaries (homological internal &optional no-odd)
  "basis for the odd-ending boundaries in bidegree (HOM,INT).
If optional arg NO-ODD is non-nil, remove terms which are all-odd
before and after computing boundaries."
  (matrix-to-set (boundary-matrix homological internal no-odd)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; boundary-matrix:
;; arguments: homological, internal, no-odd
;; return:
;; matrix for the basis of the odd-ending boundaries in bidegree (hom,int), 
;; in coordinates in terms of lambda-basis(hom,int).
;; if no-odd is non-zero, remove terms which are all-odd before and
;; after computing boundaries.
;; 
;; plan: compute boundaries of all odd-ending terms.
;;   expand results in terms of basis of Lambda, getting matrix bdries.
;;   row-reduce this matrix to get basis for the space of boundaries.
(defun boundary-matrix (homological internal &optional no-odd raw)
  "matrix for the basis of the odd-ending boundaries in bidegree (HOM,INT), 
in coordinates in terms of lambda-basis(HOM,INT).
If optional arg NO-ODD is non-nil, remove terms which are all-odd before
and after computing boundaries."
  (let ((other-basis (lambda-basis (1- homological) internal))
	mat
	poly
	(index 1))

    ;; initialize other-basis:
    (if no-odd
	(setf other-basis (strip-odds-from-poly other-basis)))
    (setf other-basis (odd-ending other-basis))

    ;; for each elt of basis, compute its boundary to get matrix mat
    (dolist (mono other-basis)
      (setf poly (coboundary-mono mono))
      (if no-odd
	  (setf poly (strip-odds-from-poly poly)))
      (unless (and (not raw) (null poly))
	(setf mat (insert-vector-in-matrix poly mat index)
	      index (1+ index))))

    (if raw
	mat
      (row-reduce mat))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; cycles:
;; arguments: int homological, int internal, int no-odd
;; return:
;; basis for the cycles in bidegree (hom,int).
;; if no-odd is non-zero, remove terms which are all-odd before and
;; after computing baundaries.
;; 
;; strategy: call cycle-matrix, and then convert to set of polys.
(defun cycles (homological internal &optional no-odd)
  "basis for the cycles in bidegree (HOM,INT).
If optional arg NO-ODD is non-nil, remove terms which are all odd
before and after computing boundaries."
  (matrix-to-set (cycle-matrix homological internal no-odd)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; cycle-matrix:
;; arguments: int homological, int internal, int no-odd, set-of-polys * basis
;; return:
;; matrix for the basis for the cycles in bidegree (hom,int),
;; in coordinates in terms of lambda-basis(hom,int)
;; if no-odd is non-zero, remove terms which are all-odd before and
;; after computing baundaries.
;;
;; strategy: take all odd-ending basis elements in bidegree (hom, int)
;;   compute their boundaries
;;   convert to coordinates in terms of the basis for degree (hom+1, int)
;;   thus getting a matrix
;;   augment by adding an identity matrix on the right side
;;   row reduce this matrix
;;   basis for kernel: the right half of each row where the left half is zero
;;      is a basis vector
(defun cycle-matrix (homological internal &optional no-odd)
  "matrix for the basis of the cycles in bidegree (HOM,INT), 
in coordinates in terms of lambda-basis(HOM,INT).
If optional arg NO-ODD is non-nil, remove terms which are all-odd before
and after computing boundaries."
  (let ((basis (lambda-basis homological internal))
	mat)
    
    ;; initialize basis
    (if no-odd
	(setf basis (strip-odds-from-poly basis)))
    (setf basis (odd-ending basis))
 
    (setf mat (boundary-matrix (1+ homological) internal no-odd t))

    (unless (or (null mat) (null basis))
      (set-to-matrix (kernel mat basis)))))

;; unused at the moment:
(defun strip-nils-from-mat (mat)
  "remove any NIL rows from MAT and row-reduce the result"
  (row-reduce
   (loop for row in (copy-alist mat)
	 unless (null (cdr row))
	 collect row)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; ext:
;; arguments: int homological, int internal, quotient
;; return:
;; if quotient is nil, return cohomology of Lambda in bidegree (hom,int)
;; if quotient is t, return cohomology of Lambda/im(theta) in bidegree (hom,int)
;;
;; strategy: compute boundaries and cycles in bidegree(hom, int)
;; for each cycle, check: is it in the span of boundaries?
;;    if yes, discard it.
;;    if no, add it to a basis for ext, and also add it to the boundaries.
(defun ext (homological internal &optional quotient)
  "Compute cohomology of Lambda in bidegree (HOM, INT).
If optional argument QUOTIENT is non-nil, compute cohomology of
Lambda/im(theta) in bidegree (HOM,INT)."
  (let (mat-b      ;; matrix of boundaries
	mat-z      ;; matrix of cycles
	mat-h      ;; matrix of answer
	row-b ;; index of first row of zeroes in mat-b = 1 + rank(mat-b)
	(row-h 0) ;; how many ext elements we've found so far
	(basis (lambda-basis homological internal))
	span)
    (unless (> homological internal)
      (if (= homological internal)
	  (list (list (make-list homological :initial-element 0)))
	(progn
	  (setf mat-b (boundary-matrix homological internal quotient)
		mat-z (boundary-matrix (1+ homological) internal quotient))
	  (if quotient (setf basis (strip-odds-from-poly basis)))
	  (setf basis (odd-ending basis))
	  (when (< (length mat-b) (- (length basis)
				     (length mat-z)))
	    (setf mat-z (cycle-matrix homological internal quotient)
		  row-b (1+ (rank mat-b)))

	    ;; now go through cycles, seeing if each one is in span of
	    ;;   boundaries + (cycles so far)
	    (dolist (matrix-row mat-z)
	      (setf span (in-span (cdr matrix-row) mat-b))
	      (if (null span) ;; not in span
		  (setf mat-b (insert-vector-in-matrix (cdr matrix-row)
						       mat-b row-b)
			mat-h (insert-vector-in-matrix (cdr matrix-row)
						       mat-h row-h)
			row-b (1+ row-b)
			row-h (1+ row-h))))
	  
	    ;; now convert mat-h to a set of polynomials
	    (matrix-to-set mat-h)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; betti:
;; arguments: int homological, int internal, quotient
;; return:
;; betti number of cohomology of Lambda (or Lambda/im(theta)
;; if quotient is t) in bidegree (hom,int).
(defun betti (homological internal &optional quotient)
  "Compute betti number (i.e., dim of cohomology) of Lambda in
bidegree (HOM, INT).  If optional argument QUOTIENT is non-nil,
compute betti number of Lambda/im(theta) in bidegree (HOM,INT)."
  (let (mat-b      ;; matrix of boundaries
	mat-z      ;; matrix of cycles
	(basis (lambda-basis homological internal)))
    (if (> homological internal)
	0
      (if (= homological internal)
	  1
	(progn
	  (setf mat-b (boundary-matrix homological internal quotient)
		mat-z (boundary-matrix (1+ homological) internal quotient))
	  (if quotient (setf basis (strip-odds-from-poly basis)))
	  (setf basis (odd-ending basis))
	  (- (length basis) (length mat-z) (length mat-b)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; theta-inverse: 
;; arguments: poly
;; return:
;; NULL if coboundary(poly) is not all odd
;; theta^{-1}(coboundary(poly)) if it is all odd
;;
;; that is, send each lambda_{2n+1} in poly to lambda_n
(defun theta-inverse (poly)
  "nil if coboundary(poly) is not all odd, and
theta^{-1}(coboundary(poly)) if it is all odd."
  (when (null (strip-odds-from-poly poly))
    (loop for mono in poly
	  collect (loop for n in mono
			collect (floor n 2)))))
  
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; boundaries-all:
;; arguments: homological internal &optional quotient
;; return:
;; basis for all of the boundaries in bidegree (hom,int), returned
;; as a list of lists: cdr of each element is the boundary, car is
;; what it's a boundary of.
(defun boundaries-all (homological internal &optional quotient)
  "basis for the boundaries in bidegree (HOM, INT).  If optional arg
QUOTIENT is non-nil, basis for (boundaries intersect Lambda/im theta).
list of lists: cdr of each element is a boundary, car is what it's
a boundary of"
  (let ((other-basis (lambda-basis (1- homological) internal))
	(mat nil)
	(mat2 nil)
	(index 1))

    (if quotient
	(setf other-basis (strip-odds-from-poly other-basis)))
    ;; compute boundaries and store in mat

    (loop for mono in other-basis do
	  (setf mat (insert-vector-in-matrix (coboundary-mono mono) mat index)
		mat2 (insert-vector-in-matrix (list mono) mat2 index)
		index (1+ index)))

    (multiple-value-bind (reduced polys) (row-reduce mat mat2)
      (loop for row in reduced
	    for poly in polys
	    unless (null (cdr row))
	    collect (cons (cdr poly) (cdr row))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; homological:
;; arguments: poly
;; return:
;; homological degree of poly.  assumes that poly is homogeneous.
(defun homological (poly)
  "homological degree of POLY"
  (length (car poly)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; internal:
;; arguments: poly
;; return:
;; internal degree of poly.  assumes that poly is homogeneous.
(defun internal (poly)
  "internal degree of POLY"
  (+ (homological poly) (apply '+ (car poly))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; boundaryp:
;; arguments: poly
;; return:
;; t if poly is a boundary, nil otherwise
(defun boundaryp (poly &optional quotient)
  "t if POLY is a boundary, nil otherwise.  If optional arg
QUOTIENT is non-nil, t if POLY is a boundary in Lambda/im(theta)"
  (let ((hom (homological poly))
	(int (internal poly))
	(poly-copy poly)
	bdry-matrix
	index-matrix
	rows)
    (unless (= hom int)
      (setf rows
	    (loop for row in (boundaries-all hom int quotient)
		  count row into i do
		  (push (cons i (cdr row)) bdry-matrix)
		  (push (cons i (car row)) index-matrix)
		  finally (return i)))
      (if quotient
	  (setf poly-copy (strip-odds-from-poly poly-copy)))
      (push (cons (1+ rows) poly-copy) bdry-matrix)
      (push (cons (1+ rows) nil) index-matrix)
      (multiple-value-bind (a b) (row-reduce bdry-matrix index-matrix)
	(if (null (assoc (1+ rows) a))
	    (cdr (assoc (1+ rows) b)))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; boundary-preimage:
;; arguments: mono
;; return:
;; list of basis elements whose boundary contains mono as a summand
(defun boundary-preimage (mono)
  "list of basis elements whose boundary contains MONO as a summand"
  (let ((hom (homological (list mono)))
	(int (internal (list mono)))
	basis)
    (setf basis (lambda-basis (1- hom) int))
    (loop for x in basis
	  if (member mono (coboundary-mono x) :test #'equal)
	  collect x)))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; bss:
;; arguments: min-hom max-hom min-int max-int
;; return:
;; nil
;; print results of computations in the BSS in the given range
;; of dimensions (hom is homological, int is internal).  For each
;; class in the ss, print its bidegree, whether a perm cycle or 
;; supporting a differential.
;;
;; do this by computing H*(quotient) and examining
;; the cohomology classes and their lifts to Lambda.
;;
;; note that the internal degree must be even for there to be
;; a differential starting in bidegree (hom, int)
(defun bss (min-hom max-hom min-int max-int
		    &key verbose reverse differentials
		    no-gc-warnings file)
  "compute what's going on in the Bockstein spectral sequence as
homological degree ranges from MIN-HOM to MAX-HOM, and internal
degree ranges from MIN-INT to MAX-INT.

If key VERBOSE is non-nil, print some extra information (like the
current bidegree).
If key REVERSE is non-nil, use opposite sorting for monomials.
If key DIFFERENTIALS is non-nil, just report on the differentials.
If key NO-GC-WARNINGS is non-nil, don't print warnings about garbage
collection -- just automatically garbage collect when necessary.
If key FILE is non-nil, send output there."
  (if (> 0 min-hom)
      (setf min-hom 0))
  (if (< min-int min-hom)
      (setf min-int min-hom))
  (if (or (> min-hom max-hom)
	  (> min-int max-int))
      (princ "empty range.  exiting...~%")
    (let ((stream (if file (open file :direction :output
				      :if-exists :append
				      :if-does-not-exist :create)
		    t)))
      (unwind-protect
	  (let (cobdry
		index
		funny
		(*global-gc-behavior* (if no-gc-warnings :auto :auto-and-warn))
		(lex-order (if reverse #'> #'<))
		(jump (if differentials 2 1)))
	    (princ #\Newline)
	    (when (and differentials (oddp min-int))
	      (setf min-int (* 2 (floor min-int 2))))
	    (do ((hom min-hom (1+ hom)))
		((> hom max-hom))
	      (do ((int min-int (+ int jump)))
		  ((> int max-int))
		(if verbose
		    (multiple-value-bind
			(sec min hour date month) (get-decoded-time)
		      (format stream "starting bidegree [~A,~A], ~A:~2,'0D:~2,'0D on ~A/~A~%"
			      hom int hour min sec month date)))
		;; every element in H*(quotient) either lifts to a
		;; cycle in Lambda, in which case it's a permanent
		;; cycle, or lifts to something whose coboundary is
		;; in the image of theta, in which case it
		;; corresponds to a differential.  So go through
		;; quotient term by term.
		(dolist (poly (ext hom int t))
		  (setf cobdry (coboundary poly)
			index 0)
		  (loop until (null (theta-inverse cobdry)) do
			(setf cobdry (theta-inverse cobdry)
			      index (1+ index)))
		  (when (and (< 0 index)
			     (boundaryp cobdry))
		    (setf index (1- index)
			  funny t))
		  ;; index is the number of times we could apply
		  ;; theta-inverse.
		  (unless reverse
		    (setf poly (nreverse poly)
			  cobdry (nreverse cobdry)))
		  (if (zerop index)
		      (unless differentials
			(progn
			  (format stream "~A, ~A: " hom int)
			  (format stream "permanent cycle: ~A~%"
				  (poly-to-string poly :no-sort t))))
		    (progn
		      (format stream "bidegree (~A, ~A): " hom int)
		      (if funny
			  (format stream "funny?"))
		      (format stream "~%  d_~A(~A) = ~A~%"
			      index
			      (poly-to-string poly :no-sort t)
			      (poly-to-string cobdry :no-sort t)))))))
	    (if file (close stream)))
	(if file (close stream))))))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;; print-betti:
;; arguments: min-hom max-hom min-int max-int
;; return:
;; nil
;; print betti numbers for cohomology of Lambda in the given range
;; of dimensions (hom is homological, int is internal).
(defun print-betti (min-hom max-hom min-int max-int
			    &key verbose reverse quotient)
  "print betti numbers for cohomology of Lambda (or cohomology of
Lambda/im(theta) if QUOTIENT is t) as homological degree ranges
from MIN-HOM to MAX-HOM, and internal degree ranges from MIN-INT
to MAX-INT.

If key QUOTIENT is non-nil, work in Lambda/im(theta) instead of Lambda.
If key VERBOSE is non-nil, print some extra information (like the
bidegrees where the betti number is zero).
If key REVERSE is non-nil, use opposite sorting for monomials."
  (if (> 0 min-hom)
      (setf min-hom 0))
  (if (< min-int min-hom)
      (setf min-int min-hom))
  (if (or (> min-hom max-hom)
	  (> min-int max-int))
      (princ "empty range.  exiting...~%")
    (let (b
	  (*global-gc-behavior* (if no-gc-warnings :auto :auto-and-warn))
	  (lex-order (if reverse #'>= #'<=)))
      (princ #\Newline)
      (format t "hom int  betti~%---------------~%")
      (do ((hom min-hom (1+ hom)))
	  ((> hom max-hom))
	(do ((int min-int (1+ int)))
	    ((> int max-int))
	  (setf b (betti hom int quotient))
	  (if verbose
	      (multiple-value-bind
		  (sec min hour date month) (get-decoded-time)
		(format t " ~A ~3D: ~A    ~A:~2,'0D:~2,'0D on ~A/~A~%"
			hom int b hour min sec month date))
	    (if (< 0 b)
		(format t " ~A ~3D: ~A~%" hom int b))))))))

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

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

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