/* Copyright (c) 1997-2005
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Berlin, Germany)
   http://www.math.tu-berlin.de/polymake,  mailto:polymake@math.tu-berlin.de

   This program is free software; you can redistribute it and/or modify it
   under the terms of the GNU General Public License as published by the
   Free Software Foundation; either version 2, or (at your option) any
   later version: http://www.gnu.org/licenses/gpl.txt.

   This program 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 General Public License for more details.
*/

#ifndef _POLYMAKE_SPARSE_MATRIX_H
#define _POLYMAKE_SPARSE_MATRIX_H "$Project: polymake $$Id: SparseMatrix.h,v 1.30 2005/11/11 10:27:58 gawrilow Exp $"

#ifndef _POLYMAKE_SPARSE2D_H
#  include <sparse2d.h>
#endif
#ifndef _POLYMAKE_GENERIC_MATRIX_H
#  include <GenericMatrix.h>
#endif
#ifndef _POLYMAKE_SPARSE_H
#  include <sparse.h>
#endif
#ifndef _POLYMAKE_SHARED_OBJECT_H
#  include <shared_object.h>
#endif

namespace pm {

template <typename _Matrix, bool _transpose=false>
class sparse_matrix_proxy_base {
public:
   typedef typename _Matrix::value_type value_type;
protected:
   _Matrix* M;
   int i, j;

   typedef typename if_else<_transpose, typename _Matrix::col_type, typename _Matrix::row_type>::result row_vector;
   typedef typename row_vector::iterator iterator;
   typedef typename row_vector::const_iterator const_iterator;
   typedef bool2type<_transpose> Transp;

   row_vector& row(False*)
   {
      return M->row(i);
   }
   row_vector& row(True*)
   {
      return M->col(i);
   }

   const row_vector& row(False*) const
   {
      return const_cast<const _Matrix*>(M)->row(i);
   }
   const row_vector& row(True*) const
   {
      return const_cast<const _Matrix*>(M)->col(i);
   }

   const value_type& get() const
   {
      const_iterator e=row((Transp*)0).find(j);
      if (!e.at_end()) return *e;
      operations::clear<value_type> zero;
      return zero();
   }

   iterator find() { return row((Transp*)0).insert(j); }

   void insert(const value_type& x) { row((Transp*)0).insert(j,x); }

   void erase() { row((Transp*)0).erase(j); }

   void erase(const iterator& pos) { row((Transp*)0).erase(pos); }

public:
   sparse_matrix_proxy_base(_Matrix* M_arg, int i_arg, int j_arg)
      : M(M_arg), i(_transpose ? j_arg : i_arg), j(_transpose ? i_arg : j_arg) { }
};

template <typename E, sparse2d::restriction_kind restriction=sparse2d::only_rows>
class RestrictedSparseMatrix
   : public matrix_methods<RestrictedSparseMatrix<E,restriction>, E> {
protected:
   typedef sparse2d::Table<E, false, restriction> table_type;
   table_type data;

   template <typename Iterator>
   void _copy(Iterator src, True*) { std::copy(src, entire(pm::rows(*this))); }

   template <typename Iterator>
   void _copy(Iterator src, False*) { std::copy(src, entire(pm::cols(*this))); }

public:
   typedef E value_type;
   typedef E& reference;
   typedef const E& const_reference;

   explicit RestrictedSparseMatrix(int n=0) : data(n) { }

   template <typename Iterator>
   RestrictedSparseMatrix(int n, Iterator src)
      : data(n)
   {
      _copy(src, (bool2type<restriction==sparse2d::only_rows>*)0);
   }

   void swap(RestrictedSparseMatrix& M)
   {
      data.swap(M.data);
   }

   void clear() { data.clear(); }

protected:
   typedef sparse_matrix_proxy_base<RestrictedSparseMatrix, restriction==sparse2d::only_cols> proxy_base;
   typedef sparse_elem_proxy<proxy_base> proxy;

   const E& _random(int i, int j, False*) const
   {
      return this->row(i)[j];
   }
   const E& _random(int i, int j, True*) const
   {
      return this->col(j)[i];
   }
public:
   proxy operator() (int i, int j)
   {
      return proxy_base(this,i,j);
   }

   const E& operator() (int i, int j) const
   {
      return _random(i, j, (bool2type<restriction==sparse2d::only_cols>*)0);
   }

private:
   template <typename RowCol, typename _Vector>
   void _append(RowCol& row_col, const _Vector& vec, int i)
   {
      for (typename Entire<_Vector>::const_iterator v=entire(vec); !v.at_end(); ++v)
	 row_col[v.index()].push_back(i,*v);
   }

   template <typename Iterator>
   void _append_rows(int n, Iterator src, True*)
   {
      int oldrows=data.rows();
      data.resize_rows(oldrows+n);
      std::copy(src, pm::rows(*this).begin()+oldrows);
   }

   template <typename Iterator>
   void _append_rows(int, Iterator src, False*)
   {
      for (int r=data.rows(); !src.at_end(); ++src, ++r)
	 _append(pm::cols(*this), *src, r);
   }

   template <typename Iterator>
   void _append_cols(int n, Iterator src, True*)
   {
      int oldcols=data.cols();
      data.resize_cols(oldcols+n);
      std::copy(src, pm::cols(*this).begin()+oldcols);
   }

   template <typename Iterator>
   void _append_cols(int, Iterator src, False*)
   {
      for (int c=data.cols(); !src.at_end(); ++src, ++c)
	 _append(pm::rows(*this), *src, c);
   }

public:
   template <typename _Matrix>
   RestrictedSparseMatrix& operator/= (const GenericMatrix<_Matrix>& m)
   {
      _append_rows(m.rows(), entire(pm::rows(m)), (bool2type<restriction==sparse2d::only_rows>*)0);
      return *this;
   }

   template <typename _Vector>
   RestrictedSparseMatrix& operator/= (const GenericVector<_Vector>& v)
   {
      _append_rows(1, entire(item2container(v.top())), (bool2type<restriction==sparse2d::only_rows>*)0);
      return *this;
   }

   template <typename _Matrix>
   RestrictedSparseMatrix& operator|= (const GenericMatrix<_Matrix>& m)
   {
      _append_cols(m.cols(), entire(pm::cols(m)), (bool2type<restriction==sparse2d::only_cols>*)0);
      return *this;
   }

   template <typename _Vector>
   RestrictedSparseMatrix& operator|= (const GenericVector<_Vector>& v)
   {
      _append_cols(1, entire(item2container(v.top())), (bool2type<restriction==sparse2d::only_cols>*)0);
      return *this;
   }

#ifdef POLY_DEBUG
   void check() const { data.check(); }
#endif
   friend class Rows<RestrictedSparseMatrix>;
   friend class Cols<RestrictedSparseMatrix>;
   friend class SparseMatrix<E>;
};

template <typename E, sparse2d::restriction_kind restriction>
struct spec_object_traits< RestrictedSparseMatrix<E, restriction> >
   : spec_object_traits<is_container> {
   static const int dimension=2;
};

template <typename E>
class SparseMatrix : public GenericMatrix< SparseMatrix<E>, E> {
protected:
   typedef sparse2d::Table<E,false> table_type;
   shared_object<table_type> data;

   template <typename Iterator>
   void init(Iterator src, True*)
   {
      const int n=this->cols();
      for (typename Entire< Rows<SparseMatrix> >::iterator r=entire(pm::rows(*this)); !r.at_end(); ++r) {
	 for (int i=0; i<n; ++i, ++src)
	    if (*src) r->push_back(i,*src);
      }
   }

   template <typename Iterator>
   void init(Iterator src, False*)
   {
      copy(src, entire(pm::rows(*this)));
   }

public:
   typedef E value_type;
   typedef E& reference;
   typedef const E& const_reference;

   SparseMatrix() { }

   SparseMatrix(int m, int n)
      : data( make_constructor(m,n,(table_type*)0) ) { }

   template <typename Iterator>
   SparseMatrix(int m, int n, Iterator src)
      : data( make_constructor(m,n,(table_type*)0) )
   {
      init(src, (bool2type<(object_traits<typename iterator_traits<Iterator>::value_type>::total_dimension ==
			object_traits<E>::total_dimension)>*)0);
   }

   SparseMatrix(const GenericMatrix<SparseMatrix>& M)
      : data(M.top().data) { }

   template <typename Matrix2>
   SparseMatrix(const GenericMatrix<Matrix2,E>& M)
      : data( make_constructor(M.rows(), M.cols(), (table_type*)0) )
   {
      copy(pm::rows(M).begin(), entire(pm::rows(*this)));
   }

   template <typename Matrix2, typename E2>
   explicit SparseMatrix(const GenericMatrix<Matrix2,E2>& M)
      : data( make_constructor(M.rows(), M.cols(), (table_type*)0) )
   {
      typename Rows<Matrix2>::const_iterator src=pm::rows(M).begin();
      for (typename Entire< Rows<SparseMatrix> >::iterator r=entire(pm::rows(*this)); !r.at_end();  ++r, ++src)
	 r->_assign(attach_converter<E>(*src), (pure_sparse*)0);
   }

   template <sparse2d::restriction_kind restriction>
   explicit SparseMatrix(RestrictedSparseMatrix<E,restriction>& M)
      : data( constructor<table_type(typename RestrictedSparseMatrix<E,restriction>::table_type&)>(M.data) ) { }

   SparseMatrix& operator= (const SparseMatrix& M)
   {
      assign(M);
      return *this;
   }

   using SparseMatrix::generic_type::operator=;

   template <sparse2d::restriction_kind restriction>
   SparseMatrix& operator= (RestrictedSparseMatrix<E,restriction>& M)
   {
      data=make_constructor(M.data, (table_type*)0);
      return *this;
   }

   void swap(SparseMatrix& M) { data.swap(M.data); }

   friend void relocate(SparseMatrix* from, SparseMatrix* to)
   {
      pm::relocate(&from->data, &to->data);
   }

   void resize(int r, int c) { data->resize(r,c); }

   void clear() { data.clear(); }

   void clear(int r, int c)
   {
      data=make_constructor(r,c,(table_type*)0);
   }

   void remove0s()
   {
      for (typename Entire< Rows<SparseMatrix> >::iterator r=entire(pm::rows(*this)); !r.at_end(); ++r)
	 r->remove0s();
   }
protected:
   typedef sparse_matrix_proxy_base<SparseMatrix> proxy_base;
   typedef sparse_elem_proxy<proxy_base> proxy;
public:
   proxy operator() (int i, int j) { return proxy_base(this,i,j); }

   const E& operator() (int i, int j) const { return this->row(i)[j]; }

   template <typename Matrix2>
   SparseMatrix& operator/= (const GenericMatrix<Matrix2,E>& m)
   {
      DIMENSION_CHECK((this->cols() != m.cols()), "append_rows - column dimensions mismatch");
      if (is_shared()) {
	 *this = const_cast<const SparseMatrix&>(*this) / m;
      } else {
	 int oldrows=this->rows();
	 data.get()->resize_rows(oldrows+m.rows());
	 copy(entire(pm::rows(m)), pm::rows(*this)+oldrows);
      }
      return *this;
   }

   template <typename Vector2>
   SparseMatrix& operator/= (const GenericVector<Vector2,E>& v)
   {
      DIMENSION_CHECK((this->cols() != v.dim()), "append_row - column dimensions mismatch");
      if (is_shared()) {
	 *this = const_cast<const SparseMatrix&>(*this) / v;
      } else {
	 int oldrows=this->rows();
	 data.get()->resize_rows(oldrows+1);
	 this->row(oldrows)=v;
      }
      return *this;
   }

   template <typename Matrix2>
   SparseMatrix& operator|= (const GenericMatrix<Matrix2,E>& m)
   {
      DIMENSION_CHECK((this->rows() != m.rows()), "append_cols - row dimensions mismatch");
      if (is_shared()) {
	 *this = const_cast<const SparseMatrix&>(*this) | m;
      } else {
	 int oldcols=this->cols();
	 data.get()->resize_cols(oldcols+m.cols());
	 copy(entire(pm::cols(m)), pm::cols(*this)+oldcols);
      }
      return *this;
   }

   template <typename Vector2>
   SparseMatrix& operator|= (const GenericVector<Vector2,E>& v)
   {
      DIMENSION_CHECK((this->rows() != v.dim()), "append_col - row dimensions mismatch");
      if (is_shared()) {
	 *this = const_cast<const SparseMatrix&>(*this) | v;
      } else {
	 int oldcols=this->cols();
	 data.get()->resize_cols(oldcols+1);
	 this->col(oldcols)=v;
      }
      return *this;
   }

   template <typename row_number_consumer, typename col_number_consumer>
   void squeeze(row_number_consumer rnc, col_number_consumer cnc) { data->squeeze(rnc,cnc); }

   template <typename row_number_consumer>
   void squeeze(row_number_consumer rnc) { data->squeeze(rnc); }

   void squeeze() { data->squeeze(); }

   template <typename Iterator>
   void permute_rows(Iterator perm) { data->permute_rows(perm, (False*)0); }

   template <typename Iterator>
   void permute_cols(Iterator perm) { data->permute_cols(perm, (False*)0); }

   template <typename Iterator>
   void permute_inv_rows(Iterator inv_perm) { data->permute_rows(inv_perm, (True*)0); }

   template <typename Iterator>
   void permute_inv_cols(Iterator inv_perm) { data->permute_cols(inv_perm, (True*)0); }

#ifdef POLY_DEBUG
   void check() const { data->check(); }
#endif
protected:
   friend class Rows<SparseMatrix>;
   friend class Cols<SparseMatrix>;
   friend class GenericMatrix<SparseMatrix>;
   friend class alias_to_shared<SparseMatrix>;

   bool is_shared() const { return data.is_shared(); }

   const shared_object<table_type>& data_member() const { return data; }

   void assign(const GenericMatrix<SparseMatrix>& M) { data=M.top().data; }

   template <typename Matrix2>
   void assign(const GenericMatrix<Matrix2>& M)
   {
      if (is_shared() || this->rows() != M.rows() || this->cols() != M.cols()) {
	 assign(SparseMatrix(M));
      } else {
	 GenericMatrix<SparseMatrix>::assign(M);
      }
   }

   template <typename E2>
   void _fill(const E2& x, False*)
   {
      if (is_shared())
	 data=make_constructor(this->rows(), this->cols(), (table_type*)0);
      if (x)
	 GenericMatrix<SparseMatrix>::_fill(x, (False*)0);
   }
};

template <typename E>
class alias_ptr< SparseMatrix<E> >
   : public alias_to_shared< SparseMatrix<E> > {
public:
   alias_ptr(const SparseMatrix<E>* ptr_arg)
      : alias_to_shared< SparseMatrix<E> >(ptr_arg) { }
};

template <typename E>
struct check_container_feature< SparseMatrix<E>, pure_sparse > : True { };

template <typename Tree>
class sparse_matrix_line
   : public sparse2d::line<Tree>,
     public GenericVector< sparse_matrix_line<Tree>, typename Tree::data_type> {
   typedef sparse2d::line<Tree> _super;
protected:
   ~sparse_matrix_line();
public:
   typedef typename Tree::data_type element_type;
   friend class SparseMatrix<element_type>;

   sparse_matrix_line& operator= (sparse_matrix_line& l)
   {
      return GenericVector<sparse_matrix_line>::operator=(l);
   }

   using GenericVector<sparse_matrix_line>::operator=;
protected:
   typedef sparse_proxy_base<sparse_matrix_line> proxy_base;
   typedef sparse_elem_proxy<proxy_base> proxy;
public:
   proxy operator[] (int i) { return proxy_base(this,i); }

   const element_type& operator[] (int i) const
   {
      typename _super::const_iterator e=this->find(i);
      if (!e.at_end()) return *e;
      operations::clear<element_type> zero;
      return zero();
   }

   int dim() const { return this->max_size(); }
};

template <typename tree_type>
struct check_container_feature<sparse_matrix_line<tree_type>, pure_sparse> : True { };

template <typename tree_type>
struct spec_object_traits< sparse_matrix_line<tree_type> >
   : spec_object_traits<is_container> {
   static const int is_resizeable=0;
};

template <typename E, sparse2d::restriction_kind restriction>
class Rows< RestrictedSparseMatrix<E, restriction> >
   : public redirected_container< Rows< RestrictedSparseMatrix<E, restriction> >,
	       cons< Container< sparse2d::Rows<E, false, restriction, sparse_matrix_line> >,
		     Hidden< RestrictedSparseMatrix<E, restriction> > >
            > {
public:
   typedef typename if_else<restriction==sparse2d::only_rows, random_access_iterator_tag, output_iterator_tag>::result
      category;
   typename Rows::container& get_container()
   {
      return reinterpret_cast<typename Rows::container&>(this->hidden().data);
   }
   const typename Rows::container& get_container() const
   {
      return reinterpret_cast<const typename Rows::container&>(this->hidden().data);
   }
   void resize(int n)
   {
      this->hidden().data.resize_rows(n);
   }
};

template <typename E, sparse2d::restriction_kind restriction>
class Cols< RestrictedSparseMatrix<E, restriction> >
   : public redirected_container< Cols< RestrictedSparseMatrix<E, restriction> >,
	       cons< Container< sparse2d::Cols<E, false, restriction, sparse_matrix_line> >,
		     Hidden< RestrictedSparseMatrix<E, restriction> > >
            > {
public:
   typedef typename if_else<restriction==sparse2d::only_cols, random_access_iterator_tag, output_iterator_tag>::result
      category;
   typename Cols::container& get_container()
   {
      return reinterpret_cast<typename Cols::container&>(this->hidden().data);
   }
   const typename Cols::container& get_container() const
   {
      return reinterpret_cast<const typename Cols::container&>(this->hidden().data);
   }
   void resize(int n)
   {
      this->hidden().data.resize_cols(n);
   }
};

template <class E>
class Rows< SparseMatrix<E> >
   : public redirected_container< Rows< SparseMatrix<E> >,
	       cons< Container< sparse2d::Rows<E, false, sparse2d::full, sparse_matrix_line> >,
		     Hidden< SparseMatrix<E> > >
            > {
public:
   typename Rows::container& get_container()
   {
      return reinterpret_cast<typename Rows::container&>(*this->hidden().data);
   }
   const typename Rows::container& get_container() const
   {
      return reinterpret_cast<const typename Rows::container&>(*this->hidden().data);
   }
   void resize(int n)
   {
      this->hidden().data->resize_rows(n);
   }
};

template <typename E>
class Cols< SparseMatrix<E> >
   : public redirected_container< Cols< SparseMatrix<E> >,
	       cons< Container< sparse2d::Cols<E, false, sparse2d::full, sparse_matrix_line> >,
		     Hidden< SparseMatrix<E> > >
            > {
public:
   typename Cols::container& get_container()
   {
      return reinterpret_cast<typename Cols::container&>(*this->hidden().data);
   }
   const typename Cols::container& get_container() const
   {
      return reinterpret_cast<const typename Cols::container&>(*this->hidden().data);
   }
   void resize(int n)
   {
      this->hidden().data->resize_cols(n);
   }
};

/// sparse matrix statistics collection
template <typename E>
class SparseMatrixStatistics {
public:
   unsigned int maxnon0s, maxrowsize, maxcolsize;
   E maxabs;

   SparseMatrixStatistics()
      : maxnon0s(0), maxrowsize(0), maxcolsize(0), maxabs(0) { }

   void gather(const SparseMatrix<E>& m)
   {
      unsigned int non0s=0;
      for (typename Entire< Rows< SparseMatrix<E> > >::const_iterator r=entire(rows(m)); !r.at_end(); ++r) {
	 if (unsigned int s=r->size()) {
	    for (typename Entire< typename SparseMatrix<E>::row_type >::const_iterator e=entire(*r); !e.at_end(); ++e)
	       maxabs=std::max(maxabs, std::abs(*e));
	    maxrowsize=std::max(maxrowsize, s);
	    non0s+=s;
	 }
      }
      maxnon0s=std::max(maxnon0s, non0s);
      for (typename Entire< Cols< SparseMatrix<E> > >::const_iterator c=entire(cols(m)); !c.at_end(); ++c) {
	 if (unsigned int s=c->size()) {
	    maxcolsize=std::max(maxcolsize, s);
	 }
      }
   }

   void gather(const Transposed< SparseMatrix<E> >& m)
   {
      gather(m.hidden());
   }

   // statistics at two various moments can be glued together
   SparseMatrixStatistics& operator+= (const SparseMatrixStatistics& s)
   {
      maxnon0s=std::max(maxnon0s,s.maxnon0s);
      maxabs=std::max(maxabs,s.maxabs);
      maxrowsize=std::max(maxrowsize,s.maxrowsize);
      maxcolsize=std::max(maxcolsize,s.maxcolsize);
      return *this;
   }

   template <typename Traits> friend
   std::basic_ostream<char, Traits>&
   operator<< (std::basic_ostream<char, Traits>& os, const SparseMatrixStatistics& s)
   {
      return os << ">>> " << s.maxnon0s << " nonzeroes,  max abs(element)=" << s.maxabs
		<< "\n>>> max row size=" << s.maxrowsize << ",  max col size=" << s.maxcolsize << endl;
   }
};

} // end namespace pm

namespace polymake {
   using pm::SparseMatrix;
   using pm::RestrictedSparseMatrix;
}

namespace std {
   template <typename E> inline
   void swap(pm::SparseMatrix<E>& M1, pm::SparseMatrix<E>& M2)
   {
      M1.swap(M2);
   }

   template <typename E, pm::sparse2d::restriction_kind restriction> inline
   void swap(pm::RestrictedSparseMatrix<E,restriction>& M1,
	     pm::RestrictedSparseMatrix<E,restriction>& M2)
   {
      M1.swap(M2);
   }

   template <typename tree_type> inline
   void swap(pm::sparse_matrix_line<tree_type>& l1, pm::sparse_matrix_line<tree_type>& l2)
   {
      l1.swap(l2);
   }
}

#endif // _POLYMAKE_SPARSE_MATRIX_H

// Local Variables:
// mode:C++
// End:
