;;; This holds the different procedures for classifying a tree.  This includes:

;;; 1. Sepkoski & Kendrick's paraphyly
;;; 2. S&K's hard cull monophyly
;;; 3. S&K's soft cull monophyly
;;; 3.5 Smitth and Patterson's (systematist) soft cull.
;;; 4. A better monophyly to match S&K's paraphyly
;;; 5. Heirarchical (nested) monophyly

;;; To do:
;;;*** the soft cull algorithm is implementing the systematic soft cull, not
;;; the S&K soft cull.  That is, if I want to make something a monophyle, it must be a
;;; taxon node.

(in-package 'USER)

(defvar *paraphyle-nodes* ()) ; array of all the paraphyletic nodes
(defvar *next-para-node* 0) ; pointer into that array.


(defun sum-vect (vector v-length)
  (let ((sum 0))
    (dotimes (i v-length sum)
      (incf sum (svref vector i)))))


(defun classify ()
  (paraphyly)
  (hard-cull *root*)
  (soft-cull *root*)
  (syst-soft-cull *root*)
  (full-monophyly *root*) ; must be called after syst-soft-cull.
  (random-taxa *root*)    ; must be called after full-monophyly.
  (match-soft)
  (match-syst-soft)
  (match-random *root*)
  (match-distrib-rand *root*)
  (match-distrib-soft *root*)
  (nested-monophyly *root*)
  (when *save-trees*
    (record-tree *root* #'snode-id
		 (format () "trees/lineage/lineage-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-paraphyle
		 (format () "trees/paraphyly/para-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-hard-monophyle
		 (format () "trees/hardcull/hard-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-soft-monophyle
		 (format () "trees/softcull/soft-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-distrib-soft
		 (format () "trees/softcull/distrib-soft-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-syst-monophyle
		 (format () "trees/systematist/syst-soft-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-full-monophyle
		 (format () "trees/monophyly/mono-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-rand-taxon
		 (format () "trees/random/random-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-distrib-rand
		 (format () "trees/random/distrib-rand-tree-~a.txt" *run-number*))
    (record-tree *root* #'snode-match-rand
		 (format () "trees/random/match-rand-tree-~a.txt" *run-number*))))
    

;;; I will need to keep track of the number of taxa (via *next-taxon* id).
(defvar *next-paraphyle* 1)
(defvar *next-syst-paraphyle* 1)
(defvar *next-syst-monophyle* 0)
(defvar *next-hard-monophyle* 0)
(defvar *next-soft-monophyle* 0)
(defvar *next-full-monophyle* 0)
(defvar *next-random-taxon* 0)
(defvar *next-match-random-taxon* 0)
(defvar *next-nested-monophyle* 0)

;-------------------------------------Paraphyly--------------------------------

;;; nodes are selected with a *paraphyly-prob* chance as ancestors.
;;; But a taxon is only accepted with probability = (n + k - 1) / (n + k)
;;; where n is the number of lineages (nodes) in the smaller of the new 
;;; taxon and the old taxon with the new one removed.

(defvar *nodes-remaining* 0)
(defvar *taxa-remaining* 0)
(defvar *trace-probability* 0)

(defvar *para-taxa-sizes* ()) ; Data structure for keeping track of the
			      ; sizes of the paraphyletic taxa.

(defun paraphyly ()
  (setf *next-paraphyle* 1)
  (setf *para-taxa-sizes* (list (1+ *next-leaf-id*)))
  (let ((max-taxa (round (* *paraphyly-prob* (1+ *next-leaf-id*))
			  PROBSCALE)))
    (loop ; keep classifying until I've got enough taxa
     ;(format t "paraphyling...~%")
     (setf *nodes-remaining* (- (1+ *next-leaf-id*) *next-paraphyle*))
     (setf *taxa-remaining* (- max-taxa *next-paraphyle*))
     (setf *trace-probability* (round (* PROBSCALE *taxa-remaining*)
				      *nodes-remaining*))
     (trace-paraphyly *root*)
     (when (>= *next-paraphyle* max-taxa) (return ()))))
  ; save taxon sizes for human examination into paraph file.
  (record-phyles  *next-paraphyle* *para-taxa-sizes* "paraphyles.dat"))

;;; Runs through the tree and assigns paraphyletic taxa.
(defun trace-paraphyly (node)
  (unless (or (null node) (= 0 *taxa-remaining*))
    (when (and (null (snode-left node))
	       (null (snode-right node)))
      (when (< (random PROBSCALE) *trace-probability*)
	(let ((ancestor (find-ancestor node (snode-leaf-id node))))
	  (when (consider-taxon ancestor)
	    (label-paraphyle ancestor *next-paraphyle* (snode-paraphyle ancestor))
	    (decf *taxa-remaining*)
	    (incf *next-paraphyle*))))
	(decf *nodes-remaining*)
	(setf *trace-probability*
	      (if (= 0 *nodes-remaining*) 0
		(round (* PROBSCALE *taxa-remaining*)
		       *nodes-remaining*))))
      (trace-paraphyly (snode-left node))
      (trace-paraphyly (snode-right node))))

(defun consider-taxon (node)
  (if (or (null (snode-parent node)) ; not the root.
	  (not (= (snode-paraphyle node) ; already a taxon.
		  (snode-paraphyle (snode-parent node)))))
      nil ; don't accept this one.
    (let* ((new-n (count-leaves node (snode-paraphyle node)))
	   (old-n (- (get-taxon-size (snode-paraphyle node)
				     *next-paraphyle*
				     *para-taxa-sizes*)
		     new-n)))
      (if (< (random PROBSCALE) (/ (* PROBSCALE
				      (+ (min new-n old-n)
					 *k* -1))
				   (+ (min new-n old-n)
				      *k*)))
	  (progn
	    (set-taxon-size (snode-paraphyle node) old-n
			    *next-paraphyle* *para-taxa-sizes*)
	    (add-taxon-size new-n *para-taxa-sizes*)
	    t)
	nil))))
  

(defun find-ancestor (node leaf-id)
  (let ((parent (snode-parent node)))
    (if (or (null parent)
	    (not (= (snode-leaf-id parent) leaf-id)))
	node
      (find-ancestor parent leaf-id))))


(defun count-leaves (node id)
  (if (null node) 0
    (if (and (null (snode-left node))
	     (null (snode-right node))
	     (= id (snode-paraphyle node)))
	1
      (+ (count-leaves (snode-left node) id)
	 (count-leaves (snode-right node) id)))))

(defun count-all-leaves (node)
  (if (null node) 0
    (if (and (null (snode-left node))
	     (null (snode-right node)))
	1
      (+ (count-all-leaves (snode-left node))
	 (count-all-leaves (snode-right node))))))

(defun count-nodes (node)
  (if (null node) 0
      (+  1 
	  (count-nodes (snode-left node))
	  (count-nodes (snode-right node)))))


(defun label-paraphyle (node label id)
  (unless (or (null node) (not (= id (snode-paraphyle node))))
    (setf (snode-paraphyle node) label)
    (label-paraphyle (snode-left node) label id)
    (label-paraphyle (snode-right node) label id)))



(defun get-taxon-size (id max-id taxa)
  (if (< id max-id) ;legal
      (nth (- max-id (1+ id)) taxa)
    nil))


;;; Replace the size of "id" taxon with a new "size" in the taxa sizes list.
(defun set-taxon-size (id size max-id taxa)
  (nthreplace (- max-id (1+ id)) 
	      taxa
	      size))


;;; add-taxon-size is a macro.
;-----------------------------No Basal Group Paraphyly-------------------------

(defun no-basal ()
  (setf *para-forbidden* (make-array (1+ *next-paraphyle*)
				     :initial-element nil))
  (setf (svref *para-forbidden* 0) t) ;exclude the basal group
  (record-match-phyles *next-paraphyle* (1- *next-paraphyle*) *para-taxa-sizes*
		 *para-forbidden* "no-basalphyles.dat"))


;-------------------------------------Hard Cull--------------------------------


(defvar *hard-taxa-indecies* ())

(defun hard-cull (node)
  (setf *paraphyle-nodes* (make-array (round *max-diversity* 2)
				      :initial-element ()))
  (setf *next-para-node* 0)
  (setf *hard-taxa-indecies* ())
  (setf *next-hard-monophyle* 0)
  (do-hard-cull node)
  (setf *next-hard-monophyle* (length *hard-taxa-indecies*))
  ; output into hard cull file.
  (let ((hard-taxa (select-taxa *hard-taxa-indecies*
			      *para-taxa-sizes*
			      *next-paraphyle*)))
    (record-phyles  *next-hard-monophyle* hard-taxa
		 "hardphyles.dat"))
)

;;; Trace through the paraphyly ids.  Whenever I find a taxon defining node,
;;; check if it is a clade and not monotypic.  If not, do nothing, leaving the 
;;; hard-monophyle field
;;; negative.  If it is a clade, copy the paraphyle id into the hard-monophyle field.

(defun do-hard-cull (node)
  (unless (or (null node) (a-leaf? node))
    (when (taxon-node? node)
      (if (clade? node (snode-paraphyle node))
	  (progn
	    (push (snode-paraphyle node) *hard-taxa-indecies*)
	    (assign-hard-clade node))
	(progn
	  (setf (svref *paraphyle-nodes* *next-para-node*) node)
	  (incf *next-para-node*))))
    (do-hard-cull (snode-left node))
    (do-hard-cull (snode-right node))))



(defun clade? (node id)
  (if (null node) t
    (and (= id (snode-paraphyle node))
	 (clade? (snode-left node) id)
	 (clade? (snode-right node) id))))


(defun assign-hard-clade (node)
  (unless (null node)
    (progn
      (setf (snode-hard-monophyle node) (snode-paraphyle node))
      (assign-hard-clade (snode-left node))
      (assign-hard-clade (snode-right node)))))


(defun select-the-taxa (indecies a-list subset max-id)
  (if (null indecies) subset
    (select-the-taxa (rest indecies) a-list 
		     (cons (get-a-taxon (first indecies) a-list max-id)
			   subset)
		     max-id)))



;-------------------------------------Soft Cull--------------------------------
;;; I have to worry about overwriting work I've already done due to the "look ahead"
;;; of find-a-monophyle.

(defvar *soft-taxa-sizes* ()) ; Data structure for keeping track of the
			      ; sizes of the soft cull taxa.
(defvar *soft-indicies* ())

(defun soft-cull (node)
  (setf *soft-taxa-sizes* ())
  (setf *soft-indicies* ())
  (setf *next-soft-monophyle* 0)
  (do-soft-cull node)
  ; output into soft cull file.
  (record-phyles *next-soft-monophyle* *soft-taxa-sizes*
		 "softphyles.dat")
)

(defun do-soft-cull (node)
  (unless (or (null node) (a-leaf? node))
    (let ((left (snode-left node))
	  (right (snode-right node)))
      (when (taxon-node? node)
	(if (clade? node (snode-paraphyle node))
	    (progn
	      (push (get-a-taxon (snode-paraphyle node)
				 *para-taxa-sizes*
				 *next-paraphyle*)
		    *soft-taxa-sizes*)
	      (assign-soft-clade node (incf *next-soft-monophyle*))
	      (push *next-soft-monophyle* *soft-indicies*))
	  (progn
	    (find-a-monophyle left (snode-paraphyle node))
	    (find-a-monophyle right (snode-paraphyle node)))))
      (do-soft-cull left)
      (do-soft-cull right))))

(defun find-a-monophyle (node id)
  (unless (or (null node) (a-leaf? node))
    (if (and (clade? node id)
	     (potential-taxon-node? node))
	(count-and-assign-soft-clade node)
      (progn
	(find-a-monophyle (snode-left node) id)
	(find-a-monophyle (snode-right node) id)))))

(defun count-and-assign-soft-clade (node)
  (push (assign-soft-clade node (incf *next-soft-monophyle*))
	*soft-taxa-sizes*)
  (push *next-soft-monophyle* *soft-indicies*))


(defun assign-soft-clade (node id)
  (if (null node) 0
    (progn
      (setf (snode-soft-monophyle node) id)
      (+ (if (a-leaf? node) 1 0)
	 (assign-soft-clade (snode-left node) id)
	 (assign-soft-clade (snode-right node) id)))))

;--------------------------------Systematist Perspective--------------------------
;                                 (Smith and Patterson)
(defvar *syst-paraphyly-prob* 0)
(defvar *syst-para-taxa-sizes* ())
(defvar *syst-indicies* ())

(defun syst-soft-cull (node)
  (setf *syst-para-taxa-sizes* ())
  (setf *syst-indicies* ())
  (setf *next-syst-monophyle* 0)
  (do-syst-soft-cull node)
  ; output into soft cull file.
  (record-phyles *next-syst-monophyle* *syst-para-taxa-sizes*
		 "syst-softphyles.dat")
)

(defun do-syst-soft-cull (node)
  (unless (or (null node) (a-leaf? node))
    (let ((left (snode-left node))
	  (right (snode-right node)))
      (when (taxon-node? node)
	(if (clade? node (snode-paraphyle node))
	    (progn
	      (push (get-a-taxon (snode-paraphyle node)
				 *para-taxa-sizes*
				 *next-paraphyle*)
		    *syst-para-taxa-sizes*)
	      (assign-syst-clade node (incf *next-syst-monophyle*))
	      (push *next-syst-monophyle* *syst-indicies*))
	  (progn
	    (find-syst-monophyle left (snode-paraphyle node))
	    (find-syst-monophyle right (snode-paraphyle node)))))
      (do-syst-soft-cull left)
      (do-syst-soft-cull right))))

(defun find-syst-monophyle (node id)
  (unless (or (null node) (a-leaf? node))
    (if (clade? node id)
	(count-and-assign-syst-clade node)
      (progn
	(find-syst-monophyle (snode-left node) id)
	(find-syst-monophyle (snode-right node) id)))))

(defun count-and-assign-syst-clade (node)
  (push (assign-syst-clade node (incf *next-syst-monophyle*))
	*syst-para-taxa-sizes*)
  (push *next-syst-monophyle* *syst-indicies*))


(defun assign-syst-clade (node id)
  (if (null node) 0
    (progn
      (setf (snode-syst-monophyle node) id)
      (+ (if (a-leaf? node) 1 0)
	 (assign-syst-clade (snode-left node) id)
	 (assign-syst-clade (snode-right node) id)))))





;;; I scale the probability of defining a taxon to the number of total nodes in
;;; the tree, rather than by the number of leaves.

;(defun syst-paraphyly ()
;  (setf *syst-paraphyly-prob* (round (/ (* *paraphyly-prob* *next-leaf-id*)
;                                        *next-id*)))
;  (setf *next-syst-paraphyle* 1)
;  (setf *syst-para-taxa-sizes* (list (1+ *next-leaf-id*)))
;  (let ((max-taxa (round (* *syst-paraphyly-prob* *next-id*)
;                          PROBSCALE)))
;    (loop ; keep classifying until I've got enough taxa
;     (setf *nodes-remaining* (- *next-id* *next-syst-paraphyle*))
;     (setf *taxa-remaining* (- max-taxa *next-syst-paraphyle*))
;     (setf *trace-probability* (round (* PROBSCALE *taxa-remaining*)
;                                      *nodes-remaining*))
;       (syst-trace-paraphyly *root*)
;       (when (>= *next-syst-paraphyle* max-taxa) (return ()))))
;  ; save taxon sizes for human examination into paraph file.
;  (record-phyles *syst-para-taxa-sizes* "syst-paraphyles.dat"))
;
;
;
;;;; Using the Systematist's definition of speciation.
;(defun syst-trace-paraphyly (node)
;  (unless (or (null node) (= 0 *taxa-remaining*))
;    (when (and (< (random PROBSCALE) *trace-probability*)
;               (syst-consider-taxon node))
;          (label-syst-paraphyle node *next-syst-paraphyle* (snode-syst-paraphyle node))
;          (decf *taxa-remaining*)
;          (incf *next-syst-paraphyle*))
;    (decf *nodes-remaining*)
;    (setf *trace-probability*
;          (if (= 0 *nodes-remaining*) 0
;            (round (* PROBSCALE *taxa-remaining*)
;                   *nodes-remaining*)))
;    (syst-trace-paraphyly (snode-left node))
;    (syst-trace-paraphyly (snode-right node))))
;
;(defun syst-consider-taxon (node) 
;  (if (or (null (snode-parent node)) ; not the root.
;          (not (= (snode-syst-paraphyle node) ; already a taxon.
;                  (snode-syst-paraphyle (snode-parent node)))))
;      nil ; don't accept this one.
;    (let* ((new-n (syst-count-leaves node (snode-syst-paraphyle node)))
;           (old-n (- (get-taxon-size (snode-syst-paraphyle node)
;                                     *next-syst-paraphyle*
;                                     *syst-para-taxa-sizes*)
;                     new-n)))
;      (if (< (random PROBSCALE) (/ (* PROBSCALE
;                                      (+ (min new-n old-n)
;                                         *k* -1))
;                                   (+ (min new-n old-n)
;                                      *k*)))
;          (progn
;            (set-taxon-size (snode-syst-paraphyle node) old-n
;                            *next-syst-paraphyle* *syst-para-taxa-sizes*)
;            (add-taxon-size new-n *syst-para-taxa-sizes*)
;            t)
;        nil))))
;
;(defun syst-count-leaves (node id)
;  (if (null node) 0
;    (if (and (null (snode-left node))
;             (null (snode-right node))
;             (= id (snode-syst-paraphyle node)))
;        1
;      (+ (count-leaves (snode-left node) id)
;         (count-leaves (snode-right node) id)))))
;
;(defun label-syst-paraphyle (node label id)
;  (unless (or (null node) (not (= id (snode-syst-paraphyle node))))
;    (setf (snode-syst-paraphyle node) label)
;    (label-syst-paraphyle (snode-left node) label id)
;    (label-syst-paraphyle (snode-right node) label id)))

;-------------------------------------Full Monophyly--------------------------------

;;; This is the same as the systematic soft cull, except that it allows monotypic taxa.
;;; It is derived directly from the syst. soft cull by just labelling the unlablled
;;; leaves as new taxa.

(defvar *full-mono-sizes* ())

(defun full-monophyly (node)
  (setf *next-full-monophyle* *next-syst-monophyle*)
  (setf *full-mono-sizes* *syst-para-taxa-sizes*)
  (label-neglected-leaves node)
  (record-phyles *next-full-monophyle* *full-mono-sizes* "monophyles.dat")
  )

(defun label-neglected-leaves (node)
  (unless (null node)
    (if (> (snode-syst-monophyle node) -1) ; if part of a monophyle...
	(copy-syst-label node)
      (if (and (null (snode-left node)) ; if leaf...
	       (null (snode-right node)))
	  (progn
	    (push 1 *full-mono-sizes*)
	    (setf (snode-full-monophyle node) (incf *next-full-monophyle*)))
	(progn ; if internal node...
	  (label-neglected-leaves (snode-left node))
	  (label-neglected-leaves (snode-right node)))))))

(defun copy-syst-label (node)
  (unless (null node) 
    (setf (snode-full-monophyle node) (snode-syst-monophyle node))
    (copy-syst-label (snode-left node))
    (copy-syst-label (snode-right node))))



;-------------------------------------Nested Monophyly--------------------------------

(defun nested-monophyly (node-list)

)

;;; Counting the diversity through time: 

;;; Keep a list of who is alive ordered by extinction times as well as an array for
;;; each classification system with a number for each taxon indicating how many
;;; lineages are extant for that taxon.  Finally a number for each system indicating
;;; the number of extant taxa on that system (which is output to the file)


;-------------------------------------Random Taxa--------------------------------

;;; The idea here is to test if it is number of taxa that is really doing the trick.
;;; So just randomly assign nodes of the tree to n taxa where  next-full-monophyly <
;;; n < number of classified nodes in full monphyly.  

;;; This should result in the ugliest taxonomy ever invented by humans.  It is weaker
;;; on sampling than the other taxonomies because there is no ancester relationship
;;; in the taxa, so you can't infer anything from a sample of an old member then a
;;; gap then a sample of a newer member.

;;; I should only classify the same number of nodes that are classified in full monophyly
;;; to avoid confounding with the coverage statistic.

;;; Algorithm:
;;; Randomly assign M nodes into N groups.  This will leave some groups empty.  So I
;;; must fix them up by shifting a node from a group with more than one to an empty group.
;;; I'll do this by just searching down from the root and whenever I encounter an
;;; individual from a 2+ group I'll reassign it to an empty group.

;;; I now give complete coverage, rather than matching the full-monophyle coverage as
;;; before.

(defvar *rand-taxa-sizes* ()) ; Data structure for keeping track of the
			      ; sizes of the soft cull taxa.
(defvar *rand-taxa-prob* 0)
(defvar *ancestors-remaining* 0)
(defvar *labelleds-remaining* 0)
(defvar *next-empty-taxon* nil)

(defun random-taxa (node)
  (setf *next-random-taxon* (round (+ *max-diversity* 1 *next-full-monophyle*)
				   2))
  (setf *rand-taxa-sizes* (make-array (1+ *next-random-taxon*)
				      :initial-element 0))
  ; note that we want to be working with #tips - 1 = *max-diversity*
  (setf *labelleds-remaining* *max-diversity*) ;(count-assigned-ancestors node))
  (setf *ancestors-remaining* *max-diversity*)
  (setf *rand-taxa-prob* (round (* PROBSCALE *labelleds-remaining*)
				*ancestors-remaining*))
				
  (randomly-assign-taxa node)
  ; output into rand taxa file.
  (record-phyles *next-random-taxon*
		 (rest (array-to-list *rand-taxa-sizes* (1+ *next-random-taxon*)))
		 "randphyles.dat")
)

(defun count-assigned-ancestors (node)
  (if (or (null node) (null (snode-left node))
	      (null (snode-right node)))
      0
    (+ (if (< (snode-full-monophyle node) 0) 0 1)
       (count-assigned-ancestors (snode-left node))
       (count-assigned-ancestors (snode-right node)))))
       

    
(defun randomly-assign-taxa (node)
  (do-random-assignments node)
  (setf *next-empty-taxon* (find-next-empty *rand-taxa-sizes*
					    1 *next-random-taxon*))
  (fill-empty-taxa node))

(defun do-random-assignments (node)
  (unless (null node)
    (if (and (null (snode-left node)) (null (snode-right node))) ; a leaf
	(let ((taxon (1+ (random *next-random-taxon*))))
	  (incf (svref *rand-taxa-sizes* taxon))
	  (setf (snode-rand-taxon node) taxon))
      (progn ; ancestor node...
	(when (< (random PROBSCALE) *rand-taxa-prob*) ; label it...
	  (let ((taxon (1+ (random *next-random-taxon*))))
	    ;(incf (svref *rand-taxa-sizes* taxon))
	    (setf (snode-rand-taxon node) taxon)
	    (decf *labelleds-remaining*)))
	(decf *ancestors-remaining*)
	(setf *rand-taxa-prob* (if (<= *ancestors-remaining* 0) 0
				 (round (* PROBSCALE *labelleds-remaining*)
					*ancestors-remaining*)))
	(do-random-assignments (snode-left node))
	(do-random-assignments (snode-right node))))))


(defun fill-empty-taxa (node)
  (unless (or (null node) (null *next-empty-taxon*))
    (if (and (null (snode-left node)) ; a leaf
	     (null (snode-right node))
	     (>= (snode-rand-taxon node) 0)
	     (> (svref *rand-taxa-sizes* (snode-rand-taxon node)) 1)) ;then reassign...
	(progn
	  (decf (svref *rand-taxa-sizes* (snode-rand-taxon node)))
	  (setf (svref *rand-taxa-sizes* *next-empty-taxon*) 1)
	  (setf *next-empty-taxon*
		(find-next-empty *rand-taxa-sizes* (1+ *next-empty-taxon*)
				 *next-random-taxon*))
	  (fill-empty-taxa (snode-left node))
	  (fill-empty-taxa (snode-right node)))
      (progn 
	(fill-empty-taxa (snode-left node))
	(fill-empty-taxa (snode-right node))))))
 	  



(defun find-next-empty (arr start end)
  (if (> start end) nil
    (if (= 0 (svref arr start)) start
      (find-next-empty arr (1+ start) end))))



;-------------------------------------Matched Soft Cull--------------------------------

;;; randomly select *next-paraphyle* - *next-soft-monophyle* numbers from the range
;;; 1 to *next-soft-monophyle* and make these the excluded group.

;;; Then rewrite the macro to exclude these groups rather than everything above
;;; *next-paraphyle*.

(defvar *indices* ())
(defvar *soft-forbidden* ()) ; array of flags indicating which soft cull clades are
					;excluded.


(defun match-soft ()
  (setf *indices* (make-array *next-soft-monophyle*
			      :initial-element 0))
  (dotimes (i *next-soft-monophyle*) ; put the indicies into the array.
    (setf (svref *indices* i) (1+ i)))
  (setf *soft-forbidden* (make-array (1+ *next-soft-monophyle*)
				     :initial-element nil))
  (dotimes (i (- *next-soft-monophyle* *next-paraphyle*))
    (let ((num (random (- *next-soft-monophyle* i))))
      (setf (svref *soft-forbidden*
		   (svref *indices* num)) t)
      (swap num (- *next-soft-monophyle* i 1) *indices*)))
  (record-match-phyles *next-paraphyle* *next-soft-monophyle* *soft-taxa-sizes*
		 *soft-indicies* *soft-forbidden* "match-softphyles.dat"))

;-------------------------------------Matched Syst Soft Cull---------------------------

;;; Same as Matched Soft Cull but select from *next-syst-monophyle*.
(defvar *syst-soft-forbidden* ()) ; array of flags indicating which soft cull clades are
					;excluded.



(defun match-syst-soft ()
  (setf *indices* (make-array *next-syst-monophyle*
			      :initial-element 0))
  (dotimes (i *next-syst-monophyle*)
    (setf (svref *indices* i) (1+ i)))
  (setf *syst-soft-forbidden* (make-array (1+ *next-syst-monophyle*)
				     :initial-element nil))
  (dotimes (i (- *next-syst-monophyle* *next-paraphyle*))
    (let ((num (random (- *next-syst-monophyle* i))))
      (setf (svref *syst-soft-forbidden*
		   (svref *indices* num)) t)
      (swap num (- *next-syst-monophyle* i 1) *indices*)))
  (record-match-phyles *next-paraphyle* *next-syst-monophyle* *syst-para-taxa-sizes*
		 *syst-indicies* *syst-soft-forbidden* "match-systphyles.dat"))

;-------------------------------------Matched Random-----------------------------------

;;; Just redo the random classification with *next-paraphyle* number of taxa.
;;; How much coverage?


;;;*** Still must do this***


(defvar *match-rand-taxa-sizes* ()) ; Data structure for keeping track of the
			      ; sizes of the random taxa.

(defun match-random (node)
  (setf *next-match-random-taxon* *next-paraphyle*)
  (setf *match-rand-taxa-sizes* (make-array (1+ *next-match-random-taxon*)
				      :initial-element 0))
  ; note that we want to be working with #tips - 1 = *max-diversity*
  (setf *labelleds-remaining* *max-diversity*) ;(count-assigned-ancestors node))
  (setf *ancestors-remaining* *max-diversity*)
  (setf *rand-taxa-prob* (round (* PROBSCALE *labelleds-remaining*)
				*ancestors-remaining*))
				
  (match-randomly-assign-taxa node)
  ; output into rand taxa file.
  (record-phyles *next-match-random-taxon*
		 (rest (array-to-list *match-rand-taxa-sizes*
				      (1+ *next-match-random-taxon*)))
		 "matchrandphyles.dat")
)

     

    
(defun match-randomly-assign-taxa (node)
  (do-match-random-assignments node)
  (setf *next-empty-taxon* (find-next-empty *match-rand-taxa-sizes*
					    1 *next-match-random-taxon*))
  (match-fill-empty-taxa node))

(defun do-match-random-assignments (node)
  (unless (null node)
    (if (and (null (snode-left node)) (null (snode-right node))) ; a leaf
	(let ((taxon (1+ (random *next-match-random-taxon*))))
	  (incf (svref *match-rand-taxa-sizes* taxon))
	  (setf (snode-match-rand node) taxon))
      (progn ; ancestor node...
	(when (< (random PROBSCALE) *rand-taxa-prob*) ; label it...
	  (let ((taxon (1+ (random *next-match-random-taxon*))))
	    (incf (svref *match-rand-taxa-sizes* taxon))
	    (setf (snode-match-rand node) taxon)
	    (decf *labelleds-remaining*)))
	(decf *ancestors-remaining*)
	(setf *rand-taxa-prob* (if (<= *ancestors-remaining* 0) 0
				 (round (* PROBSCALE *labelleds-remaining*)
					*ancestors-remaining*)))
	(do-match-random-assignments (snode-left node))
	(do-match-random-assignments (snode-right node))))))


(defun match-fill-empty-taxa (node)
  (unless (or (null node) (null *next-empty-taxon*))
    (if (and (>= (snode-match-rand node) 0)
	     (> (svref *match-rand-taxa-sizes* (snode-match-rand node)) 1)) ;then reassign...
	(progn
	  (decf (svref *match-rand-taxa-sizes* (snode-match-rand node)))
	  (setf (svref *match-rand-taxa-sizes* *next-empty-taxon*) 1)
	  (setf *next-empty-taxon*
		(find-next-empty *match-rand-taxa-sizes* (1+ *next-empty-taxon*)
				 *next-match-random-taxon*))
	  (match-fill-empty-taxa (snode-left node))
	  (match-fill-empty-taxa (snode-right node)))
      (progn 
	(match-fill-empty-taxa (snode-left node))
	(match-fill-empty-taxa (snode-right node))))))
 	  




;-------------------------------Matched Distribution Random------------------------------

(defvar *rand-taxon-counters* ())
(defvar *rand-parental-counters* ())

(defun match-distrib-rand (node)
  (setf *rand-taxon-counters* (make-array *next-paraphyle*
					  :initial-contents *para-taxa-sizes*))
  (setf *rand-parental-counters* (make-array *next-paraphyle*
					  :initial-element 0))
  (setup-parental-counters node *rand-parental-counters*)
  (assign-rand-distrib node))

(defun assign-rand-distrib (node)
  (unless (null node)
      (if (a-leaf? node)
	  (setf (snode-distrib-rand node)
		(find-and-fill-empty-slot *rand-taxon-counters* *next-paraphyle*))
	(let ((taxon (find-and-fill-empty-slot *rand-parental-counters* 
					       *next-paraphyle*)))
	  (setf (snode-distrib-rand node) taxon)
	  (assign-rand-distrib (snode-left node))
	  (assign-rand-distrib (snode-right node))))))

(defun setup-parental-counters (node counters)
  (unless (null node)
      (unless (a-leaf? node)
	  (incf (svref counters (snode-paraphyle node)))
	  (setup-parental-counters (snode-left node) counters)
	  (setup-parental-counters (snode-right node) counters))))


;;; This takes an array, randomly selects an index with a non-zero content,
;;; decrements that content and returns the index.

(defun find-and-fill-empty-slot (counters-array num-counters)
  (let ((index (random num-counters))
	(counter 0))
    (while (and (<= (svref counters-array index) 0)
		(< counter num-counters))
      (progn
	(incf counter)
	(setf index (rem (1+ index) num-counters))))
    (if (>= counter num-counters)
	nil
      (progn
	(decf (svref counters-array index))
	index))))
		

;-------------------------------Matched Distribution Soft---------------------------------

(defvar *ancestor-queue* ())
(defvar *next-clade* 0)
(defvar *distrib-soft-sizes* ()) 

(defun match-distrib-soft (node)
  (setf *ancestor-queue* (make-queue))
  (setf *distrib-soft-sizes* ()) 
  (setf *next-clade* 0)
  (enqueue *ancestor-queue* node)
  (note-leaf-numbers *root*)
  (assign-match-distrib-soft (sort *para-taxa-sizes* #'>))
  (record-phyles  *next-clade* *distrib-soft-sizes* "distrib-soft.dat"))


(defun note-leaf-numbers (node)
  (unless (null node)
      (if (a-leaf? node)
	  (progn
	    (setf (snode-num-leaves node) 1)
	    1)
	(setf (snode-num-leaves node)
	      (+ (note-leaf-numbers (snode-left node))
		 (note-leaf-numbers (snode-right node)))))))


;;; Each call of this will assign the largest remaining clade of
;;; size <= (first target-sizes).  So I need to:
;;;    1. dequeue
;;;    2. if <= target, label this clade and recurse.
;;;    3. else split it into its two children and either requeue or enqueue them.
;;;    4. repeat.

(defun assign-match-distrib-soft (target-sizes)
  (unless (null target-sizes)
    (let ((target (first target-sizes))
	  (front-guy (dequeue *ancestor-queue*)))
      (loop
       (when (<= (snode-num-leaves front-guy) target)
	 (push (snode-num-leaves front-guy) *distrib-soft-sizes*)
	 (label-clade front-guy (incf *next-clade*))
	 (return ()))
       (let ((left-child (snode-left front-guy))
	     (right-child (snode-right front-guy)))
	 (unless (null left-child)
	   (ordered-requeue *ancestor-queue* left-child #'num-leaf-pred))
	 (unless (null right-child)
	   (ordered-requeue *ancestor-queue* right-child #'num-leaf-pred))
	 (setf front-guy (dequeue *ancestor-queue*)))
       )
      (assign-match-distrib-soft (rest target-sizes)))))

(defun num-leaf-pred (node1 node2)
  (>= (snode-num-leaves node1) (snode-num-leaves node2)))

(defun label-clade (node label) 
  (unless (null node)
    (setf (snode-distrib-soft node) label)
    (label-clade (snode-left node) label)
    (label-clade (snode-right node) label)))