fullmatrix.cpp

   1 /*
   2  *  fullmatrix.cpp
   3  *  Mothur
   4  *
   5  *  Created by Sarah Westcott on 3/6/09.
   6  *  Copyright 2009 Schloss Lab UMASS Amherst. All rights reserved.
   7  *
   8  */
   9
  10 #include "fullmatrix.h"
  11
  12 /**************************************************************************/
  13 //This constructor reads a distance matrix file and stores the data in the matrix.
  14 FullMatrix::FullMatrix(ifstream& filehandle) {
  15         try{
  16                 globaldata = GlobalData::getInstance();
  17                 groupmap = globaldata->gGroupmap;
  18
  19                 string name, group;
  20                 filehandle >> numSeqs >> name;
  21
  22                 //make the matrix filled with zeros
  23                 matrix.resize(numSeqs);
  24                 for(int i = 0; i < numSeqs; i++) {
  25                         matrix[i].resize(numSeqs, 0.0);
  26                 }
  27                 group = groupmap->getGroup(name);
  28                 if(group == "not found") {      mothurOut("Error: Sequence '" + name + "' was not found in the group file, please correct."); mothurOutEndLine(); exit(1); }
  29                 index.resize(numSeqs);
  30                 index[0].seqName = name;
  31                 index[0].groupName = group;
  32
  33                 //determine if matrix is square or lower triangle
  34                 //if it is square read the distances for the first sequence
  35                 char d;
  36                 bool square;
  37                 while((d=filehandle.get()) != EOF){
  38
  39                         //is d a number meaning its square
  40                         if(isalnum(d)){
  41                                 square = true;
  42                                 filehandle.putback(d);
  43
  44                                 for(int i=0;i<numSeqs;i++){
  45                                         filehandle >> matrix[0][i];
  46                                 }
  47                                 break;
  48                         }
  49
  50                         //is d a line return meaning its lower triangle
  51                         if(d == '\n'){
  52                                 square = false;
  53                                 break;
  54                         }
  55                 }
  56
  57                 //read rest of matrix
  58                 if (square == true) { readSquareMatrix(filehandle); }
  59                 else { readLTMatrix(filehandle); }
  60
  61                 sortGroups(0, numSeqs-1);
  62
  63         }
  64         catch(exception& e) {
  65                 errorOut(e, "FullMatrix", "FullMatrix");
  66                 exit(1);
  67         }
  68 }
  69 /**************************************************************************/
  70 void FullMatrix::readSquareMatrix(ifstream& filehandle) {
  71         try {
  72
  73                 Progress* reading;
  74                 reading = new Progress("Reading matrix:     ", numSeqs * numSeqs);
  75
  76                 int count = 0;
  77
  78                 string group, name;
  79
  80                 for(int i=1;i<numSeqs;i++){
  81                         filehandle >> name;
  82
  83                         group = groupmap->getGroup(name);
  84                         index[i].seqName = name;
  85                         index[i].groupName = group;
  86
  87                         if(group == "not found") {      mothurOut("Error: Sequence '" + name + "' was not found in the group file, please correct."); mothurOutEndLine(); exit(1); }
  88
  89                         for(int j=0;j<numSeqs;j++){
  90                                 filehandle >> matrix[i][j];
  91
  92                                 count++;
  93                                 reading->update(count);
  94                         }
  95                 }
  96                 reading->finish();
  97                 delete reading;
  98         }
  99         catch(exception& e) {
 100                 errorOut(e, "FullMatrix", "readSquareMatrix");
 101                 exit(1);
 102         }
 103 }
 104 /**************************************************************************/
 105 void FullMatrix::readLTMatrix(ifstream& filehandle) {
 106         try {
 107                 Progress* reading;
 108                 reading = new Progress("Reading matrix:     ", numSeqs * (numSeqs - 1) / 2);
 109
 110                 int count = 0;
 111                 float distance;
 112
 113                 string group, name;
 114
 115                 for(int i=1;i<numSeqs;i++){
 116                         filehandle >> name;
 117
 118                         group = groupmap->getGroup(name);
 119                         index[i].seqName = name;
 120                         index[i].groupName = group;
 121
 122                         if(group == "not found") {      mothurOut("Error: Sequence '" + name + "' was not found in the group file, please correct."); mothurOutEndLine();  exit(1); }
 123
 124                         for(int j=0;j<i;j++){
 125                                 filehandle >> distance;
 126                                 matrix[i][j] = distance;  matrix[j][i] = distance;
 127                                 count++;
 128                                 reading->update(count);
 129                         }
 130                 }
 131
 132                 reading->finish();
 133                 delete reading;
 134         }
 135         catch(exception& e) {
 136                 errorOut(e, "FullMatrix", "readLTMatrix");
 137                 exit(1);
 138         }
 139 }
 140
 141 /**************************************************************************/
 142
 143 void FullMatrix::sortGroups(int low, int high){
 144         try{
 145
 146                 if (low < high) {
 147                         int i = low+1;
 148                         int j = high;
 149                         int pivot = (low+high) / 2;
 150
 151                         swapRows(low, pivot);  //puts pivot in final spot
 152
 153                         /* compare value */
 154                         //what group does this row belong to
 155                         string key = index[low].groupName;
 156
 157                         /* partition */
 158                         while(i <= j) {
 159                                 /* find member above ... */
 160                                 while((i <= high) && (index[i].groupName <= key))       {  i++;  }
 161
 162                                 /* find element below ... */
 163                                 while((j >= low) && (index[j].groupName > key))         {  j--;  }
 164
 165                                 if(i < j) {
 166                                         swapRows(i, j);
 167                                 }
 168                         }
 169
 170                         swapRows(low, j);
 171
 172                         /* recurse */
 173                         sortGroups(low, j-1);
 174                         sortGroups(j+1, high);
 175                 }
 176
 177         }
 178         catch(exception& e) {
 179                 errorOut(e, "FullMatrix", "sortGroups");
 180                 exit(1);
 181         }
 182 }
 183
 184 /**************************************************************************/
 185 void FullMatrix::swapRows(int i, int j) {
 186         try {
 187
 188                 float y;
 189                 string z, name;
 190
 191                 /* swap rows*/
 192                 for (int h = 0; h < numSeqs; h++) {
 193                         y = matrix[i][h];
 194                         matrix[i][h] = matrix[j][h];
 195                         matrix[j][h] = y;
 196                 }
 197
 198                 /* swap columns*/
 199                 for (int b = 0; b < numSeqs; b++) {
 200                         y = matrix[b][i];
 201                         matrix[b][i] = matrix[b][j];
 202                         matrix[b][j] = y;
 203                 }
 204
 205                 //swap map elements
 206                 z = index[i].groupName;
 207                 index[i].groupName = index[j].groupName;
 208                 index[j].groupName = z;
 209
 210                 name = index[i].seqName;
 211                 index[i].seqName = index[j].seqName;
 212                 index[j].seqName = name;
 213
 214
 215         }
 216         catch(exception& e) {
 217                 errorOut(e, "FullMatrix", "swapRows");
 218                 exit(1);
 219         }
 220 }
 221 /**************************************************************************/
 222
 223 float FullMatrix::get(int i, int j){    return matrix[i][j];            }
 224
 225 /**************************************************************************/
 226
 227 vector<string> FullMatrix::getGroups(){ return groups;          }
 228
 229 /**************************************************************************/
 230
 231 vector<int> FullMatrix::getSizes(){     return sizes;           }
 232
 233 /**************************************************************************/
 234
 235 int FullMatrix::getNumGroups(){ return groups.size();           }
 236
 237 /**************************************************************************/
 238
 239 int FullMatrix::getNumSeqs(){   return numSeqs;         }
 240
 241 /**************************************************************************/
 242
 243 void FullMatrix::printMatrix(ostream& out) {
 244         try{
 245                 for (int i = 0; i < numSeqs; i++) {
 246                         out << "row " << i << " group = " << index[i].groupName << " name = " << index[i].seqName << endl;
 247                         for (int j = 0; j < numSeqs; j++) {
 248                                 out << i << '\t' << j << '\t' << matrix[i][j] << endl;
 249                         }
 250                         out << endl;
 251                 }
 252
 253                 for (int i = 0; i < numSeqs; i++) {  out << i << '\t' <<  index[i].seqName << endl;  }
 254         }
 255         catch(exception& e) {
 256                 errorOut(e, "FullMatrix", "printMatrix");
 257                 exit(1);
 258         }
 259 }
 260
 261 /**************************************************************************/
 262