SparseBlockMatrix.cpp

#include "SparseBlockMatrix.h"

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

#include <iostream>
#include <iomanip>
#include <QDebug>
#include <QMapIterator>
#include <QListIterator>

#include "IString.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, matrix<double>*> it(src);
     while (it.hasNext()) {
         it.next();

         // copy matrix block from src
         matrix<double>* m = new matrix<double>(*(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
    matrix<double>* m = new matrix<double>(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, matrix<double>*> 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, matrix<double>*> 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, matrix<double>*> 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, matrix<double>*> 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::zeroBlocks() {
    QMapIterator<int, matrix<double>*> 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 << nBlocks;

    QMapIterator<int, matrix<double>*> 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() << nRows << 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) {
    int 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
      matrix<double>* 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, matrix<double>*> it(sbcm);
     while ( it.hasNext() ) {
       it.next();

       if( !it.value() )
         continue;

       // get matrix
       matrix<double>* 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, matrix<double>*> it(src);
     while (it.hasNext()) {
         it.next();

         // copy matrix block from src
         matrix<double>* m = new matrix<double>(*(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;

    matrix<double>* m = new matrix<double>(nRows,nCols);

    if ( !m )
      return false;

    m->clear();

    this->insert(nRowBlock,m);

    return true;
  }


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

    QMapIterator<int, matrix<double>*> 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, matrix<double>*> 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::zeroBlocks() {
    QMapIterator<int, matrix<double>*> 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, matrix<double>*> it(*this);
    while ( it.hasNext() ) {
      it.next();

      nrowBlock = it.key();
      matrix<double>* 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;
    }
  }

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

    QMapIterator<int, matrix<double>*> 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() << nRows << 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) {
    int 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
      matrix<double>* 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, matrix<double>*> it(sbrm);
     while ( it.hasNext() ) {
       it.next();

       if( !it.value() )
         continue;

       // get matrix
       matrix<double>* 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, matrix<double>*> 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;
  }


  matrix<double>* 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;
    }
  }


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

    stream << nBlockColumns;

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

    return stream;
  }


  QDataStream &operator>>(QDataStream &stream, SparseBlockMatrix &sparseBlockMatrix) {
    int 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;
  }
}