weighted.cpp

   1 /*
   2  *  weighted.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 2/9/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "weighted.h"
  11
  12 /**************************************************************************************************/
  13
  14 EstOutput Weighted::getValues(Tree* t) {
  15     try {
  16                 globaldata = GlobalData::getInstance();
  17                 int numGroups;
  18
  19                 //if the user has not entered specific groups to analyze then do them all
  20                 if (globaldata->Groups.size() == 0) {
  21                         numGroups = tmap->getNumGroups();
  22                 }else {
  23                         numGroups = globaldata->Groups.size();
  24                 }
  25
  26                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
  27                 int n = 1;
  28                 for (int i=1; i<numGroups; i++) {
  29                         for (int l = n; l < numGroups; l++) {
  30                                 //initialize weighted scores
  31                                 if (globaldata->Groups.size() == 0) {
  32                                         WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] = 0.0;
  33                                 }else {
  34                                         WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] = 0.0;
  35                                 }
  36                         }
  37                 }
  38
  39                 data.clear(); //clear out old values
  40
  41                 double D = 0.0000;
  42
  43                 for(int i=0;i<t->getNumLeaves();i++){
  44                         int index = i;
  45                         double sum = 0.0000;
  46
  47                         //while you aren't at root
  48                         while(t->tree[index].getParent() != -1){
  49
  50                                 if(t->tree[index].pGroups.size() != 0){
  51                                         sum += t->tree[index].getBranchLength();
  52                                 }
  53
  54                                 //old_index = you
  55                                 int old_index = index;
  56                                 //index = your parent
  57                                 index = t->tree[index].getParent();
  58
  59                                 //if you grandparent is the root
  60                                 if(t->tree[index].getParent() == -1 && t->tree[old_index].pGroups.size() != 0){
  61                                         int lc = t->tree[t->tree[index].getLChild()].pGroups.size();
  62                                         int rc = t->tree[t->tree[index].getRChild()].pGroups.size();
  63
  64
  65                                         if(lc == 0 || rc == 0){
  66                                                 sum -= t->tree[old_index].getBranchLength();
  67                                         }
  68                                 }
  69                         }
  70
  71                         if(t->tree[i].getGroup() != ""){
  72                                 sum /= (double)tmap->seqsPerGroup[t->tree[i].getGroup()];
  73                                 D += sum;
  74                         }
  75                 }
  76
  77
  78                 for(int i=0;i<t->getNumNodes();i++){
  79                         //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
  80                         n = 1;
  81                         for (int b=1; b<numGroups; b++) {
  82                                 for (int l = n; l < numGroups; l++) {
  83                                         double u;
  84                                         //the user has not entered specific groups
  85                                         if (globaldata->Groups.size() == 0) {
  86                                                 //does this node have descendants from group b-1
  87                                                 it = t->tree[i].pcount.find(tmap->namesOfGroups[b-1]);
  88                                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
  89                                                 if (it != t->tree[i].pcount.end()) {
  90                                                         u = (double) t->tree[i].pcount[tmap->namesOfGroups[b-1]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[b-1]];
  91                                                 }else { u = 0.00; }
  92
  93                                                 //does this node have descendants from group l
  94                                                 it = t->tree[i].pcount.find(tmap->namesOfGroups[l]);
  95                                                 //if it does subtract their percentage from u
  96                                                 if (it != t->tree[i].pcount.end()) {
  97                                                         u -= (double) t->tree[i].pcount[tmap->namesOfGroups[l]] / (double) tmap->seqsPerGroup[tmap->namesOfGroups[l]];
  98                                                 }
  99
 100                                                 u = abs(u) * t->tree[i].getBranchLength();
 101
 102                                                 //save groupcombs u value
 103                                                 WScore[tmap->namesOfGroups[b-1]+tmap->namesOfGroups[l]] += u;
 104
 105                                         //the user has entered specific groups
 106                                         }else {
 107                                                 //does this node have descendants from group b-1
 108                                                 it = t->tree[i].pcount.find(globaldata->Groups[b-1]);
 109                                                 //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 110                                                 if (it != t->tree[i].pcount.end()) {
 111                                                         u = (double) t->tree[i].pcount[globaldata->Groups[b-1]] / (double) tmap->seqsPerGroup[globaldata->Groups[b-1]];
 112                                                 }else { u = 0.00; }
 113
 114                                                 //does this node have descendants from group l
 115                                                 it = t->tree[i].pcount.find(globaldata->Groups[l]);
 116                                                 //if it does subtract their percentage from u
 117                                                 if (it != t->tree[i].pcount.end()) {
 118                                                         u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]];
 119                                                 }
 120
 121                                                 u = abs(u) * t->tree[i].getBranchLength();
 122
 123                                                 //save groupcombs u value
 124                                                 WScore[globaldata->Groups[b-1]+globaldata->Groups[l]] += u;
 125                                         }
 126                                 }
 127                                 n++;
 128                         }
 129                 }
 130
 131                 //calculate weighted score for each group combination
 132                 double UN;
 133                 n = 1;
 134                 for (int i=1; i<numGroups; i++) {
 135                         for (int l = n; l < numGroups; l++) {
 136                                 //the user has not entered specific groups
 137                                 if (globaldata->Groups.size() == 0) {
 138                                         UN = (WScore[tmap->namesOfGroups[i-1]+tmap->namesOfGroups[l]] / D);
 139                                 }else {//they have entered specific groups
 140                                         UN = (WScore[globaldata->Groups[i-1]+globaldata->Groups[l]] / D);
 141                                 }
 142                                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 143                                 data.push_back(UN);
 144                         }
 145                 }
 146
 147                 return data;
 148         }
 149         catch(exception& e) {
 150                 cout << "Standard Error: " << e.what() << " has occurred in the Weighted class Function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 151                 exit(1);
 152         }
 153         catch(...) {
 154                 cout << "An unknown error has occurred in the Weighted class function getValues. Please contact Pat Schloss at pschloss@microbio.umass.edu." << "\n";
 155                 exit(1);
 156         }
 157
 158 }
 159
 160 /**************************************************************************************************/