]> git.donarmstrong.com Git - mothur.git/blobdiff - weighted.cpp
sffinfo bug with flow grams right index when clipQualRight=0
[mothur.git] / weighted.cpp
index a0d593c803f691821afdb680e79be80f179b0afa..cf1291dfa91af84e4d821541b2a8f43a169d64f7 100644 (file)
 
 EstOutput Weighted::getValues(Tree* t, int p, string o) {
     try {
-               globaldata = GlobalData::getInstance();
-               
                data.clear(); //clear out old values
                int numGroups;
                vector<double> D;
                processors = p;
                outputDir = o;
+        
+        CountTable* ct = t->getCountTable();
                
-               numGroups = globaldata->Groups.size();
+               numGroups = m->getNumGroups();
                
                if (m->control_pressed) { return data; }
                
@@ -31,14 +31,14 @@ EstOutput Weighted::getValues(Tree* t, int p, string o) {
                        for (int l = 0; l < i; l++) {   
                                //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]);
+                               vector<string> groups; groups.push_back((m->getGroups())[i]); groups.push_back((m->getGroups())[l]);
                                namesOfGroupCombos.push_back(groups);
                        }
                }
                
-               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+               #if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
                        if(processors == 1){
-                               data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size());
+                               data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), ct);
                        }else{
                                int numPairs = namesOfGroupCombos.size();
                                
@@ -52,12 +52,12 @@ EstOutput Weighted::getValues(Tree* t, int p, string o) {
                                        lines.push_back(linePair(startPos, numPairsPerProcessor));
                                }
 
-                               data = createProcesses(t, namesOfGroupCombos);
+                               data = createProcesses(t, namesOfGroupCombos, ct);
                                
                                lines.clear();
                        }
                #else
-                       data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size());
+                       data = driver(t, namesOfGroupCombos, 0, namesOfGroupCombos.size(), ct);
                #endif
                
                return data;
@@ -69,9 +69,9 @@ EstOutput Weighted::getValues(Tree* t, int p, string o) {
 }
 /**************************************************************************************************/
 
-EstOutput Weighted::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos) {
+EstOutput Weighted::createProcesses(Tree* t, vector< vector<string> > namesOfGroupCombos, CountTable* ct) {
        try {
-#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux)
+#if defined (__APPLE__) || (__MACH__) || (linux) || (__linux) || (__linux__) || (__unix__) || (__unix)
                int process = 1;
                vector<int> processIDS;
                
@@ -87,9 +87,9 @@ EstOutput Weighted::createProcesses(Tree* t, vector< vector<string> > namesOfGro
                        }else if (pid == 0){
        
                                EstOutput Myresults;
-                               Myresults = driver(t, namesOfGroupCombos, lines[process].start, lines[process].num);
+                               Myresults = driver(t, namesOfGroupCombos, lines[process].start, lines[process].num, ct);
                        
-                               m->mothurOut("Merging results."); m->mothurOutEndLine();
+                               //m->mothurOut("Merging results."); m->mothurOutEndLine();
                                
                                //pass numSeqs to parent
                                ofstream out;
@@ -110,7 +110,7 @@ EstOutput Weighted::createProcesses(Tree* t, vector< vector<string> > namesOfGro
                        }
                }
        
-               results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num);
+               results = driver(t, namesOfGroupCombos, lines[0].start, lines[0].num, ct);
        
                //force parent to wait until all the processes are done
                for (int i=0;i<(processors-1);i++) { 
@@ -141,10 +141,10 @@ EstOutput Weighted::createProcesses(Tree* t, vector< vector<string> > namesOfGro
                                m->gobble(in);
                        }
                        in.close();
-                       remove(s.c_str());
+                       m->mothurRemove(s);
                }
                
-               m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine();
+               //m->mothurOut("DONE."); m->mothurOutEndLine(); m->mothurOutEndLine();
                
                return results;
 #endif         
@@ -155,7 +155,7 @@ EstOutput Weighted::createProcesses(Tree* t, vector< vector<string> > namesOfGro
        }
 }
 /**************************************************************************************************/
-EstOutput Weighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num) { 
+EstOutput Weighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos, int start, int num, CountTable* ct) { 
  try {
                EstOutput results;
                vector<double> D;
@@ -179,30 +179,28 @@ EstOutput Weighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos,
                                int numSeqsInGroupI = it->second;
                                
                                double sum = getLengthToRoot(t, t->groupNodeInfo[groupA][j], groupA, groupB);
-                               double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groupA]);
+                               double weightedSum = ((numSeqsInGroupI * sum) / (double)ct->getGroupCount(groupA));
                        
                                D[count] += weightedSum;
                        }
-
+                       
                        //adding the wieghted sums from group l
                        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[t->groupNodeInfo[groupB][j]].pcount.find(groupB);
                                int numSeqsInGroupL = it->second;
                                
                                double sum = getLengthToRoot(t, t->groupNodeInfo[groupB][j], groupA, groupB);
-                               double weightedSum = ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groupB]);
+                               double weightedSum = ((numSeqsInGroupL * sum) / (double)ct->getGroupCount(groupB));
                        
                                D[count] += weightedSum;
                        }
-                       
                        count++;
                }
-               
+        
                //calculate u for the group comb 
-               
                for (int h = start; h < (start+num); h++) {     
                        //report progress
-                       m->mothurOut("Processing combo: " + toString(h)); m->mothurOutEndLine();
+                       //m->mothurOut("Processing combo: " + toString(h)); m->mothurOutEndLine();
                                                
                        string groupA = namesOfGroupCombos[h][0]; 
                        string groupB = namesOfGroupCombos[h][1];
@@ -218,27 +216,34 @@ EstOutput Weighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos,
                                it = t->tree[i].pcount.find(groupA);
                                //if it does u = # of its descendants with a certain group / total number in tree with a certain group
                                if (it != t->tree[i].pcount.end()) {
-                                       u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA];
+                                       u = (double) t->tree[i].pcount[groupA] / (double) ct->getGroupCount(groupA);
                                }else { u = 0.00; }
                                
                                
                                //does this node have descendants from group l
                                it = t->tree[i].pcount.find(groupB);
+                               
                                //if it does subtract their percentage from u
                                if (it != t->tree[i].pcount.end()) {
-                                       u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB];
+                                       u -= (double) t->tree[i].pcount[groupB] / (double) ct->getGroupCount(groupB);
                                }
                                
-                               //if this is not the root then add it
-                               if (rootForGrouping[namesOfGroupCombos[h]].count(i) == 0) {
+                               if (includeRoot) {
                                        if (t->tree[i].getBranchLength() != -1) {
                                                u = abs(u * t->tree[i].getBranchLength());
                                                WScore[(groupA+groupB)] += u; 
                                        }
+                               }else {
+                                       //if this is not the root then add it
+                                       if (rootForGrouping[namesOfGroupCombos[h]].count(i) == 0) {
+                                               if (t->tree[i].getBranchLength() != -1) {
+                                                       u = abs(u * t->tree[i].getBranchLength());
+                                                       WScore[(groupA+groupB)] += u; 
+                                               }
+                                       }
                                }
-                               
-                               
                        }
+                       
                }
                
                /********************************************************/
@@ -262,9 +267,10 @@ EstOutput Weighted::driver(Tree* t, vector< vector<string> > namesOfGroupCombos,
 /**************************************************************************************************/
 EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) { 
  try {
-               globaldata = GlobalData::getInstance();
                
                data.clear(); //clear out old values
+     
+        CountTable* ct = t->getCountTable();
                
                if (m->control_pressed) { return data; }
                
@@ -281,7 +287,7 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                        int numSeqsInGroupI = it->second;
                        
                        double sum = getLengthToRoot(t, t->groupNodeInfo[groups[0]][j], groups[0], groups[1]);
-                       double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]);
+                       double weightedSum = ((numSeqsInGroupI * sum) / (double)ct->getGroupCount(groups[0]));
                
                        D += weightedSum;
                }
@@ -292,7 +298,7 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                        int numSeqsInGroupL = it->second;
                        
                        double sum = getLengthToRoot(t, t->groupNodeInfo[groups[1]][j], groups[0], groups[1]);
-                       double weightedSum = ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
+                       double weightedSum = ((numSeqsInGroupL * sum) / (double)ct->getGroupCount(groups[1]));
                
                        D += weightedSum;
                }
@@ -308,7 +314,7 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                        it = t->tree[i].pcount.find(groupA);
                        //if it does u = # of its descendants with a certain group / total number in tree with a certain group
                        if (it != t->tree[i].pcount.end()) {
-                               u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA];
+                               u = (double) t->tree[i].pcount[groupA] / (double) ct->getGroupCount(groupA);
                        }else { u = 0.00; }
                        
                        
@@ -316,15 +322,22 @@ EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
                        it = t->tree[i].pcount.find(groupB);
                        //if it does subtract their percentage from u
                        if (it != t->tree[i].pcount.end()) {
-                               u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB];
+                               u -= (double) t->tree[i].pcount[groupB] / (double) ct->getGroupCount(groupB);
                        }
                        
-                       //if this is not the root then add it
-                       if (rootForGrouping[groups].count(i) == 0) {
+                       if (includeRoot) {
                                if (t->tree[i].getBranchLength() != -1) {
                                        u = abs(u * t->tree[i].getBranchLength());
                                        WScore[(groupA+groupB)] += u;
                                }
+                       }else{
+                               //if this is not the root then add it
+                               if (rootForGrouping[groups].count(i) == 0) {
+                                       if (t->tree[i].getBranchLength() != -1) {
+                                               u = abs(u * t->tree[i].getBranchLength());
+                                               WScore[(groupA+groupB)] += u;
+                                       }
+                               }
                        }
                }               
                /********************************************************/
@@ -363,37 +376,43 @@ double Weighted::getLengthToRoot(Tree* t, int v, string groupA, string groupB) {
                while(t->tree[index].getParent() != -1){
 
                        if (m->control_pressed) {  return sum; }
-                               
-                       //am I the root for this grouping? if so I want to stop "early"
-                       //does my sibling have descendants from the users groups? 
-                       int parent = t->tree[index].getParent();
-                       int lc = t->tree[parent].getLChild();
-                       int rc = t->tree[parent].getRChild();
                        
-                       int sib = lc;
-                       if (lc == index) { sib = rc; }
+                       int parent = t->tree[index].getParent();
                        
-                       map<string, int>::iterator itGroup;
-                       int pcountSize = 0;
-                       itGroup = t->tree[sib].pcount.find(groupA);
-                       if (itGroup != t->tree[sib].pcount.end()) { pcountSize++;  } 
-                       itGroup = t->tree[sib].pcount.find(groupB);
-                       if (itGroup != t->tree[sib].pcount.end()) { pcountSize++;  } 
-                                                       
-                       //if yes, I am not the root so add me
-                       if (pcountSize != 0) {
-                               if (t->tree[index].getBranchLength() != -1) {
-                                       sum += abs(t->tree[index].getBranchLength()) + tempTotal;
-                                       tempTotal = 0.0;
-                               }else {
-                                       sum += tempTotal;
-                                       tempTotal = 0.0;
-                               }
-                               rootForGrouping[grouping].clear();
-                               rootForGrouping[grouping].insert(parent);
-                       }else { //if no, I may be the root so add my br to tempTotal until I am proven innocent
-                               if (t->tree[index].getBranchLength() != -1) {
-                                       tempTotal += abs(t->tree[index].getBranchLength()); 
+                       if (includeRoot) { //add everyone
+                               if(t->tree[index].getBranchLength() != -1){     sum += abs(t->tree[index].getBranchLength());   }
+                       }else {
+                               
+                               //am I the root for this grouping? if so I want to stop "early"
+                               //does my sibling have descendants from the users groups? 
+                               int lc = t->tree[parent].getLChild();
+                               int rc = t->tree[parent].getRChild();
+                               
+                               int sib = lc;
+                               if (lc == index) { sib = rc; }
+                               
+                               map<string, int>::iterator itGroup;
+                               int pcountSize = 0;
+                               itGroup = t->tree[sib].pcount.find(groupA);
+                               if (itGroup != t->tree[sib].pcount.end()) { pcountSize++;  } 
+                               itGroup = t->tree[sib].pcount.find(groupB);
+                               if (itGroup != t->tree[sib].pcount.end()) { pcountSize++;  } 
+                               
+                               //if yes, I am not the root so add me
+                               if (pcountSize != 0) {
+                                       if (t->tree[index].getBranchLength() != -1) {
+                                               sum += abs(t->tree[index].getBranchLength()) + tempTotal;
+                                               tempTotal = 0.0;
+                                       }else {
+                                               sum += tempTotal;
+                                               tempTotal = 0.0;
+                                       }
+                                       rootForGrouping[grouping].clear();
+                                       rootForGrouping[grouping].insert(parent);
+                               }else { //if no, I may be the root so add my br to tempTotal until I am proven innocent
+                                       if (t->tree[index].getBranchLength() != -1) {
+                                               tempTotal += abs(t->tree[index].getBranchLength()); 
+                                       }
                                }
                        }
                        
@@ -402,12 +421,12 @@ double Weighted::getLengthToRoot(Tree* t, int v, string groupA, string groupB) {
                
                //get all nodes above the root to add so we don't add their u values above
                index = *(rootForGrouping[grouping].begin());
+               
                while(t->tree[index].getParent() != -1){
                        int parent = t->tree[index].getParent();
                        rootForGrouping[grouping].insert(parent);
                        index = parent;
                }
-               
                return sum;
        }
        catch(exception& e) {