binaries/src/mafft/extensions/mxscarna_src/probconsRNA/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 #include "SafeVector.h"
  12 #include "nrutil.h"
  13
  14 using namespace std;
  15
  16 const float POSTERIOR_CUTOFF = 0.01;         // minimum posterior probability
  17                                              // value that is maintained in the
  18                                              // sparse matrix representation
  19
  20 typedef pair<int,float> PIF;                 // Sparse matrix entry type
  21                                              //   first --> column
  22                                              //   second --> value
  23
  24 namespace MXSCARNA {
  25 struct PIF2 {        // Sparse matrix entry type
  26     int i;
  27     int j;
  28     float prob;
  29 };
  30 }
  31
  32 /////////////////////////////////////////////////////////////////
  33 // SparseMatrix
  34 //
  35 // Class for sparse matrix computations
  36 /////////////////////////////////////////////////////////////////
  37 namespace MXSCARNA {
  38 class SparseMatrix {
  39
  40   int seq1Length, seq2Length;                     // dimensions of matrix
  41   VI rowSize;                                     // rowSize[i] = # of cells in row i
  42   SafeVector<PIF> data;                           // data values
  43   SafeVector<SafeVector<PIF>::iterator> rowPtrs;  // pointers to the beginning of each row
  44
  45  public:
  46   SafeVector<PIF2> data2;
  47   /////////////////////////////////////////////////////////////////
  48   // SparseMatrix::SparseMatrix()
  49   //
  50   // Private constructor.1
  51   /////////////////////////////////////////////////////////////////
  52   SparseMatrix() { }
  53
  54   /////////////////////////////////////////////////////////////////
  55   // SparseMatrix::SparseMatrix()
  56   //
  57   // Constructor.  Builds a sparse matrix from a posterior matrix.
  58   // Note that the expected format for the posterior matrix is as
  59   // a (seq1Length+1) x (seq2Length+1) matrix where the 0th row
  60   // and 0th column are ignored (they should contain all zeroes).
  61   /////////////////////////////////////////////////////////////////
  62
  63   SparseMatrix (int seq1Length, int seq2Length, const VF &posterior) :
  64     seq1Length (seq1Length), seq2Length (seq2Length) {
  65
  66     int numCells = 0;
  67
  68     assert (seq1Length > 0);
  69     assert (seq2Length > 0);
  70
  71     // calculate memory required; count the number of cells in the
  72     // posterior matrix above the threshold
  73     VF::const_iterator postPtr = posterior.begin();
  74     for (int i = 0; i <= seq1Length; i++){
  75       for (int j = 0; j <= seq2Length; j++){
  76         if (*(postPtr++) >= POSTERIOR_CUTOFF){
  77           assert (i != 0 && j != 0);
  78           numCells++;
  79         }
  80       }
  81     }
  82
  83     // allocate memory
  84     data.resize(numCells);
  85     rowSize.resize (seq1Length + 1); rowSize[0] = -1;
  86     rowPtrs.resize (seq1Length + 1); rowPtrs[0] = data.end();
  87
  88     // build sparse matrix
  89     postPtr = posterior.begin() + seq2Length + 1;           // note that we're skipping the first row here
  90     SafeVector<PIF>::iterator dataPtr = data.begin();
  91     for (int i = 1; i <= seq1Length; i++){
  92       postPtr++;                                            // and skipping the first column of each row
  93       rowPtrs[i] = dataPtr;
  94       for (int j = 1; j <= seq2Length; j++){
  95         if (*postPtr >= POSTERIOR_CUTOFF){
  96           dataPtr->first = j;
  97           dataPtr->second = *postPtr;
  98           dataPtr++;
  99         }
 100         postPtr++;
 101       }
 102       rowSize[i] = dataPtr - rowPtrs[i];
 103     }
 104   }
 105
 106   //////////////////////////////////////////////////////////////////////////
 107   // SparseMatrix::SetSparseMatrix()
 108   //
 109   // Constructor.
 110   //////////////////////////////////////////////////////////////////////////
 111   void SetSparseMatrix(int inseq1Length, int inseq2Length, const Trimat<float> &bppMat, float cutoff = 0.01) {
 112       seq1Length = inseq1Length;
 113       seq2Length = inseq2Length;
 114
 115       int numCells = 0;
 116
 117       assert (seq1Length > 0);
 118       assert (seq2Length > 0);
 119
 120       data.clear();
 121       rowSize.clear();
 122       rowPtrs.clear();
 123       for (int i = 1; i <= seq1Length; i++) {
 124           for (int j = i; j <= seq2Length; j++) {
 125               if (bppMat.ref(i, j) >= cutoff ) {
 126                   numCells++;
 127               }
 128           }
 129       }
 130
 131       // allocate memory
 132       data.resize(numCells);
 133       for (int i = 0; i < numCells; i++) {
 134           data[i].first  = 0;
 135           data[i].second = 0;
 136       }
 137       rowSize.resize (seq1Length + 1); rowSize[0] = -1;
 138       rowPtrs.resize (seq1Length + 1); rowPtrs[0] = data.end();
 139
 140       SafeVector<PIF>::iterator dataPtr = data.begin();
 141       for (int i = 1; i <= seq1Length; i++) {
 142           rowPtrs[i] = dataPtr;
 143           for (int j = i; j <= seq2Length; j++) {
 144               if ( bppMat.ref(i, j) >= cutoff ) {
 145                   dataPtr->first = j;
 146                   dataPtr->second = bppMat.ref(i, j);
 147                   dataPtr++;
 148               }
 149           }
 150           rowSize[i] = dataPtr - rowPtrs[i];
 151       }
 152
 153       float tmp;
 154       for(int k = 1; k <= seq1Length; k++) {
 155           for(int m = k, n = k; n <= k + 300 && m >= 1 && n <= seq2Length; m--, n++) {
 156               if ((tmp = GetValue(m, n)) > 0) {
 157                   PIF2 p;
 158                   p.i    = m;
 159                   p.j    = n;
 160                   p.prob = tmp;
 161                   data2.push_back(p);
 162               }
 163           }
 164
 165           for(int m = k, n = k + 1; n <= k + 300 && m >= 1 && n <= seq2Length; m--, n++) {
 166               if ((tmp = GetValue(m, n)) > 0) {
 167                   PIF2 p;
 168                   p.i    = m;
 169                   p.j    = n;
 170                   p.prob = tmp;
 171                   data2.push_back(p);
 172               }
 173           }
 174       }
 175   }
 176
 177   /////////////////////////////////////////////////////////////////
 178   // SparseMatrix::GetRowPtr()
 179   //
 180   // Returns the pointer to a particular row in the sparse matrix.
 181   /////////////////////////////////////////////////////////////////
 182
 183   SafeVector<PIF>::iterator GetRowPtr (int row) const {
 184     assert (row >= 1 && row <= seq1Length);
 185     return rowPtrs[row];
 186   }
 187
 188   /////////////////////////////////////////////////////////////////
 189   // SparseMatrix::GetValue()
 190   //
 191   // Returns value at a particular row, column.
 192   /////////////////////////////////////////////////////////////////
 193
 194   float GetValue (int row, int col){
 195     assert (row >= 1 && row <= seq1Length);
 196     assert (col >= 1 && col <= seq2Length);
 197     for (int i = 0; i < rowSize[row]; i++){
 198       if (rowPtrs[row][i].first == col) return rowPtrs[row][i].second;
 199     }
 200     return 0;
 201   }
 202
 203   void SetValue(int row, int col, float value) {
 204     assert (row >= 1 && row <= seq1Length);
 205     assert (col >= 1 && col <= seq2Length);
 206     for (int i = 0; i < rowSize[row]; i++){
 207       if (rowPtrs[row][i].first == col) rowPtrs[row][i].second = value;
 208     }
 209   }
 210
 211   /////////////////////////////////////////////////////////////////
 212   // SparseMatrix::GetRowSize()
 213   //
 214   // Returns the number of entries in a particular row.
 215   /////////////////////////////////////////////////////////////////
 216
 217   int GetRowSize (int row) const {
 218     assert (row >= 1 && row <= seq1Length);
 219     return rowSize[row];
 220   }
 221
 222   /////////////////////////////////////////////////////////////////
 223   // SparseMatrix::GetSeq1Length()
 224   //
 225   // Returns the first dimension of the matrix.
 226   /////////////////////////////////////////////////////////////////
 227
 228   int GetSeq1Length () const {
 229     return seq1Length;
 230   }
 231
 232   /////////////////////////////////////////////////////////////////
 233   // SparseMatrix::GetSeq2Length()
 234   //
 235   // Returns the second dimension of the matrix.
 236   /////////////////////////////////////////////////////////////////
 237
 238   int GetSeq2Length () const {
 239     return seq2Length;
 240   }
 241
 242   /////////////////////////////////////////////////////////////////
 243   // SparseMatrix::GetRowPtr
 244   //
 245   // Returns the pointer to a particular row in the sparse matrix.
 246   /////////////////////////////////////////////////////////////////
 247
 248   int GetNumCells () const {
 249     return data.size();
 250   }
 251
 252   /////////////////////////////////////////////////////////////////
 253   // SparseMatrix::Print()
 254   //
 255   // Prints out a sparse matrix.
 256   /////////////////////////////////////////////////////////////////
 257
 258   void Print (ostream &outfile) const {
 259     outfile << "Sparse Matrix:" << endl;
 260     for (int i = 1; i <= seq1Length; i++){
 261       outfile << "  " << i << ":";
 262       for (int j = 0; j < rowSize[i]; j++){
 263         outfile << " (" << rowPtrs[i][j].first << "," << rowPtrs[i][j].second << ")";
 264       }
 265       outfile << endl;
 266     }
 267   }
 268
 269   /////////////////////////////////////////////////////////////////
 270   // SparseMatrix::ComputeTranspose()
 271   //
 272   // Returns a new sparse matrix containing the transpose of the
 273   // current matrix.
 274   /////////////////////////////////////////////////////////////////
 275
 276   SparseMatrix *ComputeTranspose () const {
 277
 278     // create a new sparse matrix
 279     SparseMatrix *ret = new SparseMatrix();
 280     int numCells = data.size();
 281
 282     ret->seq1Length = seq2Length;
 283     ret->seq2Length = seq1Length;
 284
 285     // allocate memory
 286     ret->data.resize (numCells);
 287     ret->rowSize.resize (seq2Length + 1); ret->rowSize[0] = -1;
 288     ret->rowPtrs.resize (seq2Length + 1); ret->rowPtrs[0] = ret->data.end();
 289
 290     // compute row sizes
 291     for (int i = 1; i <= seq2Length; i++) ret->rowSize[i] = 0;
 292     for (int i = 0; i < numCells; i++)
 293       ret->rowSize[data[i].first]++;
 294
 295     // compute row ptrs
 296     for (int i = 1; i <= seq2Length; i++){
 297       ret->rowPtrs[i] = (i == 1) ? ret->data.begin() : ret->rowPtrs[i-1] + ret->rowSize[i-1];
 298     }
 299
 300     // now fill in data
 301     SafeVector<SafeVector<PIF>::iterator> currPtrs = ret->rowPtrs;
 302
 303     for (int i = 1; i <= seq1Length; i++){
 304       SafeVector<PIF>::iterator row = rowPtrs[i];
 305       for (int j = 0; j < rowSize[i]; j++){
 306         currPtrs[row[j].first]->first = i;
 307         currPtrs[row[j].first]->second = row[j].second;
 308         currPtrs[row[j].first]++;
 309       }
 310     }
 311
 312     return ret;
 313   }
 314
 315   /////////////////////////////////////////////////////////////////
 316   // SparseMatrix::GetPosterior()
 317   //
 318   // Return the posterior representation of the sparse matrix.
 319   /////////////////////////////////////////////////////////////////
 320
 321   VF *GetPosterior () const {
 322
 323     // create a new posterior matrix
 324     VF *posteriorPtr = new VF((seq1Length+1) * (seq2Length+1)); assert (posteriorPtr);
 325     VF &posterior = *posteriorPtr;
 326
 327     // build the posterior matrix
 328     for (int i = 0; i < (seq1Length+1) * (seq2Length+1); i++) posterior[i] = 0;
 329     for (int i = 1; i <= seq1Length; i++){
 330       VF::iterator postPtr = posterior.begin() + i * (seq2Length+1);
 331       for (int j = 0; j < rowSize[i]; j++){
 332         postPtr[rowPtrs[i][j].first] = rowPtrs[i][j].second;
 333       }
 334     }
 335
 336     return posteriorPtr;
 337   }
 338
 339 };
 340 }
 341 #endif