X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=consensuscommand.cpp;h=8c34068bba76f4347d8783b35c74af69979921ea;hb=7e354c9abb09ea3cf5b500a16cc7f6dd79ccb6f5;hp=9a90219d14c3809a687e053851440971eecb2e9b;hpb=92dde9a6d6c638fcbbd5dbaa5c79167564f90e49;p=mothur.git diff --git a/consensuscommand.cpp b/consensuscommand.cpp index 9a90219..8c34068 100644 --- a/consensuscommand.cpp +++ b/consensuscommand.cpp @@ -17,18 +17,9 @@ ConcensusCommand::ConcensusCommand(string fileroot){ abort = false; filename = fileroot; - //allow user to run help - //if(option == "help") { help(); abort = true; } - //else { - //if (option != "") { mothurOut("There are no valid parameters for the consensus command."); mothurOutEndLine(); abort = true; } - - // //no trees were read - // if (globaldata->gTree.size() == 0) { mothurOut("You must execute the read.tree command, before you may use the consensus command."); mothurOutEndLine(); abort = true; } - //else { - t = globaldata->gTree; - // } - //} + t = globaldata->gTree; + } catch(exception& e) { errorOut(e, "ConcensusCommand", "ConcensusCommand"); @@ -97,9 +88,10 @@ int ConcensusCommand::execute(){ out2 << "Species in Order: " << endl << endl; for (int i = 0; i < treeSet.size(); i++) { out2 << i+1 << ". " << treeSet[i] << endl; } - vector temp; //output sets included out2 << endl << "Sets included in the consensus tree:" << endl << endl; + + vector temp; for (it2 = nodePairsInTree.begin(); it2 != nodePairsInTree.end(); it2++) { //only output pairs not leaves if (it2->first.size() > 1) { @@ -112,6 +104,7 @@ int ConcensusCommand::execute(){ //find spot int index = findSpot(it2->first[i]); temp[index] = "*"; + //temp[index] = it2->first[i] + " "; } //output temp @@ -173,8 +166,8 @@ int ConcensusCommand::buildConcensusTree(vector nodeSet) { //terminate recursion }else if (count == numNodes) { return 0; } else { - leftChildSet = getNextAvailableSet(nodeSet); - rightChildSet = getRestSet(nodeSet, leftChildSet); + //finds best child pair + leftChildSet = getNextAvailableSet(nodeSet, rightChildSet); int left = buildConcensusTree(leftChildSet); int right = buildConcensusTree(rightChildSet); consensusTree->tree[count].setChildren(left, right); @@ -221,6 +214,7 @@ void ConcensusCommand::getSets() { } } + //add each leaf to terminate recursion in consensus //you want the leaves in there but with insignifigant sightings value so it is added last //for each leaf node get descendant info. @@ -238,33 +232,70 @@ void ConcensusCommand::getSets() { } sort(treeSet.begin(), treeSet.end()); + + + map< vector, int> nodePairsCopy = nodePairs; + + + //set initial rating on pairs to sightings + subgroup sightings + while (nodePairsCopy.size() != 0) { + + vector small = getSmallest(nodePairsCopy); + + int subgrouprate = getSubgroupRating(small); + + nodePairsInitialRate[small] = nodePairs[small] + subgrouprate; + + nodePairsCopy.erase(small); + } + } catch(exception& e) { errorOut(e, "ConcensusCommand", "getSets"); exit(1); } } +//********************************************************************************************************************** +vector ConcensusCommand::getSmallest(map< vector, int> nodes) { + try{ + vector smallest = nodes.begin()->first; + int smallsize = smallest.size(); + + for(it2 = nodes.begin(); it2 != nodes.end(); it2++) { + if(it2->first.size() < smallsize) { smallsize = it2->first.size(); smallest = it2->first; } + } + + return smallest; + } + catch(exception& e) { + errorOut(e, "ConcensusCommand", "getSmallest"); + exit(1); + } +} //********************************************************************************************************************** -vector ConcensusCommand::getNextAvailableSet(vector bigset) { +vector ConcensusCommand::getNextAvailableSet(vector bigset, vector& rest) { try { +//cout << "new call " << endl << endl << endl; vector largest; largest.clear(); - int largestSighting = -1; - - //go through the sets - for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) { - //are you a subset of bigset - if (isSubset(bigset, it2->first) == true) { - - //are you the largest. if you are the same size as current largest refer to sighting - if (it2->first.size() > largest.size()) { largest = it2->first; largestSighting = it2->second; } - else if (it2->first.size() == largest.size()) { - if (it2->second > largestSighting) { largest = it2->first; largestSighting = it2->second; } - } - - } + rest.clear(); + + //if you are just 2 groups + if (bigset.size() == 2) { + rest.push_back(bigset[0]); + largest.push_back(bigset[1]); + }else{ + rest = bestSplit[bigset][0]; + largest = bestSplit[bigset][1]; } + + //save for printing out later and for branch lengths + nodePairsInTree[rest] = nodePairs[rest]; + + //delete whatever set you return because it is no longer available + nodePairs.erase(rest); + //save for printing out later and for branch lengths nodePairsInTree[largest] = nodePairs[largest]; @@ -280,6 +311,73 @@ vector ConcensusCommand::getNextAvailableSet(vector bigset) { } } +/**********************************************************************************************************************/ +int ConcensusCommand::getSubgroupRating(vector group) { + try { + map< vector, int>::iterator ittemp; + map< vector< vector > , int >::iterator it3; + int rate = 0; + + // ***********************************************************************************// + //1. this function must be called passing it littlest sets to biggest + // since it the rating is made from your sighting plus you best splits rating + //2. it saves the top pair to use later + // ***********************************************************************************// + + + if (group.size() < 3) { return rate; } + + map< vector, int> possiblePairing; //this is all the subsets of group + + //go through the sets + for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) { + //are you a subset of bigset, then save in possiblePairings + if (isSubset(group, it2->first) == true) { possiblePairing[it2->first] = it2->second; } + } + + map< vector< vector > , int > rating; + + while (possiblePairing.size() != 0) { + + it2 = possiblePairing.begin(); + vector temprest = getRestSet(group, it2->first); + + //is the rest a set available in possiblePairings + ittemp = possiblePairing.find(temprest); + if (ittemp != possiblePairing.end()) { //if the rest is in the possible pairings then add this pair to rating map + vector< vector > temprate; + temprate.push_back(it2->first); temprate.push_back(temprest); + + rating[temprate] = (nodePairsInitialRate[it2->first] + nodePairsInitialRate[temprest]); + + //erase so you dont add 1,2 and 2,1. + possiblePairing.erase(temprest); + } + + possiblePairing.erase(it2); + } + + + it3 = rating.begin(); + rate = it3->second; + vector< vector > topPair = it3->first; + + //choose the split with the best rating + for (it3 = rating.begin(); it3 != rating.end(); it3++) { + + if (it3->second > rate) { rate = it3->second; topPair = it3->first; } + } + + bestSplit[group] = topPair; + + return rate; + } + catch(exception& e) { + errorOut(e, "ConcensusCommand", "getSubgroupRating"); + exit(1); + } +} + //********************************************************************************************************************** vector ConcensusCommand::getRestSet(vector bigset, vector subset) { try { @@ -294,12 +392,6 @@ vector ConcensusCommand::getRestSet(vector bigset, vector ConcensusCommand::getRestSet(vector bigset, vector bigset, vector subset) { try { + + if (subset.size() > bigset.size()) { return false; } + //check if each guy in suset is also in bigset for (int i = 0; i < subset.size(); i++) { bool match = false;