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                 for (int i = 0; i < numNodes; i++) {
  82                         delete f1leaves[i];
  83                         delete f2leaves[i];
  84                 }
  85
  86         //      cout << "exiting " << endl;
  87                 data[0] = Chao;
  88                 return data;
  89         }
  90         catch(exception& e) {
  91                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  92                 exit(1);
  93         }
  94         catch(...) {
  95                 cout << "An unknown error has occurred in the SharedChao1 class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
  96                 exit(1);
  97         }
  98 }
  99
 100 /***********************************************************************/
 101 //builds trees structure with n leaf nodes initialized to 0.
 102 void SharedChao1::initialTree(int n) {
 103         try {
 104                 // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
 105                 numLeaves = pow(2, (float)n) / 2;
 106                 numNodes = 2*numLeaves - 1;
 107                 int countleft = 0;
 108                 int countright = 1;
 109
 110                 f1leaves.resize(numNodes);
 111                 f2leaves.resize(numNodes);
 112
 113                 //initialize leaf values
 114                 for (int i = 0; i < numLeaves; i++) {
 115                         f1leaves[i] = new IntNode;
 116                         f1leaves[i]->lvalue = 0;
 117                         f1leaves[i]->rvalue = 0;
 118                         f1leaves[i]->left = NULL;
 119                         f1leaves[i]->right = NULL;
 120
 121                         f2leaves[i] = new IntNode;
 122                         f2leaves[i]->lvalue = 0;
 123                         f2leaves[i]->rvalue = 0;
 124                         f2leaves[i]->left = NULL;
 125                         f2leaves[i]->right = NULL;
 126                 }
 127
 128                 //set pointers to children
 129                 for (int j = numLeaves; j < numNodes; j++) {
 130                         f1leaves[j] = new IntNode;
 131                         f1leaves[j]->left = f1leaves[countleft];
 132                         f1leaves[j]->right = f1leaves[countright];
 133
 134                         f2leaves[j] = new IntNode;
 135                         f2leaves[j]->left = f2leaves[countleft];
 136                         f2leaves[j]->right =f2leaves[countright];
 137
 138                         countleft = countleft + 2;
 139                         countright = countright + 2;
 140                 }
 141
 142                 //point to root
 143                 f1root = f1leaves[numNodes-1];
 144
 145                 //point to root
 146                 f2root = f2leaves[numNodes-1];
 147
 148                 //set coeffients
 149                 setCoef(f2root, 0);
 150         }
 151         catch(exception& e) {
 152                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 153                 exit(1);
 154         }
 155         catch(...) {
 156                 cout << "An unknown error has occurred in the SharedChao1 class function initialTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 157                 exit(1);
 158         }
 159 }
 160
 161 /***********************************************************************/
 162 //take vector containing the abundance info. for a bin and updates trees.
 163 void SharedChao1::updateTree(vector<int> bin) {
 164         try {
 165                 updateBranchf1(f1root, bin, 0);
 166                 updateBranchf2(f2root, bin, 0);
 167         }
 168         catch(exception& e) {
 169                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 170                 exit(1);
 171         }
 172         catch(...) {
 173                 cout << "An unknown error has occurred in the SharedChao1 class function updateTree. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 174                 exit(1);
 175         }
 176 }
 177
 178 /***********************************************************************/
 179 void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
 180         try {
 181                 //if you have more than one group
 182                 if (index == (bin.size()-1)) {
 183                         if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
 184                         else { node->rvalue++;  }
 185                 }else {
 186                         if (bin[index] == 1) {
 187                                 //follow path as if you are 1
 188                                 updateBranchf1(node->left, bin, index+1);
 189                         }
 190                         //follow path as if you are +
 191                         updateBranchf1(node->right, bin, index+1);
 192                 }
 193         }
 194         catch(exception& e) {
 195                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 196                 exit(1);
 197         }
 198         catch(...) {
 199                 cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf1. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 200                 exit(1);
 201         }
 202 }
 203
 204 /***********************************************************************/
 205 void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
 206         try {
 207                 //if you have more than one group
 208                 if (index == (bin.size()-1)) {
 209                         if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
 210                         else { node->rvalue++;  }
 211                 }else {
 212                         if (bin[index] == 2) {
 213                                 //follow path as if you are 1
 214                                 updateBranchf2(node->left, bin, index+1);
 215                         }
 216                         //follow path as if you are +
 217                         updateBranchf2(node->right, bin, index+1);
 218                 }
 219         }
 220         catch(exception& e) {
 221                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 222                 exit(1);
 223         }
 224         catch(...) {
 225                 cout << "An unknown error has occurred in the SharedChao1 class function updateBranchf2. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 226                 exit(1);
 227         }
 228 }
 229
 230 /***********************************************************************/
 231 void SharedChao1::setCoef(IntNode* node, int coef) {
 232         try {
 233                 if (node->left != NULL) {
 234                         setCoef(node->left, coef+1);
 235                         setCoef(node->right, coef);
 236                 }else {
 237                         node->lcoef = coef+1;
 238                         node->rcoef = coef;
 239                 }
 240         }
 241         catch(exception& e) {
 242                 cout << "Standard Error: " << e.what() << " has occurred in the SharedChao1 class Function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 243                 exit(1);
 244         }
 245         catch(...) {
 246                 cout << "An unknown error has occurred in the SharedChao1 class function getCoef. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 247                 exit(1);
 248         }
 249 }
 250
 251 /***********************************************************************/
 252 //for debugging purposes
 253 void SharedChao1::printTree() {
 254
 255         cout << "F1 leaves" << endl;
 256         printBranch(f1root);
 257
 258         cout << "F2 leaves" << endl;
 259         printBranch(f2root);
 260
 261
 262 }
 263 /*****************************************************************/
 264 void SharedChao1::printBranch(IntNode* node) {
 265         try {
 266
 267                 // you are not a leaf
 268                 if (node->left != NULL) {
 269                         printBranch(node->left);
 270                         printBranch(node->right);
 271                 }else { //you are a leaf
 272                         cout << node->lvalue << endl;
 273                         cout << node->rvalue << endl;
 274                 }
 275
 276         }
 277         catch(exception& e) {
 278                 cout << "Standard Error: " << e.what() << " has occurred in the Tree class Function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 279                 exit(1);
 280         }
 281         catch(...) {
 282                 cout << "An unknown error has occurred in the Tree class function printBranch. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 283                 exit(1);
 284         }
 285 }
 286
 287 /*****************************************************************/
 288
 289
 290
 291