SparseMatrix.hh

/*
This file is part of the BIAS library (Basic ImageAlgorithmS).

Copyright (C) 2003-2009    (see file CONTACT for details)
Multimediale Systeme der Informationsverarbeitung
Institut fuer Informatik
Christian-Albrechts-Universitaet Kiel


BIAS is free software; you can redistribute it and/or modify
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation; either version 2.1 of the License, or
(at your option) any later version.

BIAS is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public License
along with BIAS; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
*/

#ifndef __SparseMatrix_hh__
#define __SparseMatrix_hh__

#include "bias_config.h"


#include "Base/Math/Matrix.hh"
#include "Base/Math/Vector.hh"
#include <map>
#include <list>
#include <iostream>

// threshold for zero elements
#define SPARSE_MATRIX_EPSILON 1e-12

namespace BIAS {

  /** @class SparseMatrix
      @brief Implementation of sparse matrix operations.
      @author beder, rework by esquivel
      @date 01/2007, rework 11/2008, 
      @note This implementation of a sparse matrix uses a flexible structure
                        that may be extended by further elements, without the need to
                        reallocate the data manually. This structure is realized using
                        a list of pairs.
   */
class BIASMathBase_EXPORT SparseMatrix {

  public:

    /** @brief Default constructor for zero sparse matrix */
  SparseMatrix(): maxRow_(0), maxCol_(0), numRow_(0), numCol_(0) {}
   
    /** @brief Default constructor for zero sparse matrix */
  SparseMatrix(unsigned int numRow, unsigned int numCol): maxRow_(0), maxCol_(0), numRow_(numRow), numCol_(numCol) {}

    /** @brief Constructor from BIAS matrix */
    SparseMatrix(BIAS::Matrix<double> const &A);

    
    /** @brief Default destructor */
    ~SparseMatrix() {}

    SparseMatrix& operator=(const SparseMatrix& S);
    
    /** @brief Returns the maximum number of rows (max non-zero entry). */
    unsigned int GetRows() const { return maxRow_+1; }

    /** @brief Returns the maximum number of columns (max non-zero entry). */
    unsigned int GetCols() const { return maxCol_+1; }

    /** @brief Returns the total number of rows (matrix size). */
    unsigned int GetRowsNum() const { return numRow_; }

    /** @brief Returns the total number of columns (matrix size). */
    unsigned int GetColsNum() const { return numCol_; }

    /** @brief Returns value of an element.
        @note Not very efficient since it has to be tested first if
              the element has been set already!
        @author esquivel 11/2008 */
    double GetElement(unsigned int row, unsigned int col) const;

    /** @brief Sets a new element or changes the value of an existing element.
        @note Existing elements are removed by setting value to 0!
        @author esquivel 11/2008 */
    void SetElement(unsigned int row, unsigned int col, double val);

    /** @brief Sets a submatrix beginning at (row,col) from a dense matrix.
        @brief NO part of the specified submatrix must have been set before!
     */
    void SetSubMatrix(unsigned int row, unsigned int col,
                      BIAS::Matrix<double> const &A,
                      unsigned int startRowInA, unsigned int startColInA,
                      unsigned int row_count, unsigned int col_count);

    /** @brief Returns transposed matrix in result (as sparse matrix). */
    void Transpose(SparseMatrix &result);

    /** @brief Multiplies this elementwise with 'multiplier'. */
    void MultiplyIP(const double multiplier);

    /** @brief Multiplies with A and returns result M*A. */
    void Multiply(SparseMatrix &A, SparseMatrix &result);

    /** @brief Multiplies with transpose of A and returns result M*A'. */
    void MultiplyWithTransposeOf(SparseMatrix &A, SparseMatrix &result);

    /** @brief Multiplies with vector x and returns result A*x. */
    void Multiply(BIAS::Vector<double> &x, BIAS::Vector<double> &result);

    /** @brief Inverts matrix and returns inverse (as dense matrix) and
               solution of linear equation system A^(-1)*x = b in result.
        @attention Valid only for square matrices! */
    int InvertAndSolve(BIAS::Vector<double> &b, BIAS::Matrix<double> &inverse,
                       BIAS::Vector<double> &result);

    /** @brief Inverts matrix and returns inverse (as sparse matrix) and
               solution of linear equation system A^(-1)*x = b in result.
        @attention Valid only for square matrices! */
    int InvertAndSolve(BIAS::Vector<double> &b, SparseMatrix &inverse,
                       BIAS::Vector<double> &result);

    /** @brief Returns solution of linear equation system A*x = b. */
    int Solve(BIAS::Vector<double> &b, BIAS::Vector<double> &result);

    /// @brief computes pseudo inverse and returns result in inverse
    int PseudoInverse(BIAS::Matrix<double>& inverse);
    
    /// @brief computes pseudo inverse and returns result in inverse
    int PseudoInverse(BIAS::SparseMatrix& inverse);
    
    /// @brief clears old matrix and reinits with given size
    void Reinit(unsigned rows, unsigned cols);
    
    /** @brief Returns inverse matrix as dense matrix.
        @attention Valid only for square matrices! */
    int Invert(BIAS::Matrix<double> &inverse);

    /** @brief Returns inverse matrix as sparse matrix.
        @attention Valid only for square matrices! */
    int Invert(SparseMatrix &inverse);

    /* @brief Appends matrix H at the bottom (MATLAB notation: A = [A; H]). */
    void AppendMatrixBottom (SparseMatrix &H);

    /* @brief Appends matrix H at the right (MATLAB notation: A = [A, H]). */
    void AppendMatrixRight (SparseMatrix &H);

    /** @brief Returns this matrix as dense BIAS matrix in M. */
    void GetAsDense(BIAS::Matrix<double> &M) const;

    /** @brief Returns submatrix as dense BIAS matrix in M.
        @todo Not efficient because of row->column mapping... */
    void GetAsDense(unsigned int row, unsigned int col, unsigned int numRows,
                    unsigned int numCols, BIAS::Matrix<double> &M) const;

    /** @brief Prints the contents of this matrix. */
    void Print(std::string const &name, std::ostream &stream=std::cout) const;

    /** @brief Returns maximum element at diagonal (at least zero).
        @param row Returns diagonal index of maximum (at least 0).
        @author koeser 04/2007 */
    double GetMaxDiagonalElement(unsigned int &index) const;
    inline double GetMaxDiagonalElement() const {
      unsigned int index; return GetMaxDiagonalElement(index);
    };

    /** @brief Returns maximum element at given row (at least zero).
        @param col Returns column index of maximum (at least 0).
        @author esquivel 11/2008 */
    double GetMaxRowElement(unsigned int row, unsigned int &col) const;
    inline double GetMaxRowElement(unsigned int row) const {
      unsigned int col; return GetMaxRowElement(row, col);
    };

    /** @brief Returns maximum element at given column (at least zero).
        @param row Returns row index of maximum (at least 0).
        @todo Not efficient because of row->column mapping...
        @author esquivel 11/2008 */
    double GetMaxColumnElement(unsigned int col, unsigned int &row) const;
    inline double GetMaxColumnElement(unsigned int col) const {
      unsigned int row; return GetMaxColumnElement(col, row);
    };

    /** @brief Runs across diagonal and adds value to diagonal elements.
        @note If diagonal elements does not exist, it is created.
        @attention Implemented only for square matrices!
        @author koeser 04/2007 */
    void AddToDiagonal(const double &value);

    /** @brief Runs across diagonal and multiplies diagonal element with value.
        @author koeser 04/2007 */
    void MultiplyDiagonalBy(const double &value);

    /** @brief Checks if this sparse matrix is still consistent, i.e.
               order of rows and columns is strictly ascending, and there
               are no empty rows, no values below epsilon and no values
               outside of the matrix boundaries given by maxRow_, maxCol_!
        @return Returns true if matrix is consistent, false otherwise.
        @author esquivel 11/2008 */
    bool CheckConsistency();

    /** @brief Return rows, columns and values for all elements in matrix.
        @param rows Returns row indices of elements.
        @param cols Returns column indices of elements.
        @param values Returns values of elements.
    */
    void GetAllElements(std::vector<unsigned int> &rows,
                        std::vector<unsigned int> &cols,
                        std::vector<double> &values) const;
   
  protected: 

    /** @brief Insert an element.
        @note The element must NOT have been set before!
     */
    void InsertElement_(unsigned int row, unsigned int col, double value);

    /** @brief Performs Gauss-Jordan algorithm on sparse matrix. */
    int GaussJordan_();

    /** @brief Maps row indices to mappings from column indices to values. */
    std::map<unsigned int, std::list<std::pair<unsigned int, double> > > rows_;

    /** @brief Current max. row and column containing non-zero elements. */
    unsigned int maxRow_, maxCol_;

    /** @brief  Indicates the current number of rows and clumns.
                              values are negative, if the dimensions of the matrix are unknown.*/
    unsigned int numRow_, numCol_;

  };

  inline std::ostream& operator<<(std::ostream& stream,
                                const BIAS::SparseMatrix& matrix)
  {
    matrix.Print("SparseMatrix", stream);
    return stream;
  }

}

#endif