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                 vector<double> D;
  19
  20                 numGroups = globaldata->Groups.size();
  21
  22                 //calculate number of comparisons i.e. with groups A,B,C = AB, AC, BC = 3;
  23                 int count = 0;
  24                 for (int i=0; i<numGroups; i++) {
  25                         for (int l = i+1; l < numGroups; l++) {
  26                                 //initialize weighted scores
  27                                 WScore[globaldata->Groups[i]+globaldata->Groups[l]] = 0.0;
  28
  29                                 vector<string> groups; groups.push_back(globaldata->Groups[i]); groups.push_back(globaldata->Groups[l]);
  30
  31                                 D.push_back(0.0000); //initialize a spot in D for each combination
  32
  33                                 /********************************************************/
  34                                 //calculate a D value for each group combo
  35                                 for(int v=0;v<t->getNumLeaves();v++){
  36                                         int index = v;
  37                                         double sum = 0.0000;
  38
  39                                         //while you aren't at root
  40                                         while(t->tree[index].getParent() != -1){
  41
  42                                                 //if you have a BL
  43                                                 if(t->tree[index].getBranchLength() != -1){
  44                                                         sum += abs(t->tree[index].getBranchLength());
  45                                                 }
  46                                                 index = t->tree[index].getParent();
  47                                         }
  48
  49                                         //get last breanch length added
  50                                         if(t->tree[index].getBranchLength() != -1){
  51                                                 sum += abs(t->tree[index].getBranchLength());
  52                                         }
  53
  54                                         //is this sum from a sequence which is in one of the users groups
  55                                         if (inUsersGroups(t->tree[v].getGroup(), globaldata->Groups) == true) {
  56                                                 //is this sum from a sequence which is in this groupCombo
  57                                                 if (inUsersGroups(t->tree[v].getGroup(), groups)) {
  58                                                         int numSeqsInGroupI, numSeqsInGroupL;
  59
  60                                                         map<string, int>::iterator it;
  61                                                         it = t->tree[v].pcount.find(groups[0]);
  62                                                         if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i
  63                                                                 numSeqsInGroupI = it->second;
  64                                                         }else{ numSeqsInGroupI = 0; }
  65
  66                                                         it = t->tree[v].pcount.find(groups[1]);
  67                                                         if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l
  68                                                                 numSeqsInGroupL = it->second;
  69                                                         }else{ numSeqsInGroupL = 0; }
  70
  71                                                         double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
  72
  73                                                         //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
  74
  75                                                         D[count] += weightedSum;
  76                                                 }
  77                                         }
  78                                 }
  79                                 /*********************************************************/
  80                                 count++;
  81                         }
  82                 }
  83
  84                 data.clear(); //clear out old values
  85
  86                 for(int i=0;i<t->getNumNodes();i++){
  87                         //calculate weighted score for each of the group comb i.e. with groups A,B,C = AB, AC, BC.
  88                         for (int b=0; b<numGroups; b++) {
  89                                 for (int l = b+1; l < numGroups; l++) {
  90                                         //calculate a u value for each combo
  91                                         double u;
  92                                         //does this node have descendants from group b-1
  93                                         it = t->tree[i].pcount.find(globaldata->Groups[b]);
  94                                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
  95                                         if (it != t->tree[i].pcount.end()) {
  96                                                 u = (double) t->tree[i].pcount[globaldata->Groups[b]] / (double) tmap->seqsPerGroup[globaldata->Groups[b]];
  97                                         }else { u = 0.00; }
  98
  99                                         //does this node have descendants from group l
 100                                         it = t->tree[i].pcount.find(globaldata->Groups[l]);
 101                                         //if it does subtract their percentage from u
 102                                         if (it != t->tree[i].pcount.end()) {
 103                                                 u -= (double) t->tree[i].pcount[globaldata->Groups[l]] / (double) tmap->seqsPerGroup[globaldata->Groups[l]];
 104                                         }
 105
 106                                         u = abs(u * t->tree[i].getBranchLength());
 107
 108                                         //save groupcombs u value
 109                                         WScore[globaldata->Groups[b]+globaldata->Groups[l]] += u;
 110                                 /*********************************************************/
 111                                 }
 112                         }
 113                 }
 114
 115                 //calculate weighted score for each group combination
 116                 double UN;
 117                 count = 0;
 118                 for (int i=0; i<numGroups; i++) {
 119                         for (int l = i+1; l < numGroups; l++) {
 120                                 UN = (WScore[globaldata->Groups[i]+globaldata->Groups[l]] / D[count]);
 121
 122                                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 123                                 data.push_back(UN);
 124                                 count++;
 125                         }
 126                 }
 127                 return data;
 128         }
 129         catch(exception& e) {
 130                 m->errorOut(e, "Weighted", "getValues");
 131                 exit(1);
 132         }
 133 }
 134
 135 /**************************************************************************************************/
 136 EstOutput Weighted::getValues(Tree* t, string groupA, string groupB) {
 137  try {
 138                 globaldata = GlobalData::getInstance();
 139
 140                 data.clear(); //clear out old values
 141
 142                 //initialize weighted score
 143                 WScore[(groupA+groupB)] = 0.0;
 144                 float D = 0.0;
 145
 146                 vector<string> groups; groups.push_back(groupA); groups.push_back(groupB);
 147
 148                 /********************************************************/
 149                 //calculate a D value for the group combo
 150                 for(int v=0;v<t->getNumLeaves();v++){
 151                         int index = v;
 152                         double sum = 0.0000;
 153
 154                         //while you aren't at root
 155                         while(t->tree[index].getParent() != -1){
 156
 157                                 //if you have a BL
 158                                 if(t->tree[index].getBranchLength() != -1){
 159                                         sum += abs(t->tree[index].getBranchLength());
 160                                 }
 161                                 index = t->tree[index].getParent();
 162                         }
 163
 164                         //get last breanch length added
 165                         if(t->tree[index].getBranchLength() != -1){
 166                                 sum += abs(t->tree[index].getBranchLength());
 167                         }
 168
 169                         if (inUsersGroups(t->tree[v].getGroup(), groups)) {
 170                                 int numSeqsInGroupI, numSeqsInGroupL;
 171
 172                                 map<string, int>::iterator it;
 173                                 it = t->tree[v].pcount.find(groups[0]);
 174                                 if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group i
 175                                         numSeqsInGroupI = it->second;
 176                                 }else{ numSeqsInGroupI = 0; }
 177
 178                                 it = t->tree[v].pcount.find(groups[1]);
 179                                 if (it != t->tree[v].pcount.end()) { //this leaf node contains seqs from group l
 180                                         numSeqsInGroupL = it->second;
 181                                 }else{ numSeqsInGroupL = 0; }
 182
 183                                 double weightedSum = ((numSeqsInGroupI * sum) / (double)tmap->seqsPerGroup[groups[0]]) + ((numSeqsInGroupL * sum) / (double)tmap->seqsPerGroup[groups[1]]);
 184
 185                                 //sum /= (double)tmap->seqsPerGroup[t->tree[v].getGroup()];
 186
 187                                 D += weightedSum;
 188                         }
 189                 }
 190                 /********************************************************/
 191
 192                 //calculate u for the group comb
 193                 for(int i=0;i<t->getNumNodes();i++){
 194                         double u;
 195                         //does this node have descendants from groupA
 196                         it = t->tree[i].pcount.find(groupA);
 197                         //if it does u = # of its descendants with a certain group / total number in tree with a certain group
 198                         if (it != t->tree[i].pcount.end()) {
 199                                 u = (double) t->tree[i].pcount[groupA] / (double) tmap->seqsPerGroup[groupA];
 200                         }else { u = 0.00; }
 201
 202
 203                         //does this node have descendants from group l
 204                         it = t->tree[i].pcount.find(groupB);
 205                         //if it does subtract their percentage from u
 206                         if (it != t->tree[i].pcount.end()) {
 207                                 u -= (double) t->tree[i].pcount[groupB] / (double) tmap->seqsPerGroup[groupB];
 208                         }
 209
 210                         u = abs(u * t->tree[i].getBranchLength());
 211
 212                         //save groupcombs u value
 213                         WScore[(groupA+groupB)] += u;
 214                 }
 215
 216                 /********************************************************/
 217
 218                 //calculate weighted score for the group combination
 219                 double UN;
 220                 UN = (WScore[(groupA+groupB)] / D);
 221
 222                 if (isnan(UN) || isinf(UN)) { UN = 0; }
 223                 data.push_back(UN);
 224
 225                 return data;
 226         }
 227         catch(exception& e) {
 228                 m->errorOut(e, "Weighted", "getValues");
 229                 exit(1);
 230         }
 231 }
 232
 233
 234
 235
 236
 237
 238
 239
 240