;;; Carlo C. Maley 9/20/96

;;; These procedures trace through a tree and gather diversity statistics across time
;;; for both lineages and taxa.  Also for both perfect sampling and imperfect sampling.


(in-package 'USER)

(format t "start ")
(defvar *time* 0)
(defvar *para-taxa* ())
(defvar *num-para-taxa* 0)
(format t "1a ")
;(defvar *unclassified-para-taxa* 0) just (svref *para-taxa* 0)
(defvar *no-basal-num-para-taxa* 0)
(defvar *hard-taxa* ())
(defvar *num-hard-taxa* -1)
(defvar *soft-taxa* ())
(format t "1b ")
(defvar *num-soft-taxa* -1)
(defvar *num-soft-taxa-match* -1)
(defvar *unclassifed-soft-match* -1)
(defvar *syst-soft-taxa* ())
(defvar *num-syst-soft-taxa* -1)
(defvar *num-syst-soft-taxa-match* -1)
(format t "1c ")
(defvar *unclassifed-syst-soft-match* -1)
(defvar *full-taxa* ())
(defvar *num-full-taxa* -1)
(defvar *rand-taxa* ())
(format t "1d ")
(defvar *num-rand-taxa* -1)
(defvar *match-rand-taxa* ())
(defvar *match-num-rand-taxa* -1)
(defvar *match-distrib-rand-taxa* ())
(defvar *match-distrib-num-rand-taxa* -1)
(defvar *match-distrib-soft-taxa* ())
(defvar *match-distrib-num-soft-taxa* -1)
(defvar *diversity* 0)
(format t "1 ")

(defvar *samp-para-taxa* ())
(defvar *samp-num-para-taxa* 0)
;(defvar *samp-unclassified-para-taxa* 0)
(defvar *samp-no-basal-num-para-taxa* 0)
(defvar *samp-hard-taxa* ())
(defvar *samp-num-hard-taxa* -1)
(defvar *samp-soft-taxa* ())
(defvar *samp-num-soft-taxa* -1)
(defvar *samp-num-soft-taxa-match* -1)
(defvar *samp-unclassifed-soft-match* -1)
(defvar *samp-syst-soft-taxa* ())
(defvar *samp-num-syst-soft-taxa* -1)
(defvar *samp-num-syst-soft-taxa-match* -1)
(defvar *samp-unclassifed-syst-soft-match* -1)
(defvar *samp-full-taxa* ())
(defvar *samp-num-full-taxa* -1)
(defvar *samp-rand-taxa* ())
(defvar *samp-num-rand-taxa* -1)
(defvar *samp-match-rand-taxa* ())
(defvar *samp-match-num-rand-taxa* -1)
(defvar *samp-match-distrib-rand-taxa* ())
(defvar *samp-match-distrib-num-rand-taxa* -1)
(defvar *samp-match-distrib-soft-taxa* ())
(defvar *samp-match-distrib-num-soft-taxa* -1)
(defvar *samp-diversity* 0)
(format t "2 ")


;;; I'm recording the presence of the root node without checking to see if it is sampled.
;;; hmm...I guess I am assuming we know the root node was extant at the start.

(defun record-taxa-diversity (mass-ext speciation extinction)
  (with-open-file
   (div-file (format () "output/diversity-~a.log" *run-number*)
		    :direction :output
		    :if-does-not-exist :create
		    :if-exists :overwrite)
   (format div-file "# Run ~a~%" *run-number*)
   (format div-file "# N = ~a  Speciation = ~a  Extinction = ~a~%" (1+ *max-diversity*)
	   speciation extinction)
   (format div-file "# K = ~a  Sampling = ~a  Proportion of taxa = ~a~%" *k*
	   (/ *sample* (* 1.0 PROBSCALE)) (/ *paraphyly-prob* (* 1.0 PROBSCALE)))
   (format div-file "# Extinctions (when, magnitude) in percents: ~a~%" mass-ext)
   (setf (snode-first-sample *root*) 0) ;assume the root lineage was there from start.
   (setf (snode-last-sample *root*) 0) ;assume the root lineage was there from start.

   (format t "Filling in the sampling gaps by inference...~%")
   
   (infer-samples *root*)
   (setf *time* 0)
   (setf *diversity* 1)
   (iter-make-iter-list)
   (iter-insert-node *root*)
   (setf *para-taxa* (make-array *next-paraphyle* :initial-element 0))
   (setf *num-para-taxa* 0)
   (setf *no-basal-num-para-taxa* 0)
   (setf *hard-taxa* (make-array *next-paraphyle* :initial-element 0))
   (setf *num-hard-taxa* 0)
   (setf *soft-taxa* (make-array (1+ *next-soft-monophyle*) :initial-element 0))
   (setf *num-soft-taxa* 0)
   (setf *num-soft-taxa-match* 0)
   (setf *unclassifed-soft-match* 0)
   (setf *syst-soft-taxa* (make-array (1+ *next-syst-monophyle*) :initial-element 0))
   (setf *num-syst-soft-taxa* 0)
   (setf *num-syst-soft-taxa-match* 0)
   (setf *unclassifed-syst-soft-match* 0)
   (setf *full-taxa* (make-array (1+ *next-full-monophyle*) :initial-element 0))
   (setf *num-full-taxa* 0)
   (setf *rand-taxa* (make-array (1+ *next-random-taxon*) :initial-element 0))
   (setf *num-rand-taxa* 0)
   (setf *match-rand-taxa* (make-array (1+ *next-match-random-taxon*) :initial-element 0))
   (setf *match-num-rand-taxa* 0)
   (setf *match-distrib-rand-taxa* (make-array *next-paraphyle*
					       :initial-element 0))
   (setf *match-distrib-num-rand-taxa* 0)
   (setf *match-distrib-soft-taxa* (make-array (1+ *next-paraphyle*)
					       :initial-element 0))
   (setf *match-distrib-num-soft-taxa* 0)
   (record-node (snode-paraphyle *root*) *para-taxa* *num-para-taxa*)
   (record-node (snode-hard-monophyle *root*) *hard-taxa* *num-hard-taxa*)
   (record-node (snode-soft-monophyle *root*) *soft-taxa* *num-soft-taxa*)
   (match-record-node (snode-soft-monophyle *root*) *soft-taxa*
		      *unclassifed-soft-match* *num-soft-taxa-match*
		      *soft-forbidden*)
   (record-node (snode-syst-monophyle *root*) *syst-soft-taxa* *num-syst-soft-taxa*)
   (match-record-node (snode-syst-monophyle *root*) *syst-soft-taxa*
		      *unclassifed-syst-soft-match* *num-syst-soft-taxa-match*
		      *syst-soft-forbidden*)
   (record-node (snode-full-monophyle *root*) *full-taxa* *num-full-taxa*)
   (record-node (snode-rand-taxon *root*) *rand-taxa* *num-rand-taxa*)
   (record-node (snode-match-rand *root*) *match-rand-taxa* *match-num-rand-taxa*)
   (record-node (snode-distrib-rand *root*) *match-distrib-rand-taxa*
		*match-distrib-num-rand-taxa*)
   (record-node (snode-distrib-soft *root*) *match-distrib-soft-taxa*
		*match-distrib-num-soft-taxa*)
   ; sampled data
   (setf *samp-diversity* 1)
   (backup-iter-make-iter-list)
   (backup-iter-insert-node *root*)
   (setf *samp-para-taxa* (make-array *next-paraphyle* :initial-element 0))
   (setf *samp-num-para-taxa* 0)
   (setf *samp-no-basal-num-para-taxa* 0)
   (setf *samp-hard-taxa* (make-array *next-paraphyle* :initial-element 0))
   (setf *samp-num-hard-taxa* 0)
   (setf *samp-soft-taxa* (make-array (1+ *next-soft-monophyle*)
				      :initial-element 0))
   (setf *samp-num-soft-taxa* 0)
   (setf *samp-num-soft-taxa-match* 0)
   (setf *samp-unclassifed-soft-match* 0)
   (setf *samp-syst-soft-taxa* (make-array (1+ *next-syst-monophyle*)
					   :initial-element 0))
   (setf *samp-num-syst-soft-taxa* 0)
   (setf *samp-num-syst-soft-taxa-match* 0)
   (setf *samp-unclassifed-syst-soft-match* 0)
   (setf *samp-full-taxa* (make-array (1+ *next-full-monophyle*) :initial-element 0))
   (setf *samp-num-full-taxa* 0)
   (setf *samp-rand-taxa* (make-array (1+ *next-random-taxon*) :initial-element 0))
   (setf *samp-num-rand-taxa* 0)
   (setf *samp-match-rand-taxa* (make-array (1+ *next-match-random-taxon*)
					    :initial-element 0))
   (setf *samp-match-num-rand-taxa* 0)
   (setf *samp-match-distrib-rand-taxa* (make-array *next-paraphyle*
					       :initial-element 0))
   (setf *samp-match-distrib-num-rand-taxa* 0)
   (setf *samp-match-distrib-soft-taxa* (make-array (1+ *next-paraphyle*)
					       :initial-element 0))
   (setf *samp-match-distrib-num-soft-taxa* 0)
   (record-node (snode-paraphyle *root*) *samp-para-taxa* *samp-num-para-taxa*)
   (record-node (snode-hard-monophyle *root*) *samp-hard-taxa* *samp-num-hard-taxa*)
   (record-node (snode-soft-monophyle *root*) *samp-soft-taxa* *samp-num-soft-taxa*)
   (match-record-node (snode-soft-monophyle *root*) *samp-soft-taxa*
		      *samp-unclassifed-soft-match* *samp-num-soft-taxa-match*
		      *soft-forbidden*)
   (record-node (snode-syst-monophyle *root*) *samp-syst-soft-taxa*
		*samp-num-syst-soft-taxa*)
   (match-record-node (snode-syst-monophyle *root*) *samp-syst-soft-taxa*
		*samp-unclassifed-syst-soft-match* *samp-num-syst-soft-taxa-match*
		*syst-soft-forbidden*)
   (record-node (snode-full-monophyle *root*) *samp-full-taxa* *samp-num-full-taxa*)
   (record-node (snode-rand-taxon *root*) *samp-rand-taxa* *samp-num-rand-taxa*)
   (record-node (snode-match-rand *root*) *samp-match-rand-taxa*
		*samp-match-num-rand-taxa*)
   (record-node (snode-distrib-rand *root*) *samp-match-distrib-rand-taxa*
		*samp-match-distrib-num-rand-taxa*)
   (record-node (snode-distrib-soft *root*) *samp-match-distrib-soft-taxa*
		*samp-match-distrib-num-soft-taxa*)

   (format t "Gathering Diversity stats...~%")
   
   (recap div-file))
  (record-classification-stats))

(format t "3 ")

;;; Loop through the extant list checking the age.  If age = *time* then remove
;;; that node and add any offspring.
(defun recap (div-file)
  (loop
   (when (> *time* *clock*) (return ()))
   (recap-step)
   (record-step-diversity div-file)))

(format t "4 ")

;;; For each classification system I'm printing the number of taxa present, the "coverage"
;;; (proportion of the known tree being counted), and the average size of a taxon.  
(defun record-step-diversity (div-file)
  (format div-file "~a~a~a~a~a~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~a~a~a~4,3F~a~4,3F~%"
	  *time* #\tab
	  *diversity* #\tab
	  *samp-diversity* #\tab 
	  
	  *num-para-taxa*  #\tab
	  1.000  #\tab 
	  (if (= 0 *num-para-taxa*) 0
	    (/ *diversity*  *num-para-taxa*)) #\tab 
	  
	  *samp-num-para-taxa*  #\tab
	  1.000  #\tab
	  (if (= 0 *samp-num-para-taxa*) 0
	    (/ *samp-diversity*  *samp-num-para-taxa*)) #\tab 
	  
	  *num-hard-taxa*  #\tab
	  (- 1.0 (/ (svref *hard-taxa*  0) (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-hard-taxa*) 0
	    (/ (- *diversity* (svref *hard-taxa*  0)) *num-hard-taxa*)) #\tab 
	  
	  *samp-num-hard-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-hard-taxa*  0)
		    (* 1.0 *samp-diversity*)))	 #\tab 
	  (if (= 0 *samp-num-hard-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-hard-taxa*  0))
	       *samp-num-hard-taxa*)) #\tab 
	  
	  *num-soft-taxa*  #\tab
	  (- 1.0 (/ (svref *soft-taxa*  0) (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-soft-taxa*) 0
	    (/ (- *diversity* (svref *soft-taxa*  0)) *num-soft-taxa*)) #\tab 
	  
	  *samp-num-soft-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-soft-taxa*  0)
		    (* 1.0 *samp-diversity*)))	 #\tab 
	  (if (= 0 *samp-num-soft-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-soft-taxa*  0))
	       *samp-num-soft-taxa*)) #\tab 
	  
	  *num-syst-soft-taxa*  #\tab
	  (- 1.0 (/ (svref *syst-soft-taxa*  0) (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-syst-soft-taxa*) 0
	    (/ (- *diversity* (svref *syst-soft-taxa*  0)) *num-syst-soft-taxa*)) #\tab 
	  
	  *samp-num-syst-soft-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-syst-soft-taxa*  0)
		    (* 1.0 *samp-diversity*))) #\tab 
	  (if (= 0 *samp-num-syst-soft-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-syst-soft-taxa*  0))
	     *samp-num-syst-soft-taxa*)) #\tab 
	  
	  *num-full-taxa*  #\tab
	  (- 1.0 (/ (svref *full-taxa*  0) (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-full-taxa*) 0
	    (/ (- *diversity* (svref *full-taxa*  0)) *num-full-taxa*)) #\tab 
	  
	  *samp-num-full-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-full-taxa*  0)
		    (* 1.0 *samp-diversity*)))	 #\tab 
	  (if (= 0 *samp-num-full-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-full-taxa*  0))
	       *samp-num-full-taxa*)) #\tab 
	  
	  *num-rand-taxa*  #\tab
	  (- 1.0 (/ (svref *rand-taxa*  0) (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-rand-taxa*) 0
	    (/ (- *diversity* (svref *rand-taxa*  0)) *num-rand-taxa*)) #\tab 
	  
	  *samp-num-rand-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-rand-taxa*  0)
		    (* 1.0 *samp-diversity*))) #\tab 
	  (if (= 0 *samp-num-rand-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-rand-taxa*  0))
	       *samp-num-rand-taxa*)) #\tab 

	  *num-soft-taxa-match*  #\tab
	  (- 1.0 (/ *unclassifed-soft-match* (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-soft-taxa-match*) 0
	    (/ (- *diversity* *unclassifed-soft-match*) *num-soft-taxa-match*)) #\tab 
	  
	  *samp-num-soft-taxa-match*  #\tab
	  (- 1.0 (/ *samp-unclassifed-soft-match*
		    (* 1.0 *samp-diversity*)))	 #\tab 
	  (if (= 0 *samp-num-soft-taxa-match*) 0
	    (/ (- *samp-diversity* *samp-unclassifed-soft-match*)
	       *samp-num-soft-taxa-match*)) #\tab 
	  
	  *num-syst-soft-taxa-match*  #\tab
	  (- 1.0 (/ *unclassifed-syst-soft-match*
		    (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *num-syst-soft-taxa-match*) 0
	    (/ (- *diversity* *unclassifed-syst-soft-match*)
	       *num-syst-soft-taxa-match*)) #\tab 
	  
	  *samp-num-syst-soft-taxa-match*  #\tab
	  (- 1.0 (/ *samp-unclassifed-syst-soft-match*
		    (* 1.0 *samp-diversity*))) #\tab 
	  (if (= 0 *samp-num-syst-soft-taxa-match*) 0
	    (/ (- *samp-diversity* *samp-unclassifed-syst-soft-match*)
	       *samp-num-syst-soft-taxa-match*)) #\tab 
	  
	  *match-num-rand-taxa*  #\tab
	  (- 1.0 (/ (svref *match-rand-taxa*  0)
		    (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *match-num-rand-taxa*) 0
	    (/ (- *diversity* (svref *match-rand-taxa*  0)) *match-num-rand-taxa*)) #\tab 
	  
	  *samp-match-num-rand-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-match-rand-taxa*  0)
		    (* 1.0 *samp-diversity*))) #\tab 
	  (if (= 0 *samp-match-num-rand-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-match-rand-taxa*  0))
	       *samp-match-num-rand-taxa*)) #\tab 

	  *no-basal-num-para-taxa*  #\tab
	  (- 1.0 (/ (svref *para-taxa*  0) (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *no-basal-num-para-taxa*) 0
	    (/ (- *diversity*  (svref *para-taxa*  0))
	       *no-basal-num-para-taxa*)) #\tab 
	  
	  *samp-no-basal-num-para-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-para-taxa*  0)
		    (* 1.0 *samp-diversity*))) #\tab
	  (if (= 0 *samp-no-basal-num-para-taxa*) 0
	    (/ (- *samp-diversity*  (svref *samp-para-taxa*  0))
	       *samp-no-basal-num-para-taxa*)) #\tab 
	  
	  *match-distrib-num-rand-taxa*  #\tab
	  (- 1.0 (/ (svref *match-distrib-rand-taxa*  0)
		    (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *match-distrib-num-rand-taxa*) 0
	    (/ (- *diversity* (svref *match-distrib-rand-taxa*  0))
	       *match-distrib-num-rand-taxa*)) #\tab 
	  
	  *samp-match-distrib-num-rand-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-match-distrib-rand-taxa*  0)
		    (* 1.0 *samp-diversity*))) #\tab 
	  (if (= 0 *samp-match-distrib-num-rand-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-match-distrib-rand-taxa*  0))
	       *samp-match-distrib-num-rand-taxa*)) #\tab 

	  *match-distrib-num-soft-taxa*  #\tab
	  (- 1.0 (/ (svref *match-distrib-soft-taxa*  0)
		    (* 1.0 *diversity*))) #\tab 
	  (if (= 0 *match-distrib-num-soft-taxa*) 0
	    (/ (- *diversity* (svref *match-distrib-soft-taxa*  0))
	       *match-distrib-num-soft-taxa*)) #\tab 
	  
	  *samp-match-distrib-num-soft-taxa*  #\tab
	  (- 1.0 (/ (svref *samp-match-distrib-soft-taxa*  0)
		    (* 1.0 *samp-diversity*))) #\tab 
	  (if (= 0 *samp-match-distrib-num-soft-taxa*) 0
	    (/ (- *samp-diversity* (svref *samp-match-distrib-soft-taxa*  0))
	       *samp-match-distrib-num-soft-taxa*)) 

	  ))

(format t "5 ")

;;; We assume that if we have the last sample of an acestor at time t and the first
;;; sample of an descendant at time t + n, then the splitting event must have have
;;; happened at time t + (n - 1).  That is, we don't believe a split happened until
;;; we see evidence of it.

(defun infer-samples (node)
  (unless (null node)
    (let ((left (snode-left node))
	  (right (snode-right node)))
      (unless (and (null left)
	       (null right))
      	(infer-samples left)
	(infer-samples right)
	(let* ((left-time (snode-first-sample left))
	       (right-time (snode-first-sample right))
	       (branch-time (get-branch-time node left-time right-time)))
	  (when (> branch-time 0)
	    (when (< (snode-first-sample node) 0)
	      (setf (snode-first-sample node) (1- branch-time)))
	    (setf (snode-last-sample node) (1- branch-time))
	    (unless (< left-time 0)
	      (setf (snode-first-sample left) branch-time))
	    (unless (< right-time 0)
	      (setf (snode-first-sample right) branch-time))))))))
	    
	  



		     
(format t "done~%")