5 * Created by Sarah Westcott on 1/8/09.
6 * Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
10 #include "sharedchao1.h"
12 /***********************************************************************/
13 EstOutput SharedChao1::getValues(vector<SharedRAbundVector*> shared){
17 int numGroups = shared.size();
18 float Chao = 0.0; float leftvalue, rightvalue;
20 // IntNode is defined in mothur.h
21 // The tree used here is a binary tree used to represent the f1+++, f+1++, f++1+, f+++1, f11++, f1+1+...
22 // combinations required to solve the chao estimator equation for any number of groups. Conceptually, think
23 // of each node as having a 1 and a + value, or for f2 values a 2 and a + value, and 2 pointers to intnodes, and 2 coeffient values.
24 // The coeffient value is how many times you chose branch 1 to get to that fvalue.
25 // If you choose left you are selecting the 1 or 2 value and right means the + value. For instance, to find
26 // the number of bins that have f1+1+ you would start at the root, go left, right, left, and select the rightvalue.
27 // the coeffient is 2. Note: we only set the coeffient in f2 values.
29 //create and initialize trees to 0.
30 initialTree(numGroups);
32 //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values
33 for (int i = 0; i < shared[0]->size(); i++) {
34 //get bin values and calc shared
35 bool sharedByAll = true;
37 for (int j = 0; j < numGroups; j++) {
38 temp.push_back(shared[j]->getAbundance(i));
39 if (temp[j] == 0) { sharedByAll = false; }
43 if (sharedByAll == true) {
44 //find f1 and f2values
50 //calculate chao1, (numleaves-1) because numleaves contains the ++ values.
52 for(int i=0;i<numLeaves;i++){
53 if (f2leaves[i]->lvalue == 0) { bias = true;}// break;}
57 for (int i = 0; i < numLeaves; i++) {
59 leftvalue = (float)(f1leaves[i]->lvalue * (f1leaves[i]->lvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * (f2leaves[i]->lvalue + 1));
60 if (i != (numLeaves-1)) {
61 rightvalue = (float)(f1leaves[i]->rvalue * (f1leaves[i]->rvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * (f2leaves[i]->rvalue + 1));
63 rightvalue = (float)(f1leaves[i]->rvalue);
65 Chao += leftvalue + rightvalue;
70 for (int i = 0; i < numLeaves; i++) {
72 leftvalue = (float)(f1leaves[i]->lvalue * f1leaves[i]->lvalue) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * f2leaves[i]->lvalue);
73 if (i != (numLeaves-1)) {
74 rightvalue = (float)(f1leaves[i]->rvalue * f1leaves[i]->rvalue) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * f2leaves[i]->rvalue);
76 rightvalue = (float)(f1leaves[i]->rvalue);
78 Chao += leftvalue + rightvalue;
82 for (int i = 0; i < numNodes; i++) {
92 errorOut(e, "SharedChao1", "getValues");
97 /***********************************************************************/
98 //builds trees structure with n leaf nodes initialized to 0.
99 void SharedChao1::initialTree(int n) {
101 // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
102 numLeaves = pow(2, (float)n) / 2;
103 numNodes = 2*numLeaves - 1;
107 f1leaves.resize(numNodes);
108 f2leaves.resize(numNodes);
110 //initialize leaf values
111 for (int i = 0; i < numLeaves; i++) {
112 f1leaves[i] = new IntNode;
113 f1leaves[i]->lvalue = 0;
114 f1leaves[i]->rvalue = 0;
115 f1leaves[i]->left = NULL;
116 f1leaves[i]->right = NULL;
118 f2leaves[i] = new IntNode;
119 f2leaves[i]->lvalue = 0;
120 f2leaves[i]->rvalue = 0;
121 f2leaves[i]->left = NULL;
122 f2leaves[i]->right = NULL;
125 //set pointers to children
126 for (int j = numLeaves; j < numNodes; j++) {
127 f1leaves[j] = new IntNode;
128 f1leaves[j]->left = f1leaves[countleft];
129 f1leaves[j]->right = f1leaves[countright];
131 f2leaves[j] = new IntNode;
132 f2leaves[j]->left = f2leaves[countleft];
133 f2leaves[j]->right =f2leaves[countright];
135 countleft = countleft + 2;
136 countright = countright + 2;
140 f1root = f1leaves[numNodes-1];
143 f2root = f2leaves[numNodes-1];
148 catch(exception& e) {
149 errorOut(e, "SharedChao1", "initialTree");
154 /***********************************************************************/
155 //take vector containing the abundance info. for a bin and updates trees.
156 void SharedChao1::updateTree(vector<int> bin) {
158 updateBranchf1(f1root, bin, 0);
159 updateBranchf2(f2root, bin, 0);
161 catch(exception& e) {
162 errorOut(e, "SharedChao1", "updateTree");
167 /***********************************************************************/
168 void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
170 //if you have more than one group
171 if (index == (bin.size()-1)) {
172 if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
173 else { node->rvalue++; }
175 if (bin[index] == 1) {
176 //follow path as if you are 1
177 updateBranchf1(node->left, bin, index+1);
179 //follow path as if you are +
180 updateBranchf1(node->right, bin, index+1);
183 catch(exception& e) {
184 errorOut(e, "SharedChao1", "updateBranchf1");
189 /***********************************************************************/
190 void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
192 //if you have more than one group
193 if (index == (bin.size()-1)) {
194 if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
195 else { node->rvalue++; }
197 if (bin[index] == 2) {
198 //follow path as if you are 1
199 updateBranchf2(node->left, bin, index+1);
201 //follow path as if you are +
202 updateBranchf2(node->right, bin, index+1);
205 catch(exception& e) {
206 errorOut(e, "SharedChao1", "updateBranchf2");
211 /***********************************************************************/
212 void SharedChao1::setCoef(IntNode* node, int coef) {
214 if (node->left != NULL) {
215 setCoef(node->left, coef+1);
216 setCoef(node->right, coef);
218 node->lcoef = coef+1;
222 catch(exception& e) {
223 errorOut(e, "SharedChao1", "setCoef");
228 /***********************************************************************/
229 //for debugging purposes
230 void SharedChao1::printTree() {
232 mothurOut("F1 leaves"); mothurOutEndLine();
235 mothurOut("F2 leaves"); mothurOutEndLine();
240 /*****************************************************************/
241 void SharedChao1::printBranch(IntNode* node) {
244 // you are not a leaf
245 if (node->left != NULL) {
246 printBranch(node->left);
247 printBranch(node->right);
248 }else { //you are a leaf
249 mothurOut(toString(node->lvalue)); mothurOutEndLine();
250 mothurOut(toString(node->rvalue)); mothurOutEndLine();
254 catch(exception& e) {
255 errorOut(e, "SharedChao1", "printBranch");
260 /*****************************************************************/