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
  33                 for (int i = 0; i < shared[0]->size(); i++) {
  34                         //get bin values and calc shared
  35                         bool sharedByAll = true;
  36                         temp.clear();
  37                         for (int j = 0; j < numGroups; j++) {
  38                                 temp.push_back(shared[j]->getAbundance(i));
  39                                 if (temp[j] == 0) { sharedByAll = false; }
  40                         }
  41
  42                         //they are shared
  43                         if (sharedByAll == true) {
  44                                 //find f1 and f2values
  45                                 updateTree(temp);
  46                         }
  47                 }
  48
  49
  50                 //calculate chao1, (numleaves-1) because numleaves contains the ++ values.
  51                 bool bias;
  52                 for(int i=0;i<numLeaves;i++){
  53                         if ((f2leaves[i]->lvalue == 0) || (f2leaves[i]->rvalue == 0))  { bias = true;  }// break;}
  54                 }
  55
  56                 if(bias){
  57                         for (int i = 0; i < numLeaves; i++) {
  58
  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));
  62                                 }else{
  63                                         //add in sobs
  64                                         rightvalue = (float)(f1leaves[i]->rvalue);
  65                                 }
  66                                 Chao += leftvalue + rightvalue;
  67                         }
  68                 }
  69                 else{
  70
  71                         for (int i = 0; i < numLeaves; i++) {
  72
  73                                 leftvalue = (float)(f1leaves[i]->lvalue * f1leaves[i]->lvalue) / (float)((pow(2, (float)f2leaves[i]->lcoef)) * f2leaves[i]->lvalue);
  74                                 if (i != (numLeaves-1)) {
  75                                         rightvalue = (float)(f1leaves[i]->rvalue * f1leaves[i]->rvalue) / (float)((pow(2, (float)f2leaves[i]->rcoef)) * f2leaves[i]->rvalue);
  76                                 }else{
  77                                         //add in sobs
  78                                         rightvalue = (float)(f1leaves[i]->rvalue);
  79                                 }
  80                                 Chao += leftvalue + rightvalue;
  81                         }
  82                 }
  83
  84                 for (int i = 0; i < numNodes; i++) {
  85                         delete f1leaves[i];
  86                         delete f2leaves[i];
  87                 }
  88
  89
  90                 data[0] = Chao;
  91                 return data;
  92         }
  93         catch(exception& e) {
  94                 errorOut(e, "SharedChao1", "getValues");
  95                 exit(1);
  96         }
  97 }
  98
  99 /***********************************************************************/
 100 //builds trees structure with n leaf nodes initialized to 0.
 101 void SharedChao1::initialTree(int n) {
 102         try {
 103                 // (2^n) / 2. Divide by 2 because each leaf node contains 2 values. One for + and one for 1 or 2.
 104                 numLeaves = pow(2, (float)n) / 2;
 105                 numNodes = 2*numLeaves - 1;
 106                 int countleft = 0;
 107                 int countright = 1;
 108
 109                 f1leaves.resize(numNodes);
 110                 f2leaves.resize(numNodes);
 111
 112                 //initialize leaf values
 113                 for (int i = 0; i < numLeaves; i++) {
 114                         f1leaves[i] = new IntNode;
 115                         f1leaves[i]->lvalue = 0;
 116                         f1leaves[i]->rvalue = 0;
 117                         f1leaves[i]->left = NULL;
 118                         f1leaves[i]->right = NULL;
 119
 120                         f2leaves[i] = new IntNode;
 121                         f2leaves[i]->lvalue = 0;
 122                         f2leaves[i]->rvalue = 0;
 123                         f2leaves[i]->left = NULL;
 124                         f2leaves[i]->right = NULL;
 125                 }
 126
 127                 //set pointers to children
 128                 for (int j = numLeaves; j < numNodes; j++) {
 129                         f1leaves[j] = new IntNode;
 130                         f1leaves[j]->left = f1leaves[countleft];
 131                         f1leaves[j]->right = f1leaves[countright];
 132
 133                         f2leaves[j] = new IntNode;
 134                         f2leaves[j]->left = f2leaves[countleft];
 135                         f2leaves[j]->right =f2leaves[countright];
 136
 137                         countleft = countleft + 2;
 138                         countright = countright + 2;
 139                 }
 140
 141                 //point to root
 142                 f1root = f1leaves[numNodes-1];
 143
 144                 //point to root
 145                 f2root = f2leaves[numNodes-1];
 146
 147                 //set coeffients
 148                 setCoef(f2root, 0);
 149         }
 150         catch(exception& e) {
 151                 errorOut(e, "SharedChao1", "initialTree");
 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                 errorOut(e, "SharedChao1", "updateTree");
 165                 exit(1);
 166         }
 167 }
 168
 169 /***********************************************************************/
 170 void SharedChao1::updateBranchf1(IntNode* node, vector<int> bin, int index) {
 171         try {
 172                 //if you have more than one group
 173                 if (index == (bin.size()-1)) {
 174                         if (bin[index] == 1) { node->lvalue++; node->rvalue++; }
 175                         else { node->rvalue++;  }
 176                 }else {
 177                         if (bin[index] == 1) {
 178                                 //follow path as if you are 1
 179                                 updateBranchf1(node->left, bin, index+1);
 180                         }
 181                         //follow path as if you are +
 182                         updateBranchf1(node->right, bin, index+1);
 183                 }
 184         }
 185         catch(exception& e) {
 186                 errorOut(e, "SharedChao1", "updateBranchf1");
 187                 exit(1);
 188         }
 189 }
 190
 191 /***********************************************************************/
 192 void SharedChao1::updateBranchf2(IntNode* node, vector<int> bin, int index) {
 193         try {
 194                 //if you have more than one group
 195                 if (index == (bin.size()-1)) {
 196                         if (bin[index] == 2) { node->lvalue++; node->rvalue++; }
 197                         else { node->rvalue++;  }
 198                 }else {
 199                         if (bin[index] == 2) {
 200                                 //follow path as if you are 1
 201                                 updateBranchf2(node->left, bin, index+1);
 202                         }
 203                         //follow path as if you are +
 204                         updateBranchf2(node->right, bin, index+1);
 205                 }
 206         }
 207         catch(exception& e) {
 208                 errorOut(e, "SharedChao1", "updateBranchf2");
 209                 exit(1);
 210         }
 211 }
 212
 213 /***********************************************************************/
 214 void SharedChao1::setCoef(IntNode* node, int coef) {
 215         try {
 216                 if (node->left != NULL) {
 217                         setCoef(node->left, coef+1);
 218                         setCoef(node->right, coef);
 219                 }else {
 220                         node->lcoef = coef+1;
 221                         node->rcoef = coef;
 222                 }
 223         }
 224         catch(exception& e) {
 225                 errorOut(e, "SharedChao1", "setCoef");
 226                 exit(1);
 227         }
 228 }
 229
 230 /***********************************************************************/
 231 //for debugging purposes
 232 void SharedChao1::printTree() {
 233
 234         mothurOut("F1 leaves"); mothurOutEndLine();
 235         printBranch(f1root);
 236
 237         mothurOut("F2 leaves"); mothurOutEndLine();
 238         printBranch(f2root);
 239
 240
 241 }
 242 /*****************************************************************/
 243 void SharedChao1::printBranch(IntNode* node) {
 244         try {
 245
 246                 // you are not a leaf
 247                 if (node->left != NULL) {
 248                         printBranch(node->left);
 249                         printBranch(node->right);
 250                 }else { //you are a leaf
 251                         mothurOut(toString(node->lvalue)); mothurOutEndLine();
 252                         mothurOut(toString(node->rvalue)); mothurOutEndLine();
 253                 }
 254
 255         }
 256         catch(exception& e) {
 257                 errorOut(e, "SharedChao1", "printBranch");
 258                 exit(1);
 259         }
 260 }
 261
 262 /*****************************************************************/
 263
 264
 265
 266