binaries/src/MSAProbs-0.9.7/MSAProbs/SparseMatrix.h

   1 /////////////////////////////////////////////////////////////////
   2 // SparseMatrix.h
   3 //
   4 // Sparse matrix computations
   5 /////////////////////////////////////////////////////////////////
   6
   7 #ifndef SPARSEMATRIX_H
   8 #define SPARSEMATRIX_H
   9
  10 #include <iostream>
  11
  12 using namespace std;
  13
  14 const float POSTERIOR_CUTOFF = 0.01;         // minimum posterior probability
  15 // value that is maintained in the
  16 // sparse matrix representation
  17
  18 typedef pair<int, float> PIF;                 // Sparse matrix entry type
  19 //   first --> column
  20 //   second --> value
  21
  22 /////////////////////////////////////////////////////////////////
  23 // SparseMatrix
  24 //
  25 // Class for sparse matrix computations
  26 /////////////////////////////////////////////////////////////////
  27
  28 class SparseMatrix {
  29
  30         int seq1Length, seq2Length;                     // dimensions of matrix
  31         VI rowSize;                              // rowSize[i] = # of cells in row i
  32         SafeVector<PIF> data;                           // data values
  33         SafeVector<SafeVector<PIF>::iterator> rowPtrs; // pointers to the beginning of each row
  34
  35         /////////////////////////////////////////////////////////////////
  36         // SparseMatrix::SparseMatrix()
  37         //
  38         // Private constructor.
  39         /////////////////////////////////////////////////////////////////
  40
  41         SparseMatrix() {
  42         }
  43
  44 public:
  45
  46         /////////////////////////////////////////////////////////////////
  47         // SparseMatrix::SparseMatrix()
  48         //
  49         // Constructor.  Builds a sparse matrix from a posterior matrix.
  50         // Note that the expected format for the posterior matrix is as
  51         // a (seq1Length+1) x (seq2Length+1) matrix where the 0th row
  52         // and 0th column are ignored (they should contain all zeroes).
  53         /////////////////////////////////////////////////////////////////
  54
  55         SparseMatrix(int seq1Length, int seq2Length, const VF &posterior) :
  56                         seq1Length(seq1Length), seq2Length(seq2Length) {
  57
  58                 int numCells = 0;
  59
  60                 assert(seq1Length > 0);
  61                 assert(seq2Length > 0);
  62
  63                 // calculate memory required; count the number of cells in the
  64                 // posterior matrix above the threshold
  65                 VF::const_iterator postPtr = posterior.begin();
  66                 for (int i = 0; i <= seq1Length; i++) {
  67                         for (int j = 0; j <= seq2Length; j++) {
  68                                 if (*(postPtr++) >= POSTERIOR_CUTOFF) {
  69                                         assert(i != 0 && j != 0);
  70                                         numCells++;
  71                                 }
  72                         }
  73                 }
  74
  75                 // allocate memory
  76                 data.resize(numCells);
  77                 rowSize.resize(seq1Length + 1);
  78                 rowSize[0] = -1;
  79                 rowPtrs.resize(seq1Length + 1);
  80                 rowPtrs[0] = data.end();
  81
  82                 // build sparse matrix
  83                 postPtr = posterior.begin() + seq2Length + 1; // note that we're skipping the first row here
  84                 SafeVector<PIF>::iterator dataPtr = data.begin();
  85                 for (int i = 1; i <= seq1Length; i++) {
  86                         postPtr++;              // and skipping the first column of each row
  87                         rowPtrs[i] = dataPtr;
  88                         for (int j = 1; j <= seq2Length; j++) {
  89                                 if (*postPtr >= POSTERIOR_CUTOFF) {
  90                                         dataPtr->first = j;
  91                                         dataPtr->second = *postPtr;
  92                                         dataPtr++;
  93                                 }
  94                                 postPtr++;
  95                         }
  96                         rowSize[i] = dataPtr - rowPtrs[i];
  97                 }
  98         }
  99
 100         /////////////////////////////////////////////////////////////////
 101         // SparseMatrix::GetRowPtr()
 102         //
 103         // Returns the pointer to a particular row in the sparse matrix.
 104         /////////////////////////////////////////////////////////////////
 105
 106         SafeVector<PIF>::iterator GetRowPtr(int row) const {
 107                 assert(row >= 1 && row <= seq1Length);
 108                 return rowPtrs[row];
 109         }
 110
 111         /////////////////////////////////////////////////////////////////
 112         // SparseMatrix::GetValue()
 113         //
 114         // Returns value at a particular row, column.
 115         /////////////////////////////////////////////////////////////////
 116
 117         float GetValue(int row, int col) {
 118                 assert(row >= 1 && row <= seq1Length);
 119                 assert(col >= 1 && col <= seq2Length);
 120                 for (int i = 0; i < rowSize[row]; i++) {
 121                         if (rowPtrs[row][i].first == col)
 122                                 return rowPtrs[row][i].second;
 123                 }
 124                 return 0;
 125         }
 126
 127         /////////////////////////////////////////////////////////////////
 128         // SparseMatrix::GetRowSize()
 129         //
 130         // Returns the number of entries in a particular row.
 131         /////////////////////////////////////////////////////////////////
 132
 133         int GetRowSize(int row) const {
 134                 assert(row >= 1 && row <= seq1Length);
 135                 return rowSize[row];
 136         }
 137
 138         /////////////////////////////////////////////////////////////////
 139         // SparseMatrix::GetSeq1Length()
 140         //
 141         // Returns the first dimension of the matrix.
 142         /////////////////////////////////////////////////////////////////
 143
 144         int GetSeq1Length() const {
 145                 return seq1Length;
 146         }
 147
 148         /////////////////////////////////////////////////////////////////
 149         // SparseMatrix::GetSeq2Length()
 150         //
 151         // Returns the second dimension of the matrix.
 152         /////////////////////////////////////////////////////////////////
 153
 154         int GetSeq2Length() const {
 155                 return seq2Length;
 156         }
 157
 158         /////////////////////////////////////////////////////////////////
 159         // SparseMatrix::GetRowPtr
 160         //
 161         // Returns the pointer to a particular row in the sparse matrix.
 162         /////////////////////////////////////////////////////////////////
 163
 164         int GetNumCells() const {
 165                 return data.size();
 166         }
 167
 168         /////////////////////////////////////////////////////////////////
 169         // SparseMatrix::Print()
 170         //
 171         // Prints out a sparse matrix.
 172         /////////////////////////////////////////////////////////////////
 173
 174         void Print(ostream &outfile) const {
 175                 outfile << "Sparse Matrix:" << endl;
 176                 for (int i = 1; i <= seq1Length; i++) {
 177                         outfile << "  " << i << ":";
 178                         for (int j = 0; j < rowSize[i]; j++) {
 179                                 outfile << " (" << rowPtrs[i][j].first << ","
 180                                                 << rowPtrs[i][j].second << ")";
 181                         }
 182                         outfile << endl;
 183                 }
 184         }
 185
 186         /////////////////////////////////////////////////////////////////
 187         // SparseMatrix::ComputeTranspose()
 188         //
 189         // Returns a new sparse matrix containing the transpose of the
 190         // current matrix.
 191         /////////////////////////////////////////////////////////////////
 192
 193         SparseMatrix *ComputeTranspose() const {
 194
 195                 // create a new sparse matrix
 196                 SparseMatrix *ret = new SparseMatrix();
 197                 int numCells = data.size();
 198
 199                 ret->seq1Length = seq2Length;
 200                 ret->seq2Length = seq1Length;
 201
 202                 // allocate memory
 203                 ret->data.resize(numCells);
 204                 ret->rowSize.resize(seq2Length + 1);
 205                 ret->rowSize[0] = -1;
 206                 ret->rowPtrs.resize(seq2Length + 1);
 207                 ret->rowPtrs[0] = ret->data.end();
 208
 209                 // compute row sizes
 210                 for (int i = 1; i <= seq2Length; i++)
 211                         ret->rowSize[i] = 0;
 212                 for (int i = 0; i < numCells; i++)
 213                         ret->rowSize[data[i].first]++;
 214
 215                 // compute row ptrs
 216                 for (int i = 1; i <= seq2Length; i++) {
 217                         ret->rowPtrs[i] =
 218                                         (i == 1) ?
 219                                                         ret->data.begin() :
 220                                                         ret->rowPtrs[i - 1] + ret->rowSize[i - 1];
 221                 }
 222
 223                 // now fill in data
 224                 SafeVector<SafeVector<PIF>::iterator> currPtrs = ret->rowPtrs;
 225
 226                 for (int i = 1; i <= seq1Length; i++) {
 227                         SafeVector<PIF>::iterator row = rowPtrs[i];
 228                         for (int j = 0; j < rowSize[i]; j++) {
 229                                 currPtrs[row[j].first]->first = i;
 230                                 currPtrs[row[j].first]->second = row[j].second;
 231                                 currPtrs[row[j].first]++;
 232                         }
 233                 }
 234
 235                 return ret;
 236         }
 237
 238         /////////////////////////////////////////////////////////////////
 239         // SparseMatrix::GetPosterior()
 240         //
 241         // Return the posterior representation of the sparse matrix.
 242         /////////////////////////////////////////////////////////////////
 243
 244         VF *GetPosterior() const {
 245
 246                 // create a new posterior matrix
 247                 VF *posteriorPtr = new VF((seq1Length + 1) * (seq2Length + 1));
 248                 assert(posteriorPtr);
 249                 VF &posterior = *posteriorPtr;
 250
 251                 // build the posterior matrix
 252                 for (int i = 0; i < (seq1Length + 1) * (seq2Length + 1); i++)
 253                         posterior[i] = 0;
 254                 for (int i = 1; i <= seq1Length; i++) {
 255                         VF::iterator postPtr = posterior.begin() + i * (seq2Length + 1);
 256                         for (int j = 0; j < rowSize[i]; j++) {
 257                                 postPtr[rowPtrs[i][j].first] = rowPtrs[i][j].second;
 258                         }
 259                 }
 260
 261                 return posteriorPtr;
 262         }
 263
 264 };
 265
 266 #endif