#! /bin/sh
#
#  sort_PAUP_trees  tree_file  tree_lengths
#
#     tree_file     -- PHYLIP user tree file (first line is number of trees)
#     tree_lengths  -- PAUP display buffer after running  "Tree Scores" command
#                     with length of all trees
#
if test $# -ne 2; then
      echo "" >&2
	echo "Usage:  sort_PAUP_trees  tree_file  file_with_scores" >&2
      echo "" >&2
      echo "        sorted trees appear on standard output" >&2
      echo "" >&2
	exit
elif test ! -f "$1"; then
	echo "Unable to find file of trees: $1" >&2
	exit
elif test ! -f "$2"; then
	echo "Unable to find file of scores: $2" >&2
	exit
fi

#  Read list of tree scores from PAUP log file and print on one line:

gawk '
BEGIN {state = 0; ntree = 0; prefix = "";}

NF == 0 {next;}

/^Lengths of trees in memory/ {state = 1; next;}

state == 1 {
	if ($0 !~ /^Tree #/) exit;
	nt1 = $3;
	nt2 = $NF;
	if (nt1 != ntree + 1) {
		print "Tree numbering problem.  Expected " ntree+1 ", but found " nt1> "/dev/stderr";
		print $0 > "/dev/stderr";
		exit;
		}
	state = 2;
	next;
	}

state == 2 {
	if ($1 != "Length")  {
		print "Expected tree lengths.  Found:" > "/dev/stderr";
		print $0 > "/dev/stderr";
		exit;
		}

	for (i = 2; i <= NF; i++) {printf("%s%d", prefix, $i); prefix = "\t"; ntree++;}

	if (ntree != nt2) {
		print "Tree numbering problem.  Expected " nt2 " scores, but found " ntree > "/dev/stderr";
		print $0 > "/dev/stderr";
		exit;
		}
	state = 1;
	next;
	}

END {printf("\n");}' "$2" | \

#  print tree length in front of tree, sort and split

gawk '
NF == 0 {next; }

NR == 1 {for (i = 1; i <= NF; i++) len[i] = $i; printf("0\t%d\n", NF); next;}

{ printf("%d\t%s\n", len[++n], $0); next;
	}' - "$1" | sort +0n | sed 's/^[^	]*	//'

#  do the protein maximum likelihood in the order of most parsimonious to least

echo ""  >&2
echo "To use these topologies, you might use a command like:"  >&2
echo ""  >&2
echo "   protml -I -u  file.phy  file.tpl > file.protml"  >&2
echo ""  >&2
echo "   where: file.phy is a PHYLIP format protein sequence file"  >&2
echo "          file.tpl is the file containing the trees from this script"  >&2
echo "          file.protml is the file to save the protml results into"  >&2
echo ""  >&2
echo "   Note that '-I' is a capital '-i', NOT a lower case '-L'"  >&2
echo ""  >&2

# cat anal_2.BLOSUM45.protml | grep 'ln L' | sed 's/^.*ln L:  */	/' | grep -n '.' | sort +1nr