/* Copyright (c) 1997-2010
   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: core_point_algo.cc 9891 2010-11-10 18:01:52Z gawrilow $
*/

#include <Integer.h>
#include <Rational.h>
#include <Matrix.h>
#include <polymake.h>


namespace polymake {

  namespace {
    inline Integer round_down(const Rational rat) {
      const Integer truncated_rat(rat);
      return (truncated_rat!=rat && rat<0) ? truncated_rat-1 : truncated_rat;
    }
    
    inline bool feasible(const Vector<Rational>& point, const Matrix<Rational>& Inequalities) {
      for(Entire< Rows< Matrix<Rational> > >::const_iterator i=entire(rows(Inequalities)); !i.at_end(); ++i)
	if (point*(*i)<0) return false;
      return true;
    }

  }

namespace polytope {

  //Maximization with utility vector c=(1,...,1)
  //Assumption: the LP is feasible and bounded in the direction of c
  perl::ListReturn find_transitive_lp_sol(const Matrix<Rational>& Inequalities){
    int N_cols=Inequalities.cols();
    int dim=N_cols-1;
    const Rational mininf=-999999999; //FIXME: What is appropriate here?
    const Rational maxinf=999999999; //FIXME: What is appropriate here?
    Rational lower_bound=mininf;
    Rational upper_bound=maxinf;
    bool max_bounded=true;
    bool feasible=true;
    Rational a=maxinf;
    for(Entire< Rows< Matrix<Rational> > >::const_iterator i=entire(rows(Inequalities)); !i.at_end(); ++i){
      Rational sum=0;
      for(int j=1;j<N_cols;j++){ //add up all coefficients of one ineq (=project to fixed space spanned by (1,...,1))
	sum+=(*i)[j];
      }
      Rational constraint_val=-(*i)[0]/sum;
      Rational bound=dim*constraint_val;
      if(sum>0 && bound>lower_bound){ //new lower bound
	lower_bound=bound;
      } else if(sum<0 &&bound<upper_bound){ //new upper bound
	a=constraint_val;
	upper_bound=bound;
      }  
    }
    if(upper_bound<lower_bound){
      feasible=false;
    } else if(upper_bound==maxinf){
      max_bounded=false;
    }
    Vector<Rational> optLPsolution(1|same_element_vector<Rational>(a,dim));
    Rational optLPvalue=upper_bound;
    perl::ListReturn result;
    result << optLPsolution
	   << optLPvalue
	   << feasible
	   << max_bounded;
    return result;
  }


  perl::ListReturn core_point_algo(perl::Object p, const Rational optLPvalue, perl::OptionSet options)
  {
    const Matrix<Rational> Inequalities=p.give("FACETS|INEQUALITIES");
    const int n=p.give("AMBIENT_DIM");
    const Integer d(round_down(optLPvalue/n));
    const Integer::div_t this_division(round_down(optLPvalue),n);
    const Integer r_start(this_division.rem);
    
    const bool verbose = options["verbose"];
    
    if (verbose)
      cout << "dimension=" << n << ", "
	   << "LP approximation=" << optLPvalue << ", "
	   << "d=" << d << endl;

    Integer r(r_start);
    bool ILP_not_feasible(true);
    Vector<Rational> optILPsolution;
    Rational optILPvalue;
    while (r>=0 && ILP_not_feasible) {
      if (verbose)
	cout << "trying r=" << r << endl;
      const Vector<Rational> core_point(1|same_element_vector<Rational>(d+1,r)|same_element_vector<Rational>(d,n-r));
      if (feasible(core_point,Inequalities)) {
	optILPsolution = core_point;
	optILPvalue=(d+1)*r+d*(n-r);
	ILP_not_feasible=false;
      } else
	--r;
    }
    
    perl::ListReturn result;
    
    if (ILP_not_feasible) {
      if (verbose)
	cout << "ILP infeasible" << endl;
    } else {
      if (verbose)
	cout << "optimal ILP solution=" << optILPsolution << ", "
	     << "value=" << optILPvalue << endl;
      result << optILPsolution
	     << optILPvalue;
    }
    
    return result;
  }

UserFunction4perl("# @category Optimization"
		  "# Algorithm to solve highly symmetric integer linear programs (ILP)."
		  "# It is required that the group of the ILP induces the alternating or symmetric group"
		  "# on the set of coordinate directions."
		  "# The linear objective function is the vector (0,1,1,..,1)."
		  "# "
		  "# @param Polytope p"
		  "# @param Rational optLPvalue optimal value of LP approximation"
		  "# @option Bool verbose"
		  "# @return perl::ListReturn (optimal solution, optimal value) or empty",
		  &core_point_algo, "core_point_algo(Polytope, $; {verbose => undef})");

UserFunction4perl("# @category Optimization"
		  "# Algorithm to solve symmetric linear programs (LP) of the form"
		  "# max c<sup>t</sup>x , c=(0,1,1,..,1)"
		  "# subject to the inequality system given by //Inequalities//."
		  "# It is required that the symmetry group of the LP acts transitively"
		  "# on the coordinate directions."
		  "# "
		  "# @param Matrix Inequalities the inequalities describing the feasible region"
		  "# @return perl::ListReturn (optLPsolution,optLPvalue,feasible,max_bounded)",
		  &find_transitive_lp_sol, "find_transitive_lp_sol(Matrix)");
} }
