]> git.donarmstrong.com Git - mothur.git/blobdiff - indicatorcommand.cpp
fixed indicator grouping problem
[mothur.git] / indicatorcommand.cpp
index 21c70324611794f0574eb07cd97f2d10980a8bef..01627c0ec1b9aa5d7971a564f131924b179ffc94 100644 (file)
@@ -301,8 +301,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                ofstream out;
                m->openOutputFile(outputFileName, out);
                out.setf(ios::fixed, ios::floatfield); out.setf(ios::showpoint);
-               out << "Node\tOTU#\tIndVal" << endl;
-               
+                               
                string treeOutputDir = outputDir;
                if (outputDir == "") {  treeOutputDir += m->hasPath(treefile);  }
                string outputTreeFileName = treeOutputDir + m->getRootName(m->getSimpleName(treefile)) + "indicator.tre";
@@ -319,7 +318,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                
                //you need the distances to leaf to decide grouping below
                //this will also set branch lengths if the tree does not include them
-               map<int, float> distToLeaf = getLengthToLeaf(T);
+               map<int, float> distToRoot = getDistToRoot(T);
                        
                //for each node
                for (int i = T->getNumLeaves(); i < T->getNumNodes(); i++) {
@@ -335,18 +334,10 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                        if (sharedfile != "") {
                                vector< vector<SharedRAbundVector*> > groupings;
                                
-                               /*groupings.resize(1);
-                               groupings[0].push_back(lookup[0]);
-                               groupings[0].push_back(lookup[1]);
-                               groupings[0].push_back(lookup[2]);
-                               groupings[0].push_back(lookup[3]);
-                               groupings[0].push_back(lookup[4]);*/
-                               
                                //get nodes that will be a valid grouping
                                //you are valid if you are not one of my descendants
-                               //AND your distToLeaf is <= mine
-                               //AND your distToLeaf is >= my smallest childs
-                               //AND you were not added as part of a larger groupings
+                               //AND your distToRoot is >= mine
+                               //AND you were not added as part of a larger grouping. Largest nodes are added first.
                                
                                set<string> groupsAlreadyAdded;
                                //create a grouping with my grouping
@@ -366,7 +357,9 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                                
                                
                                for (int j = (T->getNumNodes()-1); j >= 0; j--) {
-                                       if ((descendantNodes[i].count(j) == 0) && (distToLeaf[j] <= distToLeaf[i]) && ((distToLeaf[j] >= distToLeaf[T->tree[i].getLChild()]) || (distToLeaf[j] >= distToLeaf[T->tree[i].getRChild()]))) {
+       
+                                       
+                                       if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) { 
                                                vector<SharedRAbundVector*> subset;
                                                int count = 0;
                                                int doneCount = nodeToDescendants[j].size();
@@ -385,10 +378,10 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                                        }
                                }
                                
-                               if (groupsAlreadyAdded.size() != lookup.size()) { cout << i << '\t' << groupsAlreadyAdded.size() << '\t' << lookup.size() << endl; m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
-                               for (int k = 0; k < lookup.size(); k++) {
-                                       if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; }
-                               }
+                               if (groupsAlreadyAdded.size() != lookup.size()) {  m->mothurOut("[ERROR]: could not make proper groupings."); m->mothurOutEndLine(); }
+                               //for (int k = 0; k < lookup.size(); k++) {
+                               //      if (groupsAlreadyAdded.count(lookup[k]->getGroup()) == 0) { cout << lookup[k]->getGroup() << endl; }
+                               //}
                                
                                indicatorValues = getValues(groupings);
                                
@@ -397,9 +390,8 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                                
                                //get nodes that will be a valid grouping
                                //you are valid if you are not one of my descendants
-                               //AND your distToLeaf is <= mine
-                               //AND your distToLeaf is >= my smallest childs
-                               //AND you were not added as part of a larger grouping
+                               //AND your distToRoot is >= mine
+                               //AND you were not added as part of a larger grouping. Largest nodes are added first.
                                
                                set<string> groupsAlreadyAdded;
                                //create a grouping with my grouping
@@ -418,7 +410,7 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                                if (subset.size() != 0) { groupings.push_back(subset); }
                                
                                for (int j = (T->getNumNodes()-1); j >= 0; j--) {
-                                       if ((descendantNodes[i].count(j) == 0) && (distToLeaf[j] <= distToLeaf[i]) && ((distToLeaf[j] >= distToLeaf[T->tree[i].getLChild()]) || (distToLeaf[j] >= distToLeaf[T->tree[i].getRChild()]))) {
+                                       if ((descendantNodes[i].count(j) == 0) && (distToRoot[j] >= distToRoot[i])) {
                                                vector<SharedRAbundFloatVector*> subset;
                                                int count = 0;
                                                int doneCount = nodeToDescendants[j].size();
@@ -444,15 +436,18 @@ int IndicatorCommand::GetIndicatorSpecies(Tree*& T){
                        
                        if (m->control_pressed) { out.close(); return 0; }
                        
+                       
                        /******************************************************/
                        //output indicator values to table form + label tree  //
                        /*****************************************************/
+                       out << (i+1) << '\t';
                        for (int j = 0; j < indicatorValues.size(); j++) {
                                
                                if (m->control_pressed) { out.close(); return 0; }
                                
-                               out << (i+1) << '\t' << (j+1) << '\t' << indicatorValues[j] << endl;
+                               out << indicatorValues[j] << '\t';
                        }
+                       out << endl;
                        
                        T->tree[i].setLabel((i+1));
                        
@@ -581,19 +576,23 @@ vector<float> IndicatorCommand::getValues(vector< vector<SharedRAbundVector*> >&
        }
 }
 //**********************************************************************************************************************
-//you need the distances to leaf to decide groupings
+//you need the distances to root to decide groupings
 //this will also set branch lengths if the tree does not include them
-map<int, float> IndicatorCommand::getLengthToLeaf(Tree*& T){
+map<int, float> IndicatorCommand::getDistToRoot(Tree*& T){
        try {
                map<int, float> dists;
                
-               for (int i = 0; i < T->getNumNodes(); i++) {
-                       
-                       int lc = T->tree[i].getLChild();
-                       int rc = T->tree[i].getRChild();
-                                
-                       //if you have no branch length, set it then calc
-                       if (T->tree[i].getBranchLength() <= 0.0) {
+               bool hasBranchLengths = false;
+               for (int i = 0; i < T->getNumNodes(); i++) { 
+                       if (T->tree[i].getBranchLength() > 0.0) {  hasBranchLengths = true; break; }
+               }
+               
+               //set branchlengths if needed
+               if (!hasBranchLengths) { 
+                       for (int i = 0; i < T->getNumNodes(); i++) {
+                               
+                               int lc = T->tree[i].getLChild();
+                               int rc = T->tree[i].getRChild();
                                
                                if (lc == -1) { // you are a leaf
                                        //if you are a leaf set you priliminary length to 1.0, this may adjust later
@@ -605,28 +604,32 @@ map<int, float> IndicatorCommand::getLengthToLeaf(Tree*& T){
                                        float rdist = dists[rc];
                                        
                                        float greater = ldist;
-                                       if (rdist > greater) { greater = rdist; dists[i] = ldist; }
-                                       else { dists[i] = rdist; }
+                                       if (rdist > greater) { greater = rdist; dists[i] = ldist + 1.0;}
+                                       else { dists[i] = rdist + 1.0; }
+                                       
                                        
                                        //branch length = difference + 1
                                        T->tree[lc].setBranchLength((abs(ldist-greater) + 1.0));
                                        T->tree[rc].setBranchLength((abs(rdist-greater) + 1.0));
                                }
-                               
-                       }else{
-                               if (lc == -1) { dists[i] = T->tree[i].getBranchLength(); }
-                               else { //smaller of my two children distances plus my branch length
-                                       //look at your children's length to leaf
-                                       float ldist = dists[lc];
-                                       float rdist = dists[rc];
-                                       
-                                       float smaller = ldist;
-                                       if (rdist < smaller) { smaller = rdist;  }
-                                       
-                                       dists[i] = smaller + T->tree[i].getBranchLength(); 
+                       }
+               }
+                       
+               dists.clear();
+               
+               for (int i = 0; i < T->getNumNodes(); i++) {
+                       
+                       double sum = 0.0;
+                       int index = i;
+                       
+                       while(T->tree[index].getParent() != -1){
+                               if (T->tree[index].getBranchLength() != -1) {
+                                       sum += abs(T->tree[index].getBranchLength()); 
                                }
+                               index = T->tree[index].getParent();
                        }
                        
+                       dists[i] = sum;
                }
                
                return dists;