SparseBlockMatrix.cpp

#include "SparseBlockMatrix.h"

// std lib
#include <iostream>
#include <iomanip>

// qt lib
#include <QDataStream>
#include <QDebug>
#include <QMapIterator>
#include <QListIterator>

// boost lib
#include <boost/numeric/ublas/matrix_sparse.hpp>
#include <boost/numeric/ublas/matrix_proxy.hpp>
#include <boost/numeric/ublas/io.hpp>

// Isis lib
#include "IString.h"
#include "LinearAlgebra.h"

using namespace boost::numeric::ublas;

namespace Isis {

  SparseBlockColumnMatrix::~SparseBlockColumnMatrix() {
    wipe();
  }


  void SparseBlockColumnMatrix::wipe() {
    qDeleteAll(values());
    clear();
  }


  SparseBlockColumnMatrix::SparseBlockColumnMatrix(const SparseBlockColumnMatrix& src) {
    copy(src);
  }


  void SparseBlockColumnMatrix::copy(const SparseBlockColumnMatrix& src) {
    // handi-wipe
    wipe();

    // copy matrix blocks from src
    QMapIterator<int, LinearAlgebra::Matrix *> it(src);
    while ( it.hasNext() ) {
      it.next();

      // copy matrix block from src
      LinearAlgebra::Matrix *m = new LinearAlgebra::Matrix(*(it.value()));

      // insert matrix into map
      this->insert(it.key(),m);
    }
  }


  SparseBlockColumnMatrix&
      SparseBlockColumnMatrix::operator=(const SparseBlockColumnMatrix& src) {
    if ( this == &src )
      return *this;

    copy(src);

    return *this;
  }


  bool SparseBlockColumnMatrix::insertMatrixBlock(int nColumnBlock, int nRows, int nCols) {
    // check if matrix already exists at the key "nColumnBlock"
    if ( this->contains(nColumnBlock) )
      return true;

    // allocate matrix block with nRows and nCols
    LinearAlgebra::Matrix *m = new LinearAlgebra::Matrix(nRows,nCols);

    if ( !m )
      return false;

    // zero matrix elements
    m->clear();

    // insert matrix into map
    this->insert(nColumnBlock,m);

    return true;
  }


  int SparseBlockColumnMatrix::numberOfElements() {
    int nElements = 0;

    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      nElements += it.value()->size1()*it.value()->size2();
    }

    return nElements;
  }


  int SparseBlockColumnMatrix::numberOfColumns() {

    int nColumns = 0;

    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      nColumns = it.value()->size2();
      break;
    }

    return nColumns;
  }


  int SparseBlockColumnMatrix::numberOfRows() {

    // iterate to last block (the diagonal one)
    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;
    }

    int nRows = it.value()->size1();

    return nRows;
  }


  void SparseBlockColumnMatrix::print(std::ostream& outstream) {
    if ( size() == 0 ) {
      outstream << "Empty SparseBlockColumnMatrix..." << std::endl;
      return;
    }

    outstream << "Printing SparseBlockColumnMatrix..." << std::endl;
    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      if( it.value() )
        outstream << it.key() << std::endl << *(it.value()) << std::endl
                  << std::endl;
      else
        outstream << "NULL block pointer at row[" << IString(it.key())
                  << "]!" << std::endl;
    }
  }


  void SparseBlockColumnMatrix::printClean(std::ostream& outstream) {
    if ( size() == 0 ) {
      outstream << "Empty SparseBlockColumnMatrix..." << std::endl;
      return;
    }

    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      LinearAlgebra::Matrix * m = it.value();

      int rows = m->size1();
      int cols = m->size2();

      for ( int i = 0; i < rows; i++ ) {
        for ( int j = 0; j < cols; j++ ) {
          double d = m->at_element(i,j);
          if ( j == cols-1 )
            outstream << std::setprecision(12) << d << std::endl;
          else
            outstream << std::setprecision(12) << d << ",";
        }
      }

    }
    outstream << std::endl;
  }


  void SparseBlockColumnMatrix::zeroBlocks() {
    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();
      it.value()->clear();
    }
  }


  QDataStream &operator<<(QDataStream &stream, const SparseBlockColumnMatrix &sbcm) {
    // write number of blocks in this column
    int nBlocks = sbcm.size();
    stream << (qint32)nBlocks;

    QMapIterator<int, LinearAlgebra::Matrix *> it(sbcm);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      int nRows = it.value()->size1();
      int nCols = it.value()->size2();

      // write block number (key); rows (size1); and columns (size2)
      stream << it.key() << (qint32)nRows << (qint32)nCols;

      double* data = &it.value()->data()[0];

      // write raw matrix data
      stream.writeRawData((const char*)data, nRows*nCols*sizeof(double));
    }

    return stream;
  }


  QDataStream &operator>>(QDataStream &stream, SparseBlockColumnMatrix &sbcm) {
    qint32 nBlocks, nBlockNumber, nRows, nCols;
    int i, r, c;

    stream >> nBlocks;

    for ( i = 0; i < nBlocks; i++ ) {
      // read block number (key); rows (size1); and columns (size2)
      stream >> nBlockNumber >> nRows >> nCols;

      double data[nRows*nCols];

      // read raw matrix data
      stream.readRawData((char*)data, nRows*nCols*sizeof(double));

      // insert matrix at correct key
      sbcm.insertMatrixBlock(nBlockNumber, nRows, nCols);

      // get matrix
      LinearAlgebra::Matrix *matrix = sbcm[nBlockNumber];

      // fill with data
      for ( r = 0; r < nRows; r++ ) {
        for ( c = 0; c < nCols; c++ ) {
          int nLocation = r*nRows + c;
          (*matrix)(r,c) = data[nLocation];
        }
      }
    }

    return stream;
  }


  QDebug operator<<(QDebug dbg, const SparseBlockColumnMatrix &sbcm) {
    dbg.space() << "New Block" << endl;

    QMapIterator<int, LinearAlgebra::Matrix *> it(sbcm);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      // get matrix
      LinearAlgebra::Matrix *matrix = it.value();

      // matrix rows, columns
      int nRows = matrix->size1();
      int nCols = matrix->size2();

      dbg.nospace() << qSetFieldWidth(4);
      dbg.nospace() << qSetRealNumberPrecision(8);

      for ( int r = 0; r < nRows; r++ ) {
        for ( int c = 0; c < nCols; c++ ) {
          dbg.space() << (*matrix)(r,c);
        }
        dbg.space() << endl;
      }

      dbg.space() << endl;
    }

    return dbg;
  }


  // SparseBlockRowMatrix methods

  SparseBlockRowMatrix::~SparseBlockRowMatrix() {
    wipe();
  }


  void SparseBlockRowMatrix::wipe() {
    qDeleteAll(values());
    clear();
  }


  SparseBlockRowMatrix::SparseBlockRowMatrix(const SparseBlockRowMatrix& src) {
    copy(src);
  }


  void SparseBlockRowMatrix::copy(const SparseBlockRowMatrix& src) {
    // handi-wipe
    wipe();

    // copy matrix blocks from src
    QMapIterator<int, LinearAlgebra::Matrix *> it(src);
    while ( it.hasNext() ) {
      it.next();

      // copy matrix block from src
      LinearAlgebra::Matrix *m = new LinearAlgebra::Matrix(*(it.value()));

      // insert matrix into map
      this->insert(it.key(),m);
    }
  }


  SparseBlockRowMatrix&
      SparseBlockRowMatrix::operator=(const SparseBlockRowMatrix& src) {
    if ( this == &src )
      return *this;

    copy(src);

    return *this;
  }


  bool SparseBlockRowMatrix::insertMatrixBlock(int nRowBlock, int nRows, int nCols) {
    if ( this->contains(nRowBlock) )
      return false;

    LinearAlgebra::Matrix *m = new LinearAlgebra::Matrix(nRows,nCols);

    if ( !m )
      return false;

    m->clear();

    this->insert(nRowBlock,m);

    return true;
  }


  int SparseBlockRowMatrix::numberOfElements() {
    int nElements = 0;

    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      nElements += it.value()->size1()*it.value()->size2();
    }

    return nElements;
  }


  void SparseBlockRowMatrix::print(std::ostream& outstream) {
    if ( size() == 0 ) {
      outstream << "Empty SparseBlockRowMatrix..." << std::endl;
      return;
    }

    outstream << "Printing SparseBlockRowMatrix..." << std::endl;
    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      if( it.value() )
        outstream << it.key() << std::endl << *(it.value()) << std::endl
                  << std::endl;
      else
        outstream << "NULL block pointer at column[" << IString(it.key())
                  << "]!" << std::endl;
    }
  }


  void SparseBlockRowMatrix::printClean(std::ostream& outstream) {
    if ( size() == 0 ) {
      outstream << "Empty SparseBlockRowMatrix..." << std::endl;
      return;
    }

    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      LinearAlgebra::Matrix *m = it.value();

      int rows = m->size1();
      int cols = m->size2();

      for ( int i = 0; i < rows; i++ ) {
        for ( int j = 0; j < cols; j++ ) {
          double d = m->at_element(i,j);
          if ( j == cols-1 )
            outstream << std::setprecision(9) << d << std::endl;
          else
            outstream << std::setprecision(9) << d << ",";
        }
      }

    }
    outstream << std::endl;
  }


  void SparseBlockRowMatrix::zeroBlocks() {
    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();
      it.value()->clear();
    }
  }


  void SparseBlockRowMatrix::copyToBoost(compressed_matrix<double>& B) {
    B.clear();

    int ncols, nstart, nend, nrowBlock;
    range rRow = range(0,3);
    range rCol;

    QMapIterator<int, LinearAlgebra::Matrix *> it(*this);
    while ( it.hasNext() ) {
      it.next();

      nrowBlock = it.key();
      LinearAlgebra::Matrix *m = it.value();

      ncols = m->size2();

      nstart = nrowBlock*ncols;
      nend = nstart + ncols;

      rCol = range(nstart,nend);

      matrix_range<compressed_matrix<double> > m1 (B, rRow, rCol);

      m1 = *m;
    }
  }


    int SparseBlockRowMatrix::getLeadingColumnsForBlock(int nblockColumn) {

      if ( nblockColumn == 0 )
        return 0;

      int nLeadingColumnsElements = 0;

      int nCol = 0;

      while ( nCol < nblockColumn ) {
        if ( !(*this)[nCol] ) {
          nCol++;
          continue;
        }

        int ncolumns = (*this)[nCol]->size2();

        if ( ncolumns == -1 )
          continue;

        nLeadingColumnsElements += ncolumns;

        nCol++;
      }

      return nLeadingColumnsElements;
    }


  QDataStream &operator<<(QDataStream &stream, const SparseBlockRowMatrix &sbrm) {
    // write number of blocks in this column
    int nBlocks = sbrm.size();
    stream << (qint32)nBlocks;

    QMapIterator<int, LinearAlgebra::Matrix *> it(sbrm);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      int nRows = it.value()->size1();
      int nCols = it.value()->size2();

      // write block number (key); rows (size1); and columns (size2)
      stream << it.key() << (qint32)nRows << (qint32)nCols;

      double* data = &it.value()->data()[0];

      // write raw matrix data
      stream.writeRawData((const char*)data, nRows*nCols*sizeof(double));
    }

    return stream;
  }


  QDataStream &operator>>(QDataStream &stream, SparseBlockRowMatrix &sbrm) {
    qint32 nBlocks, nBlockNumber, nRows, nCols;
    int i, r, c;

    stream >> nBlocks;

    for ( i = 0; i < nBlocks; i++ ) {
      // read block number (key); rows (size1); and columns (size2)
      stream >> nBlockNumber >> nRows >> nCols;

      double data[nRows*nCols];

      // read raw matrix data
      stream.readRawData((char*)data, nRows*nCols*sizeof(double));

      // insert matrix at correct key
      sbrm.insertMatrixBlock(nBlockNumber, nRows, nCols);

      // get matrix
      LinearAlgebra::Matrix *matrix = sbrm[nBlockNumber];

      // fill with data
      for ( r = 0; r < nRows; r++ ) {
        for ( c = 0; c < nCols; c++ ) {
          int nLocation = r*nRows + c;
          (*matrix)(r,c) = data[nLocation];
        }
      }
    }

    return stream;
  }


  QDebug operator<<(QDebug dbg, const SparseBlockRowMatrix &sbrm) {
    dbg.space() << "New Block" << endl;

    QMapIterator<int, LinearAlgebra::Matrix *> it(sbrm);
    while ( it.hasNext() ) {
      it.next();

      if( !it.value() )
        continue;

      // get matrix
      LinearAlgebra::Matrix *matrix = it.value();

      // matrix rows, columns
      int nRows = matrix->size1();
      int nCols = matrix->size2();

      dbg.nospace() << qSetFieldWidth(4);
      dbg.nospace() << qSetRealNumberPrecision(8);

      for ( int r = 0; r < nRows; r++ ) {
        for ( int c = 0; c < nCols; c++ ) {
          dbg.space() << (*matrix)(r,c);
        }
        dbg.space() << endl;
      }
      dbg.space() << endl;
    }

    return dbg;
  }


  // SparseBlockMatrix methods

  SparseBlockMatrix::~SparseBlockMatrix() {
    wipe();
  }


  void SparseBlockMatrix::wipe() {
    qDeleteAll(*this);
    clear();
  }


  SparseBlockMatrix::SparseBlockMatrix(const SparseBlockMatrix& src) {
    copy(src);
  }


  void SparseBlockMatrix::copy(const SparseBlockMatrix& src) {
    // handi-wipe
    wipe();

    // copy all SparseBlockColumnMatrix objects from src
    for( int i = 0; i < src.size(); i++ ) {
      append( new SparseBlockColumnMatrix(*(src.at(i))));
    }
  }


  SparseBlockMatrix& SparseBlockMatrix::operator=(const SparseBlockMatrix& src) {
    if ( this == &src )
      return *this;

    copy(src);

    return *this;
  }


  bool SparseBlockMatrix::setNumberOfColumns( int n ) {

    for( int i = 0; i < n; i++ )
      append( new SparseBlockColumnMatrix() );

    return true;
  }


  bool SparseBlockMatrix::insertMatrixBlock(int nColumnBlock, int nRowBlock, int nRows, int nCols) {
    return (*this)[nColumnBlock]->insertMatrixBlock(nRowBlock, nRows, nCols);
  }


  int SparseBlockMatrix::numberOfBlocks() {
    int nBlocks = 0;

    for( int i = 0; i < size(); i++ ) {
      if ( !(*this)[i] )
        continue;

      nBlocks += (*this)[i]->size();
    }

    return nBlocks;
  }


  int SparseBlockMatrix::numberOfDiagonalBlocks() {
    int ndiagBlocks = 0;

    for( int i = 0; i < size(); i++ ) {
      SparseBlockColumnMatrix* column = at(i);

      if ( !column )
        continue;

      QMapIterator<int, LinearAlgebra::Matrix *> it(*column);
      while ( it.hasNext() ) {
        it.next();

        if( it.key() == i ) {
          ndiagBlocks++;
          break;
        }
      }
    }

    return ndiagBlocks;
  }


  int SparseBlockMatrix::numberOfOffDiagonalBlocks() {
    return (numberOfBlocks() - numberOfDiagonalBlocks());
  }


  int SparseBlockMatrix::numberOfElements() {
    int nElements = 0;

    for( int i = 0; i < size(); i++ ) {
      if ( !(*this)[i] )
        continue;

      nElements += (*this)[i]->numberOfElements();
    }

    return nElements;
  }


  LinearAlgebra::Matrix *SparseBlockMatrix::getBlock(int column, int row) {
    return (*(*this)[column])[row];
  }


  void SparseBlockMatrix::zeroBlocks() {
    for ( int i = 0; i < size(); i++ )
      (*this)[i]->zeroBlocks();
  }


  void SparseBlockMatrix::print(std::ostream& outstream) {
    if ( size() == 0 ) {
      outstream << "Empty SparseBlockMatrix..." << std::endl;
      return;
    }

    outstream << "Printing SparseBlockMatrix..." << std::endl;
    for( int i = 0; i < size(); i++ ) {
      SparseBlockColumnMatrix* column = at(i);

      if ( column )
        column->print(outstream);
      else
        outstream << "NULL column pointer at column[" << IString(i)
                  << "]!" << std::endl;
    }
  }


  void SparseBlockMatrix::printClean(std::ostream& outstream) {
    if ( size() == 0 ) {
      outstream << "Empty SparseBlockMatrix..." << std::endl;
      return;
    }

    for( int i = 0; i < size(); i++ ) {
      SparseBlockColumnMatrix* column = at(i);

      if ( column )
        column->printClean(outstream);
      else
        outstream << "NULL column pointer at column[" << IString(i)
                  << "]!" << std::endl;
    }
  }


  int SparseBlockMatrix::getLeadingColumnsForBlock(int nblockColumn) {

    if ( nblockColumn == 0 )
      return 0;

    int nLeadingColumnsElements = 0;

    int nCol = 0;

    while ( nCol < nblockColumn ) {
      if ( !(*this)[nCol] )
        continue;

      int ncolumns = (*this)[nCol]->numberOfColumns();

      if ( ncolumns == -1 )
        continue;

      nLeadingColumnsElements += ncolumns;

      nCol++;
    }

    return nLeadingColumnsElements;
  }


  int SparseBlockMatrix::getLeadingRowsForBlock(int nblockRow) {

    if ( nblockRow == 0 )
      return 0;

    int i = 0;
    int nLeadingRows = 0;

    while ( i < nblockRow ) {
      SparseBlockColumnMatrix* column = at(i);

      if ( !column )
        continue;

      QMapIterator<int, LinearAlgebra::Matrix *> it(*column);
      // iterate to last element in column
      while ( it.hasNext() ) {
        it.next();

        if( it.key() == i )
          nLeadingRows += it.value()->size1();
      }
      i++;
    }

    return nLeadingRows;
  }


  QDataStream &operator<<(QDataStream &stream, const SparseBlockMatrix &sparseBlockMatrix) {
    int nBlockColumns = sparseBlockMatrix.size();

    stream << (qint32)nBlockColumns;

    for ( int i =0; i < nBlockColumns; i++ )
      stream << *sparseBlockMatrix.at(i);

    return stream;
  }


  QDataStream &operator>>(QDataStream &stream, SparseBlockMatrix &sparseBlockMatrix) {
    qint32 nBlockColumns;

    // read and set number of block columns
    stream >> nBlockColumns;
    sparseBlockMatrix.setNumberOfColumns(nBlockColumns);

    for ( int i =0; i < nBlockColumns; i++ )
      stream >> *sparseBlockMatrix.at(i);

    return stream;
  }

  QDebug operator<<(QDebug dbg, const SparseBlockMatrix &m) {
    int nBlockColumns = m.size();

    for ( int i =0; i < nBlockColumns; i++ )
      dbg << *m.at(i);

    return dbg;
  }
}