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