--- /dev/null
+/////////////////////////////////////////////////////////////////
+// SparseMatrix.h
+//
+// Sparse matrix computations
+/////////////////////////////////////////////////////////////////
+
+#ifndef SPARSEMATRIX_H
+#define SPARSEMATRIX_H
+
+#include <iostream>
+
+using namespace std;
+
+const float POSTERIOR_CUTOFF = 0.01; // minimum posterior probability
+// value that is maintained in the
+// sparse matrix representation
+
+typedef pair<int, float> PIF; // Sparse matrix entry type
+// first --> column
+// second --> value
+
+/////////////////////////////////////////////////////////////////
+// SparseMatrix
+//
+// Class for sparse matrix computations
+/////////////////////////////////////////////////////////////////
+
+class SparseMatrix {
+
+ int seq1Length, seq2Length; // dimensions of matrix
+ VI rowSize; // rowSize[i] = # of cells in row i
+ SafeVector<PIF> data; // data values
+ SafeVector<SafeVector<PIF>::iterator> rowPtrs; // pointers to the beginning of each row
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::SparseMatrix()
+ //
+ // Private constructor.
+ /////////////////////////////////////////////////////////////////
+
+ SparseMatrix() {
+ }
+
+public:
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::SparseMatrix()
+ //
+ // Constructor. Builds a sparse matrix from a posterior matrix.
+ // Note that the expected format for the posterior matrix is as
+ // a (seq1Length+1) x (seq2Length+1) matrix where the 0th row
+ // and 0th column are ignored (they should contain all zeroes).
+ /////////////////////////////////////////////////////////////////
+
+ SparseMatrix(int seq1Length, int seq2Length, const VF &posterior) :
+ seq1Length(seq1Length), seq2Length(seq2Length) {
+
+ int numCells = 0;
+
+ assert(seq1Length > 0);
+ assert(seq2Length > 0);
+
+ // calculate memory required; count the number of cells in the
+ // posterior matrix above the threshold
+ VF::const_iterator postPtr = posterior.begin();
+ for (int i = 0; i <= seq1Length; i++) {
+ for (int j = 0; j <= seq2Length; j++) {
+ if (*(postPtr++) >= POSTERIOR_CUTOFF) {
+ assert(i != 0 && j != 0);
+ numCells++;
+ }
+ }
+ }
+
+ // allocate memory
+ data.resize(numCells);
+ rowSize.resize(seq1Length + 1);
+ rowSize[0] = -1;
+ rowPtrs.resize(seq1Length + 1);
+ rowPtrs[0] = data.end();
+
+ // build sparse matrix
+ postPtr = posterior.begin() + seq2Length + 1; // note that we're skipping the first row here
+ SafeVector<PIF>::iterator dataPtr = data.begin();
+ for (int i = 1; i <= seq1Length; i++) {
+ postPtr++; // and skipping the first column of each row
+ rowPtrs[i] = dataPtr;
+ for (int j = 1; j <= seq2Length; j++) {
+ if (*postPtr >= POSTERIOR_CUTOFF) {
+ dataPtr->first = j;
+ dataPtr->second = *postPtr;
+ dataPtr++;
+ }
+ postPtr++;
+ }
+ rowSize[i] = dataPtr - rowPtrs[i];
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetRowPtr()
+ //
+ // Returns the pointer to a particular row in the sparse matrix.
+ /////////////////////////////////////////////////////////////////
+
+ SafeVector<PIF>::iterator GetRowPtr(int row) const {
+ assert(row >= 1 && row <= seq1Length);
+ return rowPtrs[row];
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetValue()
+ //
+ // Returns value at a particular row, column.
+ /////////////////////////////////////////////////////////////////
+
+ float GetValue(int row, int col) {
+ assert(row >= 1 && row <= seq1Length);
+ assert(col >= 1 && col <= seq2Length);
+ for (int i = 0; i < rowSize[row]; i++) {
+ if (rowPtrs[row][i].first == col)
+ return rowPtrs[row][i].second;
+ }
+ return 0;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetRowSize()
+ //
+ // Returns the number of entries in a particular row.
+ /////////////////////////////////////////////////////////////////
+
+ int GetRowSize(int row) const {
+ assert(row >= 1 && row <= seq1Length);
+ return rowSize[row];
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetSeq1Length()
+ //
+ // Returns the first dimension of the matrix.
+ /////////////////////////////////////////////////////////////////
+
+ int GetSeq1Length() const {
+ return seq1Length;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetSeq2Length()
+ //
+ // Returns the second dimension of the matrix.
+ /////////////////////////////////////////////////////////////////
+
+ int GetSeq2Length() const {
+ return seq2Length;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetRowPtr
+ //
+ // Returns the pointer to a particular row in the sparse matrix.
+ /////////////////////////////////////////////////////////////////
+
+ int GetNumCells() const {
+ return data.size();
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::Print()
+ //
+ // Prints out a sparse matrix.
+ /////////////////////////////////////////////////////////////////
+
+ void Print(ostream &outfile) const {
+ outfile << "Sparse Matrix:" << endl;
+ for (int i = 1; i <= seq1Length; i++) {
+ outfile << " " << i << ":";
+ for (int j = 0; j < rowSize[i]; j++) {
+ outfile << " (" << rowPtrs[i][j].first << ","
+ << rowPtrs[i][j].second << ")";
+ }
+ outfile << endl;
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::ComputeTranspose()
+ //
+ // Returns a new sparse matrix containing the transpose of the
+ // current matrix.
+ /////////////////////////////////////////////////////////////////
+
+ SparseMatrix *ComputeTranspose() const {
+
+ // create a new sparse matrix
+ SparseMatrix *ret = new SparseMatrix();
+ int numCells = data.size();
+
+ ret->seq1Length = seq2Length;
+ ret->seq2Length = seq1Length;
+
+ // allocate memory
+ ret->data.resize(numCells);
+ ret->rowSize.resize(seq2Length + 1);
+ ret->rowSize[0] = -1;
+ ret->rowPtrs.resize(seq2Length + 1);
+ ret->rowPtrs[0] = ret->data.end();
+
+ // compute row sizes
+ for (int i = 1; i <= seq2Length; i++)
+ ret->rowSize[i] = 0;
+ for (int i = 0; i < numCells; i++)
+ ret->rowSize[data[i].first]++;
+
+ // compute row ptrs
+ for (int i = 1; i <= seq2Length; i++) {
+ ret->rowPtrs[i] =
+ (i == 1) ?
+ ret->data.begin() :
+ ret->rowPtrs[i - 1] + ret->rowSize[i - 1];
+ }
+
+ // now fill in data
+ SafeVector<SafeVector<PIF>::iterator> currPtrs = ret->rowPtrs;
+
+ for (int i = 1; i <= seq1Length; i++) {
+ SafeVector<PIF>::iterator row = rowPtrs[i];
+ for (int j = 0; j < rowSize[i]; j++) {
+ currPtrs[row[j].first]->first = i;
+ currPtrs[row[j].first]->second = row[j].second;
+ currPtrs[row[j].first]++;
+ }
+ }
+
+ return ret;
+ }
+
+ /////////////////////////////////////////////////////////////////
+ // SparseMatrix::GetPosterior()
+ //
+ // Return the posterior representation of the sparse matrix.
+ /////////////////////////////////////////////////////////////////
+
+ VF *GetPosterior() const {
+
+ // create a new posterior matrix
+ VF *posteriorPtr = new VF((seq1Length + 1) * (seq2Length + 1));
+ assert(posteriorPtr);
+ VF &posterior = *posteriorPtr;
+
+ // build the posterior matrix
+ for (int i = 0; i < (seq1Length + 1) * (seq2Length + 1); i++)
+ posterior[i] = 0;
+ for (int i = 1; i <= seq1Length; i++) {
+ VF::iterator postPtr = posterior.begin() + i * (seq2Length + 1);
+ for (int j = 0; j < rowSize[i]; j++) {
+ postPtr[rowPtrs[i][j].first] = rowPtrs[i][j].second;
+ }
+ }
+
+ return posteriorPtr;
+ }
+
+};
+
+#endif
git://source.jalview.org / jabaws.git / blobdiff