fix weighted unifrac bug in findIndex and randomize labels

diff --git a/Mothur.xcodeproj/project.pbxproj b/Mothur.xcodeproj/project.pbxproj
index 62d929a1ba847330373bd88396c2904869525f37..483badd56a1de0838f8337547f0c3038bdd1322f 100644 (file)
--- a/Mothur.xcodeproj/project.pbxproj
+++ b/Mothur.xcodeproj/project.pbxproj
@@ -740,6 +740,7 @@
                                GCC_OPTIMIZATION_LEVEL = 3;
                                GCC_WARN_ABOUT_RETURN_TYPE = YES;
                                GCC_WARN_UNUSED_VARIABLE = YES;
+                               INSTALL_MODE_FLAG = "a-w,a+rX";
                                PREBINDING = NO;
                                SDKROOT = "$(DEVELOPER_SDK_DIR)/MacOSX10.5.sdk";
                        };
@@ -755,6 +756,7 @@
                                GCC_OPTIMIZATION_LEVEL = 0;
                                GCC_WARN_ABOUT_RETURN_TYPE = YES;
                                GCC_WARN_UNUSED_VARIABLE = YES;
+                               INSTALL_MODE_FLAG = "a-w,a+rX";
                                PREBINDING = NO;
                                SDKROOT = "$(DEVELOPER_SDK_DIR)/MacOSX10.5.sdk";
                        };

diff --git a/parsimony.cpp b/parsimony.cpp
index da8412ce71d590fe57366fb453126db00ff0f715..3840514523002b545b7cf0f2fbb848b4c0ad0dd5 100644 (file)
--- a/parsimony.cpp
+++ b/parsimony.cpp
@@ -10,6 +10,7 @@
 #include "parsimony.h"
 
 /**************************************************************************************************/
+
 EstOutput Parsimony::getValues(Tree* t) {
        try {
                globaldata = GlobalData::getInstance();
@@ -38,12 +39,8 @@ EstOutput Parsimony::getValues(Tree* t) {
                        else if(iSize > rcSize || iSize > lcSize){
                                score++;
                        }
-                       cout << i << ' ' << score << ": ";
-                       t->tree[i].printNode();
                } 
                
-               //string hold;
-               //cin >> hold;
                
                data[0] = score;
                

diff --git a/tree.cpp b/tree.cpp
index 2689b037f6a494474ce1e008ba5efb1d250f0791..3b3b6840a6fd1a427258e7b84d0afe3bf181bd7e 100644 (file)
--- a/tree.cpp
+++ b/tree.cpp
@@ -354,33 +354,27 @@ void Tree::randomLabels() {
 
 void Tree::randomLabels(string groupA, string groupB) {
        try {
-               for(int i = 0; i < numLeaves; i++) {
-                       int z;
-                       //get random index to switch with
-                       z = int((float)(i+1) * (float)(rand()) / ((float)RAND_MAX+1.0));        
-                       
-                       //you only want to randomize the nodes that are from a group the user wants analyzed, so
-                       //if either of the leaf nodes you are about to switch are not in the users groups then you don't want to switch them.
-                       if (((tree[z].getGroup() == groupA) || (tree[z].getGroup() == groupB)) && ((tree[i].getGroup() == groupA) || (tree[i].getGroup() == groupB))) {
-                               //switches node i and node z's info.
-                               map<string,int> lib_hold = tree[z].pGroups;
-                               tree[z].pGroups = (tree[i].pGroups);
-                               tree[i].pGroups = (lib_hold);
-                               
-                               string zgroup = tree[z].getGroup();
-                               tree[z].setGroup(tree[i].getGroup());
-                               tree[i].setGroup(zgroup);
-                               
-                               string zname = tree[z].getName();
-                               tree[z].setName(tree[i].getName());
-                               tree[i].setName(zname);
+               int numSeqsA = globaldata->gTreemap->seqsPerGroup[groupA];
+               int numSeqsB = globaldata->gTreemap->seqsPerGroup[groupB];
+
+               vector<string> randomGroups(numSeqsA+numSeqsB, groupA);
+               for(int i=numSeqsA;i<randomGroups.size();i++){
+                       randomGroups[i] = groupB;
+               }
+               random_shuffle(randomGroups.begin(), randomGroups.end());
                                
-                               map<string,int> gcount_hold = tree[z].pcount;
-                               tree[z].pcount = (tree[i].pcount);
-                               tree[i].pcount = (gcount_hold);
+               int randomCounter = 0;                          
+               for(int i=0;i<numLeaves;i++){
+                       if(tree[i].getGroup() == groupA || tree[i].getGroup() == groupB){
+                               tree[i].setGroup(randomGroups[randomCounter]);
+                               tree[i].pcount.clear();
+                               tree[i].pcount[randomGroups[randomCounter]] = 1;
+                               tree[i].pGroups.clear();
+                               tree[i].pGroups[randomGroups[randomCounter]] = 1;
+                               randomCounter++;
                        }
                }
-       }
+       }               
        catch(exception& e) {
                cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function randomLabels. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
                exit(1);
@@ -528,7 +522,7 @@ void Tree::printBranch(int node) {
                        printBranch(tree[node].getRChild());
                        out << ")";
                }else { //you are a leaf
-                       out << tree[node].getName() << ":" << tree[node].getBranchLength();
+                       out << tree[node].getGroup() << ":" << tree[node].getBranchLength();
                }
                
        }

diff --git a/unifracweightedcommand.cpp b/unifracweightedcommand.cpp
index 8fc5b3bcedd56c76200dfcb8889374906c9f28b3..a99459deeef74b6e99eab76c43b7a76787366b20 100644 (file)
--- a/unifracweightedcommand.cpp
+++ b/unifracweightedcommand.cpp
@@ -72,10 +72,10 @@ int UnifracWeightedCommand::execute() {
                        
                        //get scores for random trees
                        for (int j = 0; j < iters; j++) {
-                               int n = 1;
+//                             int n = 1;
                                int count = 0;
                                for (int r=0; r<numGroups; r++) { 
-                                       for (int l = n; l < numGroups; l++) {
+                                       for (int l = r+1; l < numGroups; l++) {
                                                //copy T[i]'s info.
                                                randT->getCopy(T[i]);
                                                 
@@ -90,18 +90,18 @@ int UnifracWeightedCommand::execute() {
                                                        //get wscore of random tree
                                                        randomData = weighted->getValues(randT, tmap->namesOfGroups[r], tmap->namesOfGroups[l]);
                                                }
+//                                             randT->createNewickFile("hold"+toString(r)+toString(l)+toString(j));
 
                                                //save scores
                                                rScores[count].push_back(randomData[0]);
                                                validScores[count][randomData[0]] = randomData[0];
                                                count++;
                                        }
-                                       n++;
+//                                     n++;
                                }
                        }
 
                        removeValidScoresDuplicates(); 
-                       
                        //find the signifigance of the score for summary file
                        for (int f = 0; f < numComp; f++) {
                                //sort random scores
@@ -111,7 +111,7 @@ int UnifracWeightedCommand::execute() {
                                //so if you have 1000 random trees the index returned is 100 
                                //then there are 900 trees with a score greater then you. 
                                //giving you a signifigance of 0.900
-                               int index = findIndex(userData[f]);    if (index == -1) { cout << "error in UnifracWeightedCommand" << endl; exit(1); } //error code
+                               int index = findIndex(userData[f], f);    if (index == -1) { cout << "error in UnifracWeightedCommand" << endl; exit(1); } //error code
                        
                                //the signifigance is the number of trees with the users score or higher 
                                WScoreSig.push_back((iters-index)/(float)iters);
@@ -228,15 +228,12 @@ void UnifracWeightedCommand::removeValidScoresDuplicates() {
 }
 
 /***********************************************************/
-int UnifracWeightedCommand::findIndex(float score) {
+int UnifracWeightedCommand::findIndex(float score, int index) {
        try{
-               for (int e = 0; e < numComp; e++) {
-                       for (int i = 0; i < rScores[e].size(); i++) {
-//cout << rScores[e][i] << " number " << i << endl;
-                               if (rScores[e][i] >= score) { return i; }
-                       }
+               for (int i = 0; i < rScores[index].size(); i++) {
+                       if (rScores[index][i] >= score) {       return i;       }
                }
-               return -1;
+               return rScores[index].size();
        }
        catch(exception& e) {
                cout << "Standard Error: " << e.what() << " has occurred in the UnifracWeightedCommand class Function findIndex. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";

diff --git a/unifracweightedcommand.h b/unifracweightedcommand.h
index d309566e45b0d58fb2631eeb1a0ad1816907f437..e8daaefccbead969b54f2996f83ae378fc407a1c 100644 (file)
--- a/unifracweightedcommand.h
+++ b/unifracweightedcommand.h
@@ -47,7 +47,7 @@ class UnifracWeightedCommand : public Command {
                void printWSummaryFile();
        //      void printWeightedFile();  
                void removeValidScoresDuplicates();
-               int findIndex(float);
+               int findIndex(float, int);
                void setGroups(); 
 };