diff --git a/src/anromero/persistent-homology.lisp b/src/anromero/persistent-homology.lisp new file mode 100644 index 0000000..7296199 --- /dev/null +++ b/src/anromero/persistent-homology.lisp @@ -0,0 +1,520 @@ +;; PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY +;; PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY +;; PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY PERSISTENT-HOMOLOGY + +(IN-PACKAGE #:cat) + +(provide "persistent-homology") + +(defun eff-prst-hmlg-basis-dvs-quotient (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (the list + (let* ((degr (+ p q)) + (fltr-p-basis (fltrd-basis fltrcm degr p)) + ;;(p-basis-l (length fltr-p-basis)) + (num-mtrx (spsq-num-mtrx fltrcm (1+ p) p q)) + (den-mtrx (spsq-den-mtrx fltrcm r p q)) + (cmbn-list nil) + (divs nil)) + (declare + (list fltr-p-basis cmbn-list divs) + ;;(type fixnum degr p-basis-l) + (type matrix num-mtrx den-mtrx) + ) + (progn + (if (= 0 (array-total-size num-mtrx)) + (progn + (setq cmbn-list nil) + (setq divs nil)) + (if (= 0 (array-total-size den-mtrx)) + (progn + (setq cmbn-list (gnrt-mtrx-to-list num-mtrx)) + (setq divs (make-list (length cmbn-list) :initial-element 0))) + (let ((quotient-list (mtrx-quotient num-mtrx den-mtrx))) + (declare (list quotient-list)) + (setq cmbn-list (first quotient-list) + divs (second quotient-list))))) + (list + (mapcar + #'(lambda (int-list) + (declare (list int-list)) + (the cmbn + (let* ((cmbn (cmbn degr)) + (cmpr (cmpr fltrcm))) + (declare + (type cmbn cmbn) + (type cmprf cmpr)) + (do ((mark1 int-list (cdr mark1)) + (mark2 fltr-p-basis (cdr mark2))) + ((or (endp mark1) (endp mark2))) + (declare (type list mark1 mark2)) + (if (not (= 0 (car mark1))) + (setq cmbn (2cmbn-add cmpr cmbn + (cmbn degr (car mark1) (car mark2)))))) + cmbn))) + cmbn-list) + divs))))) + + +(defun eff-prst-hmlg-basis-dvs-total (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (the list + (let* ((degr (+ p q)) + (fltr-p-basis (fltrd-basis fltrcm degr p)) + (p-basis-l (length fltr-p-basis)) + ;;(fltr-p-1-list (fltr-basis-gnrt-list fltrcm degr (1- p))) + (z1-list (fltr-chcm-z-gnrt-list fltrcm (1+ p) p q)) + (z1-mtrx (gnrt-list-to-mtrx z1-list)) + (z2-list (fltr-chcm-z-gnrt-list fltrcm (1- r) (1- (+ p r)) (+ 2 (- q r)))) + (z2-mtrx (gnrt-list-to-mtrx z2-list)) + ;;(fltr-p-1-mtrx (gnrt-list-to-mtrx fltr-p-1-list)) + (dffr-mtrx (flcc-dffr-mtrx fltrcm (1+ degr) (1- (+ p r)))) + (nil-mtrx + #-ACLPC (make-array (list 0 0) + :element-type 'fixnum + :initial-element 0) + #+ACLPC (make-array (list 0 0) + :element-type 'fixnum)) + (dffr-z2-mtrx + (if (or (= 0 (array-total-size z2-mtrx)) (= 0 (array-total-size dffr-mtrx))) + nil-mtrx + (mtrx-prdc dffr-mtrx z2-mtrx))) + (bnd-z2-line-n p-basis-l) + (bnd-z2-column-n (column-number dffr-z2-mtrx)) + (bnd-z2-mtrx + (if (or (= 0 bnd-z2-line-n) (= 0 bnd-z2-column-n)) + nil-mtrx + (submatrix dffr-z2-mtrx 0 (1- bnd-z2-line-n) 0 (1- bnd-z2-column-n)))) + (num-mtrx z1-mtrx) + (den-mtrx bnd-z2-mtrx) + (cmbn-list nil) + (divs nil)) + (declare + (list fltr-p-basis cmbn-list divs) + (type fixnum degr p-basis-l) + (type matrix num-mtrx den-mtrx)) + (progn + (if (= 0 (array-total-size num-mtrx)) + (progn + (setq cmbn-list nil) + (setq divs nil)) + (if (= 0 (array-total-size den-mtrx)) + (progn + (setq cmbn-list (gnrt-mtrx-to-list num-mtrx)) + (setq divs (make-list (length cmbn-list) :initial-element 0))) + (let ((quotient-list (mtrx-quotient num-mtrx den-mtrx))) + (declare (list quotient-list)) + (setq cmbn-list (first quotient-list) + divs (second quotient-list))))) + (list + (mapcar + #'(lambda (int-list) + (declare (list int-list)) + (the cmbn + (let* ((cmbn (cmbn degr)) + (cmpr (cmpr fltrcm))) + (declare + (type cmbn cmbn) + (type cmprf cmpr)) + (do ((mark1 int-list (cdr mark1)) + (mark2 fltr-p-basis (cdr mark2))) + ((or (endp mark1) (endp mark2))) + (declare (type list mark1 mark2)) + (if (not (= 0 (car mark1))) + (setq cmbn (2cmbn-add cmpr cmbn + (cmbn degr (car mark1) (car mark2)))))) + cmbn))) + cmbn-list) + divs))))) + + +(defun eff-prst-hmlg-basis-dvs-bar (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (the list + (let* ((degr (+ p q)) + (fltr-p-basis (fltrd-basis fltrcm degr p)) + ;;(p-basis-l (length fltr-p-basis)) + (num-mtrx (spsq-den-mtrx fltrcm r p q)) + (den-mtrx (spsq-den-mtrx fltrcm (- r 1) p q)) + (cmbn-list nil) + (divs nil)) + (declare + (list fltr-p-basis cmbn-list divs) + ;;(type fixnum degr p-basis-l) + (type matrix num-mtrx den-mtrx)) + (progn + (if (= 0 (array-total-size num-mtrx)) + (progn + (setq cmbn-list nil) + (setq divs nil)) + ;; (if (= 0 (array-total-size den-mtrx)) + ;; (progn + ;; (setq cmbn-list (gnrt-mtrx-to-list num-mtrx)) + ;; (setq divs (make-list (length cmbn-list) :initial-element 0))) + (let ((quotient-list (mtrx-quotient num-mtrx den-mtrx))) + (declare (list quotient-list)) + (setq cmbn-list (first quotient-list) + divs (second quotient-list)))) ;;) + (list + (mapcar + #'(lambda (int-list) + (declare (list int-list)) + (the cmbn + (let* ((cmbn (cmbn degr)) + (cmpr (cmpr fltrcm))) + (declare + (type cmbn cmbn) + (type cmprf cmpr)) + (do ((mark1 int-list (cdr mark1)) + (mark2 fltr-p-basis (cdr mark2))) + ((or (endp mark1) (endp mark2))) + (declare (type list mark1 mark2)) + (if (not (= 0 (car mark1))) + (setq cmbn (2cmbn-add cmpr cmbn + (cmbn degr (car mark1) (car mark2)))))) + cmbn))) + cmbn-list) + divs))))) + + + + + +(defun eff-prst-hmlg-gnrts-total (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((basis-dvs (eff-prst-hmlg-basis-dvs-total fltrcm r p q)) + (basis (first basis-dvs)) + (dvs (second basis-dvs)) + (rslt nil)) + (declare (type list basis-dvs basis dvs rslt)) + (do ((mark1 basis (cdr mark1)) + (mark2 dvs (cdr mark2))) + ((endp mark1)) + (if (not (eql 1 (car mark2))) + (setf rslt (append rslt (list (car mark1)))))) + rslt)) + +(defun eff-prst-hmlg-gnrts-quotient (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((basis-dvs (eff-prst-hmlg-basis-dvs-quotient fltrcm r p q)) + (basis (first basis-dvs)) + (dvs (second basis-dvs)) + (rslt nil)) + (declare (type list basis-dvs basis dvs rslt)) + (do ((mark1 basis (cdr mark1)) + (mark2 dvs (cdr mark2))) + ((endp mark1)) + (if (not (eql 1 (car mark2))) + (setf rslt (append rslt (list (car mark1)))))) + rslt)) + +(defun eff-prst-hmlg-gnrts-bar (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((basis-dvs (eff-prst-hmlg-basis-dvs-bar fltrcm r p q)) + (basis (first basis-dvs)) + (dvs (second basis-dvs)) + (rslt nil)) + (declare (type list basis-dvs basis dvs rslt)) + (do ((mark1 basis (cdr mark1)) + (mark2 dvs (cdr mark2))) + ((endp mark1)) + (if (not (eql 1 (car mark2))) + (setf rslt (append rslt (list (car mark1)))))) + rslt)) + + +(defun eff-prst-hmlg-group-total (fltrcm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((basis-dvs (eff-prst-hmlg-basis-dvs-total fltrcm r p q)) + ;; (cmbns (first basis-dvs)) + (divs (second basis-dvs))) + (declare (list basis-dvs divs)) + + (dolist (item divs) + (declare (type fixnum item)) + (if (not (eql 1 item)) + (progn + (format t "~2%Component Z") + (unless (zerop item) + (format t "/~DZ" item)) + ))) + (terpri))) + +(defun eff-prst-hmlg-group-quotient (fltrcm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((basis-dvs (eff-prst-hmlg-basis-dvs-quotient fltrcm r p q)) + ;;(cmbns (first basis-dvs)) + (divs (second basis-dvs))) + (declare (list basis-dvs divs)) + + (dolist (item divs) + (declare (type fixnum item)) + (if (not (eql 1 item)) + (progn + (format t "~2%Component Z") + (unless (zerop item) + (format t "/~DZ" item)) + ))) + (terpri))) + +(defun eff-prst-hmlg-group-bar (fltrcm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((basis-dvs (eff-prst-hmlg-basis-dvs-bar fltrcm r p q)) + ;; (cmbns (first basis-dvs)) + (divs (second basis-dvs))) + (declare (list basis-dvs divs)) + + (dolist (item divs) + (declare (type fixnum item)) + (if (not (eql 1 item)) + (progn + (format t "~2%Component Z") + (unless (zerop item) + (format t "/~DZ" item)) + ))) + (terpri))) + +(defun prst-hmlg-group-total-rpq (fltrcm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((hmtp-eq (efhm fltrcm)) + (eff-chcm (rbcc hmtp-eq)) + ;;(degr (+ p q)) + ) + (declare + (type homotopy-equivalence hmtp-eq) + (type chain-complex eff-chcm) + ;;(type fixnum degr) + ) + (progn + (if (not (typep eff-chcm 'filtered-chain-complex)) + (translate-filtration hmtp-eq)) + (eff-prst-hmlg-group-total eff-chcm r p q)))) + +(defun prst-hmlg-group-quotient-rpq (fltrcm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((hmtp-eq (efhm fltrcm)) + (eff-chcm (rbcc hmtp-eq)) + ;;(degr (+ p q)) + ) + (declare + (type homotopy-equivalence hmtp-eq) + (type chain-complex eff-chcm) + ;; (type fixnum degr) + ) + (progn + (if (not (typep eff-chcm 'filtered-chain-complex)) + (translate-filtration hmtp-eq)) + (eff-prst-hmlg-group-quotient eff-chcm r p q)))) + +(defun prst-hmlg-group-bar-rpq (fltrcm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((hmtp-eq (efhm fltrcm)) + (eff-chcm (rbcc hmtp-eq)) + ;;(degr (+ p q)) + ) + (declare + (type homotopy-equivalence hmtp-eq) + (type chain-complex eff-chcm) + ;;(type fixnum degr) + ) + (progn + (if (not (typep eff-chcm 'filtered-chain-complex)) + (translate-filtration hmtp-eq)) + (eff-prst-hmlg-group-bar eff-chcm r p q)))) + +(defun prst-hmlg-gnrts-total-rpq (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((hmtp-eq (efhm fltrcm)) + (eff-chcm (rbcc hmtp-eq)) + ;;(degr (+ p q)) + ) + (declare + (type homotopy-equivalence hmtp-eq) + (type chain-complex eff-chcm) + ;;(type fixnum degr) + ) + (progn + (if (not (typep eff-chcm 'filtered-chain-complex)) + (translate-filtration hmtp-eq)) + (eff-prst-hmlg-gnrts-total eff-chcm r p q)))) + +(defun prst-hmlg-gnrts-quotient-rpq (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((hmtp-eq (efhm fltrcm)) + (eff-chcm (rbcc hmtp-eq)) + ;;(degr (+ p q)) + ) + (declare + (type homotopy-equivalence hmtp-eq) + (type chain-complex eff-chcm) + ;;(type fixnum degr) + ) + (progn + (if (not (typep eff-chcm 'filtered-chain-complex)) + (translate-filtration hmtp-eq)) + (eff-prst-hmlg-gnrts-quotient eff-chcm r p q)))) + +(defun prst-hmlg-gnrts-bar-rpq (FltrCm r p q) + (declare + (type filtered-chain-complex fltrcm) + (type fixnum r p q)) + (let* ((hmtp-eq (efhm fltrcm)) + (eff-chcm (rbcc hmtp-eq)) + ;;(degr (+ p q)) + ) + (declare + (type homotopy-equivalence hmtp-eq) + (type chain-complex eff-chcm) + ;;(type fixnum degr) + ) + (progn + (if (not (typep eff-chcm 'filtered-chain-complex)) + (translate-filtration hmtp-eq)) + (eff-prst-hmlg-gnrts-bar eff-chcm r p q)))) + + (defun total-prst-hmlg-gnrts (fltrcm i j n) + (prst-hmlg-gnrts-total-rpq fltrcm (1+ (- j i)) i (- n i) )) + +(defun total-prst-hmlg-group (fltrcm i j n) + (format t "Persistent Homology H^{~D,~D}_~D" i j n) + (prst-hmlg-group-total-rpq fltrcm (1+ (- j i)) i (- n i) )) + + (defun prst-hmlg-gnrts-quotient (fltrcm i j n) + (prst-hmlg-gnrts-quotient-rpq fltrcm (1+ (- j i)) i (- n i) )) + +(defun prst-hmlg-group-quotient (fltrcm i j n) + (format t "Persistent Homology H^{~D,~D}_~D" i j n) + (prst-hmlg-group-quotient-rpq fltrcm (1+ (- j i)) i (- n i) )) + + + (defun prst-hmlg-gnrts (fltrcm i j n) + (prst-hmlg-gnrts-bar-rpq fltrcm (1+ (- j i)) i (- n i) )) + +(defun prst-hmlg-group (fltrcm i j n) + (format t "Persistent Homology B^{~D,~D}_~D" i j n) + (prst-hmlg-group-bar-rpq fltrcm (1+ (- j i)) i (- n i) )) + + + + + + + +;;; EXAMPLES + + +#| + +(defun contains (l e) + (let ((rslt nil)) + (declare (list rslt)) + (dotimes (i (length l)) + (if (eq e (nth i l)) + (setf rslt 't))) + rslt)) + + +(setf triangle + (build-finite-ss + '(a b c + 1 ab (b a) ac (c a) bc (c b) + 2 abc (bc ac ab) ))) + + + +(setf triangle-flin + #'(lambda (degr gnrt) + (declare + (type fixnum degr) + (type gnrt gnrt)) + (let ( + (l1 (list 'a)) + (l2 (list 'b)) + (l3 (list 'c)) + (l4 (list 'ab)) + (l5 (list 'bc)) + (l6 (list 'ac)) + (l7 (list 'abc)) + (rslt -1)) + (declare (list l1 l2 l3 l4 l5 l6 l7) + (type fixnum rslt)) + (if (contains l1 gnrt) (setf rslt 1) + (if (contains l2 gnrt) (setf rslt 2) + (if (contains l3 gnrt) (setf rslt 3) + (if (contains l4 gnrt) (setf rslt 4) + (if (contains l5 gnrt) (setf rslt 5) + (if (contains l6 gnrt) (setf rslt 6) + (if (contains l7 gnrt) (setf rslt 7) + + + ))))))) + rslt))) + + +(setf triangle (change-chcm-to-flcc triangle triangle-flin 'triangle-flin)) + + +(dotimes (p 8) + (print (spsq-group triangle 1 p (- 0 p))) + (print (spsq-group triangle 1 p (- 1 p))) + (print (spsq-group triangle 1 p (- 2 p))) + (print (spsq-group triangle 1 p (- 3 p)))) + + +(dotimes (p 8) + (print (spsq-group triangle 2 p (- 0 p))) + (print (spsq-group triangle 2 p (- 1 p))) + (print (spsq-group triangle 2 p (- 2 p))) + (print (spsq-group triangle 2 p (- 3 p)))) + + +(dotimes (n 3) + (dotimes (p 8) + (print (prst-hmlg-group triangle p p n)) + (print (prst-hmlg-group triangle p (1+ p) n)) + (print (prst-hmlg-group triangle p (+ 2 p) n)) + (print (prst-hmlg-group triangle p (+ 3 p) n)) + (print (prst-hmlg-group triangle p (+ 4 p) n)) + )) + +(dotimes (n 3) + (dotimes (p 8) + (print (total-prst-hmlg-group triangle p p n)) + (print (total-prst-hmlg-group triangle p (1+ p) n)) + (print (total-prst-hmlg-group triangle p (+ 2 p) n)) + (print (total-prst-hmlg-group triangle p (+ 3 p) n)) + (print (total-prst-hmlg-group triangle p (+ 4 p) n)) + )) + + +(print (prst-hmlg-group triangle 1 8 0)) + +|# \ No newline at end of file