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, GroupMap* g, bool s) : groupmap(g), sim(s) {
  15         try{
  16                 m = MothurOut::getInstance();
  17
  18                 string name, group;
  19
  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") {      m->mothurOut("Error: Sequence '" + name + "' was not found in the group file, please correct."); m->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                                         if (sim) {  matrix[0][i] = 1.0 - 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                 filehandle.close();
  63
  64                 if (!m->control_pressed) { sortGroups(0, numSeqs-1); }
  65
  66         }
  67         catch(exception& e) {
  68                 m->errorOut(e, "FullMatrix", "FullMatrix");
  69                 exit(1);
  70         }
  71 }
  72 /**************************************************************************/
  73 int FullMatrix::readSquareMatrix(ifstream& filehandle) {
  74         try {
  75
  76                 Progress* reading;
  77                 reading = new Progress("Reading matrix:     ", numSeqs * numSeqs);
  78
  79                 int count = 0;
  80
  81                 string group, name;
  82
  83                 for(int i=1;i<numSeqs;i++){
  84                         filehandle >> name;
  85
  86                         group = groupmap->getGroup(name);
  87                         index[i].seqName = name;
  88                         index[i].groupName = group;
  89
  90                         if(group == "not found") {      m->mothurOut("Error: Sequence '" + name + "' was not found in the group file, please correct."); m->mothurOutEndLine(); exit(1); }
  91
  92                         for(int j=0;j<numSeqs;j++){
  93                                 if (m->control_pressed) { delete reading;  return 0; }
  94
  95                                 filehandle >> matrix[i][j];
  96                                 if (sim) {  matrix[i][j] = 1.0 - matrix[i][j];  }
  97
  98                                 count++;
  99                                 reading->update(count);
 100                         }
 101                 }
 102
 103                 if (m->control_pressed) { delete reading;  return 0; }
 104
 105                 reading->finish();
 106                 delete reading;
 107
 108                 return 0;
 109         }
 110         catch(exception& e) {
 111                 m->errorOut(e, "FullMatrix", "readSquareMatrix");
 112                 exit(1);
 113         }
 114 }
 115 /**************************************************************************/
 116 int FullMatrix::readLTMatrix(ifstream& filehandle) {
 117         try {
 118
 119                 Progress* reading;
 120                 reading = new Progress("Reading matrix:     ", numSeqs * (numSeqs - 1) / 2);
 121
 122                 int count = 0;
 123                 float distance;
 124
 125                 string group, name;
 126
 127                 for(int i=1;i<numSeqs;i++){
 128                         filehandle >> name;
 129
 130                         group = groupmap->getGroup(name);
 131                         index[i].seqName = name;
 132                         index[i].groupName = group;
 133
 134                         if(group == "not found") {      m->mothurOut("Error: Sequence '" + name + "' was not found in the group file, please correct."); m->mothurOutEndLine();  exit(1); }
 135
 136                         for(int j=0;j<i;j++){
 137                                 if (m->control_pressed) { delete reading;  return 0; }
 138
 139                                 filehandle >> distance;
 140                                 if (sim) {  distance = 1.0 - distance;  }
 141
 142                                 matrix[i][j] = distance;  matrix[j][i] = distance;
 143
 144                                 count++;
 145                                 reading->update(count);
 146                         }
 147                 }
 148
 149                 if (m->control_pressed) { delete reading;  return 0; }
 150
 151                 reading->finish();
 152                 delete reading;
 153
 154                 return 0;
 155         }
 156         catch(exception& e) {
 157                 m->errorOut(e, "FullMatrix", "readLTMatrix");
 158                 exit(1);
 159         }
 160 }
 161
 162 /**************************************************************************/
 163
 164 void FullMatrix::sortGroups(int low, int high){
 165         try{
 166
 167                 if (low < high) {
 168                         int i = low+1;
 169                         int j = high;
 170                         int pivot = (low+high) / 2;
 171
 172                         swapRows(low, pivot);  //puts pivot in final spot
 173
 174                         /* compare value */
 175                         //what group does this row belong to
 176                         string key = index[low].groupName;
 177
 178                         /* partition */
 179                         while(i <= j) {
 180                                 /* find member above ... */
 181                                 while((i <= high) && (index[i].groupName <= key))       {  i++;  }
 182
 183                                 /* find element below ... */
 184                                 while((j >= low) && (index[j].groupName > key))         {  j--;  }
 185
 186                                 if(i < j) {
 187                                         swapRows(i, j);
 188                                 }
 189                         }
 190
 191                         swapRows(low, j);
 192
 193                         /* recurse */
 194                         sortGroups(low, j-1);
 195                         sortGroups(j+1, high);
 196                 }
 197
 198         }
 199         catch(exception& e) {
 200                 m->errorOut(e, "FullMatrix", "sortGroups");
 201                 exit(1);
 202         }
 203 }
 204
 205 /**************************************************************************/
 206 void FullMatrix::swapRows(int i, int j) {
 207         try {
 208
 209                 float y;
 210                 string z, name;
 211
 212                 /* swap rows*/
 213                 for (int h = 0; h < numSeqs; h++) {
 214                         y = matrix[i][h];
 215                         matrix[i][h] = matrix[j][h];
 216                         matrix[j][h] = y;
 217                 }
 218
 219                 /* swap columns*/
 220                 for (int b = 0; b < numSeqs; b++) {
 221                         y = matrix[b][i];
 222                         matrix[b][i] = matrix[b][j];
 223                         matrix[b][j] = y;
 224                 }
 225
 226                 //swap map elements
 227                 z = index[i].groupName;
 228                 index[i].groupName = index[j].groupName;
 229                 index[j].groupName = z;
 230
 231                 name = index[i].seqName;
 232                 index[i].seqName = index[j].seqName;
 233                 index[j].seqName = name;
 234
 235
 236         }
 237         catch(exception& e) {
 238                 m->errorOut(e, "FullMatrix", "swapRows");
 239                 exit(1);
 240         }
 241 }
 242 /**************************************************************************/
 243
 244 float FullMatrix::get(int i, int j){    return matrix[i][j];            }
 245
 246 /**************************************************************************/
 247
 248 vector<string> FullMatrix::getGroups(){ return groups;          }
 249
 250 /**************************************************************************/
 251
 252 vector<int> FullMatrix::getSizes(){     return sizes;           }
 253
 254 /**************************************************************************/
 255
 256 int FullMatrix::getNumGroups(){ return groups.size();           }
 257
 258 /**************************************************************************/
 259
 260 int FullMatrix::getNumSeqs(){   return numSeqs;         }
 261
 262 /**************************************************************************/
 263
 264 void FullMatrix::printMatrix(ostream& out) {
 265         try{
 266                 for (int i = 0; i < numSeqs; i++) {
 267                         out << "row " << i << " group = " << index[i].groupName << " name = " << index[i].seqName << endl;
 268                         for (int j = 0; j < numSeqs; j++) {
 269                                 out << i << '\t' << j << '\t' << matrix[i][j] << endl;
 270                         }
 271                         out << endl;
 272                 }
 273
 274                 for (int i = 0; i < numSeqs; i++) {  out << i << '\t' <<  index[i].seqName << endl;  }
 275         }
 276         catch(exception& e) {
 277                 m->errorOut(e, "FullMatrix", "printMatrix");
 278                 exit(1);
 279         }
 280 }
 281
 282 /**************************************************************************/
 283