/* Copyright (c) 1997-2006
   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_INCIDENCE_MATRIX_H
#define _POLYMAKE_INCIDENCE_MATRIX_H "$Project: polymake $$Id: IncidenceMatrix.h,v 1.32 2006/03/09 16:21:35 gawrilow Exp $"

#ifndef _POLYMAKE_GENERIC_INCIDENCE_MATRIX_H
#  include <GenericIncidenceMatrix.h>
#endif
#ifndef _POLYMAKE_SPARSE2D_H
#  include <sparse2d.h>
#endif
#ifndef _POLYMAKE_SET_H
#  include <Set.h>
#endif

namespace pm {

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

   typedef typename if_else<_transpose, typename _Matrix::col_type, typename _Matrix::row_type>::result row_set;
   typedef bool2type<_transpose> Transp;

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

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

   bool get() const
   {
      return row((Transp*)0).exists(j);
   }

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

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

   void toggle()
   {
      row((Transp*)0).toggle(j);
   }
public:
   typedef bool value_type;

   incidence_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 Tree>
class incidence_line
   : public modified_tree< incidence_line<Tree>,
	       cons< Hidden< sparse2d::line<Tree> >,
		     Accessor< BuildUnaryIt<operations::index2element> > >
            >,
     public GenericMutableSet< incidence_line<Tree>, int, operations::cmp<> > {
protected:
   ~incidence_line();
public:
   using modified_tree<incidence_line>::get_comparator_object;

   incidence_line& operator= (const incidence_line& l)
   {
      return incidence_line::generic_mutable_type::operator=(l);
   }
   using incidence_line::generic_mutable_type::operator=;

   template <class _Set>
   incidence_line& operator= (const Complement<_Set, int, operations::cmp<> >& cl)
   {
      return incidence_line::generic_mutable_type::operator=(sequence(0, this->dim()) * cl);
   }

   bool contains (int e) const { return this->exists(e); }
};

template <typename Tree>
struct check_container_feature<incidence_line<Tree>, sparse_compatible> : True { };

template <sparse2d::restriction_kind restriction=sparse2d::only_rows>
class RestrictedIncidenceMatrix
   : public matrix_methods<RestrictedIncidenceMatrix<restriction>, bool> {
protected:
   typedef sparse2d::Table<nothing, 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 bool value_type;
   typedef bool reference;
   typedef const bool const_reference;
   static const bool prefer_concat=false;

   RestrictedIncidenceMatrix() : data(nothing()) { }
   explicit RestrictedIncidenceMatrix(int n) : data(n) { }

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

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

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

protected:
   typedef incidence_proxy_base<RestrictedIncidenceMatrix, restriction==sparse2d::only_cols> proxy_base;
   typedef sparse_elem_proxy<proxy_base> proxy;

   bool _random(int i, int j, True*) const
   {
      return this->row(i).exists(j);
   }
   bool _random(int i, int j, False*) const
   {
      return this->col(j).exists(i);
   }
public:
   proxy operator() (int i, int j)
   {
      return proxy_base(this,i,j);
   }

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

   bool exists(int i, int j) const
   {
      return _random(i, j, (bool2type<restriction==sparse2d::only_rows>*)0);
   }
private:
   template <typename RowCol, typename _Set>
   void _append(RowCol& row_col, const _Set& set, int i)
   {
      for (typename Entire<_Set>::const_iterator s=entire(set); !s.at_end(); ++s)
	 row_col[*s].push_back(i);
   }

   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>
   RestrictedIncidenceMatrix& operator/= (const GenericIncidenceMatrix<_Matrix>& m)
   {
      _append_rows(m.rows(), entire(pm::rows(m)), (bool2type<restriction==sparse2d::only_rows>*)0);
      return *this;
   }

   template <typename _Set>
   RestrictedIncidenceMatrix& operator/= (const GenericSet<_Set, int, operations::cmp<> >& s)
   {
      _append_rows(1, entire(item2container(s.top())), (bool2type<restriction==sparse2d::only_rows>*)0);
      return *this;
   }

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

   template <typename _Set>
   RestrictedIncidenceMatrix& operator|= (const GenericSet<_Set, int, operations::cmp<> >& s)
   {
      _append_cols(1, entire(item2container(s.top())), (bool2type<restriction==sparse2d::only_cols>*)0);
      return *this;
   }

#ifdef POLY_DEBUG
   void check() const { data.check(); }
#endif
   friend class Rows<RestrictedIncidenceMatrix>;
   friend class Cols<RestrictedIncidenceMatrix>;
   template <bool> friend class IncidenceMatrix;
};

template <sparse2d::restriction_kind restriction>
class Rows< RestrictedIncidenceMatrix<restriction> >
   : public redirected_container< Rows< RestrictedIncidenceMatrix<restriction> >,
	       cons< Container< sparse2d::Rows<nothing, false, restriction, incidence_line> >,
                     Hidden< RestrictedIncidenceMatrix<restriction> > >
            > {
   typedef redirected_container<Rows> _super;
public:
   typename _super::container& get_container()
   {
      return reinterpret_cast<typename _super::container&>(this->hidden().data);
   }
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(this->hidden().data);
   }

   int size() const
   {
      return this->hidden().data.rows();
   }
   bool empty() const
   {
      return !size();
   }
   void resize(int n)
   {
      this->hidden().data.resize_rows(n);
   }
};

template <sparse2d::restriction_kind restriction>
class Cols< RestrictedIncidenceMatrix<restriction> >
   : public redirected_container< Cols< RestrictedIncidenceMatrix<restriction> >,
	       cons< Container< sparse2d::Cols<nothing, false, restriction, incidence_line> >,
                     Hidden< RestrictedIncidenceMatrix<restriction> > >
            > {
   typedef redirected_container<Cols> _super;
public:
   typename _super::container& get_container()
   {
      return reinterpret_cast<typename _super::container&>(this->hidden().data);
   }
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(this->hidden().data);
   }

   int size() const
   {
      return this->hidden().data.cols();
   }
   bool empty() const {
      return !size();
   }
   void resize(int n)
   {
      this->hidden().data.resize_cols(n);
   }
};

template <sparse2d::restriction_kind restriction>
struct spec_object_traits< RestrictedIncidenceMatrix<restriction> >
   : spec_object_traits<is_container> {
   static const int dimension=2;

   typedef typename if_else<restriction==sparse2d::only_rows,
			    Rows< RestrictedIncidenceMatrix<restriction> >,
			    Cols< RestrictedIncidenceMatrix<restriction> > >::result serialized;

   static serialized& serialize(RestrictedIncidenceMatrix<restriction>& M)
   {
      return reinterpret_cast<serialized&>(M);
   }
   static const serialized& serialize(const RestrictedIncidenceMatrix<restriction>& M)
   {
      return reinterpret_cast<const serialized&>(M);
   }
};

template <bool symmetric>
class IncidenceMatrix : public GenericIncidenceMatrix< IncidenceMatrix<symmetric>, bool2type<symmetric> > {
protected:
   typedef sparse2d::Table<nothing,symmetric> table_type;
   shared_object<table_type> data;

   // initialize from a dense boolean sequence in row order
   template <typename Iterator>
   void _init(Iterator src, True*, True*)
   {
      const int n=this->cols();
      for (typename Entire< Rows<IncidenceMatrix> >::iterator r_i=entire(pm::rows(*this));
	   !r_i.at_end(); ++r_i) {
	 int i=0;
	 if (symmetric) {
	    i=r_i.index();
	    std::advance(src,i);
	 }
	 for (; i<n; ++i, ++src)
	    if (*src) r_i->push_back(i);
      }
   }

   // initialize rowwise from a sequence of sets
   template <typename Iterator>
   void _init(Iterator src, False*, True*)
   {
      std::copy(src, entire(pm::rows(*this)));
   }

public:
   IncidenceMatrix() { }

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

   template <typename Iterator>
   IncidenceMatrix(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==0)>*)0, (True*)0);
   }

   IncidenceMatrix(const GenericIncidenceMatrix<IncidenceMatrix>& M)
      : data(M.top().data) { }

   template <typename _Matrix>
   IncidenceMatrix(const GenericIncidenceMatrix<_Matrix>& M)
      : data( make_constructor(M.rows(), M.cols(),(table_type*)0) )
   {
      _init(pm::rows(M).begin(), (False*)0, (bool2type<!symmetric || _Matrix::is_symmetric>*)0);
   }

   template <sparse2d::restriction_kind restriction>
   explicit IncidenceMatrix(RestrictedIncidenceMatrix<restriction>& M)
      : data( make_constructor(M.data,(table_type*)0) ) { }

private:
   template <typename Input>
   void input(GenericInput<Input>& is, False*)
   {
      RestrictedIncidenceMatrix<> Mt;
      is.top() >> Mt;
      *this=Mt;
   }

   template <typename Input>
   void input(GenericInput<Input>& is, True*)
   {
      is.top() >> pm::rows(*this);
   }

   template <typename Input>
   void input(GenericInput<Input>& is)
   {
      input(is, (bool2type<symmetric>*)0);
   }
public:
   template <typename Input> friend
   Input& operator>> (GenericInput<Input>& is, IncidenceMatrix& M)
   {
      M.input(is);
      return is.top();
   }

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

   using IncidenceMatrix::generic_type::operator=;

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

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

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

   void resize(int m, int n) { data->resize(m,n); }

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

protected:
   typedef sparse_elem_proxy< incidence_proxy_base<IncidenceMatrix> > proxy;
public:
   proxy operator() (int i, int j)
   {
      return incidence_proxy_base<IncidenceMatrix>(this,i,j);
   }

   bool operator() (int i, int j) const
   {
      return exists(i,j);
   }

   bool exists(int i, int j) const
   {
      return data->row(i).exists(j);
   }

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

   template <typename _Set>
   IncidenceMatrix& operator/= (const GenericSet<_Set, int, operations::cmp<> >& s)
   {
      DIMENSION_CHECK((!s.top().empty() && s.top().back() >= this->cols()), "append_row - set too long");
      if (is_shared()) {
	 *this = const_cast<const IncidenceMatrix&>(*this) / s;
      } else {
	 int oldrows=this->rows();
	 data.get()->resize_rows(oldrows+1);
	 this->row(oldrows)=s;
      }
      return *this;
   }

   template <typename _Set>
   IncidenceMatrix& operator/= (const Complement<_Set, int, operations::cmp<> >& s)
   {
      DIMENSION_CHECK((!s.base().empty() && s.base().back() >= this->cols()), "append_row - set too long");
      if (is_shared()) {
	 *this = const_cast<const IncidenceMatrix&>(*this) / (sequence(0,this->cols()) * s);
      } else {
	 int oldrows=this->rows();
	 data.get()->resize_rows(oldrows+1);
	 this->row(oldrows)=sequence(0,this->cols()) * s;
      }
      return *this;
   }

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

   template <typename _Set>
   IncidenceMatrix& operator |= (const GenericSet<_Set, int, operations::cmp<> >& s)
   {
      DIMENSION_CHECK((!s.top().empty() && s.top().back() >= this->rows()), "append_col - set too long");
      if (is_shared()) {
	 *this = const_cast<const IncidenceMatrix&>(*this) | s;
      } else {
	 int oldcols=this->cols();
	 data.get()->resize_cols(oldcols+1);
	 this->col(oldcols)=s;
      }
      return *this;
   }

   template <typename _Set>
   IncidenceMatrix& operator|= (const Complement<_Set, int, operations::cmp<> >& s)
   {
      DIMENSION_CHECK((!s.base().empty() && s.base().back() >= this->rows()), "append_col - set too long");
      if (is_shared()) {
	 *this = const_cast<const IncidenceMatrix&>(*this) | (sequence(0,this->rows()) * s);
      } else {
	 int oldcols=this->cols();
	 data.get()->resize_cols(oldcols+1);
	 this->col(oldcols)=sequence(0,this->rows()) * s;
      }
      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<IncidenceMatrix>;
   friend class Cols<IncidenceMatrix>;
   friend class GenericIncidenceMatrix<IncidenceMatrix>;
   friend class alias_to_shared<IncidenceMatrix>;

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

   template <typename Matrix>
   void assign(const GenericIncidenceMatrix<Matrix>& M, True*)
   {
      if (is_shared() || this->rows() != M.rows() || this->cols() != M.cols()) {
	 assign(IncidenceMatrix(M));
      } else {
	 GenericIncidenceMatrix<IncidenceMatrix>::assign(M, (True*)0);
      }
   }

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

template <bool symmetric>
class alias_ptr< IncidenceMatrix<symmetric> >
   : public alias_to_shared< IncidenceMatrix<symmetric> > {
public:
   alias_ptr(const IncidenceMatrix<symmetric>* ptr_arg)
      : alias_to_shared< IncidenceMatrix<symmetric> >(ptr_arg) { }
};

template <bool symmetric>
class Rows< IncidenceMatrix<symmetric> >
   : public redirected_container< Rows< IncidenceMatrix<symmetric> >,
	       cons< Container< sparse2d::Rows<nothing, symmetric, sparse2d::full, incidence_line> >,
                     Hidden< IncidenceMatrix<symmetric> > >
            > {
   typedef redirected_container<Rows> _super;
public:
   typename _super::container& get_container()
   {
      return reinterpret_cast<typename _super::container&>(*this->hidden().data);
   }
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(*this->hidden().data);
   }
   void resize(int n)
   {
      this->hidden().data->resize_rows(n);
   }
};

template <bool symmetric>
class Cols< IncidenceMatrix<symmetric> >
   : public redirected_container< Cols< IncidenceMatrix<symmetric> >,
	       cons< Container< sparse2d::Cols<nothing, symmetric, sparse2d::full, incidence_line> >,
                     Hidden< IncidenceMatrix<symmetric> > >
            > {
   typedef redirected_container<Cols> _super;
public:
   typename _super::container& get_container()
   {
      return reinterpret_cast<typename _super::container&>(*this->hidden().data);
   }
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(*this->hidden().data);
   }
   void resize(int n)
   {
      this->hidden().data->resize_cols(n);
   }
};

// convolution of two incidence relations
template <typename _Matrix1, typename _Matrix2> inline
IncidenceMatrix<>
convolute(const GenericIncidenceMatrix<_Matrix1>& M1, const GenericIncidenceMatrix<_Matrix2>& M2)
{
   DIMENSION_CHECK((M1.cols() != M2.rows()), "convolute - dimensions mismatch");

   IncidenceMatrix<> result(M1.rows(), M2.cols());
   typename Rows<_Matrix1>::const_iterator r1=rows(M1).begin();
   for (typename Entire< Rows< IncidenceMatrix<> > >::iterator dst=entire(rows(result));
	!dst.at_end();  ++dst, ++r1)
      accumulate_in(entire(rows(M2.minor(*r1,All))), BuildBinary<operations::add>(), *dst);

   return result;
}

} // end namespace pm

namespace polymake {
   using pm::IncidenceMatrix;
   using pm::RestrictedIncidenceMatrix;
}

namespace std {
   template <bool symmetric> inline
   void swap(pm::IncidenceMatrix<symmetric>& M1, pm::IncidenceMatrix<symmetric>& M2) { M1.swap(M2); }

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

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

#endif // _POLYMAKE_INCIDENCE_MATRIX_H

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