// -*- 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