sharedchao1.cpp

   1 /*
   2  *  sharedchao1.cpp
   3  *  Dotur
   4  *
   5  *  Created by Sarah Westcott on 1/8/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "sharedchao1.h"
  11
  12 /***********************************************************************/
  13 EstOutput SharedChao1::getValues(vector<SharedRAbundVector*> shared){
  14         try {
  15                 data.resize(1,0);
  16                 vector<int> temp;
  17                 int numGroups = shared.size();
  18                 float Chao = 0.0; float leftvalue, rightvalue;
  19
  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.
  28
  29                 //create and initialize trees to 0.
  30                 initialTree(numGroups);
  31
  32                 for (int i = 0; i < shared[0]->getNumBins(); i++) {
  33                         //get bin values and calc shared
  34                         bool sharedByAll = true;
  35                         temp.clear();
  36                         for (int j = 0; j < numGroups; j++) {
  37                                 temp.push_back(shared[j]->getAbundance(i));
  38                                 if (temp[j] == 0) { sharedByAll = false; }
  39                         }
  40
  41                         //they are shared
  42                         if (sharedByAll == true) {
  43                                 //find f1 and f2values
  44                                 updateTree(temp);
  45                         }
  46                 }
  47
  48
  49                 //calculate chao1, (numleaves-1) because numleaves contains the ++ values.
  50                 bool bias = false;
  51                 for(int i=0;i<numLeaves;i++){
  52                         if (f2leaves[i]->lvalue == 0 || f2leaves[i]->rvalue == 0) { bias = true;}// break;}
  53                 }
  54
  55                 if(bias){
  56                         for (int i = 0; i < numLeaves; i++) {
  57
  58                                 leftvalue = (float)(f1leaves[i]->lvalue * (f1leaves[i]->lvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * (f2leaves[i]->lvalue + 1));
  59                                 if (i != (numLeaves-1)) {
  60                                         rightvalue = (float)(f1leaves[i]->rvalue * (f1leaves[i]->rvalue - 1)) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * (f2leaves[i]->rvalue + 1));
  61                                 }else{
  62                                         //add in sobs
  63                                         rightvalue = (float)(f1leaves[i]->rvalue);
  64                                 }
  65                                 Chao += leftvalue + rightvalue;
  66                         }
  67                 }
  68                 else{
  69
  70                         for (int i = 0; i < numLeaves; i++) {
  71
  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);
  75                                 }else{
  76                                         //add in sobs
  77                                         rightvalue = (float)(f1leaves[i]->rvalue);
  78                                 }
  79                                 Chao += leftvalue + rightvalue;
  80                         }
  81                 }
  82
  83                 for (int i = 0; i < numNodes; i++) {
  84                         delete f1leaves[i];
  85                         delete f2leaves[i];
  86                 }
  87
  88
  89                 data[0] = Chao;
  90                 return data;
  91         }
  92         catch(exception& e) {
  93                 m->errorOut(e, "SharedChao1", "getValues");
  94                 exit(1);
  95         }
  96 }
  97
  98 /***********************************************************************/
  99 //builds trees structure with n leaf nodes initialized to 0.
 100 void SharedChao1::initialTree(int n) {
 101         try {
 102                 // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
 103                 numLeaves = pow(2, (float)n) / 2;
 104                 numNodes = 2*numLeaves - 1;
 105                 int countleft = 0;
 106                 int countright = 1;
 107
 108                 f1leaves.resize(numNodes);
 109                 f2leaves.resize(numNodes);
 110
 111                 //initialize leaf values
 112                 for (int i = 0; i < numLeaves; i++) {
 113                         f1leaves[i] = new IntNode(0, 0, NULL, NULL);
 114                         f2leaves[i] = new IntNode(0, 0, NULL, NULL);
 115                 }
 116
 117                 //set pointers to children
 118                 for (int j = numLeaves; j < numNodes; j++) {
 119                         f1leaves[j] = new IntNode();
 120                         f1leaves[j]->left = f1leaves[countleft];
 121                         f1leaves[j]->right = f1leaves[countright];
 122
 123                         f2leaves[j] = new IntNode();
 124                         f2leaves[j]->left = f2leaves[countleft];
 125                         f2leaves[j]->right =f2leaves[countright];
 126
 127                         countleft = countleft + 2;
 128                         countright = countright + 2;
 129                 }
 130
 131                 //point to root
 132                 f1root = f1leaves[numNodes-1];
 133
 134                 //point to root
 135                 f2root = f2leaves[numNodes-1];
 136
 137                 //set coeffients
 138                 setCoef(f2root, 0);
 139         }
 140         catch(exception& e) {
 141                 if ((toString(e.what()) == "vector::_M_fill_insert") || (toString(e.what()) == "St9bad_alloc")) { m->mothurOut("You are using " + toString(n) + " groups which creates 2^" + toString(n+1) + " nodes. Try reducing the number of groups you selected. "); m->mothurOutEndLine(); exit(1); }
 142                 m->errorOut(e, "SharedChao1", "initialTree");
 143                 exit(1);
 144         }
 145 }
 146
 147 /***********************************************************************/
 148 //take vector containing the abundance info. for a bin and updates trees.
 149 void SharedChao1::updateTree(vector<int> bin) {
 150         try {
 151                 updateBranchf1(f1root, bin, 0);
 152                 updateBranchf2(f2root, bin, 0);
 153         }
 154         catch(exception& e) {
 155                 m->errorOut(e, "SharedChao1", "updateTree");
 156                 exit(1);
 157         }
 158 }
 159
 160 /***********************************************************************/
 161 void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
 162         try {
 163                 //if you have more than one group
 164                 if (index == (bin.size()-1)) {
 165                         if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
 166                         else { node->rvalue++;  }
 167                 }else {
 168                         if (bin[index] == 1) {
 169                                 //follow path as if you are 1
 170                                 updateBranchf1(node->left, bin, index+1);
 171                         }
 172                         //follow path as if you are +
 173                         updateBranchf1(node->right, bin, index+1);
 174                 }
 175         }
 176         catch(exception& e) {
 177                 m->errorOut(e, "SharedChao1", "updateBranchf1");
 178                 exit(1);
 179         }
 180 }
 181
 182 /***********************************************************************/
 183 void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
 184         try {
 185                 //if you have more than one group
 186                 if (index == (bin.size()-1)) {
 187                         if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
 188                         else { node->rvalue++;  }
 189                 }else {
 190                         if (bin[index] == 2) {
 191                                 //follow path as if you are 1
 192                                 updateBranchf2(node->left, bin, index+1);
 193                         }
 194                         //follow path as if you are +
 195                         updateBranchf2(node->right, bin, index+1);
 196                 }
 197         }
 198         catch(exception& e) {
 199                 m->errorOut(e, "SharedChao1", "updateBranchf2");
 200                 exit(1);
 201         }
 202 }
 203
 204 /***********************************************************************/
 205 void SharedChao1::setCoef(IntNode* node, int coef) {
 206         try {
 207                 if (node->left != NULL) {
 208                         setCoef(node->left, coef+1);
 209                         setCoef(node->right, coef);
 210                 }else {
 211                         node->lcoef = coef+1;
 212                         node->rcoef = coef;
 213                 }
 214         }
 215         catch(exception& e) {
 216                 m->errorOut(e, "SharedChao1", "setCoef");
 217                 exit(1);
 218         }
 219 }
 220
 221 /***********************************************************************/
 222 //for debugging purposes
 223 void SharedChao1::printTree() {
 224
 225         m->mothurOut("F1 leaves"); m->mothurOutEndLine();
 226         printBranch(f1root);
 227
 228         m->mothurOut("F2 leaves"); m->mothurOutEndLine();
 229         printBranch(f2root);
 230
 231
 232 }
 233 /*****************************************************************/
 234 void SharedChao1::printBranch(IntNode* node) {
 235         try {
 236
 237                 // you are not a leaf
 238                 if (node->left != NULL) {
 239                         printBranch(node->left);
 240                         printBranch(node->right);
 241                 }else { //you are a leaf
 242                         m->mothurOut(toString(node->lvalue)); m->mothurOutEndLine();
 243                         m->mothurOut(toString(node->rvalue)); m->mothurOutEndLine();
 244                 }
 245
 246         }
 247         catch(exception& e) {
 248                 m->errorOut(e, "SharedChao1", "printBranch");
 249                 exit(1);
 250         }
 251 }
 252
 253 /*****************************************************************/
 254
 255
 256
 257