#ifndef _utl_SparseMatrix_h_
#define _utl_SparseMatrix_h_

#include <map>
#include <numeric>

#include <boost/lambda/lambda.hpp>


namespace utl {

  /**
    \class SparseMatrix SparseMatrix.h "utl/SparseMatrix.h"
    \brief Sparse container class for matrix data

    Importatant notes:
    <ul>
    <li>
      the operator[i][j] access is not creating entries. It should be used for
      reading of the whole matrix.
    <li>
      the operator(i, j) access creates the corresponding entry on the fly and
      initializes it to GetInitializer() <b>before</b> the reference to the
      element is returned
    </ul>

    Please refer to the testSparseMatrix.cc and test SparseMatrixVectorOps.cc
    files for detailed usage cases.

    \author Darko Veberic
    \version $Id: SparseMatrix.h 23590 2013-05-06 17:59:58Z paul $
    \date 5 Jun 2008
    \ingroup math
  */

  template<typename T = double>
  class SparseMatrix {
  public:
    typedef int IndexType;

    /// row proxy
    class Row {
    public:
      Row(const SparseMatrix& m, const IndexType row)
        : fSparseMatrix(m), fRow(row) { }

      const T&
      operator[](const IndexType column)
        const
      { return fSparseMatrix(fRow, column); }

    private:
      const SparseMatrix& fSparseMatrix;
      const IndexType fRow;
    };

    /// index class
    class Index {
    public:
      Index(const IndexType row, const IndexType column)
        : fRow(row), fColumn(column) { }

      IndexType GetRow() const { return fRow; }

      IndexType GetColumn() const { return fColumn; }

      bool operator<(const Index& i) const
      { return fRow < i.fRow || (fRow == i.fRow && fColumn < i.fColumn); }

      bool operator==(const Index& i) const
      { return fRow == i.fRow && fColumn == i.fColumn; }

    private:
      const IndexType fRow;
      const IndexType fColumn;
    };

  private:
    typedef typename SparseMatrix<T>::Index IndexKeyType;
    typedef std::map<IndexKeyType, T> MapType;
    typedef typename MapType::value_type MapValueType;
    typedef typename MapType::iterator MapIterator;
    typedef typename MapType::const_iterator MapConstIterator;

  public:
    /// iterator class
    template<class IteratorType, typename ValueType>
    class IteratorTransformer {
    public:
      IteratorTransformer(const IteratorType& it)
        : fIterator(it) { }

      IndexType GetRow() const { return fIterator->first.GetRow(); }

      IndexType GetColumn() const { return fIterator->first.GetColumn(); }

      ValueType& operator()() const { return fIterator->second; }

      ValueType& operator*() const { return (*this)(); }

      IteratorTransformer& operator++() { ++fIterator; return *this; }

      bool operator==(const IteratorTransformer& it) const { return fIterator == it.fIterator; }

      bool operator!=(const IteratorTransformer& it) const { return fIterator != it.fIterator; }

    private:
      IteratorType fIterator;
    };

    SparseMatrix(const T& init = T())
      : fInitializer(init) { }

    template<typename U>
    explicit
    SparseMatrix(const SparseMatrix<U>& m)
      : fInitializer(m.GetInitializer())
    {
      for (typename SparseMatrix<U>::ConstIterator it = m.SparseBegin();
           it != m.SparseEnd(); ++it)
        (*this)(it.GetRow(), it.GetColumn()) = it();
    }

    unsigned int GetNonEmptySize() const { return fMatrix.size(); }

    const T& GetInitializer() const { return fInitializer; }

    void Clear() { fMatrix.clear(); }

    /// noncreational [] access trough proxy
    typename SparseMatrix::Row
    operator[](const IndexType row)
      const
    { return typename SparseMatrix::Row(*this, row); }

    /// creational () access
    T&
    operator()(const IndexType row, const IndexType column)
    {
      return
        fMatrix.insert(MapValueType(IndexKeyType(row, column), fInitializer)).first->second;
    }

    /// noncreational () access for const
    const T&
    operator()(const IndexType row, const IndexType column)
      const
    {
      const MapConstIterator it = fMatrix.find(IndexKeyType(row, column));
      return (it == fMatrix.end()) ? fInitializer : it->second;
    }

    bool IsEmpty() const { return fMatrix.empty(); }

    bool
    IsEmpty(const IndexType row, const IndexType column)
      const
    {
      const MapConstIterator it = fMatrix.find(IndexKeyType(row, column));
      return it == fMatrix.end();
    }

    unsigned int
    Reclaim()
    {
      unsigned int erased = 0;
      for (MapIterator it = fMatrix.begin();
           it != fMatrix.end(); )
        if (it->second != fInitializer)
          ++it;
        else {
          Erase(it++);
          ++erased;
        }
      return erased;
    }

    /// be sure to call Reclaim() first
    bool operator==(const SparseMatrix& m) const
    { return fMatrix == m.fMatrix && fInitializer == m.fInitializer; }

    /// be sure to call Reclaim() first
    bool operator!=(const SparseMatrix& m) const { return !operator==(m); }

    template<typename U>
    bool
    operator==(const SparseMatrix<U>& m)
      const
    {
      if (fInitializer != m.GetInitializer())
        return false;

      ConstIterator it1 = SparseBegin();
      typename SparseMatrix<U>::ConstIterator it2 = m.SparseBegin();
      while (it1 != SparseEnd() && it2 != m.SparseEnd()) {
        if (it1.GetRow() != it2.GetRow() ||
            it1.GetColumn() != it2.GetColumn() ||
            it1() != it2())
          return false;
        ++it1;
        ++it2;
      }

      if (it1 != SparseEnd() || it2 != m.SparseEnd())
        return false;

      return true;
    }

    template<typename U>
    bool operator!=(const SparseMatrix<U>& m) const { return !operator==(m); }

    typedef IteratorTransformer<MapIterator, T> Iterator;
    typedef IteratorTransformer<MapConstIterator, const T> ConstIterator;

    Iterator SparseBegin() { return Iterator(fMatrix.begin()); }

    ConstIterator SparseBegin() const { return ConstIterator(fMatrix.begin()); }

    Iterator SparseEnd() { return Iterator(fMatrix.end()); }

    ConstIterator SparseEnd() const { return ConstIterator(fMatrix.end()); }

    Iterator Find(const IndexType row, const IndexType column)
    { return Iterator(fMatrix.find(Index(row, column))); }

    ConstIterator Find(const IndexType row, const IndexType column) const
    { return ConstIterator(fMatrix.find(Index(row, column))); }

    /// multiplication with scalar
    template<typename U>
    SparseMatrix&
    operator*=(const U& value)
    {
      for (MapIterator it = fMatrix.begin(); it != fMatrix.end(); ++it)
        it->second *= value;
      return *this;
    }

    /// matrix addition
    template<typename U>
    SparseMatrix&
    operator+=(const SparseMatrix<U>& m)
    {
      for (typename SparseMatrix<U>::ConstIterator it = m.SparseBegin();
           it != m.SparseEnd(); ++it)
        (*this)(it.GetRow(), it.GetColumn()) += it();
      return *this;
    }

    /// matrix subtraction
    template<typename U>
    SparseMatrix&
    operator-=(const SparseMatrix<U>& m)
    {
      for (typename SparseMatrix<U>::ConstIterator it = m.SparseBegin();
           it != m.SparseEnd(); ++it)
        (*this)(it.GetRow(), it.GetColumn()) -= it();
      return *this;
    }

    /// in-situ transpose
    void
    Transpose()
    {
      MapType t;
      for (MapConstIterator it = fMatrix.begin(); it != fMatrix.end(); ++it) {
        const IndexKeyType& i = it->first;
        t.insert(t.begin(), std::make_pair(IndexKeyType(i.GetColumn(), i.GetRow()), it->second));
      }
      fMatrix.swap(t);
    }

  private:
    void Erase(const MapIterator it) { fMatrix.erase(it); }

    MapType fMatrix;
    const T fInitializer;

    template<typename U> friend class SparseMatrix;
  };


  /** general matrix multiplication
    In case matrix elements are not simple scalar values, the corresponding
    operator* will be used. For example, imagine matrix elements are
    vectors and the element * element multiplication is equal to scalar
    product. The required dot product functor is then used.
  */
  template<typename T, typename U>
  typename boost::lambda::return_type_2<
    boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
    SparseMatrix<T>,
    SparseMatrix<U>
  >::type
  operator*(const SparseMatrix<T>& a, const SparseMatrix<U>& b)
  {
    typedef typename boost::lambda::return_type_2<
      boost::lambda::arithmetic_action<boost::lambda::multiply_action>,
      SparseMatrix<T>,
      SparseMatrix<U>
    >::type ResultMatrix;

    ResultMatrix c;

    for (typename SparseMatrix<U>::ConstIterator bIt = b.SparseBegin();
         bIt != b.SparseEnd(); ++bIt) {
      const typename SparseMatrix<U>::IndexType k = bIt.GetRow();
      for (typename SparseMatrix<T>::ConstIterator aIt = a.SparseBegin();
           aIt != a.SparseEnd(); ++aIt) {
        const typename SparseMatrix<T>::IndexType kk = aIt.GetColumn();
        if (kk == k) {
          const typename SparseMatrix<T>::IndexType i = aIt.GetRow();
          const typename SparseMatrix<U>::IndexType j = bIt.GetColumn();
          c(i, j) += aIt() * bIt();
        }
      }
    }

    return c;
  }


  // output stuff

  template<class Stream, typename T>
  Stream&
  operator<<(Stream& s, const SparseMatrix<T>& m)
  {
    s << '(';
    typename SparseMatrix<T>::ConstIterator it = m.SparseBegin();
    if (it != m.SparseEnd()) {
      s << '{' << it.GetRow() << ',' << it.GetColumn() << ", " << it() << '}';
      ++it;
      while (it != m.SparseEnd())
        s << " {" << it.GetRow() << ',' << it.GetColumn() << ", " << it() << '}';
    }
    return s << ')';
  }


  // output functions assume the lowest index starts with 0, which is not
  // necessarily true

  template<class Stream, typename T, typename SparseType>
  Stream&
  OutputSparse(Stream& s,
               const SparseMatrix<T>& m,
               const typename SparseMatrix<T>::IndexType rows,
               const typename SparseMatrix<T>::IndexType columns,
               const char separator, const SparseType& sparse)
  {
    typedef typename SparseMatrix<T>::IndexType Index;
    for (Index i = 0; i < rows; ++i) {
      if (m.IsEmpty(i, 0))
        s << sparse;
      else
        s << m[i][0];
      for (Index j = 1; j < columns; ++j) {
        s << separator;
        if (m.IsEmpty(i, j))
          s << sparse;
        else
          s << m[i][j];
      }
      s << '\n';
    }
    return s;
  }


  template<class Stream, typename T>
  inline
  Stream&
  OutputSparse(Stream& s,
               const SparseMatrix<T>& m,
               const typename SparseMatrix<T>::IndexType rows,
               const typename SparseMatrix<T>::IndexType columns)
  {
    return OutputSparse(s, m, rows, columns, ' ', '.');
  }


  template<class Stream, typename T>
  inline
  Stream&
  Output(Stream& s,
         const SparseMatrix<T>& m,
         const typename SparseMatrix<T>::IndexType rows,
         const typename SparseMatrix<T>::IndexType columns,
         const char separator = ' ')
  {
    return OutputSparse(s, m, rows, columns, separator, m.GetInitializer());
  }


}


namespace boost {

  namespace lambda {

    /// specialization for custom class SMatrix<>
    template<typename T, typename U>
    class plain_return_type_2<
      arithmetic_action<multiply_action>,
      utl::SparseMatrix<T>,
      utl::SparseMatrix<U>
    > {
    private:
      // find type of T * U
      typedef typename
        return_type_2<
          arithmetic_action<multiply_action>,
          T,
          U
        >::type res_type;

    public:
      typedef utl::SparseMatrix<res_type> type;
    };

  }

}


#include <utl/SparseMatrixVectorOp.h>


#endif