hcluster.cpp

   1 /*
   2  *  hcluster.cpp
   3  *  Mothur
   4  *
   5  *  Created by westcott on 10/13/09.
   6  *  Copyright 2009 Schloss Lab. All rights reserved.
   7  *
   8  */
   9
  10 #include "hcluster.h"
  11 #include "rabundvector.hpp"
  12 #include "listvector.hpp"
  13 #include "sparsematrix.hpp"
  14
  15 /***********************************************************************/
  16
  17 HCluster::HCluster(RAbundVector* rav, ListVector* lv) :  rabund(rav), list(lv){
  18         try {
  19
  20                 numSeqs = list->getNumSeqs();
  21
  22                 //initialize cluster array
  23                 for (int i = 0; i < numSeqs; i++) {
  24                         clusterNode temp(1, -1, i);
  25                         clusterArray.push_back(temp);
  26                 }
  27
  28         }
  29         catch(exception& e) {
  30                 errorOut(e, "HCluster", "HCluster");
  31                 exit(1);
  32         }
  33 }
  34 /***********************************************************************/
  35
  36 void HCluster::clusterBins(){
  37         try {
  38                 //cout << smallCol << '\t' << smallRow << '\t' << smallDist << '\t' << rabund->get(clusterArray[smallRow].smallChild) << '\t' << rabund->get(clusterArray[smallCol].smallChild);
  39
  40                 rabund->set(clusterArray[smallCol].smallChild, rabund->get(clusterArray[smallRow].smallChild)+rabund->get(clusterArray[smallCol].smallChild));
  41                 rabund->set(clusterArray[smallRow].smallChild, 0);
  42                 rabund->setLabel(toString(smallDist));
  43
  44                 //cout << '\t' << rabund->get(clusterArray[smallRow].smallChild) << '\t' << rabund->get(clusterArray[smallCol].smallChild) << endl;
  45         }
  46         catch(exception& e) {
  47                 errorOut(e, "HCluster", "clusterBins");
  48                 exit(1);
  49         }
  50
  51
  52 }
  53
  54 /***********************************************************************/
  55
  56 void HCluster::clusterNames(){
  57         try {
  58                 ///cout << smallCol << '\t' << smallRow << '\t' << smallDist << '\t' << list->get(clusterArray[smallRow].smallChild) << '\t' << list->get(clusterArray[smallCol].smallChild);
  59
  60                 list->set(clusterArray[smallCol].smallChild, list->get(clusterArray[smallRow].smallChild)+','+list->get(clusterArray[smallCol].smallChild));
  61                 list->set(clusterArray[smallRow].smallChild, "");
  62                 list->setLabel(toString(smallDist));
  63
  64                 //cout << '\t' << list->get(clusterArray[smallRow].smallChild) << '\t' << list->get(clusterArray[smallCol].smallChild) << endl;
  65
  66     }
  67         catch(exception& e) {
  68                 errorOut(e, "HCluster", "clusterNames");
  69                 exit(1);
  70         }
  71
  72 }
  73 /***********************************************************************/
  74 int HCluster::getUpmostParent(int node){
  75         try {
  76
  77                 while (clusterArray[node].parent != -1) {
  78                         node = clusterArray[node].parent;
  79                 }
  80
  81                 return node;
  82         }
  83         catch(exception& e) {
  84                 errorOut(e, "HCluster", "getUpmostParent");
  85                 exit(1);
  86         }
  87 }
  88 /***********************************************************************/
  89 void HCluster::printInfo(){
  90         try {
  91
  92                 cout << "link table" << endl;
  93                 for (it = activeLinks.begin(); it!= activeLinks.end(); it++) {
  94                         cout << it->first << " = " << it->second << endl;
  95                 }
  96                 cout << endl;
  97                 for (int i = 0; i < linkTable.size(); i++) {
  98                         for (it = linkTable[i].begin(); it != linkTable[i].end(); it++) {
  99                                 cout << it->first << '\t' << it->second << '\t' << '\t';
 100                         }
 101                         cout << endl;
 102                 }
 103                 cout << endl << "clusterArray" << endl;
 104
 105                 for (int i = 0; i < clusterArray.size(); i++) {
 106                         cout << i << '\t' << clusterArray[i].numSeq << '\t' << clusterArray[i].parent << '\t' << clusterArray[i].smallChild << endl;
 107                 }
 108                 cout << endl;
 109
 110
 111         }
 112         catch(exception& e) {
 113                 errorOut(e, "HCluster", "getUpmostParent");
 114                 exit(1);
 115         }
 116 }
 117 /***********************************************************************/
 118 int HCluster::makeActive() {
 119         try {
 120
 121                 int linkValue = 1;
 122 //cout << "active - here" << endl;
 123                 it = activeLinks.find(smallRow);
 124                 it2 = activeLinks.find(smallCol);
 125
 126                 if ((it == activeLinks.end()) && (it2 == activeLinks.end())) { //both are not active so add them
 127                         int size = linkTable.size();
 128                         map<int, int> temp; map<int, int> temp2;
 129
 130                         //add link to eachother
 131                         temp[smallRow] = 1;                                                     //         1    2
 132                         temp2[smallCol] = 1;                                            // 1   0        1
 133                                                                                                                 // 2   1        0
 134                         linkTable.push_back(temp);
 135                         linkTable.push_back(temp2);
 136
 137                         //add to activeLinks
 138                         activeLinks[smallRow] = size;
 139                         activeLinks[smallCol] = size+1;
 140 //cout << "active - here1" << endl;
 141                 }else if ((it != activeLinks.end()) && (it2 == activeLinks.end())) {  //smallRow is active, smallCol is not
 142                          int size = linkTable.size();
 143                          int alreadyActiveRow = it->second;
 144                          map<int, int> temp;
 145
 146                         //add link to eachother
 147                         temp[smallRow] = 1;                                                     //         6    2       3       5
 148                         linkTable.push_back(temp);                                      // 6   0        1       2       0
 149                         linkTable[alreadyActiveRow][smallCol] = 1;      // 2   1        0       1       1
 150                                                                                                                 // 3   2        1       0       0
 151                                                                                                                 // 5   0    1   0   0
 152                         //add to activeLinks
 153                         activeLinks[smallCol] = size;
 154 //cout << "active - here2" << endl;
 155                 }else if ((it == activeLinks.end()) && (it2 != activeLinks.end())) {  //smallCol is active, smallRow is not
 156                          int size = linkTable.size();
 157                          int alreadyActiveCol = it2->second;
 158                          map<int, int> temp;
 159
 160                         //add link to eachother
 161                         temp[smallCol] = 1;                                                     //         6    2       3       5
 162                         linkTable.push_back(temp);                                      // 6   0        1       2       0
 163                         linkTable[alreadyActiveCol][smallRow] = 1;      // 2   1        0       1       1
 164                                                                                                                 // 3   2        1       0       0
 165                                                                                                                 // 5   0    1   0   0
 166                         //add to activeLinks
 167                         activeLinks[smallRow] = size;
 168 //cout << "active - here3" << endl;
 169                 }else { //both are active so add one
 170                         int row = it->second;
 171                         int col = it2->second;
 172 //cout << row << '\t' << col << endl;
 173
 174                         linkTable[row][smallCol]++;
 175                         linkTable[col][smallRow]++;
 176                         linkValue = linkTable[row][smallCol];
 177 //cout << "active - here4" << endl;
 178                 }
 179
 180                 return linkValue;
 181         }
 182         catch(exception& e) {
 183                 errorOut(e, "HCluster", "makeActive");
 184                 exit(1);
 185         }
 186 }
 187 /***********************************************************************/
 188 void HCluster::updateArrayandLinkTable() {
 189         try {
 190                 //if cluster was made update clusterArray and linkTable
 191                         int size = clusterArray.size();
 192
 193                         //add new node
 194                         clusterNode temp(clusterArray[smallRow].numSeq + clusterArray[smallCol].numSeq, -1, clusterArray[smallCol].smallChild);
 195                         clusterArray.push_back(temp);
 196
 197                         //update child nodes
 198                         clusterArray[smallRow].parent = size;
 199                         clusterArray[smallCol].parent = size;
 200
 201                         //update linkTable by merging clustered rows and columns
 202                         int rowSpot = activeLinks[smallRow];
 203                         int colSpot = activeLinks[smallCol];
 204         //cout << "here" << endl;
 205                         //fix old rows
 206                         for (int i = 0; i < linkTable.size(); i++) {
 207                                 //check if they are in map
 208                                 it = linkTable[i].find(smallRow);
 209                                 it2 = linkTable[i].find(smallCol);
 210
 211                                 if ((it!=linkTable[i].end()) && (it2!=linkTable[i].end())) { //they are both there
 212                                         linkTable[i][size] = linkTable[i][smallRow]+linkTable[i][smallCol];
 213                                         linkTable[i].erase(smallCol); //delete col row
 214                                         linkTable[i].erase(smallRow); //delete col row
 215                                 }else if ((it==linkTable[i].end()) && (it2!=linkTable[i].end())) { //only col
 216                                         linkTable[i][size] = linkTable[i][smallCol];
 217                                         linkTable[i].erase(smallCol); //delete col
 218                                 }else if ((it!=linkTable[i].end()) && (it2==linkTable[i].end())) { //only row
 219                                         linkTable[i][size] = linkTable[i][smallRow];
 220                                         linkTable[i].erase(smallRow); //delete col
 221                                 }
 222                         }
 223         //printInfo();
 224 //cout << "here2" << endl;
 225                         //merge their values
 226                         for (it = linkTable[rowSpot].begin(); it != linkTable[rowSpot].end(); it++) {
 227                                 it2 = linkTable[colSpot].find(it->first);  //does the col also have this
 228
 229                                 if (it2 == linkTable[colSpot].end()) { //not there so add it
 230                                         linkTable[colSpot][it->first] = it->second;
 231                                 }else { //merge them
 232                                         linkTable[colSpot][it->first] = it->second+it2->second;
 233                                 }
 234                         }
 235 //cout << "here3" << endl;
 236                         linkTable[colSpot].erase(size);
 237                         linkTable.erase(linkTable.begin()+rowSpot);  //delete col
 238         //printInfo();
 239                         //update activerows
 240                         activeLinks.erase(smallRow);
 241                         activeLinks.erase(smallCol);
 242
 243                         if(rowSpot>colSpot) {   activeLinks[size] = colSpot;    }
 244                         else{   activeLinks[size] = colSpot-1;  }
 245
 246                         //adjust everybody elses spot since you deleted - time vs. space
 247                         for (it = activeLinks.begin(); it != activeLinks.end(); it++) {
 248                                 if (it->second > rowSpot) {  activeLinks[it->first]--;  }
 249                         }
 250
 251
 252 //cout << "here4" << endl;
 253
 254         }
 255         catch(exception& e) {
 256                 errorOut(e, "HCluster", "updateArrayandLinkTable");
 257                 exit(1);
 258         }
 259 }
 260 /***********************************************************************/
 261 bool HCluster::update(int row, int col, float distance){
 262         try {
 263
 264                 smallRow = row;
 265                 smallCol = col;
 266                 smallDist = distance;
 267                 bool clustered = false;
 268
 269                 //find upmost parent of row and col
 270                 smallRow = getUpmostParent(smallRow);
 271                 smallCol = getUpmostParent(smallCol);
 272         //cout << "smallRow = " << smallRow << " smallCol = " << smallCol << endl;
 273
 274                 //are they active in the link table
 275                 int linkValue = makeActive(); //after this point this nodes info is active in linkTable
 276         //printInfo();
 277 //cout << "linkValue = " << linkValue << " times = " << (clusterArray[smallRow].numSeq * clusterArray[smallCol].numSeq) << endl;
 278                 //can we cluster???
 279                 if (linkValue == (clusterArray[smallRow].numSeq * clusterArray[smallCol].numSeq)) {
 280                         updateArrayandLinkTable();
 281                         clusterBins();
 282                         clusterNames();
 283                         clustered = true;
 284                         //printInfo();
 285                 }
 286
 287
 288                 return clustered;
 289         }
 290         catch(exception& e) {
 291                 errorOut(e, "HCluster", "update");
 292                 exit(1);
 293         }
 294 }
 295
 296
 297 /***********************************************************************/
 298
 299