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