/* Copyright (c) 1997-2009
   Ewgenij Gawrilow, Michael Joswig (Technische Universitaet Darmstadt, Germany)
   http://www.polymake.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.
--------------------------------------------------------------------------------
   $Project: polymake $$Id: linsolver.h 9298 2009-09-02 23:10:56Z gawrilow $
*/

#ifndef _POLYMAKE_LINSOLVER_H
#define _POLYMAKE_LINSOLVER_H

#include <polymake.h>
#include <Matrix.h>
#include <Vector.h>
#include <linalg_exceptions.h>

namespace polymake { namespace polytope {

class not_pointed : public linalg_error {
protected:
   int lin_dim;
public:
   not_pointed(int dim_arg)
      : linalg_error("polyhedron not pointed"),
	lin_dim(dim_arg) {}

   int lineality_dim() const { return lin_dim; }
};

class unbounded : public linalg_error {
public:
   unbounded() : linalg_error("unbounded linear program") {}
};

template <typename Solver>
void ch_primal(perl::Object& p, Solver& solver)
{
   typedef typename Solver::coord_type coord_type;
   const Matrix<coord_type> Points=p.give("VERTICES | POINTS");
   typename Solver::matrix_pair F=solver.enumerate_facets(Points);
   p.take("FACETS") << F.first;
   p.take("AFFINE_HULL") << F.second;
}

template <typename Solver>
void count_facets(perl::Object& p, Solver& solver)
{
   typedef typename Solver::coord_type coord_type;
   const Matrix<coord_type> Points=p.give("VERTICES | POINTS");
   p.take("N_FACETS") << solver.count_facets(Points);
}

template <typename Solver>
void ch_dual(perl::Object& p, Solver& solver)
{
   typedef typename Solver::coord_type coord_type;
   const Matrix<coord_type> H=p.give("FACETS | INEQUALITIES"),
                            EQ=p.lookup("AFFINE_HULL | EQUATIONS");

   // recall that all output properties lacking after the rule completion automatically get undefined values
   if (H.rows() || EQ.rows()) {
      try {
	 const Matrix<coord_type> V=solver.enumerate_vertices(H,EQ);
	 p.take("VERTICES") << V;
	 p.take("FEASIBLE") << true;
	 p.take("POINTED") << true;
	 p.take("LINEALITY_DIM") << 0;
      }
      catch (const not_pointed& ex) {
	 p.take("FEASIBLE") << true;
	 p.take("POINTED") << false;
	 p.take("LINEALITY_DIM") << ex.lineality_dim();
      }
      catch (infeasible) {
	 p.take("FEASIBLE") << false;
      }
   } else {
      // empty polytope
      p.take("FEASIBLE") << false;
   }
}

template <typename Solver>
void count_vertices(perl::Object& p, Solver& solver, bool only_bounded=false)
{
   typedef typename Solver::coord_type coord_type;
   const Matrix<coord_type> H=p.give("FACETS | INEQUALITIES"),
                            EQ=p.lookup("AFFINE_HULL | EQUATIONS");
   if (H.rows() || EQ.rows()) {
      try {
	 std::pair<long,long> n_vertices=solver.count_vertices(H,EQ,only_bounded);
	 if (!only_bounded) p.take("N_VERTICES") << n_vertices.first;
	 p.take("N_BOUNDED_VERTICES") << n_vertices.second;
	 p.take("FEASIBLE") << true;
	 p.take("POINTED") << true;
	 p.take("LINEALITY_DIM") << 0;
      }
      catch (const not_pointed& ex) {
	 if (!only_bounded) p.take("N_VERTICES") << 0;
	 p.take("N_BOUNDED_VERTICES") << 0;
	 p.take("FEASIBLE") << true;
	 p.take("LINEALITY_DIM") << ex.lineality_dim();
      }
      catch (infeasible) {
	 if (!only_bounded) p.take("N_VERTICES") << 0;
	 p.take("N_BOUNDED_VERTICES") << 0;
	 p.take("FEASIBLE") << false;
      }
   } else {
      // empty polytope
      if (!only_bounded) p.take("N_VERTICES") << 0;
      p.take("N_BOUNDED_VERTICES") << 0;
      p.take("FEASIBLE") << false;
   }
}

} }

#endif // _POLYMAKE_LINSOLVER_H

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