SpMat.h

//  
//  Declarations/template-bodies for sparse matrix class SpMat
//
//  SpMat.h
//
//  Implements bare-bones sparse matrix class. 
//  Main considerations has been efficiency when constructing
//  from Compressed Column format, when multiplying with vector,
//  transposing and multiplying with a vector and when concatenating.
//  Other operations which have not been prioritised such as 
//  for example inserting elements in a random order may be
//  a bit slow.  
//
//
// Jesper Andersson, FMRIB Image Analysis Group
//
// Copyright (C) 2007 University of Oxford 
//

#ifndef SpMat_h
#define SpMat_h

#include <vector>
#include <fstream>
#include <iomanip>
#include <boost/shared_ptr.hpp>
#include "newmat.h"
#include "cg.h"
#include "bicg.h"
#include "miscmaths.h"

namespace MISCMATHS {

class SpMatException: public std::exception
{
private:
  std::string m_msg;
public:
  SpMatException(const std::string& msg) throw(): m_msg(msg) {}

  virtual const char * what() const throw() {
    return string("SpMat::" + m_msg).c_str();
  }

  ~SpMatException() throw() {}
};

enum MatrixType {UNKNOWN, ASYM, SYM, SYM_POSDEF};

template<class T>
class Preconditioner;

template<class T>
class Accumulator;

//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
//
// Class SpMat:
// Interface includes:
// Multiplication with scalar:  A*=s, B=s*A, B=A*s, A and B SpMat
// Multiplication with vector:  b=A*x, A SpMat, b and x ColumnVector
// Transpose and mul with vector: b=A.trans_mult(x), A SpMat, b and x ColumnVector
// Multiplication with sparse matrix: C=A*B, A, B and C SpMat
// Addition with sparse matrix: A+=B, C=A+B, A, B and C SpMat
// Horisontal concatenation: A|=B, C=A|B, A, B and C SpMat
// Vertical concatenation: A&=B, C=A&B, A, B and C SpMat
//
// Multiplications and addition with NEWMAT matrices are
// accomplished through type-conversions. For example
// A = B*SpMat(C), A and B SpMat, C NEWMAT
// A = B.AsNewmat()*C, B SpMat, A and C NEWMAT
//
// Important implementation detail:
// _nz or .NZ() isn't strictly speaking the # of non-zero elements,
// but rather the number of elements that has an explicit 
// representation, where that representation may in principle
// be 0. This is in contrast to e.g. Matlab which chooses
// not to represent an element when its value is zero. I have
// chosen this variant because of my main use of the class where 
// it is very convenient if e.g. my Hessian and the Gibbs form
// of membrane energy has the same sparsity pattern.
// For most users this is of no consequence and they will 
// never explicitly represent a zero.
//
//@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

template<class T>
class SpMat
{
public:
  SpMat() : _m(0), _n(0), _nz(0), _ri(0), _val(0), _pw(false) {}
  SpMat(unsigned int m, unsigned int n) : _m(m), _n(n), _nz(0), _ri(n), _val(n), _pw(false) {}
  SpMat(unsigned int m, unsigned int n, const unsigned int *irp, const unsigned int *jcp, const double *sp);
  SpMat(const NEWMAT::GeneralMatrix& M);
  SpMat(const std::string& fname);
  ~SpMat() {}

  unsigned int Nrows() const {return(_m);}
  unsigned int Ncols() const {return(_n);}
  unsigned int NZ() const {return(_nz);}

  NEWMAT::ReturnMatrix AsNEWMAT() const;
  void Save(const std::string&   fname,
            unsigned int         precision) const;
  void Save(const std::string&   fname) const {Save(fname,8);}
  void Print(const std::string&  fname,
             unsigned int        precision) const;
  void Print(const std::string&  fname) const {Print(fname,8);}
  void Print(unsigned int        precision) const {Print(std::string(""),precision);}
  void Print() const {Print(8);}
  void WarningsOn() {_pw=true;}
  void WarningsOff() {_pw=false;}
  bool IsSorted() const { return(is_sorted()); }  // Returns true if all _ri arrays are sorted. For debugging.

  T Peek(unsigned int r, unsigned int c) const;
  T operator()(unsigned int r, unsigned int c) const {return(Peek(r,c));}    // Read-only

  void Set(unsigned int r, unsigned int c, const T& v) {here(r,c) = v;}             // Set a single value
  void SetColumn(unsigned int c, const NEWMAT::ColumnVector& col, double eps=0.0);  // Set a whole column (obliterating what was there before)
  SpMat<T>& MultiplyColumns(const NEWMAT::Matrix& vals);                             // Multiply each column of a matrix by a value
  void AddTo(unsigned int r, unsigned int c, const T& v) {here(r,c) += v;}          // Add value to a single (possibly existing) value

  SpMat<T>& operator+=(const SpMat& M) 
  {
    if (same_sparsity(M)) return(add_same_sparsity_mat_to_me(M,1));
    else return(add_diff_sparsity_mat_to_me(M,1));
  }
  SpMat<T>& operator-=(const SpMat& M) 
  {
    if (same_sparsity(M)) return(add_same_sparsity_mat_to_me(M,-1)); 
    else return(add_diff_sparsity_mat_to_me(M,-1));
  }
  const NEWMAT::ReturnMatrix operator*(const NEWMAT::ColumnVector& x) const;       // Multiplication with column vector
  const NEWMAT::ReturnMatrix trans_mult(const NEWMAT::ColumnVector& x) const;      // Multiplication of transpose with column vector
  const NEWMAT::ReturnMatrix TransMult(const NEWMAT::ColumnVector& x) const {
    return(trans_mult(x));                                                         // Duplication for compatibility with IML++
  }
  const SpMat<T> TransMult(const SpMat<T>& B) const { return(this->t()*B); }       // Multiplication of transpose(*this) with sparse matrix B
  SpMat<T>& operator*=(double s);                                                  // Multiplication of self with scalar
  SpMat<T> operator-(const SpMat<T>& M) const {return(SpMat<T>(M) *= -1.0);}       // Unary minus
  SpMat<T>& operator|=(const SpMat<T>& rh);                                        // Hor concat to right
  SpMat<T>& operator&=(const SpMat<T>& bh);                                        // Vert concat below

  const SpMat<T> TransMultSelf() const {return(this->t() * *this);}                 // Returns transpose(*this)*(*this)

  const SpMat<T> t() const;                                                        // Returns transpose(*this). 
  
  friend class Accumulator<T>; 

  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
  //
  // Class ColumnIterator
  //
  // This implements a const forward iterator for a given column
  //
  //@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

  // template<class TT>
  class ColumnIterator
  {
  public:
    ColumnIterator(const SpMat<T>& mat, unsigned int col, bool end=false) {
      if (end) { _ri_it = mat._ri[col-1].end(); _val_it = mat._val[col-1].end(); }
      else { _ri_it = mat._ri[col-1].begin(); _val_it = mat._val[col-1].begin(); }
    }
    ~ColumnIterator() {}
    T operator*() const { return(*_val_it); }
    unsigned int Row() const { return((*_ri_it)+1); }
    bool operator==(const ColumnIterator& rhs) const { return(_val_it == rhs._val_it); }
    bool operator!=(const ColumnIterator& rhs) const { return(!(*this == rhs)); }
    // Prefix increment. Use whenever possible
    ColumnIterator& operator++() { ++_ri_it; ++_val_it; return(*this); }
    // Postfix increment. Avoid.
    ColumnIterator& operator++(int dummy) { ColumnIterator clone(*this); ++_ri_it; ++_val_it; return(clone); }
  private:
    typename std::vector<T>::const_iterator    _val_it;
    std::vector<unsigned int>::const_iterator  _ri_it;
  };

  ColumnIterator begin(unsigned int col) const { if (col>_n) throw SpMatException("ColumnIterator: col out of range"); return(ColumnIterator(*this,col)); }
  ColumnIterator end(unsigned int col) const { if (col>_n) throw SpMatException("ColumnIterator: col out of range"); return(ColumnIterator(*this,col,true)); }

  template<class TT>
  friend const SpMat<TT> operator*(const SpMat<TT>& lh, const SpMat<TT>& rh);      // Multiplication of two sparse matrices

  NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&           b,          // Solve for x in b=(*this)*x
                                 MatrixType                            type = UNKNOWN,
                                 double                                tol = 1e-4, 
                                 unsigned int                          miter = 200,
				 boost::shared_ptr<Preconditioner<T> > C = boost::shared_ptr<Preconditioner<T> >()) const;

  NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
			         MatrixType                             type,
			         double                                 tol,
			         unsigned int                           miter,
				 const NEWMAT::ColumnVector&            x_init) const;

  NEWMAT::ReturnMatrix SolveForx(const NEWMAT::ColumnVector&            b,
			         MatrixType                             type,
			         double                                 tol,