5 * Created by Sarah Westcott on 4/29/09.
6 * Copyright 2009 Schloss Lab UMASS AMherst. All rights reserved.
10 #include "consensus.h"
12 //**********************************************************************************************************************
13 Tree* Consensus::getTree(vector<Tree*>& t){
15 numNodes = t[0]->getNumNodes();
16 numLeaves = t[0]->getNumLeaves();
19 //get the possible pairings
22 if (m->control_pressed) { return 0; }
24 consensusTree = new Tree(t[0]->getCountTable());
26 it2 = nodePairs.find(treeSet);
28 nodePairsInTree[treeSet] = it2->second;
30 //erase treeset because you are adding it
31 nodePairs.erase(treeSet);
33 //set count to numLeaves;
36 buildConsensusTree(treeSet);
38 if (m->control_pressed) { delete consensusTree; return 0; }
40 consensusTree->assembleTree();
42 if (m->control_pressed) { delete consensusTree; return 0; }
49 m->errorOut(e, "Consensus", "execute");
53 //**********************************************************************************************************************
54 int Consensus::printSetsInfo() {
56 //open file for pairing not included in the tree
57 string notIncluded = "cons.pairs";
59 m->openOutputFile(notIncluded, out2);
61 //output species in order
62 out2 << "Species in Order: " << endl << endl;
63 for (int i = 0; i < treeSet.size(); i++) { out2 << i+1 << ". " << treeSet[i] << endl; }
65 //output sets included
66 out2 << endl << "Sets included in the consensus tree:" << endl << endl;
68 if (m->control_pressed) { return 0; }
71 for (it2 = nodePairsInTree.begin(); it2 != nodePairsInTree.end(); it2++) {
73 if (m->control_pressed) { return 0; }
75 //only output pairs not leaves
76 if (it2->first.size() > 1) {
78 //initialize temp to all "."
79 temp.resize(treeSet.size(), ".");
81 //set the spot in temp that represents it2->first[i] to a "*"
82 for (int i = 0; i < it2->first.size(); i++) {
84 int index = findSpot(it2->first[i]);
86 //temp[index] = it2->first[i] + " ";
90 for (int j = 0; j < temp.size(); j++) {
93 out2 << '\t' << it2->second << endl;
97 //output sets not included
98 out2 << endl << "Sets NOT included in the consensus tree:" << endl << endl;
99 for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) {
101 if (m->control_pressed) { return 0; }
104 //initialize temp to all "."
105 temp.resize(treeSet.size(), ".");
107 //set the spot in temp that represents it2->first[i] to a "*"
108 for (int i = 0; i < it2->first.size(); i++) {
110 int index = findSpot(it2->first[i]);
115 for (int j = 0; j < temp.size(); j++) {
118 out2 << '\t' << it2->second << endl;
123 catch(exception& e) {
124 m->errorOut(e, "Consensus", "printSetsInfo");
128 //**********************************************************************************************************************
129 int Consensus::buildConsensusTree(vector<string> nodeSet) {
131 vector<string> leftChildSet;
132 vector<string> rightChildSet;
134 if (m->control_pressed) { return 1; }
136 //if you are at a leaf
137 if (nodeSet.size() == 1) {
138 //return the vector index of the leaf you are at
139 return consensusTree->getIndex(nodeSet[0]);
140 //terminate recursion
141 }else if (count == numNodes) { return 0; }
143 //finds best child pair
144 leftChildSet = getNextAvailableSet(nodeSet, rightChildSet);
145 int left = buildConsensusTree(leftChildSet);
146 int right = buildConsensusTree(rightChildSet);
147 consensusTree->tree[count].setChildren(left, right);
148 consensusTree->tree[count].setLabel((nodePairsInTree[nodeSet]/(float)numTrees));
149 consensusTree->tree[left].setParent(count);
150 consensusTree->tree[right].setParent(count);
156 catch(exception& e) {
157 m->errorOut(e, "Consensus", "buildConcensusTree");
162 //**********************************************************************************************************************
163 int Consensus::getSets(vector<Tree*>& t) {
168 //for each tree add the possible pairs you find
169 for (int i = 0; i < t.size(); i++) {
171 //for each non-leaf node get descendant info.
172 for (int j = numLeaves; j < numNodes; j++) {
174 if (m->control_pressed) { return 1; }
177 //go through pcounts and pull out descendants
178 for (it = t[i]->tree[j].pcount.begin(); it != t[i]->tree[j].pcount.end(); it++) {
179 temp.push_back(it->first);
183 sort(temp.begin(), temp.end());
185 it2 = nodePairs.find(temp);
186 if (it2 != nodePairs.end()) {
195 //add each leaf to terminate recursion in consensus
196 //you want the leaves in there but with insignifigant sightings value so it is added last
197 //for each leaf node get descendant info.
198 for (int j = 0; j < numLeaves; j++) {
200 if (m->control_pressed) { return 1; }
202 //only need the first one since leaves have no descendants but themselves
203 it = t[0]->tree[j].pcount.begin();
204 temp.clear(); temp.push_back(it->first);
207 treeSet.push_back(it->first);
209 //add leaf to list but with sighting value less then all non leaf pairs
213 sort(treeSet.begin(), treeSet.end());
216 map< vector<string>, int> nodePairsCopy = nodePairs;
219 //set initial rating on pairs to sightings + subgroup sightings
220 while (nodePairsCopy.size() != 0) {
221 if (m->control_pressed) { return 1; }
223 vector<string> smallOne = getSmallest(nodePairsCopy);
225 int subgrouprate = getSubgroupRating(smallOne);
227 nodePairsInitialRate[smallOne] = nodePairs[smallOne] + subgrouprate;
229 nodePairsCopy.erase(smallOne);
234 catch(exception& e) {
235 m->errorOut(e, "Consensus", "getSets");
239 //**********************************************************************************************************************
240 vector<string> Consensus::getSmallest(map< vector<string>, int> nodes) {
242 vector<string> smallest = nodes.begin()->first;
243 int smallsize = smallest.size();
245 for(it2 = nodes.begin(); it2 != nodes.end(); it2++) {
246 if(it2->first.size() < smallsize) { smallsize = it2->first.size(); smallest = it2->first; }
251 catch(exception& e) {
252 m->errorOut(e, "Consensus", "getSmallest");
257 //**********************************************************************************************************************
258 vector<string> Consensus::getNextAvailableSet(vector<string> bigset, vector<string>& rest) {
260 //cout << "new call " << endl << endl << endl;
261 vector<string> largest; largest.clear();
264 //if you are just 2 groups
265 if (bigset.size() == 2) {
266 rest.push_back(bigset[0]);
267 largest.push_back(bigset[1]);
269 rest = bestSplit[bigset][0];
270 largest = bestSplit[bigset][1];
274 //save for printing out later and for branch lengths
275 nodePairsInTree[rest] = nodePairs[rest];
277 //delete whatever set you return because it is no longer available
278 nodePairs.erase(rest);
280 //save for printing out later and for branch lengths
281 nodePairsInTree[largest] = nodePairs[largest];
283 //delete whatever set you return because it is no longer available
284 nodePairs.erase(largest);
289 catch(exception& e) {
290 m->errorOut(e, "Consensus", "getNextAvailableSet");
295 /**********************************************************************************************************************/
296 int Consensus::getSubgroupRating(vector<string> group) {
298 map< vector<string>, int>::iterator ittemp;
299 map< vector< vector<string> > , int >::iterator it3;
302 // ***********************************************************************************//
303 //1. this function must be called passing it littlest sets to biggest
304 // since it the rating is made from your sighting plus you best splits rating
305 //2. it saves the top pair to use later
306 // ***********************************************************************************//
309 if (group.size() < 3) { return rate; }
311 map< vector<string>, int> possiblePairing; //this is all the subsets of group
313 //go through the sets
314 for (it2 = nodePairs.begin(); it2 != nodePairs.end(); it2++) {
315 //are you a subset of bigset, then save in possiblePairings
316 if (isSubset(group, it2->first) == true) { possiblePairing[it2->first] = it2->second; }
319 map< vector< vector<string> > , int > rating;
321 while (possiblePairing.size() != 0) {
323 it2 = possiblePairing.begin();
324 vector<string> temprest = getRestSet(group, it2->first);
326 //is the rest a set available in possiblePairings
327 ittemp = possiblePairing.find(temprest);
328 if (ittemp != possiblePairing.end()) { //if the rest is in the possible pairings then add this pair to rating map
329 vector< vector<string> > temprate;
330 temprate.push_back(it2->first); temprate.push_back(temprest);
332 rating[temprate] = (nodePairsInitialRate[it2->first] + nodePairsInitialRate[temprest]);
334 //erase so you dont add 1,2 and 2,1.
335 possiblePairing.erase(temprest);
338 possiblePairing.erase(it2);
342 it3 = rating.begin();
344 vector< vector<string> > topPair = it3->first;
346 //choose the split with the best rating
347 for (it3 = rating.begin(); it3 != rating.end(); it3++) {
349 if (it3->second > rate) { rate = it3->second; topPair = it3->first; }
352 bestSplit[group] = topPair;
356 catch(exception& e) {
357 m->errorOut(e, "Consensus", "getSubgroupRating");
362 //**********************************************************************************************************************
363 vector<string> Consensus::getRestSet(vector<string> bigset, vector<string> subset) {
367 for (int i = 0; i < bigset.size(); i++) {
368 bool inSubset = false;
369 for (int j = 0; j < subset.size(); j++) {
370 if (bigset[i] == subset[j]) { inSubset = true; break; }
373 //its not in the subset so put it in the rest
374 if (inSubset == false) { rest.push_back(bigset[i]); }
380 catch(exception& e) {
381 m->errorOut(e, "Consensus", "getRestSet");
386 //**********************************************************************************************************************
387 bool Consensus::isSubset(vector<string> bigset, vector<string> subset) {
391 if (subset.size() > bigset.size()) { return false; }
393 //check if each guy in suset is also in bigset
394 for (int i = 0; i < subset.size(); i++) {
396 for (int j = 0; j < bigset.size(); j++) {
397 if (subset[i] == bigset[j]) { match = true; break; }
400 //you have a guy in subset that had no match in bigset
401 if (match == false) { return false; }
407 catch(exception& e) {
408 m->errorOut(e, "Consensus", "isSubset");
412 //**********************************************************************************************************************
413 int Consensus::findSpot(string node) {
417 //check if each guy in suset is also in bigset
418 for (int i = 0; i < treeSet.size(); i++) {
419 if (treeSet[i] == node) { spot = i; break; }
425 catch(exception& e) {
426 m->errorOut(e, "Consensus", "findSpot");
430 //**********************************************************************************************************************