]> git.donarmstrong.com Git - mothur.git/blobdiff - weighted.cpp
sped up unifrac weighted.
[mothur.git] / weighted.cpp
index 581078085d4e573ad100dbc48f1ec87d64558efa..c8b64e20446d363af2a6413964a25a1ad7b88dd4 100644 (file)
@@ -18,6 +18,10 @@ EstOutput Weighted::getValues(Tree* t) {
                vector<double> D;
                
                numGroups = globaldata->Groups.size();
+       
+               vector<double> sums = getBranchLengthSums(t);
+               
+               if (m->control_pressed) { return data; }
                
                //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
                int count = 0;
@@ -26,35 +30,58 @@ EstOutput Weighted::getValues(Tree* t) {
                                //initialize weighted scores
                                WScore[globaldata->Groups[i]+globaldata->Groups[l]] = 0.0;
                                
+                               vector<string> groups; groups.push_back(globaldata->Groups[i]); groups.push_back(globaldata->Groups[l]);
+                               
                                D.push_back(0.0000); //initialize a spot in D for each combination
+                       
+                               //adding the wieghted sums from group i
+                               for (int j = 0; j < t->groupNodeInfo[groups[0]].size(); j++) { //the leaf nodes that have seqs from group i
+                                       map<string, int>::iterator it = t->tree[t->groupNodeInfo[groups[0]][j]].pcount.find(groups[0]);
+                                       int numSeqsInGroupI = it->second;
+                                       
+                                       double weightedSum = ((numSeqsInGroupI * sums[t->groupNodeInfo[groups[0]][j]]) / (double)tmap->seqsPerGroup[groups[0]]);
+                               
+                                       D[count] += weightedSum;
+                               }
                                
-                               /********************************************************/
+                               //adding the wieghted sums from group l
+                               for (int j = 0; j < t->groupNodeInfo[groups[1]].size(); j++) { //the leaf nodes that have seqs from group l
+                                       map<string, int>::iterator it = t->tree[t->groupNodeInfo[groups[1]][j]].pcount.find(groups[1]);
+                                       int numSeqsInGroupL = it->second;
+                                       
+                                       double weightedSum = ((numSeqsInGroupL * sums[t->groupNodeInfo[groups[1]][j]]) / (double)tmap->seqsPerGroup[groups[1]]);
+                               
+                                       D[count] += weightedSum;
+                               }
+
+                               /********************************************************
                                //calculate a D value for each group combo
                                for(int v=0;v<t->getNumLeaves();v++){
-                                       int index = v;
-                                       double sum = 0.0000;
-               
-                                       //while you aren't at root
-                                       while(t->tree[index].getParent() != -1){
-                                                       
-                                               //if you have a BL
-                                               if(t->tree[index].getBranchLength() != -1){
-                                                       sum += abs(t->tree[index].getBranchLength());
-                                               }
-                                               index = t->tree[index].getParent();
-                                       }
-                                               
-                                       //get last breanch length added
-                                       if(t->tree[index].getBranchLength() != -1){
-                                               sum += abs(t->tree[index].getBranchLength());
-                                       }
-                                               
+                               
+                                       if (m->control_pressed) { return data; }
+                                       
                                        //is this sum from a sequence which is in one of the users groups
-                                       if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
+                                       if (m->inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
                                                //is this sum from a sequence which is in this groupCombo
-                                               if ((t->tree[v].getGroup() == globaldata->Groups[i]) || (t->tree[v].getGroup() == globaldata->Groups[l])) {
-                                                       sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
-                                                       D[count] += sum; 
+                                               if (m->inUsersGroups(t->tree[v].getGroup(), groups)) {
+                                                       int numSeqsInGroupI, numSeqsInGroupL;
+                                                       
+                                                       map<string, int>::iterator it;
+                                                       it = t->tree[v].pcount.find(groups[0]);
+                                                       if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i
+                                                               numSeqsInGroupI = it->second;
+                                                       }else{ numSeqsInGroupI = 0; }
+                                                       
+                                                       it = t->tree[v].pcount.find(groups[1]);
+                                                       if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l
+                                                               numSeqsInGroupL = it->second;
+                                                       }else{ numSeqsInGroupL = 0; }
+                                                       
+                                                       double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
+
+                                                       //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
+                                                       
+                                                       D[count] += weightedSum; 
                                                }
                                        }
                                }
@@ -69,6 +96,9 @@ EstOutput Weighted::getValues(Tree* t) {
                        //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
                        for (int b=0; b<numGroups; b++) { 
                                for (int l = b+1; l < numGroups; l++) {
+                               
+                                       if (m->control_pressed) { return data; }
+                                       
                                        //calculate a u value for each combo
                                        double u;
                                        //does this node have descendants from group b-1
@@ -109,14 +139,9 @@ EstOutput Weighted::getValues(Tree* t) {
                return data;
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               m->errorOut(e, "Weighted", "getValues");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
-               exit(1);
-       }
-
 }
 
 /**************************************************************************************************/
@@ -126,14 +151,42 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                
                data.clear(); //clear out old values
                
+               vector<double> sums = getBranchLengthSums(t);
+               
+               if (m->control_pressed) { return data; }
+               
                //initialize weighted score
                WScore[(groupA+groupB)] = 0.0;
-               float D = 0.0;
+               double D = 0.0;
                
-                                               
-               /********************************************************/
+               vector<string> groups; groups.push_back(groupA); groups.push_back(groupB);
+               
+               //adding the wieghted sums from groupA
+               for (int j = 0; j < t->groupNodeInfo[groupA].size(); j++) { //the leaf nodes that have seqs from group i
+                       map<string, int>::iterator it = t->tree[j].pcount.find(groupA);
+                       int numSeqsInGroupA = it->second;
+                       
+                       double weightedSum = ((numSeqsInGroupA * sums[j]) / (double)tmap->seqsPerGroup[groupA]);
+               
+                       D += weightedSum;
+               }
+               
+               //adding the wieghted sums from groupB
+               for (int j = 0; j < t->groupNodeInfo[groupB].size(); j++) { //the leaf nodes that have seqs from group l
+                       map<string, int>::iterator it = t->tree[j].pcount.find(groupB);
+                       int numSeqsInGroupB = it->second;
+                       
+                       double weightedSum = ((numSeqsInGroupB * sums[j]) / (double)tmap->seqsPerGroup[groupB]);
+               
+                       D += weightedSum;
+               }
+                               
+                                                                                                               
+               /********************************************************
                //calculate a D value for the group combo
                for(int v=0;v<t->getNumLeaves();v++){
+                       if (m->control_pressed) { return data; }
+                       
                        int index = v;
                        double sum = 0.0000;
                
@@ -152,15 +205,34 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                                sum += abs(t->tree[index].getBranchLength());
                        }
                                                
-                       if ((t->tree[v].getGroup() == groupA) || (t->tree[v].getGroup() == groupB)) {
-                               sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
-                               D += sum; 
+                       if (m->inUsersGroups(t->tree[v].getGroup(), groups)) {
+                               int numSeqsInGroupI, numSeqsInGroupL;
+                                                       
+                               map<string, int>::iterator it;
+                               it = t->tree[v].pcount.find(groups[0]);
+                               if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i
+                                       numSeqsInGroupI = it->second;
+                               }else{ numSeqsInGroupI = 0; }
+                               
+                               it = t->tree[v].pcount.find(groups[1]);
+                               if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l
+                                       numSeqsInGroupL = it->second;
+                               }else{ numSeqsInGroupL = 0; }
+                               
+                               double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
+                               
+                               //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
+                               
+                               D += weightedSum; 
                        }
                }
                /********************************************************/                              
                
                //calculate u for the group comb 
                for(int i=0;i<t->getNumNodes();i++){
+               
+                       if (m->control_pressed) { return data; }
+                       
                        double u;
                        //does this node have descendants from groupA
                        it = t->tree[i].pcount.find(groupA);
@@ -195,21 +267,49 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                return data; 
        }
        catch(exception& e) {
-               cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+               m->errorOut(e, "Weighted", "getValues");
                exit(1);
        }
-       catch(...) {
-               cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
+}
+/**************************************************************************************************/
+vector<double> Weighted::getBranchLengthSums(Tree* t) { 
+       try {
+               
+               vector<double> sums;
+               
+               for(int v=0;v<t->getNumLeaves();v++){
+                       
+                       if (m->control_pressed) { return sums; }
+                               
+                       int index = v;
+                       double sum = 0.0000;
+       
+                       //while you aren't at root
+                       while(t->tree[index].getParent() != -1){
+                                               
+                               //if you have a BL
+                               if(t->tree[index].getBranchLength() != -1){
+                                       sum += abs(t->tree[index].getBranchLength());
+                               }
+                               index = t->tree[index].getParent();
+                       }
+                                       
+                       //get last breanch length added
+                       if(t->tree[index].getBranchLength() != -1){
+                               sum += abs(t->tree[index].getBranchLength());
+                       }
+                       
+                       sums.push_back(sum);
+               }
+               
+               return sums;
+       }
+       catch(exception& e) {
+               m->errorOut(e, "Weighted", "getBranchLengthSums");
                exit(1);
        }
-
 }
-
-
-
-
-
-
+/**************************************************************************************************/