// -*- c++ -*-

 #ifndef LEMON_LP_SOLVER_WRAPPER_H

 #define LEMON_LP_SOLVER_WRAPPER_H

 ///\ingroup misc

 ///\file

 ///\brief Dijkstra algorithm.

 // #include <stdio.h>

 #include <stdlib.h>

 #include <iostream>

 #include <map>

 // #include <stdio>

 //#include <stdlib>

 extern "C" {

 #include "glpk.h"

 }

 #include <iostream>

 #include <vector>

 #include <string>

 #include <list>

 #include <memory>

 #include <utility>

 //#include <sage_graph.h>

 //#include <lemon/list_graph.h>

 //#include <lemon/graph_wrapper.h>

 #include <lemon/invalid.h>

 //#include <bfs_dfs.h>

 //#include <stp.h>

 //#include <lemon/max_flow.h>

 //#include <augmenting_flow.h>

 //#include <iter_map.h>

 using std::cout;

 using std::cin;

 using std::endl;

 namespace lemon {

   /// \addtogroup misc

   /// @{

   /// \brief A partitioned vector with iterable classes.

   ///

   /// This class implements a container in which the data is stored in an 

   /// stl vector, the range is partitioned into sets and each set is 

   /// doubly linked in a list. 

   /// That is, each class is iterable by lemon iterators, and any member of 

   /// the vector can bo moved to an other class.

   template <typename T>

   class IterablePartition {

   protected:

     struct Node {

       T data;

       int prev; //invalid az -1

       int next; 

     };

     std::vector<Node> nodes;

     struct Tip {

       int first;

       int last;

     };

     std::vector<Tip> tips;

   public:

     /// The classes are indexed by integers from \c 0 to \c classNum()-1.

     int classNum() const { return tips.size(); }

     /// This lemon style iterator iterates through a class. 

     class ClassIt;

     /// Constructor. The number of classes is to be given which is fixed 

     /// over the life of the container. 

     /// The partition classes are indexed from 0 to class_num-1. 

     IterablePartition(int class_num) { 

       for (int i=0; i<class_num; ++i) {

 	Tip t;

 	t.first=t.last=-1;

 	tips.push_back(t);

       }

     }

   protected:

     void befuz(ClassIt it, int class_id) {

       if (tips[class_id].first==-1) {

 	if (tips[class_id].last==-1) {

 	  nodes[it.i].prev=nodes[it.i].next=-1;

 	  tips[class_id].first=tips[class_id].last=it.i;

 	}

       } else {

 	nodes[it.i].prev=tips[class_id].last;

 	nodes[it.i].next=-1;

 	nodes[tips[class_id].last].next=it.i;

 	tips[class_id].last=it.i;

       }

     }

     void kifuz(ClassIt it, int class_id) {

       if (tips[class_id].first==it.i) {

 	if (tips[class_id].last==it.i) {

 	  tips[class_id].first=tips[class_id].last=-1;

 	} else {

 	  tips[class_id].first=nodes[it.i].next;

 	  nodes[nodes[it.i].next].prev=-1;

 	}

       } else {

 	if (tips[class_id].last==it.i) {

 	  tips[class_id].last=nodes[it.i].prev;

 	  nodes[nodes[it.i].prev].next=-1;

 	} else {

 	  nodes[nodes[it.i].next].prev=nodes[it.i].prev;

 	  nodes[nodes[it.i].prev].next=nodes[it.i].next;

 	}

       }

     }

   public:

     /// A new element with data \c t is pushed into the vector and into class 

     /// \c class_id.

     ClassIt push_back(const T& t, int class_id) { 

       Node n;

       n.data=t;

       nodes.push_back(n);

       int i=nodes.size()-1;

       befuz(i, class_id);

       return i;

     }

     /// A member is moved to an other class.

     void set(ClassIt it, int old_class_id, int new_class_id) {

       kifuz(it.i, old_class_id);

       befuz(it.i, new_class_id);

     }

     /// Returns the data pointed by \c it.

     T& operator[](ClassIt it) { return nodes[it.i].data; }

     /// Returns the data pointed by \c it.

     const T& operator[](ClassIt it) const { return nodes[it.i].data; }

     ///.

     class ClassIt {

       friend class IterablePartition;

     protected:

       int i;

     public:

       /// Default constructor.

       ClassIt() { }

       /// This constructor constructs an iterator which points

       /// to the member of th container indexed by the integer _i.

       ClassIt(const int& _i) : i(_i) { }

       /// Invalid constructor.

       ClassIt(const Invalid&) : i(-1) { }

     };

     /// First member of class \c class_id.

     ClassIt& first(ClassIt& it, int class_id) const {

       it.i=tips[class_id].first;

       return it;

     }

     /// Next member.

     ClassIt& next(ClassIt& it) const {

       it.i=nodes[it.i].next;

       return it;

     }

     /// True iff the iterator is valid.

     bool valid(const ClassIt& it) const { return it.i!=-1; }

   };

   template <typename _Col, typename _Value>

   class Expr;

   template <typename _Col, typename _Value>

   class SmallExpr {

     template <typename _C, typename _V> 

     friend class Expr;

   protected:

     _Col col;

     _Value value;

   public:

     SmallExpr(_Col _col) : col(_col), value(1) { 

     }

     SmallExpr& operator *= (_Value _value) {

       value*=_value;

       return *this;

     }

     //    template <typename _C, typename _V>

     //    friend SmallExpr<_C, _V> operator* (_V _value, 

     //					const SmallExpr<_C, _V>& expr);

     template <typename _C, typename _V>

     friend std::ostream& operator<<(std::ostream& os, 

 				    const SmallExpr<_C, _V>& expr);

   };

   template <typename _Col, typename _Value>

   SmallExpr<_Col, _Value> 

   operator* (_Value value, 

 	     const SmallExpr<_Col, _Value>& expr) {

     SmallExpr<_Col, _Value> tmp;

     tmp=expr;

     tmp*=value;

     return tmp;

   }

   template <typename _Col, typename _Value>

   std::ostream& operator<<(std::ostream& os, 

 			   const SmallExpr<_Col, _Value>& expr) {

     os << expr.value << "*" << expr.col;

     return os;

   }

   template <typename _Col, typename _Value>

   class Expr {

   protected:

     typedef 

     typename std::map<_Col, _Value> Data; 

     Data data;

   public:

     Expr() { }

     Expr(SmallExpr<_Col, _Value> expr) { 

       data.insert(std::make_pair(expr.col, expr.value));

     }

 //     Expr(_Col col) { 

 //       data.insert(std::make_pair(col, 1));

 //     }

     Expr& operator*=(_Value _value) {

       for (typename Data::iterator i=data.begin(); 

 	   i!=data.end(); ++i) {

 	(*i).second *= _value;

       }

       return *this;

     }

     Expr& operator+=(SmallExpr<_Col, _Value> expr) {

       typename Data::iterator i=data.find(expr.col);

       if (i==data.end()) {

 	data.insert(std::make_pair(expr.col, expr.value));

       } else {

 	(*i).second+=expr.value;

       }

       return *this;

     }

     //    template <typename _C, typename _V> 

     //    friend Expr<_C, _V> operator*(_V _value, const Expr<_C, _V>& expr);

     template <typename _C, typename _V> 

     friend std::ostream& operator<<(std::ostream& os, 

 				    const Expr<_C, _V>& expr);

   };

   template <typename _Col, typename _Value>

   Expr<_Col, _Value> operator*(_Value _value, 

 			       const Expr<_Col, _Value>& expr) {

     Expr<_Col, _Value> tmp;

     tmp=expr;

     tmp*=_value;

     return tmp;

   }

   template <typename _Col, typename _Value>

   std::ostream& operator<<(std::ostream& os, 

 			   const Expr<_Col, _Value>& expr) {

     for (typename Expr<_Col, _Value>::Data::const_iterator i=

 	   expr.data.begin(); 

 	 i!=expr.data.end(); ++i) {

       os << (*i).second << "*" << (*i).first << " ";

     }

     return os;

   }

   /*! \e

    */

   template <typename _Value>

   class LPSolverBase {

   public:

     /// \e

     typedef _Value Value;

     /// \e

     typedef IterablePartition<int>::ClassIt RowIt;

     /// \e

     typedef IterablePartition<int>::ClassIt ColIt;

   public:

     /// \e

     IterablePartition<int> row_iter_map;

     /// \e

     IterablePartition<int> col_iter_map;

     /// \e

     const int VALID_CLASS;

     /// \e

     const int INVALID_CLASS;

   public:

     /// \e

     LPSolverBase() : row_iter_map(2), 

 		     col_iter_map(2), 

 		     VALID_CLASS(0), INVALID_CLASS(1) { }

     /// \e

     virtual ~LPSolverBase() { }

     //MATRIX INDEPEDENT MANIPULATING FUNCTIONS

   public:

     /// \e

     virtual void setMinimize() = 0;

     /// \e

     virtual void setMaximize() = 0;

     //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS

   protected:

     /// \e

     virtual int _addRow() = 0;

     /// \e

     virtual int _addCol() = 0;

     /// \e

     virtual void _setRowCoeffs(int i, 

 			       std::vector<std::pair<int, double> > coeffs) = 0;

     /// \e

     virtual void _setColCoeffs(int i, 

 			       std::vector<std::pair<int, double> > coeffs) = 0;

     /// \e

     virtual void _eraseCol(int i) = 0;

     /// \e

     virtual void _eraseRow(int i) = 0;

   public:

     /// \e

     enum Bound { FREE, LOWER, UPPER, DOUBLE, FIXED };

   protected:

     /// \e

     virtual void _setColBounds(int i, Bound bound, 

 			       _Value lo, _Value up) = 0; 

     /// \e

     virtual void _setRowBounds(int i, Bound bound, 

 			       _Value lo, _Value up) = 0; 

     /// \e

     virtual void _setObjCoef(int i, _Value obj_coef) = 0;

     /// \e

     virtual _Value _getObjCoef(int i) = 0;

     //LOW LEVEL, SOLUTION RETRIEVING FUNCTIONS

   protected:

     virtual _Value _getPrimal(int i) = 0;

     //HIGH LEVEL INTERFACE, MATRIX MANIPULATING FUNTIONS

   public:

     /// \e

     RowIt addRow() {

       int i=_addRow(); 

       RowIt row_it;

       row_iter_map.first(row_it, INVALID_CLASS);

       if (row_iter_map.valid(row_it)) { //van hasznalhato hely

 	row_iter_map.set(row_it, INVALID_CLASS, VALID_CLASS);

 	row_iter_map[row_it]=i;

       } else { //a cucc vegere kell inzertalni mert nincs szabad hely

 	row_it=row_iter_map.push_back(i, VALID_CLASS);

       }

       return row_it;

     }

     /// \e

     ColIt addCol() {

       int i=_addCol();  

       ColIt col_it;

       col_iter_map.first(col_it, INVALID_CLASS);

       if (col_iter_map.valid(col_it)) { //van hasznalhato hely

 	col_iter_map.set(col_it, INVALID_CLASS, VALID_CLASS);

 	col_iter_map[col_it]=i;

       } else { //a cucc vegere kell inzertalni mert nincs szabad hely

 	col_it=col_iter_map.push_back(i, VALID_CLASS);

       }

       return col_it;

     }

     /// \e

     template <typename Begin, typename End>

     void setRowCoeffs(RowIt row_it, Begin begin, End end) {

       std::vector<std::pair<int, double> > coeffs;

       for ( ; begin!=end; ++begin) {

 	coeffs.push_back(std::

 			 make_pair(col_iter_map[begin->first], begin->second));

       }

       _setRowCoeffs(row_iter_map[row_it], coeffs);

     }

     /// \e

     template <typename Begin, typename End>

     void setColCoeffs(ColIt col_it, Begin begin, End end) {

       std::vector<std::pair<int, double> > coeffs;

       for ( ; begin!=end; ++begin) {

 	coeffs.push_back(std::

 			 make_pair(row_iter_map[begin->first], begin->second));

       }

       _setColCoeffs(col_iter_map[col_it], coeffs);

     }

     /// \e

     void eraseCol(const ColIt& col_it) {

       col_iter_map.set(col_it, VALID_CLASS, INVALID_CLASS);

       int cols[2];

       cols[1]=col_iter_map[col_it];

       _eraseCol(cols[1]);

       col_iter_map[col_it]=0; //glpk specifikus, de kell ez??

       ColIt it;

       for (col_iter_map.first(it, VALID_CLASS); 

 	   col_iter_map.valid(it); col_iter_map.next(it)) {

 	if (col_iter_map[it]>cols[1]) --col_iter_map[it];

       }

     }

     /// \e

     void eraseRow(const RowIt& row_it) {

       row_iter_map.set(row_it, VALID_CLASS, INVALID_CLASS);

       int rows[2];

       rows[1]=row_iter_map[row_it];

       _eraseRow(rows[1]);

       row_iter_map[row_it]=0; //glpk specifikus, de kell ez??

       RowIt it;

       for (row_iter_map.first(it, VALID_CLASS); 

 	   row_iter_map.valid(it); row_iter_map.next(it)) {

 	if (row_iter_map[it]>rows[1]) --row_iter_map[it];

       }

     }

     /// \e

     void setColBounds(const ColIt& col_it, Bound bound, 

 		      _Value lo, _Value up) {

       _setColBounds(col_iter_map[col_it], bound, lo, up);

     }

     /// \e

     void setRowBounds(const RowIt& row_it, Bound bound, 

 		      _Value lo, _Value up) {

       _setRowBounds(row_iter_map[row_it], bound, lo, up);

     }

     /// \e

     void setObjCoef(const ColIt& col_it, _Value obj_coef) {

       _setObjCoef(col_iter_map[col_it], obj_coef);

     }

     /// \e

     _Value getObjCoef(const ColIt& col_it) {

       return _getObjCoef(col_iter_map[col_it]);

     }

     //SOLVER FUNCTIONS

     /// \e

     virtual void solveSimplex() = 0;

     /// \e

     virtual void solvePrimalSimplex() = 0;

     /// \e

     virtual void solveDualSimplex() = 0;

     /// \e

     //HIGH LEVEL, SOLUTION RETRIEVING FUNCTIONS

   public:

     _Value getPrimal(const ColIt& col_it) {

       return _getPrimal(col_iter_map[col_it]);

     }

     /// \e

     virtual _Value getObjVal() = 0;

     //OTHER FUNCTIONS

     /// \e

     virtual int rowNum() const = 0;

     /// \e

     virtual int colNum() const = 0;

     /// \e

     virtual int warmUp() = 0;

     /// \e

     virtual void printWarmUpStatus(int i) = 0;

     /// \e

     virtual int getPrimalStatus() = 0;

     /// \e

     virtual void printPrimalStatus(int i) = 0;

     /// \e

     virtual int getDualStatus() = 0;

     /// \e

     virtual void printDualStatus(int i) = 0;

     /// Returns the status of the slack variable assigned to row \c row_it.

     virtual int getRowStat(const RowIt& row_it) = 0;

     /// \e

     virtual void printRowStatus(int i) = 0;

     /// Returns the status of the variable assigned to column \c col_it.

     virtual int getColStat(const ColIt& col_it) = 0;

     /// \e

     virtual void printColStatus(int i) = 0;

   };

   /// \brief Wrappers for LP solvers

   /// 

   /// This class implements a lemon wrapper for glpk.

   /// Later other LP-solvers will be wrapped into lemon.

   /// The aim of this class is to give a general surface to different 

   /// solvers, i.e. it makes possible to write algorithms using LP's, 

   /// in which the solver can be changed to an other one easily.

   class LPGLPK : public LPSolverBase<double> {

   public:

     typedef LPSolverBase<double> Parent;

   public:

     /// \e

     LPX* lp;

   public:

     /// \e

     LPGLPK() : Parent(), 

 			lp(lpx_create_prob()) {

       lpx_set_int_parm(lp, LPX_K_DUAL, 1);

     }

     /// \e

     ~LPGLPK() {

       lpx_delete_prob(lp);

     }

     //MATRIX INDEPEDENT MANIPULATING FUNCTIONS

     /// \e

     void setMinimize() { 

       lpx_set_obj_dir(lp, LPX_MIN);

     }

     /// \e

     void setMaximize() { 

       lpx_set_obj_dir(lp, LPX_MAX);

     }

     //LOW LEVEL INTERFACE, MATRIX MANIPULATING FUNCTIONS

   protected:

     /// \e

     int _addCol() { 

       return lpx_add_cols(lp, 1);

     }

     /// \e

     int _addRow() { 

       return lpx_add_rows(lp, 1);

     }

     /// \e

     virtual void _setRowCoeffs(int i, 

 			       std::vector<std::pair<int, double> > coeffs) {

       int mem_length=1+colNum();

       int* indices = new int[mem_length];

       double* doubles = new double[mem_length];

       int length=0;

       for (std::vector<std::pair<int, double> >::

 	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {

 	++length;

 	indices[length]=it->first;

 	doubles[length]=it->second;

 // 	std::cout << "  " << indices[length] << " " 

 // 		  << doubles[length] << std::endl;

       }

 //      std::cout << i << " " << length << std::endl;

       lpx_set_mat_row(lp, i, length, indices, doubles);

       delete [] indices;

       delete [] doubles;

     }

     /// \e

     virtual void _setColCoeffs(int i, 

 			       std::vector<std::pair<int, double> > coeffs) {

       int mem_length=1+rowNum();

       int* indices = new int[mem_length];

       double* doubles = new double[mem_length];

       int length=0;

       for (std::vector<std::pair<int, double> >::

 	     const_iterator it=coeffs.begin(); it!=coeffs.end(); ++it) {

 	++length;

 	indices[length]=it->first;

 	doubles[length]=it->second;

       }

       lpx_set_mat_col(lp, i, length, indices, doubles);

       delete [] indices;

       delete [] doubles;

     }

     /// \e

     virtual void _eraseCol(int i) {

       int cols[2];

       cols[1]=i;

       lpx_del_cols(lp, 1, cols);

     }

     virtual void _eraseRow(int i) {

       int rows[2];

       rows[1]=i;

       lpx_del_rows(lp, 1, rows);

     }

     virtual void _setColBounds(int i, Bound bound, 

 			       double lo, double up) {

       switch (bound) {

       case FREE:

 	lpx_set_col_bnds(lp, i, LPX_FR, lo, up);

 	break;

       case LOWER:

 	lpx_set_col_bnds(lp, i, LPX_LO, lo, up);

 	break;

       case UPPER:

 	lpx_set_col_bnds(lp, i, LPX_UP, lo, up);

 	break;

       case DOUBLE:

 	lpx_set_col_bnds(lp, i, LPX_DB, lo, up);

 	break;

       case FIXED:

 	lpx_set_col_bnds(lp, i, LPX_FX, lo, up);

 	break;

       }

     } 

     virtual void _setRowBounds(int i, Bound bound, 

 			       double lo, double up) {

       switch (bound) {

       case FREE:

 	lpx_set_row_bnds(lp, i, LPX_FR, lo, up);

 	break;

       case LOWER:

 	lpx_set_row_bnds(lp, i, LPX_LO, lo, up);

 	break;

       case UPPER:

 	lpx_set_row_bnds(lp, i, LPX_UP, lo, up);

 	break;

       case DOUBLE:

 	lpx_set_row_bnds(lp, i, LPX_DB, lo, up);

 	break;

       case FIXED:

 	lpx_set_row_bnds(lp, i, LPX_FX, lo, up);

 	break;

       }

     } 

   protected:

     /// \e

     virtual double _getObjCoef(int i) { 

       return lpx_get_obj_coef(lp, i);

     }

     /// \e

     virtual void _setObjCoef(int i, double obj_coef) { 

       lpx_set_obj_coef(lp, i, obj_coef);

     }

   public:

     /// \e

     void solveSimplex() { lpx_simplex(lp); }

     /// \e

     void solvePrimalSimplex() { lpx_simplex(lp); }

     /// \e

     void solveDualSimplex() { lpx_simplex(lp); }

     /// \e

   protected:

     virtual double _getPrimal(int i) {

       return lpx_get_col_prim(lp, i);

     }

   public:

     /// \e

     double getObjVal() { return lpx_get_obj_val(lp); }

     /// \e

     int rowNum() const { return lpx_get_num_rows(lp); }

     /// \e

     int colNum() const { return lpx_get_num_cols(lp); }

     /// \e

     int warmUp() { return lpx_warm_up(lp); }

     /// \e

     void printWarmUpStatus(int i) {

       switch (i) {

       case LPX_E_OK: cout << "LPX_E_OK" << endl; break;

       case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;	

       case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;

       case LPX_E_SING: cout << "LPX_E_SING" << endl; break;

       }

     }

     /// \e

     int getPrimalStatus() { return lpx_get_prim_stat(lp); }

     /// \e

     void printPrimalStatus(int i) {

       switch (i) {

       case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;

       case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;	

       case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;

       case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;

       }

     }

     /// \e

     int getDualStatus() { return lpx_get_dual_stat(lp); }

     /// \e

     void printDualStatus(int i) {

       switch (i) {

       case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;

       case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;	

       case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;

       case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;

       }

     }

     /// Returns the status of the slack variable assigned to row \c row_it.

     int getRowStat(const RowIt& row_it) { 

       return lpx_get_row_stat(lp, row_iter_map[row_it]); 

     }

     /// \e

     void printRowStatus(int i) {

       switch (i) {

       case LPX_BS: cout << "LPX_BS" << endl; break;

       case LPX_NL: cout << "LPX_NL" << endl; break;	

       case LPX_NU: cout << "LPX_NU" << endl; break;

       case LPX_NF: cout << "LPX_NF" << endl; break;

       case LPX_NS: cout << "LPX_NS" << endl; break;

       }

     }

     /// Returns the status of the variable assigned to column \c col_it.

     int getColStat(const ColIt& col_it) { 

       return lpx_get_col_stat(lp, col_iter_map[col_it]); 

     }

     /// \e

     void printColStatus(int i) {

       switch (i) {

       case LPX_BS: cout << "LPX_BS" << endl; break;

       case LPX_NL: cout << "LPX_NL" << endl; break;	

       case LPX_NU: cout << "LPX_NU" << endl; break;

       case LPX_NF: cout << "LPX_NF" << endl; break;

       case LPX_NS: cout << "LPX_NS" << endl; break;

       }

     }

   };

   /// @}

 } //namespace lemon

 #endif //LEMON_LP_SOLVER_WRAPPER_H