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