SparseBlockMatrix.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#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() {
    m_startColumn = 0;
  }
 
 
  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);
    }
 
    m_startColumn = src.startColumn();
  }
 
 
  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;
  }
 
 
  void SparseBlockColumnMatrix::setStartColumn(int nStartColumn) {
    m_startColumn = nStartColumn;
  }
 
 
  int SparseBlockColumnMatrix::startColumn() const {
    return m_startColumn;
  }
 
 
  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" << Qt::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() << Qt::endl;
      }
 
      dbg.space() << Qt::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" << Qt::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() << Qt::endl;
      }
      dbg.space() << Qt::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;
  }
}