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