/* Copyright (c) 1997-2007
   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_MATRIX_H
#define _POLYMAKE_MATRIX_H "$Project: polymake $$Id: Matrix.h 7563 2007-01-16 12:07:34Z gawrilow $"

#ifndef _POLYMAKE_GENERIC_MATRIX_H
#  include <GenericMatrix.h>
#endif
#ifndef _POLYMAKE_SHARED_OBJECT_H
#  include <shared_object.h>
#endif

namespace pm {

template <typename E>
class Matrix : public GenericMatrix<Matrix<E>, E> {
protected:
   int dimr, dimc;
   shared_array<E> data;
public:
   typedef E value_type;
   typedef E& reference;
   typedef const E& const_reference;
   static const bool prefer_concat=true;

   Matrix() : dimr(0), dimc(0) { }

   Matrix(int r, int c) : dimr(r), dimc(c), data(r*c) { }

   Matrix(int r, int c, typename function_argument<E>::type init)
      : dimr(r), dimc(c), data(r*c, constant(init).begin()) { }

   /** Matrix with given dimensions.  Elements are initialized from an input sequence.
       Elements are assumed to come in the row order.
   */
   template <typename Iterator>
   Matrix(int r, int c, Iterator src)
      : dimr(r), dimc(c), data(r*c, construct_matching_iterator<E>(src, (dense*)0)) { }

   /// Copy of a disguised Matrix object.
   Matrix(const GenericMatrix<Matrix>& M)
      : dimr(M.rows()), dimc(M.cols()), data(M.top().data) { }

   /// Copy of an abstract matrix of the same element type.
   template <typename Matrix2>
   Matrix(const GenericMatrix<Matrix2,E>& M)
      : dimr(M.rows()), dimc(M.cols()),
	data(dimr*dimc, ensure(concat_rows(M), (dense*)0).begin()) { }

   /// Copy of an abstract matrix with element conversion.
   template <typename Matrix2, typename E2>
   explicit Matrix(const GenericMatrix<Matrix2,E2>& M)
      : dimr(M.rows()), dimc(M.cols()),
	data(dimr*dimc, ensure(attach_converter<E>(concat_rows(M)), (dense*)0).begin()) { }

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

   using Matrix::generic_type::operator=;

   void swap(Matrix& M)
   {
      data.swap(M.data);
      std::swap(dimr,M.dimr);
      std::swap(dimc,M.dimc);
   }

   friend void relocate(Matrix* from, Matrix* to)
   {
      pm::relocate(&from->data, &to->data);
      to->dimr=from->dimr;
      to->dimc=from->dimc;
   }

   /// Resize to new dimensions, added elements initialized with default constructor.
   void resize(int r, int c)
   {
      if (c==dimc) {
	 data.resize(r*c);  dimr=r;
      } else if (c<dimc && r<=dimr) {
	 *this=this->minor(sequence(0,r),sequence(0,c));
      } else {
	 Matrix M(r,c);
	 if (c<dimc)
	    M.minor(sequence(0,dimr),All)=this->minor(All,sequence(0,c));
	 else
	    M.minor(sequence(0,std::min(dimr,r)), sequence(0,dimc))=this->minor(sequence(0,std::min(dimr,r)),All);
	 *this=M;
      }
   }

   template <typename E2>
   void assign(int r, int c, const E2& x)
   {
      data.assign(r*c, constant(x).begin());
      dimr=r; dimc=c;
   }

   void clear() { dimr=dimc=0; data.clear(); }
   void clear(int r, int c) { dimr=r; dimc=c; data.resize(r*c); }

   int rows() const { return dimr; }
   int cols() const { return dimc; }

   E& operator() (int i, int j)
   {
      INDEX_CHECK((i<0 || i>=dimr), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=dimc), "operator() - column index out of range");
      return (*data)[i*dimc+j];
   }

   const E& operator() (int i, int j) const
   {
      INDEX_CHECK((i<0 || i>=dimr), "operator() - row index out of range");
      INDEX_CHECK((j<0 || j>=dimc), "operator() - column index out of range");
      return (*data)[i*dimc+j];
   }

   template <typename Matrix2>
   Matrix& operator/= (const GenericMatrix<Matrix2,E>& M)
   {
      const int add_rows=M.rows();
      if (add_rows) {
	 DIMENSION_CHECK((data.size() && dimc != M.cols()), "operator/= - column dimensions mismatch");
	 if (data.has_aliases()) {
	    *this = const_cast<const Matrix&>(*this) / M;
	 } else {
	    if (!data.size()) {
	       dimr=0; dimc=M.cols();
	    }
	    data.append(add_rows*dimc, ensure(concat_rows(M), (dense*)0).begin());
	    dimr+=add_rows;
	 }
      }
      return *this;
   }

   template <typename Vector2>
   Matrix& operator/= (const GenericVector<Vector2,E>& v)
   {
      DIMENSION_CHECK((data.size() && dimc != v.dim()), "operator/= - column dimensions mismatch");
      if (data.has_aliases()) {
	 *this = const_cast<const Matrix&>(*this) / v;
      } else {
	 if (!data.size()) {
	    dimr=0; dimc=v.dim();
	 }
	 data.append(v.dim(), ensure(v.top(), (dense*)0).begin());
	 ++dimr;
      }
      return *this;
   }

   template <typename Matrix2>
   Matrix& operator|= (const GenericMatrix<Matrix2,E>& M)
   {
      const int add_cols=M.cols();
      if (add_cols) {
	 DIMENSION_CHECK((data.size() && dimr != M.rows()), "operator|= - row dimensions mismatch");
	 if (data.has_aliases()) {
	    *this = const_cast<const Matrix&>(*this) | M;
	 } else {
	    if (!data.size()) {
	       dimr=M.rows(); dimc=0;
	    }
	    data.weave(dimr*add_cols, dimc, pm::rows(M).begin());
	    dimc+=add_cols;
	 }
      }
      return *this;
   }

   template <typename Vector2>
   Matrix& operator|= (const GenericVector<Vector2,E>& v)
   {
      DIMENSION_CHECK((data.size() && dimr != v.dim()), "operator|= - row dimensions mismatch");
      if (data.has_aliases()) {
	 *this = const_cast<const Matrix&>(*this) | v;
      } else {
	 if (!data.size()) {
	    dimr=v.dim(); dimc=0;
	 }
	 data.weave(dimr, dimc, ensure(v.top(), (dense*)0).begin());
	 ++dimc;
      }
      return *this;
   }
protected:
   friend class ConcatRows<Matrix>;
   friend class GenericMatrix<Matrix>;
   friend class alias_to_shared<Matrix>;

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

   const shared_array<E>& data_member() const { return data; }

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

   template <typename Matrix2>
   void assign(const GenericMatrix<Matrix2>& M)
   {
      const int r=M.rows(), c=M.cols();
      data.assign(r*c, ensure(concat_rows(M), (dense*)0).begin());
      dimr=r; dimc=c;
   }

   template <typename E2>
   void _fill(const E2& x, True)
   {
      data.assign(data.size(), constant(x).begin());
   }
};

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

template <typename E>
class ConcatRows< Matrix<E> >
   : public GenericVector< ConcatRows< Matrix<E> >, E> {
public:
   Matrix<E>& hidden() { return *reinterpret_cast<Matrix<E>*>(this); }
   const Matrix<E>& hidden() const { return *reinterpret_cast<const Matrix<E>*>(this); }

   typedef E value_type;
   typedef E& reference;
   typedef const E& const_reference;
   typedef E* iterator;
   typedef const E* const_iterator;
   typedef std::reverse_iterator<iterator> reverse_iterator;
   typedef std::reverse_iterator<const_iterator> const_reverse_iterator;

   int size() const { return hidden().data.size(); }

   iterator begin() { return *hidden().data; }
   iterator end() { return *hidden().data+size(); }
   const_iterator begin() const { return *hidden().data; }
   const_iterator end() const { return *hidden().data+size(); }

   reverse_iterator rbegin() { return reverse_iterator(end()); }
   reverse_iterator rend() { return reverse_iterator(begin()); }
   const_reverse_iterator rbegin() const { return const_reverse_iterator(end()); }
   const_reverse_iterator rend() const { return const_reverse_iterator(begin()); }

   E& operator[] (int i) { return begin()[i]; }
   const E& operator[] (int i) const { return begin()[i]; }

   ConcatRows& operator= (const ConcatRows& m)
   {
      return ConcatRows::generic_type::operator=(m);
   }
   using ConcatRows::generic_type::operator=;
};

template <typename E>
class Matrix_array
   : public array_as_container<E,1>,
     public GenericVector<Matrix_array<E>, E> { };

template <typename ElemRef, bool for_Rows>
class matrix_line_factory {
   typedef typename deref<ElemRef>::type E;
   const Matrix<E> *_matrix;

   template <typename,bool> friend class matrix_line_factory;
public:
   typedef ElemRef argument_type;
   typedef typename inherit_ref<Matrix_array<E>, ElemRef>::type array;
   typedef typename inherit_const<IndexedSlice<array, Series<int,for_Rows> >, ElemRef>::type result_type;

   matrix_line_factory(const Matrix<E> *matrix_arg=0) : _matrix(matrix_arg) { }

   matrix_line_factory(const matrix_line_factory<E&, for_Rows>& f) : _matrix(f._matrix) { }

   result_type operator() (ElemRef start) const
   {
      return result_type(reinterpret_cast<array>(start),
			 Series<int,for_Rows>(0, for_Rows ? _matrix->cols() : _matrix->rows(), for_Rows ? 1 : _matrix->cols()));
   }
};
template <bool for_Rows> struct Build_matrix_line_factory { };

template <bool for_Rows, typename Iterator, typename ElemRef>
struct unary_op_builder<Build_matrix_line_factory<for_Rows>, Iterator, ElemRef> {
   typedef matrix_line_factory<ElemRef, for_Rows> operation;

   static operation create(const Build_matrix_line_factory<for_Rows>& M)
   {
      return reinterpret_cast<const Matrix<typename deref<ElemRef>::type>*>(&M);
   }
   static const operation& create(const operation& op) { return op; }

   typedef typename if_else<attrib<ElemRef>::is_const,
			    const matrix_line_factory<typename attrib<ElemRef>::minus_const, for_Rows>, operation>::type
      alt_operation;
   static operation create(alt_operation& op) { return operation(op); }
};

template <typename E, bool for_Rows>
class Matrix_first_line
   : public indexed_subset_impl< Matrix_first_line<E,for_Rows>,
				 list( Container1< ConcatRows< Matrix<E> > >,
				       Container2< Series<int,!for_Rows> >,
				       Hidden< Matrix<E> > ) > {
   typedef indexed_subset_impl<Matrix_first_line> _super;
protected:
   ~Matrix_first_line();
public:
   typename _super::container1& get_container1() { return concat_rows(this->hidden()); }
   const typename _super::container1& get_container1() const { return concat_rows(this->hidden()); }

   typename _super::container2 get_container2() const
   {
      return typename _super::container2(0, for_Rows ? this->hidden().rows() : this->hidden().cols(),
					    for_Rows ? this->hidden().cols() : 1);
   }
};

template <typename E>
class Rows< Matrix<E> >
   : public modified_container_impl< Rows< Matrix<E> >,
				     list( Container< Matrix_first_line<E, true> >,
					   Operation< Build_matrix_line_factory<true> >,
					   MasqueradedTop ) > {
   typedef modified_container_impl<Rows> _super;
protected:
   ~Rows();
public:
   const typename _super::operation& get_operation() const
   {
      return reinterpret_cast<const typename _super::operation&>(this->hidden());
   }
   void resize(int n)
   {
      this->hidden().resize(n, this->hidden().cols());
   }
};

template <class E>
class Cols< Matrix<E> >
   : public modified_container_impl< Cols< Matrix<E> >,
				     list( Container< Matrix_first_line<E, false> >,
					   Operation< Build_matrix_line_factory<false> >,
					   MasqueradedTop ) > {
   typedef modified_container_impl<Cols> _super;
protected:
   ~Cols();
public:
   const typename _super::operation& get_operation() const
   {
      return reinterpret_cast<const typename _super::operation&>(this->hidden());
   }
   void resize(int n)
   {
      this->hidden().resize(this->hidden().rows(), n);
   }
};

} // end namespace pm

namespace polymake {
   using pm::Matrix;
}

namespace std {
   template <typename E> inline
   void swap(pm::Matrix<E>& m1, pm::Matrix<E>& m2) { m1.swap(m2); }
}

#endif // _POLYMAKE_MATRIX_H

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