/* 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_CDD_INTERFACE_IMPL_H
#define _POLYMAKE_CDD_INTERFACE_IMPL_H "$Project: polymake $$Id: cdd_interface_impl.h,v 1.2 2004/06/22 16:23:15 gawrilow Exp $"

#include <cdd_interface.h>

// cddlib headers:
extern "C" {
#include "setoper.h"
#include "cdd.h"
}

namespace polymake { namespace polytope { namespace cdd_interface {

template <typename Coord> struct traits;

#ifdef GMPRATIONAL
template <>
struct traits<Rational> {
   static const dd_NumberType number_type=dd_Rational;

   // PTL -> cddlib
   static void store(mytype *dst, const Rational *src)
   {
      dd_set(*dst, src->get_rep());
   }

   // cddlib -> PTL
   static void fetch(Rational* dst, const mytype* src)
   {
      dst->set(*src);
   }
};
#else
template <>
struct traits<double> {
   static const dd_NumberType number_type=dd_Real;

   static void store(mytype *dst, const double *src)
   {
      dd_set_d(*dst, *src);
   }

   static void fetch(double *dst, const mytype* src)
   {
      *dst=dd_get_d(*const_cast<mytype*>(src));		// cddlib is not const-clean
   }
};
#endif

template <typename Coord> class cdd_polyhedron;
template <typename Coord> class cdd_lp;
template <typename Coord> class cdd_lp_sol;
template <typename Coord> class cdd_matrix;

template <typename Coord>
class cdd_vector : public traits<Coord> {
private:
   int dim;
   dd_Arow ptr;
public:
   explicit cdd_vector(int dim_arg) : dim(dim_arg)
   {
      dd_InitializeArow(dim, &ptr);
   }

   ~cdd_vector() { dd_FreeArow(dim, ptr); }

   Vector<Coord> get(int start_at=0) const;

   friend class cdd_matrix<Coord>;
   friend class cdd_polyhedron<Coord>;
   friend class cdd_lp<Coord>;
   friend class solver<Coord>;
};

template <typename Coord>
Vector<Coord> cdd_vector<Coord>::get(int start_at) const
{
   Vector<Coord> ret(dim-start_at);
   typename Vector<Coord>::iterator dst=ret.begin();
   for (mytype *src=ptr+start_at, *end=ptr+dim;  src!=end;  ++src, ++dst)
      fetch(dst, src);
   return ret;
}

template <typename Coord>
class cdd_matrix : public traits<Coord> {
private:
   dd_MatrixPtr ptr;
public:
   // build from points and rays
   cdd_matrix(const Matrix<Coord>& P);
   // build from inequalities and equations
   cdd_matrix(const Matrix<Coord>& I, const Matrix<Coord>& E);
   // build from computed convex hull
   cdd_matrix(const cdd_polyhedron<Coord>&, bool primal);

   ~cdd_matrix() { dd_FreeMatrix(ptr); }

   void add_objective(const Vector<Coord>& obj, bool maximize);

   typename solver<Coord>::matrix_pair facets_and_ah() const;

   Matrix<Coord> vertices() const;

   ListMatrix< Vector<Coord> > vertex_normals(Bitset& Vertices)
      throw (std::bad_alloc, std::runtime_error);

   friend class cdd_polyhedron<Coord>;
   friend class cdd_lp<Coord>;
   friend class solver<Coord>;
};

extern dd_ErrorType err;

template <typename Coord>
class cdd_polyhedron {
private:
   dd_PolyhedraPtr ptr;
public:
   cdd_polyhedron(const cdd_matrix<Coord>& M);
   ~cdd_polyhedron() { dd_FreePolyhedra(ptr); }

   friend class cdd_matrix<Coord>;
   friend class solver<Coord>;
};

template <typename Coord>
class cdd_lp : public traits<Coord> {
private:
   dd_LPPtr ptr;
public:
   cdd_lp(const cdd_matrix<Coord>& M);
   ~cdd_lp() { dd_FreeLPData(ptr); }

   Vector<Coord> optimal_vertex() const;
   friend class cdd_lp_sol<Coord>;
};

template <typename Coord>
class cdd_lp_sol {
private:
   dd_LPSolutionPtr ptr;
public:
   cdd_lp_sol(const cdd_lp<Coord>& L);
   ~cdd_lp_sol() { dd_FreeLPSolution(ptr); }

   Coord optimal_value() const;
};

template <typename Coord>
cdd_matrix<Coord>::cdd_matrix(const Matrix<Coord>& P)
   : ptr(dd_CreateMatrix(P.rows(), P.cols()))
{
   // get size of the input matrix
   int m = P.rows();
   int n = P.cols();

   ptr->representation = dd_Generator;    // Input type: points
   ptr->numbtype = this->number_type;
  
   // copy data: polymake -> cdd
   typename pm::Serialized< Matrix<Coord> >::const_iterator p=serialize(P).begin();
   for (mytype **r=ptr->matrix, **rend=r+m;  r!=rend;  ++r)
      for (mytype *c=*r, *cend=c+n;  c!=cend;  ++c, ++p)
	 store(c,p);
}

template <typename Coord>
cdd_matrix<Coord>::cdd_matrix(const Matrix<Coord>& I, const Matrix<Coord>& E)
   : ptr(dd_CreateMatrix(I.rows() + E.rows(), I.cols()))
{
   // get size of the input matrix
   int mi = I.rows();
   int me = E.rows();
   int n  = I.cols();

   ptr->representation = dd_Inequality;   // Input type: inequalities
   ptr->numbtype = this->number_type;
  
   typename pm::Serialized< Matrix<Coord> >::const_iterator p=serialize(I).begin();
   mytype **r, **rend;
   for (r=ptr->matrix, rend=r+mi;  r!=rend;  ++r)
      for (mytype *c=*r, *cend=c+n;  c!=cend;  ++c, ++p)
	 store(c,p);

   p=serialize(E).begin();
   ++mi; // cddlib starts counting at 1
   for (rend+=me;  r!=rend;  ++r, ++mi) {
      for (mytype *c=*r, *cend=c+n;  c!=cend;  ++c, ++p)
	 store(c,p);
      set_addelem(ptr->linset, mi);
   }
}

template <typename Coord>
void cdd_matrix<Coord>::add_objective(const Vector<Coord>& obj, bool maximize)
{
   typename Vector<Coord>::const_iterator o=obj.begin();
   for (mytype *r=ptr->rowvec, *rend=r+obj.size();  r!=rend;  ++r, ++o)
      store(r,o);
   ptr->objective= maximize ? dd_LPmax : dd_LPmin;
}

template <typename Coord>
cdd_polyhedron<Coord>::cdd_polyhedron(const cdd_matrix<Coord>& M)
   : ptr(dd_DDMatrix2Poly(M.ptr, &err))
{
   if (err != dd_NoError) {
      std::ostringstream err_msg;
      err_msg << "Error in dd_DDMatrix2Poly: " << err << endl;
      throw std::runtime_error(err_msg.str());
   }
}

template <typename Coord>
cdd_lp<Coord>::cdd_lp(const cdd_matrix<Coord>& M)
   : ptr(dd_Matrix2LP(M.ptr, &err))
{
   if (err != dd_NoError) {
      std::ostringstream err_msg;
      err_msg << "Error in dd_Matrix2LP: " << err << endl;
      throw std::runtime_error(err_msg.str());
   }
   
   if (! dd_LPSolve(ptr, dd_DualSimplex, &err)) {
      std::ostringstream err_msg;
      err_msg << "Error in dd_LPSolve: " << err << endl;
      throw std::runtime_error(err_msg.str());
   }
}

template <typename Coord>
cdd_matrix<Coord>::cdd_matrix(const cdd_polyhedron<Coord>& P, bool primal)
   : ptr(primal ? dd_CopyInequalities(P.ptr) : dd_CopyGenerators(P.ptr)) { }

template <typename Coord>
Matrix<Coord> cdd_matrix<Coord>::vertices() const
{
   if (set_card(ptr->linset)) throw not_pointed();

   dd_rowrange m = ptr->rowsize;
   dd_colrange n = ptr->colsize;
   if (m<=0) throw infeasible();

   Matrix<Coord> V(m,n);

   typename pm::Serialized< Matrix<Coord> >::iterator v=serialize(V).begin();
   for (mytype **r=ptr->matrix, **rend=r+m;  r!=rend;  ++r) {
      for (mytype *c=*r, *cend=c+n;  c!=cend;  ++c, ++v)
	 fetch(v,c);
   }
   return V;
}

template <typename Coord>
typename solver<Coord>::matrix_pair
cdd_matrix<Coord>::facets_and_ah() const
{
   dd_rowrange m = ptr->rowsize;
   dd_rowrange n = ptr->colsize;
   long linsize = set_card(ptr->linset);

   Matrix<Coord> F(m - linsize, n);
   Matrix<Coord> AH(linsize, n);

   typename pm::Serialized< Matrix<Coord> >::iterator f=serialize(F).begin(), ah=serialize(AH).begin();
   int index=1;		// cddlib starts counting at 1
   for (mytype **r=ptr->matrix, **rend=r+m;  r!=rend;  ++r, ++index) {
      if (set_member(index, ptr->linset))
	 // the current row is an equation of the lineality space
	 for (mytype *c=*r, *cend=c+n;  c!=cend;  ++c, ++ah)
	    fetch(ah,c);
      else
	 // the current row defines a facet
	 for (mytype *c=*r, *cend=c+n;  c!=cend;  ++c, ++f)
	    fetch(f,c);
   }

   if (AH.rows()+1 == AH.cols()) F.resize(0,0);		// cddlib does not correctly treate the special case of a single point
   return typename solver<Coord>::matrix_pair(F,AH);
}


template <typename Coord>
ListMatrix< Vector<Coord> >
cdd_matrix<Coord>::vertex_normals(Bitset& Vertices)
   throw (std::bad_alloc, std::runtime_error)
{
   ListMatrix< Vector<Coord> > VN;
   typename Rows< ListMatrix< Vector<Coord> > >::iterator vn_front=rows(VN).begin();
   cdd_vector<Coord> cert(ptr->colsize+1);
   for (int i=ptr->rowsize; i>=1; --i) {
      const bool is_redundant=dd_Redundant(ptr, i, cert.ptr, &err);
      if (err != dd_NoError) {
	 std::ostringstream err_msg;
	 err_msg << "Error in dd_Redundant: " << err << endl;
	 throw std::runtime_error(err_msg.str());
      }
      if (is_redundant) {
	 dd_MatrixRowRemove(&ptr, i);
      } else {
	 Vertices += i-1;
	 Vector<Coord> vertex_normal=cert.get(1);
	 vertex_normal[0]=0; vertex_normal.negate();
	 vn_front=VN.insert_row(vn_front, vertex_normal);
      }
   }
   return VN;
}

template <typename Coord>
cdd_lp_sol<Coord>::cdd_lp_sol(const cdd_lp<Coord>& L)
   : ptr(dd_CopyLPSolution(L.ptr))
{
   switch (ptr->LPS) {
   case dd_Optimal: 
       break;
   case dd_StrucInconsistent:   // if the problem is inconsistent then mark that there is no valid point
   case dd_Inconsistent:       
      throw infeasible();
   case dd_Unbounded:      // if the problem is unbounded then mark that there is no optimal solution
      throw unbounded();
   case dd_StrucDualInconsistent:  // if the dual problem is either unbounded or
   case dd_DualUnbounded:          // inconsistent we don't know whether the primal 
   case dd_DualInconsistent:       // is inconsistent or unbounded
      throw std::runtime_error("cannot handle lp solution: problem is either inconsistent or unbounded");
   default:  // all other cases are handled here ....
      std::ostringstream err_msg;
      err_msg << "cannot handle lp solution: cdd returned: " << ptr->LPS;
      throw std::runtime_error(err_msg.str());
   }
}

template <typename Coord>
Vector<Coord> cdd_lp<Coord>::optimal_vertex() const
{
   const int n=ptr->d;
   Vector<Coord> Opt(n);
   typename Vector<Coord>::iterator o=Opt.begin();
   for (mytype *s=ptr->sol, *send=s+n;  s!=send;  ++s, ++o)
      fetch(o,s);
   return Opt;
}

template <typename Coord>
struct solver<Coord>::initializer {
   initializer() { dd_set_global_constants(); }
};

template <typename Coord>
solver<Coord>::solver() {
   static initializer init;
}

template <typename Coord>
typename solver<Coord>::matrix_pair
solver<Coord>::enumerate_facets(const Matrix<Coord>& Points)
   throw (std::bad_alloc, std::runtime_error)
{
   cdd_matrix<Coord> IN(Points);
   cdd_polyhedron<Coord> P(IN);
   return cdd_matrix<Coord>(P, true).facets_and_ah();
}

template <typename Coord>
Matrix<Coord>
solver<Coord>::enumerate_vertices(const Matrix<Coord>& Inequalities, const Matrix<Coord>& Equations)
   throw (std::bad_alloc, std::runtime_error)
{
   cdd_matrix<Coord> IN(Inequalities, Equations);
   cdd_polyhedron<Coord> P(IN);
   return cdd_matrix<Coord>(P, false).vertices();
}

template <typename Coord>
typename solver<Coord>::non_redundant
solver<Coord>::find_vertices_among_points(const Matrix<Coord>& Points)
   throw (std::bad_alloc, std::runtime_error)
{
   cdd_matrix<Coord> IN(Points);
   Bitset Vertices(Points.rows());
   return non_redundant(Vertices, IN.vertex_normals(Vertices));
}

template <typename Coord>
typename solver<Coord>::lp_solution
solver<Coord>::solve_lp(const Matrix<Coord>& Inequalities, const Matrix<Coord>& Equations,
			const Vector<Coord>& Objective, bool maximize)
   throw (std::bad_alloc, std::runtime_error)
{
   cdd_matrix<Coord> IN(Inequalities, Equations);
   IN.add_objective(Objective, maximize);
   cdd_lp<Coord> LP(IN);
   cdd_lp_sol<Coord> Sol(LP);
   return lp_solution(Sol.optimal_value(), LP.optimal_vertex());
}

} } }

#endif // _POLYMAKE_CDD_INTERFACE_IMPL_H

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