/* Copyright (c) 1997-2004
   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,v 1.48 2004/09/13 18:09:36 gawrilow Exp $"

#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_serialized=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, (cons<dense, reference_relaxed>*)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(serialize(M), (cons<dense, reference_relaxed>*)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>(serialize(M)), (cons<dense, reference_relaxed>*)0).begin()) { }

protected:
   template <typename RowCursor>
   void read(RowCursor row_cursor);

public:
   template <typename Input> friend
   Input& operator>> (GenericInput<Input>& is, Matrix& M)
   {
      M.read(is.top().begin_list(&pm::rows(M)));
      return is.top();
   }

   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(); }

   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)
   {
      DIMENSION_CHECK((data.size() && dimc != M.cols()), "operator/= - column dimensions mismatch");
      if (data.has_aliases()) {
	 *this = const_cast<const Matrix&>(*this) / M;
      } else {
	 const int r=M.rows();
	 if (!data.size()) {
	    dimr=0; dimc=M.cols();
	 }
	 data.append(r*dimc, ensure(serialize(M), (cons<dense, reference_relaxed>*)0).begin());
	 dimr+=r;
      }
      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(), (cons<dense, reference_relaxed>*)0).begin());
	 ++dimr;
      }
      return *this;
   }

   template <typename Matrix2>
   Matrix& operator|= (const GenericMatrix<Matrix2,E>& M)
   {
      DIMENSION_CHECK((data.size() && dimr != M.rows()), "operator|= - row dimensions mismatch");
      if (data.has_aliases()) {
	 *this = const_cast<const Matrix&>(*this) | M;
      } else {
	 const int c=M.cols();
	 if (!data.size()) {
	    dimr=M.rows(); dimc=0;
	 }
	 data.weave(dimr*c, dimc, rows(M).begin());
	 dimc+=c;
      }
      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(), (cons<dense, reference_relaxed>*)0).begin());
	 ++dimc;
      }
      return *this;
   }
protected:
   friend class Serialized<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(serialize(M), (cons<dense, reference_relaxed>*)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>
template <typename RowCursor>
void Matrix<E>::read(RowCursor row_cursor)
{
   dimr=row_cursor.size();
   if (dimr) {
      dimc=row_cursor.begin_list((typename Matrix::row_type*)0).size();
      data.resize(dimc*dimr);
      for (typename Entire< Rows<Matrix> >::iterator r=entire(pm::rows(*this)); !r.at_end(); ++r)
	 row_cursor.retrieve_array(*r);
      
   } else {
      dimc=0;
      data.clear();
   }
}

template <typename E>
class Serialized< Matrix<E> >
   : public GenericVector< Serialized< 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]; }

   Serialized& operator= (const Serialized& m)
   {
      return Serialized::generic_type::operator=(m);
   }
   using Serialized::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() (argument_type 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>::result
      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>,
	       cons< Container1< Serialized< Matrix<E> > >,
	       cons< 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 serialize(this->hidden()); }
   const typename _super::container1& get_container1() const { return serialize(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> >,
	       cons< Container< Matrix_first_line<E, true> >,
               cons< Accessor< Build_matrix_line_factory<true> >,
                     Hidden< Matrix<E> > > >
            > {
   typedef modified_container_impl<Rows> _super;
protected:
   ~Rows();
public:
   typename _super::container& get_container()
   {
      return reinterpret_cast<typename _super::container&>(this->hidden());
   }
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(this->hidden());
   }
   const typename _super::accessor& get_accessor() const
   {
      return reinterpret_cast<const typename _super::accessor&>(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> >,
               cons< Container< Matrix_first_line<E, false> >,
               cons< Accessor< Build_matrix_line_factory<false> >,
		     Hidden< Matrix<E> > > >
            > {
   typedef modified_container_impl<Cols> _super;
protected:
   ~Cols();
public:
   typename _super::container& get_container()
   {
      return reinterpret_cast<typename _super::container&>(this->hidden());
   }
   const typename _super::container& get_container() const
   {
      return reinterpret_cast<const typename _super::container&>(this->hidden());
   }
   const typename _super::accessor& get_accessor() const
   {
      return reinterpret_cast<const typename _super::accessor&>(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:
