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