fixed indicator grouping problem

diff --git a/indicatorcommand.cpp b/indicatorcommand.cpp
index 9f96f35e95c5b5f74d58690e8d485cb0987114a8..01627c0ec1b9aa5d7971a564f131924b179ffc94 100644 (file)
--- a/indicatorcommand.cpp
+++ b/indicatorcommand.cpp
@@ -318,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++) {
@@ -336,9 +336,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 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
@@ -358,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();
@@ -377,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);
                                
@@ -389,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
@@ -410,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();
@@ -576,9 +576,9 @@ 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;
                
@@ -587,13 +587,13 @@ map<int, float> IndicatorCommand::getLengthToLeaf(Tree*& T){
                        if (T->tree[i].getBranchLength() > 0.0) {  hasBranchLengths = true; break; }
                }
                
-               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 (!hasBranchLengths) {
+               //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
                                        T->tree[i].setBranchLength(1.0);
@@ -604,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 + 1.0; }
+                                       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;

diff --git a/indicatorcommand.h b/indicatorcommand.h
index 1a41ffd33c1e0a175ed6560b300ba9469fa20cb8..0ec100aecba915fb304e72146d06cbe2376ce8aa 100644 (file)
--- a/indicatorcommand.h
+++ b/indicatorcommand.h
@@ -47,7 +47,7 @@ private:
        set<string> getDescendantList(Tree*&, int, map<int, set<string> >, map<int, set<int> >&);
        vector<float> getValues(vector< vector<SharedRAbundVector*> >&);
        vector<float> getValues(vector< vector<SharedRAbundFloatVector*> >&);
-       map<int, float> getLengthToLeaf(Tree*&);
+       map<int, float> getDistToRoot(Tree*&);
        
 };
 

diff --git a/mothur b/mothur
index a9d33c7e5e04e7e7dd2c83ec43b70a54e02cdb7c..835aa5ccd580ae1232f8b123075f7e3052aeb06e 100755 (executable)
Binary files a/mothur and b/mothur differ

diff --git a/trimseqscommand.cpp b/trimseqscommand.cpp
index fce4763dcb891ef56a57f24644b9574ffa51cd1a..3c0d378f709042cc3e1c2ca2b87e730a62dbf359 100644 (file)
--- a/trimseqscommand.cpp
+++ b/trimseqscommand.cpp
@@ -360,12 +360,12 @@ int TrimSeqsCommand::execute(){
                #endif
                
                if (m->control_pressed) {  return 0; }                  
-               cout << "done with driver " << endl;                                                                    
+                                                                               
                for(int i=0;i<fastaFileNames.size();i++){
                        cout << fastaFileNames[i] << endl;
                        
-                       if (m->isBlank(fastaFileNames[i])) { cout << fastaFileNames[i] << " was blank" << endl; remove(fastaFileNames[i].c_str());      }
-                       else if (filesToRemove.count(fastaFileNames[i]) > 0) { cout << fastaFileNames[i] << " was on the remove list" << endl; remove(fastaFileNames[i].c_str()); }
+                       if (m->isBlank(fastaFileNames[i])) {  remove(fastaFileNames[i].c_str());        }
+                       else if (filesToRemove.count(fastaFileNames[i]) > 0) { remove(fastaFileNames[i].c_str()); }
                        else {
                                ifstream inFASTA;
                                string seqName;
@@ -382,7 +382,7 @@ int TrimSeqsCommand::execute(){
                                                if(itCombo->second == i){       thisGroup = itCombo->first;     combos.erase(itCombo);  break;  }
                                        }
                                }else{ thisGroup = groupVector[i]; }
-                               cout << thisGroup << '\t' << i  << '\t' << comboStarts << endl; 
+                                       
                                while(!inFASTA.eof()){
                                        if(inFASTA.get() == '>'){
                                                inFASTA >> seqName;
@@ -394,12 +394,12 @@ int TrimSeqsCommand::execute(){
                                inFASTA.close();
                        }
                }
-       cout << "done with fastaFileNames " << endl;            
+               
                if(qFileName != ""){
                        for(int i=0;i<qualFileNames.size();i++){
                                cout << qualFileNames[i] << endl;
-                               if (m->isBlank(qualFileNames[i])) { cout << qualFileNames[i] << " was blank" << endl; remove(qualFileNames[i].c_str()); }
-                               else if (filesToRemove.count(qualFileNames[i]) > 0) { cout << qualFileNames[i] << " was on the remove list" << endl; remove(qualFileNames[i].c_str()); }
+                               if (m->isBlank(qualFileNames[i])) {  remove(qualFileNames[i].c_str());  }
+                               else if (filesToRemove.count(qualFileNames[i]) > 0) {  remove(qualFileNames[i].c_str()); }
                                else {
                                        ifstream inQual;
                                        string seqName;
@@ -419,7 +419,7 @@ int TrimSeqsCommand::execute(){
                                }
                        }
                }
-               cout << "done with qualFileNames " << endl;
+               
                
                if (m->control_pressed) { 
                        for (int i = 0; i < outputNames.size(); i++) {  remove(outputNames[i].c_str()); }
@@ -460,8 +460,6 @@ int TrimSeqsCommand::driverCreateTrim(string filename, string qFileName, string
                }
                
                ofstream outGroups;
-               //vector<ofstream*> fastaFileNames;
-               //vector<ofstream*> qualFileNames;
                
                if (oligoFile != "") {          
                        m->openOutputFile(groupFile, outGroups);   
@@ -753,7 +751,7 @@ int TrimSeqsCommand::createProcessesCreateTrim(string filename, string qFileName
                                }
                        }
                        
-                       if (allFiles) { m->mothurOut("Done with allfile"); m->mothurOutEndLine(); }
+                       if (allFiles) { m->mothurOut("Done with allfiles"); m->mothurOutEndLine(); }
                }
        
                return exitCommand;