Refactor TipIndex

cmd/booster.go

@@ -3,6 +3,7 @@ package cmd
 import (
 	"fmt"
 	goio "io"
+	"os"
 	"time"
 
 	"github.com/evolbioinfo/gotree/io"
@@ -22,6 +23,7 @@ var boosterCmd = &cobra.Command{
 		var rawtree *tree.Tree
 		var boottreefile goio.Closer
 		var boottreechan <-chan tree.Trees
+		var f *os.File
 
 		writeLogBooster()
 		if refTree, err = readTree(supportIntree); err != nil {

cmd/edgetrees.go

@@ -74,7 +74,7 @@ The resulting trees are star trees to which we added one biparition. All branch
 					if !edgeS.e.Right().Tip() {
 						var edgeOut *os.File
 						var err2 error
-						var bitsetindex uint
+						var bitsetindex int
 
 						if outtreefile == "stdout" || outtreefile == "-" {
 							if edgeOut, err2 = openWriteFile("stdout"); err2 != nil {
@@ -98,7 +98,7 @@ The resulting trees are star trees to which we added one biparition. All branch
 									err = err2
 									return
 								}
-								if edgeS.e.TipPresent(bitsetindex) {
+								if edgeS.e.TipPresent(uint(bitsetindex)) {
 									if leftnb > 0 {
 										leftbuffer.WriteRune(',')
 									}

cmd/splits.go

@@ -47,12 +47,12 @@ Then: One line per branch, and 0/1
 			}
 
 			f.WriteString("Tree\t")
-			names := t.Tree.SortedTips()
-			for i := len(names) - 1; i >= 0; i-- {
-				if i < len(names)-1 {
+			tips := t.Tree.SortedTips()
+			for i := len(tips) - 1; i >= 0; i-- {
+				if i < len(tips)-1 {
 					f.WriteString("|")
 				}
-				f.WriteString(names[i])
+				f.WriteString(tips[i].Name())
 			}
 			f.WriteString("\n")
 			for _, e := range t.Tree.Edges() {

draw/radial.go

@@ -74,7 +74,7 @@ func (layout *radialLayout) constructNode(t *tree.Tree, node *tree.Node, prev *t
 		for num, child := range node.Neigh() {
 			if child != prev {
 
-				numT, _ := node.Edges()[num].NumTipsRight()
+				numT := node.Edges()[num].NumTipsRight()
 				leafCounts = append(leafCounts, numT)
 				sumLeafCount += numT
 				i++

support/booster.go

@@ -4,7 +4,6 @@ import (
 	"fmt"
 	"math"
 	"os"
-	"sync"
 
 	"github.com/evolbioinfo/gotree/io"
 	"github.com/evolbioinfo/gotree/tree"
@@ -14,21 +13,29 @@ import (
 // If "absent" is true, then we know that the ref branch is not present in the bootstrap tree (it is faster to compute then), and we stop if dist == 1
 // Else: we do not know, we do the full postorder traversal, and speciestoadd && speciestoremove are filled
 func MinTransferDist(refedge *tree.Edge, reftree, boottree *tree.Tree, ntips int, bootedges []*tree.Edge, absent bool) (dist int, minedge *tree.Edge, speciestoadd, speciestoremove []*tree.Node) {
-	numbootedges := len(bootedges)
-	ones := make([]int, numbootedges)
+	ones := make([]int, len(bootedges))
 	p, _ := refedge.TopoDepth()
 	dist = p - 1
+	speciestoadd = make([]*tree.Node, 0, 10)
+	speciestoremove = make([]*tree.Node, 0, 10)
+
+	// If ref edge is a terminal edge
+	if p == 1 {
+		tip := refedge.Right()
+		boottip, _ := boottree.TipNode(tip.Name())
+		if boottip != nil {
+			minedge = boottip.Edges()[0]
+		}
+		return
+	}
 
 	stop := false
 	minTransferDistRecur(reftree, ntips, boottree.Root(), nil, nil, refedge, p, ones, &dist, &minedge, absent, &stop)
 
-	speciestoadd = make([]*tree.Node, 0, 10)
-	speciestoremove = make([]*tree.Node, 0, 10)
-
 	if !absent {
 		// computing species to move
 		/////////////////////////////////////////
-		n_subtree, _ := minedge.NumTipsRight()
+		n_subtree := minedge.NumTipsRight()
 		ones_subtree := ones[minedge.Id()]
 		zeros_subtree := n_subtree - ones_subtree
 		ones_total := ntips - p
@@ -62,7 +69,7 @@ func speciesToMoveRecursive(bootedge *tree.Edge, cur, prev *tree.Node, edge *tre
 	}
 
 	if edge != nil {
-		subtreesize, _ := edge.NumTipsRight()
+		subtreesize := edge.NumTipsRight()
 		if (want_ones_now && ones[edge.Id()] == subtreesize) || (!want_ones_now && ones[edge.Id()] == 0) {
 			return
 		}
@@ -81,9 +88,9 @@ func minTransferDistRecur(refTree *tree.Tree, ntips int, cur, prev *tree.Node, c
 	}
 	curOnes := 0
 	if cur.Tip() {
-		tipIndex, _ := refTree.TipIndex(cur.Name())
-		light := refEdge.TipPresent(tipIndex)
-		if r, _ := refEdge.NumTipsRight(); r > ntips/2 {
+		tipIndex := cur.TipIndex()
+		light := refEdge.TipPresent(uint(tipIndex))
+		if r := refEdge.NumTipsRight(); r > ntips/2 {
 			light = !light
 		}
 		if !light {
@@ -104,7 +111,7 @@ func minTransferDistRecur(refTree *tree.Tree, ntips int, cur, prev *tree.Node, c
 
 	if curEdge != nil {
 		ones[curEdge.Id()] = curOnes
-		r, _ := curEdge.NumTipsRight()
+		r := curEdge.NumTipsRight()
 		zero := r - curOnes
 		d := p - zero + curOnes
 		if d > ntips/2 {
@@ -144,14 +151,14 @@ func Booster(reftree *tree.Tree, boottrees <-chan tree.Trees, cpu int,
 	var mindepth int = int(math.Ceil(1.0/distcutoff + 1.0)) // For taxa move computation
 
 	if computeavgtaxa {
-		movedspecies = make([]float64, len(edges)+1)
-		movedspeciestmp = make([]int, len(edges)+1)
+		movedspecies = make([]float64, len(tips))
+		movedspeciestmp = make([]int, len(tips))
 	}
 
 	if computeperbranchtaxa {
 		movedperbranch = make([][]int, len(edges))
 		for i, _ := range edges {
-			movedperbranch[i] = make([]int, len(edges)+1)
+			movedperbranch[i] = make([]int, len(tips))
 		}
 	}
 
@@ -176,7 +183,6 @@ func Booster(reftree *tree.Tree, boottrees <-chan tree.Trees, cpu int,
 			fmt.Fprintf(os.Stderr, "CPU : %02d - Bootstrap tree %d\r", cpu, boot.Id)
 			bootedges := boot.Tree.Edges()
 			bootedgeindex := tree.NewEdgeIndex(uint64(len(bootedges)*2), 0.75)
-
 			for i, e := range bootedges {
 				e.SetId(i)
 				if !e.Right().Tip() {
@@ -187,41 +193,43 @@ func Booster(reftree *tree.Tree, boottrees <-chan tree.Trees, cpu int,
 				}
 				bootedgeindex.PutEdgeValue(e, i, e.Length())
 			}
-
-			var wg sync.WaitGroup
-			wg.Add(len(edges))
-			edgechan := make(chan *tree.Edge, cpu)
-			for _, tmpe := range edges {
-				edgechan <- tmpe
-				go func() {
-					e := <-edgechan
-					if p, _ := e.TopoDepth(); p > 1 {
-						if _, ok := bootedgeindex.Value(e); ok {
-							if p >= mindepth {
-								nbranchclose++
-							}
-							e.IncrementSupport(0.0)
-						} else if p == 2 {
-							e.IncrementSupport(1.0)
-						} else {
-							dist, minedge, sptoadd, sptoremove := MinTransferDist(e, reftree, boot.Tree, len(tips), bootedges, !(computeavgtaxa || computeperbranchtaxa))
-							//dist, edge, sptoadd, sptoremove := MinTransferDist(e, reftree, boot.Tree, len(tips), bootedges)
-							e.IncrementSupport(float64(dist))
+			//var wg sync.WaitGroup
+			//wg.Add(len(edges))
+			//edgechan := make(chan *tree.Edge, cpu)
+			for _, e := range edges {
+				//edgechan <- tmpe
+				//go func() {
+				//e := <-edgechan
+				//e := tmpe
+				if p, _ := e.TopoDepth(); p > 1 {
+					if _, ok := bootedgeindex.Value(e); ok {
+						if p >= mindepth {
+							nbranchclose++
+						}
+						e.IncrementSupport(0.0)
+					} else if p == 2 {
+						e.IncrementSupport(1.0)
+					} else {
+						dist, minedge, sptoadd, sptoremove := MinTransferDist(e, reftree, boot.Tree, len(tips), bootedges, !(computeavgtaxa || computeperbranchtaxa))
+						//dist, edge, sptoadd, sptoremove := MinTransferDist(e, reftree, boot.Tree, len(tips), bootedges)
+						e.IncrementSupport(float64(dist))
+						if computeavgtaxa || computeperbranchtaxa {
 							UpdateTaxaMoveArrays(e, minedge, dist, p,
 								movedspeciestmp, movedperbranch, &nbranchclose,
 								sptoadd, sptoremove, distcutoff, mindepth,
 								computeavgtaxa, computeperbranchtaxa)
 						}
 					}
-					wg.Done()
-				}()
+				}
+				//wg.Done()
+				//}()
 			}
-			wg.Wait()
+			//wg.Wait()
 		}
 		if computeavgtaxa || computeperbranchtaxa {
 			for _, t := range tips {
-				movedspecies[t.Id()] += float64(movedspeciestmp[t.Id()]) / float64(nbranchclose)
-				movedspeciestmp[t.Id()] = 0
+				movedspecies[t.TipIndex()] += float64(movedspeciestmp[t.TipIndex()]) / float64(nbranchclose)
+				movedspeciestmp[t.TipIndex()] = 0
 			}
 		}
 		nboot++
@@ -307,20 +315,20 @@ func UpdateTaxaMoveArrays(ref, boot *tree.Edge, dist, p int,
 
 	if computeavgtaxa && norm <= distcutoff && p >= mindepth {
 		for _, t := range speciestoadd {
-			moved_species[t.Id()]++
+			moved_species[t.TipIndex()]++
 		}
 		for _, t := range speciestoremove {
-			moved_species[t.Id()]++
+			moved_species[t.TipIndex()]++
 		}
 		(*nb_branches_close)++
 	}
 	if computeperbranchtaxa {
 		for _, t := range speciestoadd {
-			moved_species_per_branch[ref.Id()][t.Id()]++
+			moved_species_per_branch[ref.Id()][t.TipIndex()]++
 		}
 
 		for _, t := range speciestoremove {
-			moved_species_per_branch[ref.Id()][t.Id()]++
+			moved_species_per_branch[ref.Id()][t.TipIndex()]++
 		}
 	}
 }

tests/consensus_test.go

@@ -144,7 +144,7 @@ func TestConsensus2(t *testing.T) {
 	}
 
 	go func() {
-		trees <- tree.Trees{randtree1, 1, nil}
+		trees <- tree.Trees{Tree: randtree1, Id: 1, Err: nil}
 		close(trees)
 	}()
 	consens, err := tree.Consensus(trees, 0.5)
@@ -168,9 +168,9 @@ func TestConsensus2(t *testing.T) {
 	}
 
 	go func() {
-		trees2 <- tree.Trees{randtree1, 1, nil}
-		trees2 <- tree.Trees{randtree2, 2, nil}
-		trees2 <- tree.Trees{randtree3, 3, nil}
+		trees2 <- tree.Trees{Tree: randtree1, Id: 1, Err: nil}
+		trees2 <- tree.Trees{Tree: randtree2, Id: 2, Err: nil}
+		trees2 <- tree.Trees{Tree: randtree3, Id: 3, Err: nil}
 		close(trees2)
 	}()
 	consens, err = tree.Consensus(trees2, 1)

tree/algo.go

@@ -329,7 +329,7 @@ func Consensus(trees <-chan Trees, cutoff float64) (*Tree, error) {
 			if idx, err := startree.TipIndex(n); err != nil {
 				return nil, err
 			} else {
-				if bs.key.Bitset().Test(idx) {
+				if bs.key.Bitset().Test(uint(idx)) {
 					names = append(names, n)
 				}
 			}

tree/edge.go

@@ -35,7 +35,8 @@ const (
 	NIL_SUPPORT = -1.0
 	NIL_LENGTH  = -1.0
 	NIL_PVALUE  = -1.0
-	NIL_ID      = -1.0
+	NIL_ID      = -1
+	NIL_TIPID   = 0
 )
 
 /* Edge functions */
@@ -136,9 +137,6 @@ func (e *Edge) Bitset() *bitset.BitSet {
 //
 // Returns an error if not initialized.
 func (e *Edge) Id() int {
-	if e.id == NIL_ID {
-		io.ExitWithMessage(errors.New("Id has not been set"))
-	}
 	return e.id
 }
 
@@ -258,23 +256,17 @@ func (e *Edge) TipPresent(id uint) bool {
 // Number of tips on the right side of the bipartition
 // Used by "TopoDepth" function for example.
 //
-// Bitsets must be initialized, otherwise returns an error.
-func (e *Edge) NumTipsRight() (int, error) {
-	if e.bitset == nil {
-		return -1, errors.New("Cannot count right tips, subtree sizes not initialized")
-	}
-	return e.ntaxright, nil
+// Must Be initialized with ComputeEdgeHashes.
+func (e *Edge) NumTipsRight() int {
+	return e.ntaxright
 }
 
 // Number of tips on the left side of the bipartition
 // Used by "TopoDepth" function for example.
 //
-// Bitsets must be initialized, otherwise returns an error.
-func (e *Edge) NumTipsLeft() (int, error) {
-	if e.ntaxleft == 0 {
-		return -1, errors.New("Cannot count left tips, subtree sizes not initialized ")
-	}
-	return e.ntaxleft, nil
+// Must be initialized with ComputeEdgeHashes.
+func (e *Edge) NumTipsLeft() int {
+	return e.ntaxleft
 }
 
 // Return the given edge in the array of edges comparing bitsets fields

tree/node.go

@@ -15,7 +15,8 @@ type Node struct {
 	neigh   []*Node  // neighbors array
 	br      []*Edge  // Branches array (same order than neigh)
 	depth   int      // Depth of the node
-	id      int      // this field is used at discretion of the user to store information
+	id      int      // This field is used at discretion of the user to store information
+	tipid   int      // This is used by TipIndex to match tip names of different trees
 }
 
 // Uninitialized depth is coded as -1
@@ -87,6 +88,12 @@ func (n *Node) SetId(id int) {
 	n.id = id
 }
 
+// Returns the Tip Id of the node. TipId==NIL_TIPIDINDEX means that
+// it is not a tip or tipindex has not been initialized (use ReinitIndexes)
+func (n *Node) TipIndex() int {
+	return n.tipid
+}
+
 // Returns the depth of the node. Returns an error if the depth has
 // not been set yet.
 func (n *Node) Depth() (int, error) {

tree/quartets.go

@@ -2,6 +2,7 @@ package tree
 
 import (
 	"errors"
+
 	"github.com/evolbioinfo/gotree/hashmap"
 	"github.com/evolbioinfo/gotree/io"
 )
@@ -238,7 +239,7 @@ func postOrderQuartetSet(t *Tree, n *Node, prev *Node, right [][]uint) []uint {
 		if err != nil {
 			io.ExitWithMessage(err)
 		}
-		output = append(output, taxindex)
+		output = append(output, uint(taxindex))
 	} else {
 		for _, next := range n.Neigh() {
 			if next != prev {