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
  17                 vector<int> temp;
  18                 int numGroups = shared.size();
  19                 float Chao = 0.0; float leftvalue, rightvalue;
  20
  21                 // IntNode is defined in mothur.h
  22                 // The tree used here is a binary tree used to represent the f1+++, f+1++, f++1+, f+++1, f11++, f1+1+...
  23                 // combinations required to solve the chao estimator equation for any number of groups.  Conceptually, think
  24                 // 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.
  25                 // The coeffient value is how many times you chose branch 1 to get to that fvalue.
  26                 // If you choose left you are selecting the 1 or 2 value and right means the + value.  For instance, to find
  27                 // the number of bins that have f1+1+ you would start at the root, go left, right, left, and select the rightvalue.
  28                 // the coeffient is 2.  Note: we only set the coeffient in f2 values.
  29
  30                 //create and initialize trees to 0.
  31                 initialTree(numGroups);
  32
  33                 //loop through vectors calculating the f11, f1A, f2A, f1B, f2B, S12 values
  34                 for (int i = 0; i < shared[0]->size(); i++) {
  35                         //get bin values and calc shared
  36                         bool sharedByAll = true;
  37                         temp.clear();
  38                         for (int j = 0; j < numGroups; j++) {
  39                                 temp.push_back(shared[j]->getAbundance(i));
  40                                 if (temp[j] == 0) { sharedByAll = false; }
  41                         }
  42
  43                         //they are shared
  44                         if (sharedByAll == true) {
  45
  46                                 // cout << "temp = ";
  47                                 // for (int h = 0; h < temp.size(); h++) { cout << temp[h] << " "; } cout << endl;
  48                                 //find f1 and f2values
  49                                 updateTree(temp);
  50                         }
  51                 }
  52
  53
  54                 //cout << "Entering " << endl;
  55                 //calculate chao1, (numleaves-1) because numleaves contains the ++ values.
  56                 for (int i = 0; i < numLeaves; i++) {
  57                         //divide by zero error
  58                         if (f2leaves[i]->lvalue == 0) { f2leaves[i]->lvalue++; }
  59                         if (f2leaves[i]->rvalue == 0) { f2leaves[i]->rvalue++; }
  60
  61                         //cout << "f2 leaves ";
  62                         //cout << f2leaves[i]->lvalue << f2leaves[i]->rvalue << endl;
  63
  64                 //      cout << "f2 leaf coef ";
  65                         //cout << f2leaves[i]->lcoef << '\t' << f2leaves[i]->rcoef << endl;
  66
  67                 //      cout << "f1 leaves ";
  68                 //      cout << f1leaves[i]->lvalue << f1leaves[i]->rvalue << endl;
  69
  70
  71                         leftvalue = (float)(f1leaves[i]->lvalue * f1leaves[i]->lvalue) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * f2leaves[i]->lvalue);
  72                         if (i != (numLeaves-1)) {
  73                                 rightvalue = (float)(f1leaves[i]->rvalue * f1leaves[i]->rvalue) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * f2leaves[i]->rvalue);
  74                         }else{
  75                                 rightvalue = (float)(f1leaves[i]->rvalue);
  76                         }
  77                         //cout << "left = " << leftvalue << " right = " << rightvalue << endl;
  78                         Chao += leftvalue + rightvalue;
  79                 }
  80
  81         //      cout << "exiting " << endl;
  82                 data[0] = Chao;
  83                 return data;
  84         }
  85         catch(exception& e) {
  86                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  87                 exit(1);
  88         }
  89         catch(...) {
  90                 cout << "An unknown error has occurred in the SharedChao1 class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  91                 exit(1);
  92         }
  93 }
  94
  95 /***********************************************************************/
  96 //builds trees structure with n leaf nodes initialized to 0.
  97 void SharedChao1::initialTree(int n) {
  98         try {
  99                 // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
 100                 numLeaves = pow(2, (float)n) / 2;
 101                 numNodes = 2*numLeaves - 1;
 102                 int countleft = 0;
 103                 int countright = 1;
 104
 105                 f1leaves.resize(numNodes);
 106                 f2leaves.resize(numNodes);
 107
 108                 //initialize leaf values
 109                 for (int i = 0; i < numLeaves; i++) {
 110                         f1leaves[i] = new IntNode;
 111                         f1leaves[i]->lvalue = 0;
 112                         f1leaves[i]->rvalue = 0;
 113                         f1leaves[i]->left = NULL;
 114                         f1leaves[i]->right = NULL;
 115
 116                         f2leaves[i] = new IntNode;
 117                         f2leaves[i]->lvalue = 0;
 118                         f2leaves[i]->rvalue = 0;
 119                         f2leaves[i]->left = NULL;
 120                         f2leaves[i]->right = NULL;
 121                 }
 122
 123                 //set pointers to children
 124                 for (int j = numLeaves; j < numNodes; j++) {
 125                         f1leaves[j] = new IntNode;
 126                         f1leaves[j]->left = f1leaves[countleft];
 127                         f1leaves[j]->right = f1leaves[countright];
 128
 129                         f2leaves[j] = new IntNode;
 130                         f2leaves[j]->left = f2leaves[countleft];
 131                         f2leaves[j]->right =f2leaves[countright];
 132
 133                         countleft = countleft + 2;
 134                         countright = countright + 2;
 135                 }
 136
 137                 //point to root
 138                 f1root = f1leaves[numNodes-1];
 139
 140                 //point to root
 141                 f2root = f2leaves[numNodes-1];
 142
 143                 //set coeffients
 144                 setCoef(f2root, 0);
 145         }
 146         catch(exception& e) {
 147                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 148                 exit(1);
 149         }
 150         catch(...) {
 151                 cout << "An unknown error has occurred in the SharedChao1 class function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 152                 exit(1);
 153         }
 154 }
 155
 156 /***********************************************************************/
 157 //take vector containing the abundance info. for a bin and updates trees.
 158 void SharedChao1::updateTree(vector<int> bin) {
 159         try {
 160                 updateBranchf1(f1root, bin, 0);
 161                 updateBranchf2(f2root, bin, 0);
 162         }
 163         catch(exception& e) {
 164                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 165                 exit(1);
 166         }
 167         catch(...) {
 168                 cout << "An unknown error has occurred in the SharedChao1 class function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 169                 exit(1);
 170         }
 171 }
 172
 173 /***********************************************************************/
 174 void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
 175         try {
 176                 //if you have more than one group
 177                 if (index == (bin.size()-1)) {
 178                         if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
 179                         else { node->rvalue++;  }
 180                 }else {
 181                         if (bin[index] == 1) {
 182                                 //follow path as if you are 1
 183                                 updateBranchf1(node->left, bin, index+1);
 184                         }
 185                         //follow path as if you are +
 186                         updateBranchf1(node->right, bin, index+1);
 187                 }
 188         }
 189         catch(exception& e) {
 190                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 191                 exit(1);
 192         }
 193         catch(...) {
 194                 cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 195                 exit(1);
 196         }
 197 }
 198
 199 /***********************************************************************/
 200 void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
 201         try {
 202                 //if you have more than one group
 203                 if (index == (bin.size()-1)) {
 204                         if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
 205                         else { node->rvalue++;  }
 206                 }else {
 207                         if (bin[index] == 2) {
 208                                 //follow path as if you are 1
 209                                 updateBranchf2(node->left, bin, index+1);
 210                         }
 211                         //follow path as if you are +
 212                         updateBranchf2(node->right, bin, index+1);
 213                 }
 214         }
 215         catch(exception& e) {
 216                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 217                 exit(1);
 218         }
 219         catch(...) {
 220                 cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 221                 exit(1);
 222         }
 223 }
 224
 225 /***********************************************************************/
 226 void SharedChao1::setCoef(IntNode* node, int coef) {
 227         try {
 228                 if (node->left != NULL) {
 229                         setCoef(node->left, coef+1);
 230                         setCoef(node->right, coef);
 231                 }else {
 232                         node->lcoef = coef+1;
 233                         node->rcoef = coef;
 234                 }
 235         }
 236         catch(exception& e) {
 237                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 238                 exit(1);
 239         }
 240         catch(...) {
 241                 cout << "An unknown error has occurred in the SharedChao1 class function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 242                 exit(1);
 243         }
 244 }
 245
 246 /***********************************************************************/
 247 //for debugging purposes
 248 void SharedChao1::printTree() {
 249
 250         cout << "F1 leaves" << endl;
 251         printBranch(f1root);
 252
 253         cout << "F2 leaves" << endl;
 254         printBranch(f2root);
 255
 256
 257 }
 258 /*****************************************************************/
 259 void SharedChao1::printBranch(IntNode* node) {
 260         try {
 261
 262                 // you are not a leaf
 263                 if (node->left != NULL) {
 264                         printBranch(node->left);
 265                         printBranch(node->right);
 266                 }else { //you are a leaf
 267                         cout << node->lvalue << endl;
 268                         cout << node->rvalue << endl;
 269                 }
 270
 271         }
 272         catch(exception& e) {
 273                 cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 274                 exit(1);
 275         }
 276         catch(...) {
 277                 cout << "An unknown error has occurred in the Tree class function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 278                 exit(1);
 279         }
 280 }
 281
 282 /*****************************************************************/
 283
 284
 285
 286