Added Sackin index

cmd/caterpillartree.go

@@ -25,7 +25,7 @@ func caterpilarTree(nbtrees int, nbtips int, output string, seed int64, rooted b
 	}
 
 	for i := 0; i < nbtrees; i++ {
-		t, err = tree.RandomCaterpilarBinaryTree(nbtips, rooted)
+		t, err = tree.RandomCaterpillarBinaryTree(nbtips, rooted)
 		if err != nil {
 			return err
 		}

cmd/stats.go

@@ -22,7 +22,7 @@ For example:
 	Run: func(cmd *cobra.Command, args []string) {
 		/* Dividing trees */
 		f := openWriteFile(outtreefile)
-		f.WriteString("tree\tnodes\ttips\tedges\tmeanbrlen\tsumbrlen\tmeansupport\tmediansupport\trooted\tnbcherries\tcolless\n")
+		f.WriteString("tree\tnodes\ttips\tedges\tmeanbrlen\tsumbrlen\tmeansupport\tmediansupport\trooted\tnbcherries\tcolless\tsackin\n")
 		treefile, treechan := readTrees(intreefile)
 		defer treefile.Close()
 		for t := range treechan {
@@ -39,12 +39,13 @@ For example:
 			f.WriteString(fmt.Sprintf("\t%.8f", t.Tree.MeanSupport()))
 			f.WriteString(fmt.Sprintf("\t%.8f", t.Tree.MedianSupport()))
 			if t.Tree.Rooted() {
-				f.WriteString(fmt.Sprintf("\trooted\n"))
+				f.WriteString(fmt.Sprintf("\trooted"))
 			} else {
 				f.WriteString(fmt.Sprintf("\tunrooted"))
 			}
 			f.WriteString(fmt.Sprintf("\t%d", t.Tree.NbCherries()))
-			f.WriteString(fmt.Sprintf("\t%d\n", t.Tree.CollessIndex()))
+			f.WriteString(fmt.Sprintf("\t%d", t.Tree.CollessIndex()))
+			f.WriteString(fmt.Sprintf("\t%d\n", t.Tree.SackinIndex()))
 		}
 		f.Close()
 	},

test.sh

@@ -438,8 +438,8 @@ rm -f expected result clade
 
 echo "->gotree stats"
 cat > expected <<EOF
-tree	nodes	tips	edges	meanbrlen	sumbrlen	meansupport	mediansupport	rooted	nbcherries	colless
-0	18	10	17	0.10486138	1.78264354	NaN	NaN	unrooted	3	7
+tree	nodes	tips	edges	meanbrlen	sumbrlen	meansupport	mediansupport	rooted	nbcherries	colless	sackin
+0	18	10	17	0.10486138	1.78264354	NaN	NaN	unrooted	3	7	35
 EOF
 gotree generate yuletree -s 10 | gotree stats > result
 diff expected result

tests/tree_test.go

@@ -231,3 +231,72 @@ func TestColless4(t *testing.T) {
 		t.Error(fmt.Sprintf("%d cherries are found instead of %d", cidx, expected))
 	}
 }
+
+// Test counting the Sackin Index on balanced rooted and unrooted trees
+func TestSackin1(t *testing.T) {
+	// rooted, depth 7: 128 tips : Sackin should be 7*128=896
+	tr, err := tree.RandomBalancedBinaryTree(7, true)
+	if err != nil {
+		t.Error(err)
+	}
+	expected := 7 * 128
+	sidx := tr.SackinIndex()
+	if sidx != expected {
+		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
+	}
+
+	// unrooted, depth 7: 128 tips : Sackin should be 7*128=896
+	tr, err = tree.RandomBalancedBinaryTree(7, false)
+	if err != nil {
+		t.Error(err)
+	}
+	expected = 7 * 128
+	sidx = tr.SackinIndex()
+	if sidx != expected {
+		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
+	}
+
+	// rooted caterpillar, 128 tips : Sackin should be sum(1:127)+127=8255
+	tr, err = tree.RandomCaterpillarBinaryTree(128, true)
+	if err != nil {
+		t.Error(err)
+	}
+	expected = 8255
+	sidx = tr.SackinIndex()
+	if sidx != expected {
+		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
+	}
+
+	// unrooted caterpillar, 128 tips : Sackin should be (sum(2:64)+64)*2=4286
+	tr, err = tree.RandomCaterpillarBinaryTree(128, false)
+	if err != nil {
+		t.Error(err)
+	}
+	expected = 4286
+	sidx = tr.SackinIndex()
+	if sidx != expected {
+		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
+	}
+
+	treeString := "(((1,2),(3,4)),5);"
+	expected = 13
+	tr, err = newick.NewParser(strings.NewReader(treeString)).Parse()
+	if err != nil {
+		t.Error(err)
+	}
+	sidx = tr.SackinIndex()
+	if sidx != expected {
+		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
+	}
+
+	treeString = "((1,2),(3,4),5);"
+	expected = 12
+	tr, err = newick.NewParser(strings.NewReader(treeString)).Parse()
+	if err != nil {
+		t.Error(err)
+	}
+	sidx = tr.SackinIndex()
+	if sidx != expected {
+		t.Error(fmt.Sprintf("Sackin index found : %d instead of %d expected", sidx, expected))
+	}
+}

tree/stats.go

@@ -149,3 +149,35 @@ func collessIndexRecur(n *Node, prev *Node) (colless, tips int) {
 	colless += (maxtips - mintips)
 	return
 }
+
+// This functions computes the Sackin index of a rooted tree
+// As the sum of all tip depths.
+// If the tree is unrooted, then it takes as starting point the deepest
+// edge of the tree.
+// No problems with multifurcations.
+func (t *Tree) SackinIndex() (sackin int) {
+	sackin = 0
+	if !t.Rooted() {
+		edge := t.DeepestEdge()
+		leftSackin := sackinIndexRecur(edge.Left(), edge.Right(), 1)
+		rightSackin := sackinIndexRecur(edge.Right(), edge.Left(), 1)
+		sackin = leftSackin + rightSackin
+	} else {
+		sackin = sackinIndexRecur(t.Root(), nil, 0)
+	}
+	return
+}
+
+// Returns the sum of the depths of all tips down the current node
+func sackinIndexRecur(n *Node, prev *Node, depth int) (sackin int) {
+	if n.Tip() {
+		return depth
+	}
+	sackin = 0
+	for _, c := range n.Neigh() {
+		if c != prev {
+			sackin += sackinIndexRecur(c, n, depth+1)
+		}
+	}
+	return
+}

tree/treegen.go

@@ -178,7 +178,7 @@ func RandomYuleBinaryTree(nbtips int, rooted bool) (*Tree, error) {
 //nbtips : Number of tips of the random binary tree to create
 //rooted: if true, generates a rooted tree
 //branch length: follow an exponential distribution with param lambda=1/0.1
-func RandomCaterpilarBinaryTree(nbtips int, rooted bool) (*Tree, error) {
+func RandomCaterpillarBinaryTree(nbtips int, rooted bool) (*Tree, error) {
 	t := NewTree()
 	if nbtips < 2 {
 		return nil, errors.New("Cannot create an unrooted random binary tree with less than 2 tips")