/* 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_LINALG_H
#define _POLYMAKE_LINALG_H "$Project: polymake $$Id: linalg.h,v 1.29 2004/09/20 16:32:28 gawrilow Exp $"

#ifndef _POLYMAKE_MATRIX_H
#  include <Matrix.h>
#endif
#ifndef _POLYMAKE_SPARSE_MATRIX_H
#  include <SparseMatrix.h>
#endif
#ifndef _POLYMAKE_SPARSE_VECTOR_H
#  include <SparseVector.h>
#endif
#ifndef _POLYMAKE_LIST_MATRIX_H
#  include <ListMatrix.h>
#endif
#ifndef _POLYMAKE_SET_H
#  include <Set.h>
#endif
#include <linalg_exceptions.h>
#include <dense_linalg.h>
#include <sparse_linalg.h>

namespace pm { namespace operations {

template <typename OpRef,
          typename Discr=typename object_traits<typename deref<OpRef>::type>::generic_tag>
struct normalize_impl;

template <typename OpRef,
          typename Discr=typename object_traits<typename deref<OpRef>::type>::generic_tag>
struct canonicalize_impl;

template <typename OpRef>
struct normalize_impl<OpRef, is_vector>
   : div_impl<OpRef, typename deref<OpRef>::type::element_type, false, cons<is_vector, is_scalar> > {

   typedef typename deref<OpRef>::type argument_type;
   typedef div_impl<OpRef, typename argument_type::element_type, false, cons<is_vector, is_scalar> > _super;

   typename _super::result_type operator() (typename function_argument<OpRef>::type v) const
   {
      return _super::operator()(v, sqrt(sqr(v)));
   }

   typename function_argument<OpRef>::lvalue
   assign(typename function_argument<OpRef>::lvalue v) const
   {
      return v /= sqrt(sqr(v));
   }
};

template <typename OpRef>
struct canonicalize_impl<OpRef, is_vector>
   : div_impl<OpRef, typename deref<OpRef>::type::element_type, false, cons<is_vector, is_scalar> > {

   typedef typename deref<OpRef>::type argument_type;
protected:
   typedef typename argument_type::element_type scalar_type;
   typedef div_impl<OpRef, scalar_type, false, cons<is_vector, is_scalar> > _super;
public:

   typename _super::result_type operator() (typename function_argument<OpRef>::type v) const
   {
      typename Entire<argument_type>::const_iterator it=std::find_if(entire(v), conv<scalar_type, bool>());
      return _super::operator()(v, it.at_end() ? scalar_type(1) : std::abs(*it));
   }

   typename function_argument<OpRef>::lvalue
   assign(typename function_argument<OpRef>::lvalue v) const
   {
      typename Entire<argument_type>::iterator it=std::find_if(entire(v), conv<scalar_type, bool>());
      if (!it.at_end()) {
	 const scalar_type leading=std::abs(*it);
	 if (leading != 1)
	    do
	       *it /= leading;
	    while (!(++it).at_end());
      }
      return v;
   }
};

template <typename OpRef>
struct normalize : normalize_impl<OpRef> { };

template <typename OpRef>
struct canonicalize : canonicalize_impl<OpRef> { };

} } // end namespace pm::operations

namespace polymake { namespace operations {
   typedef BuildUnary<pm::operations::normalize> normalize;
   typedef BuildUnary<pm::operations::canonicalize> canonicalize;
} }

namespace pm {

/// Divide each vector in a sequence thru its length (L2-norm)
template <typename Iterator> inline
void normalize(Iterator dst)
{
   perform_assign(dst, polymake::operations::normalize());
}

template <typename _Matrix> inline
typename _Matrix::persistent_type
normalized(const GenericMatrix<_Matrix>& A)
{
   return typename _Matrix::persistent_type(A.rows(), A.cols(),
					    attach_operation(rows(A), polymake::operations::normalize()).begin());
}

/// Divide each vector in a sequence thru the abs of its first nonzero element
template <typename Iterator>
void canonicalize(Iterator dst)
{
   perform_assign(dst, polymake::operations::canonicalize());
}

/// Compute the determinant of a matrix using the Gauss elimination method
template <typename _Matrix> inline
typename _Matrix::element_type det(const GenericMatrix<_Matrix>& A)
{
   return det(typename _Matrix::persistent_type(A));
}

/// Reduce a vector with a given matrix using the Gauss elimination method
template <typename _Matrix, typename _Vector, typename E> inline
typename _Vector::persistent_type
reduce(GenericMatrix<_Matrix, E>& A, const GenericVector<_Vector, E>& V)
{
   return reduce(typename _Matrix::persistent_type(A), typename _Vector::persistent_type(V));
}

/** Compute the inverse matrix $A^-1$ using the Gauss elimination method.
    @exception degenerate_matrix if det(A)==0
*/
template <typename _Matrix> inline
typename _Matrix::persistent_type
inv(const GenericMatrix<_Matrix>& A)
{
   return inv(typename _Matrix::persistent_type(A));
}

/** Solve the linear system A*x==B
    @return x
    @exception degenerate_matrix if det(A) == 0
    @exception infeasible if rank(A) != rank(A|B)
*/
template <typename _Matrix, typename _Vector, typename E> inline
Vector<E>
lin_solve(const GenericMatrix<_Matrix, E>& A, const GenericVector<_Vector, E>& B)
{
   return lin_solve(typename _Matrix::persistent_type(A), Vector<E>(B));
}

/** Compute the basis of the subspaces spanned by a sequence of vectors and orthogonal one.
    @param v input iterator over the vectors
    @param row_basis_consumer output iterator consuming the indices of basis rows (=vectors)
    @param col_basis_consumer output iterator consuming the indices of basis columns (=coordinates)
*/
template <typename VectorIterator, typename RowBasisOutputIterator, typename ColBasisOutputIterator, typename E>
void
null_space(VectorIterator v, RowBasisOutputIterator row_basis_consumer, ColBasisOutputIterator col_basis_consumer,
	   ListMatrix< SparseVector<E> >& H)
{
   typedef ListMatrix< SparseVector<E> > AH_matrix;
   conv<E,bool> non_zero;

   for (int i=0; H.rows()>0 && !v.at_end(); ++v, ++i) {
      typename Rows<AH_matrix>::iterator h=rows(H).begin(), h_end=rows(H).end();
      E pivot=E();
      for (; h!=h_end && !non_zero(pivot=(*h)*(*v)); ++h) ;
      if (h==h_end) continue;

      typename Rows<AH_matrix>::iterator h2=h;
      for (++h2; h2!=h_end; ++h2) {
	 const E x=(*h2)*(*v);
	 if (non_zero(x)) *h2 -= (x/pivot)*(*h);
      }
      *row_basis_consumer++ = i;
      if (!derived_from_instance<ColBasisOutputIterator,black_hole>::answer)
	 *col_basis_consumer++ = h->rbegin().index();
      H.delete_row(h);
   }
}

template <typename _Matrix, typename E> inline
typename _Matrix::persistent_type
null_space(const GenericMatrix<_Matrix, E>& M)
{
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(M.cols());
   null_space(entire(rows(M)), black_hole<int>(), black_hole<int>(), H);
   return H;
}

template <typename _Vector, typename E> inline
ListMatrix< SparseVector<E> >
null_space(const GenericVector<_Vector, E>& V)
{
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(V.dim());
   null_space(entire(item2container(V.top())), black_hole<int>(), black_hole<int>(), H);
   return H;
}

/// @param req_sign expected sign of det( null_space(V) / V )
template <typename _Vector, typename E>
ListMatrix< SparseVector<E> >
null_space_oriented(const GenericVector<_Vector, E>& V, int req_sign)
{
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(V.dim());
   null_space(entire(item2container(V.top())), black_hole<int>(), black_hole<int>(), H);
   typename ensure_features<_Vector, pure_sparse>::const_iterator v_pivot=ensure(V.top(), (pure_sparse*)0).begin();
   if (v_pivot.at_end() && req_sign)
      throw infeasible("null_space_oriented: zero vector has no orientation");
   if ((sign(*v_pivot)==req_sign) == ((v_pivot.index()+V.dim()+1)%2))
      rows(H).back().negate();
   return H;
}

template <typename VectorIterator> inline
Set<int>
basis_vectors(VectorIterator v)
{
   typedef typename iterator_traits<VectorIterator>::value_type::element_type E;
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(v.at_end() ? 0 : v->dim());
   Set<int> b;
   null_space(v, std::back_inserter(b), black_hole<int>(), H);
   return b;
}

template <typename _Matrix, typename E> inline
Set<int>
basis_rows(const GenericMatrix<_Matrix,E>& M)
{
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(M.cols());
   Set<int> b;
   null_space(entire(rows(M)), std::back_inserter(b), black_hole<int>(), H);
   return b;
}

template <typename _Matrix> inline
Set<int>
basis_cols(const GenericMatrix<_Matrix>& M)
{
   return basis_rows(T(M));
}

template <typename _Matrix, typename E> inline
std::pair< Set<int>, Set<int> >
basis(const GenericMatrix<_Matrix, E>& M)
{
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(M.cols());
   Set<int> br, bc;
   null_space(entire(rows(M)), std::back_inserter(br), std::back_inserter(bc), H);
   return std::make_pair(br,bc);
}

template <typename _Matrix, typename E> inline
int rank(const GenericMatrix<_Matrix, E>& M)
{
   if (M.rows() <= M.cols()) {
      ListMatrix< SparseVector<E> > H=unit_matrix<E>(M.rows());
      null_space(entire(cols(M)), black_hole<int>(), black_hole<int>(), H);
      return M.rows()-H.rows();
   }
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(M.cols());
   null_space(entire(rows(M)), black_hole<int>(), black_hole<int>(), H);
   return M.cols()-H.rows();
}

template <typename _Matrix> inline
Set<int>
basis_rows(const GenericMatrix<_Matrix, double>& M)
{
   ListMatrix< SparseVector<double> > H=unit_matrix<double>(M.cols());
   Set<int> b;
   null_space(entire(attach_operation(rows(M), polymake::operations::normalize())),
	      std::back_inserter(b), black_hole<int>(), H);
   return b;
}

template <typename _Matrix> inline
std::pair< Set<int>, Set<int> >
basis(const GenericMatrix<_Matrix, double>& M)
{
   ListMatrix< SparseVector<double> > H=unit_matrix<double>(M.cols());
   Set<int> br, bc;
   null_space(entire(attach_operation(rows(M), polymake::operations::normalize())),
	      std::back_inserter(br), std::back_inserter(bc), H);
   return std::make_pair(br,bc);
}

template <typename _Matrix> inline
int rank(const GenericMatrix<_Matrix, double>& M)
{
   if (M.rows() <= M.cols()) {
      ListMatrix< SparseVector<double> > H=unit_matrix<double>(M.rows());
      null_space(entire(attach_operation(cols(M), polymake::operations::normalize())),
		 black_hole<int>(), black_hole<int>(), H);
      return M.rows()-H.rows();
   }
   ListMatrix< SparseVector<double> > H=unit_matrix<double>(M.cols());
   null_space(entire(attach_operation(rows(M), polymake::operations::normalize())),
	      black_hole<int>(), black_hole<int>(), H);
   return M.cols()-H.rows();
}

/// The same as basis(), but ignoring the first column of the matrix.
template <typename _Matrix, typename E> inline
std::pair< Set<int>, Set<int> >
basis_affine(const GenericMatrix<_Matrix, E>& M)
{
   int ad=M.cols()-1;
   ListMatrix< SparseVector<E> > H=unit_matrix<E>(ad);
   Set<int> br, bc;
   null_space(entire(rows(M.minor(All,range(1,ad)))), std::back_inserter(br),
	      make_output_transform_iterator(std::back_inserter(bc), operations::fix2<int, operations::add<int,int> >(1)),
	      H);
   return std::make_pair(br,bc);
}

/// Project u to v.
template <typename E, typename _Vector1, typename _Vector2> inline
typename _Vector2::persistent_type
proj(const GenericVector<_Vector1,E>& u, const GenericVector<_Vector2,E>& v)
{
   return (u*v)/sqr(v)*v;
}

template <typename VectorIterator, typename OutputIterator>
void orthogonalize(VectorIterator v, OutputIterator sqr_consumer)
{
   typedef typename iterator_traits<VectorIterator>::value_type vector_type;
   typedef typename vector_type::element_type E;
   // the first vector will not be modified
   while (!v.at_end()) {
      const E s=sqr(*v);
      if (s) {
	 VectorIterator v2=v;
	 for (++v2; !v2.at_end(); ++v2) {
	    const E x=(*v2) * (*v);
	    if (x) *v2 -= (x/s) * (*v);
	 }
      }
      *sqr_consumer++ = s;
      ++v;
   }
}

template <typename VectorIterator, typename OutputIterator>
void orthogonalize_affine(VectorIterator v, OutputIterator sqr_consumer)
{
   typedef typename iterator_traits<VectorIterator>::value_type vector_type;
   typedef typename vector_type::element_type E;
   // the first vector will not be modified
   while (!v.at_end()) {
      const E s=sqr(v->slice(1));
      if (s) {
	 VectorIterator v2=v;
	 for (++v2; !v2.at_end(); ++v2) {
	    const E x=v2->slice(1) * v->slice(1);
	    if (x) *v2 -= (x/s) * (*v);
	 }
      }
      *sqr_consumer++ = s;
      ++v;
   }
}

/// Apply the Gram-Schmidt orthogonalization to the vector sequence.
template <typename VectorIterator> inline
void orthogonalize(VectorIterator v)
{
   orthogonalize(v, black_hole<typename iterator_traits<VectorIterator>::value_type::element_type>());
}

/** The same as orthogonalize(.), but making the affine parts of the resulting vectors
    (without 0-th coordinate) orthogonal
*/
template <typename VectorIterator> inline
void orthogonalize_affine(VectorIterator v)
{
   orthogonalize_affine(v, black_hole<typename iterator_traits<VectorIterator>::value_type::element_type>());
}

/// Find row indices of all far points (that is, having zero in the first column).
template <typename _Matrix> inline
Set<int>
far_points(const GenericMatrix<_Matrix>& M)
{
   return indices(attach_selector(M.col(0), polymake::operations::logical_not()));
}

/// Find indices of rows orthogonal to the given vector
template <typename E, typename _Matrix, typename _Vector> inline
Set<int>
orthogonal_rows(const GenericMatrix<_Matrix,E>& M, const GenericVector<_Vector,E>& v)
{
   return indices(attach_selector(attach_operation(rows(M), constant(v), polymake::operations::mul()),
				  polymake::operations::logical_not()));
}

} // end namespace pm

namespace polymake {
   using pm::null_space;
   using pm::basis_vectors;
   using pm::orthogonalize;
   using pm::orthogonalize_affine;
   using pm::normalize;
   using pm::canonicalize;
}

#endif // _POLYMAKE_LINALG_H

// Local Variables:
// mode:C++
// c-basis-offset:3
// End:
