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

// #include <stdio>

//#include <stdlib>

extern "C" {

#include "glpk.h"

}

#include <iostream>

#include <vector>

#include <string>

#include <list>

#include <memory>

#include <utility>

//#include <lemon/list_graph.h>

#include <lemon/invalid.h>

#include <expression.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) { }

      friend bool operator<(const ClassIt& x, const ClassIt& y);

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

				      const ClassIt& it);

    };

    friend bool operator<(const ClassIt& x, const ClassIt& y) {

      return (x.i < y.i);

    }

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

				    const ClassIt& it) {

      os << it.i;

      return os;

    }

    /// 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; }

  };

  /*! \e

    \todo kellenene uj iterable structure bele, mert ez nem az igazi

    \todo A[x,y]-t cserel. Jobboldal, baloldal csere.

    \todo LEKERDEZESEK!!!

    \todo DOKSI!!!! Doxygen group!!!

    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.

    \nosubgrouping

  */

  template <typename _Value>

  class LPSolverBase {

    /*! @name Uncategorized functions and types (public members)

    */

    //@{

  public:

    //UNCATEGORIZED

    /// \e

    typedef _Value Value;

    /// \e

    typedef IterablePartition<int>::ClassIt Row;

    /// \e

    typedef IterablePartition<int>::ClassIt Col;

  public:

    /// \e

    IterablePartition<int> row_iter_map;

    /// \e

    IterablePartition<int> col_iter_map;

    /// \e

    std::vector<Row> int_row_map;

    /// \e

    std::vector<Col> int_col_map;

    /// \e

    const int VALID_CLASS;

    /// \e

    const int INVALID_CLASS;

    /// \e 

    static const _Value INF;

  public:

    /// \e

    LPSolverBase() : row_iter_map(2), 

		     col_iter_map(2), 

		     VALID_CLASS(0), INVALID_CLASS(1) { }

    /// \e

    virtual ~LPSolverBase() { }

    //@}

    /*! @name Medium level interface (public members)

      These functions appear in the low level and also in the high level 

      interfaces thus these each of these functions have to be implemented 

      only once in the different interfaces.

      This means that these functions have to be reimplemented for all of the 

      different lp solvers. These are basic functions, and have the same 

      parameter lists in the low and high level interfaces. 

    */

    //@{

  public:

    //UNCATEGORIZED FUNCTIONS

    /// \e

    virtual void setMinimize() = 0;

    /// \e

    virtual void setMaximize() = 0;

    //SOLVER FUNCTIONS

    /// \e

    virtual void solveSimplex() = 0;

    /// \e

    virtual void solvePrimalSimplex() = 0;

    /// \e

    virtual void solveDualSimplex() = 0;

    /// \e

    //SOLUTION RETRIEVING

    /// \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.

    virtual int getRowStat(const Row& row) = 0;

    /// \e

    virtual void printRowStatus(int i) = 0;

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

    virtual int getColStat(const Col& col) = 0;

    /// \e

    virtual void printColStatus(int i) = 0;

    //@}

    /*! @name Low level interface (protected members)

      Problem manipulating functions in the low level interface

    */

    //@{

  protected:

    //MATRIX MANIPULATING FUNCTIONS

    /// \e

    virtual int _addCol() = 0;

    /// \e

    virtual int _addRow() = 0;

    /// \e

    virtual void _eraseCol(int i) = 0;

    /// \e

    virtual void _eraseRow(int i) = 0;

    /// \e

    virtual void _setRowCoeffs(int i, 

			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;

    /// \e

    /// This routine modifies \c coeffs only by the \c push_back method.

    virtual void _getRowCoeffs(int i, 

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

    virtual void _setColCoeffs(int i, 

			       const std::vector<std::pair<int, _Value> >& coeffs) = 0;

    /// \e

    /// This routine modifies \c coeffs only by the \c push_back method.

    virtual void _getColCoeffs(int i, 

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

  public:

    /// \e

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

  protected:

    /// \e

    /// The lower bound of a variable (column) have to be given by an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or -INF.

    virtual void _setColLowerBound(int i, _Value value) = 0;

    /// \e

    /// The lower bound of a variable (column) is an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or -INF.

    virtual _Value _getColLowerBound(int i) = 0;

    /// \e

    /// The upper bound of a variable (column) have to be given by an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or INF.

    virtual void _setColUpperBound(int i, _Value value) = 0;

    /// \e

    /// The upper bound of a variable (column) is an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or INF.

    virtual _Value _getColUpperBound(int i) = 0;

    /// \e

    /// The lower bound of a linear expression (row) have to be given by an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or -INF.

    virtual void _setRowLowerBound(int i, _Value value) = 0;

    /// \e

    /// The lower bound of a linear expression (row) is an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or -INF.

    virtual _Value _getRowLowerBound(int i) = 0;

    /// \e

    /// The upper bound of a linear expression (row) have to be given by an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or INF.

    virtual void _setRowUpperBound(int i, _Value value) = 0;

    /// \e

    /// The upper bound of a linear expression (row) is an 

    /// extended number of type _Value, i.e. a finite number of type 

    /// _Value or INF.

    virtual _Value _getRowUpperBound(int i) = 0;

    /// \e

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

    /// \e

    virtual _Value _getObjCoef(int i) = 0;

    //SOLUTION RETRIEVING

    /// \e

    virtual _Value _getPrimal(int i) = 0;

    //@}

    /*! @name High level interface (public members)

      Problem manipulating functions in the high level interface

    */

    //@{

  public:

    //MATRIX MANIPULATING FUNCTIONS

    /// \e

    Col addCol() {

      int i=_addCol();  

      Col col;

      col_iter_map.first(col, INVALID_CLASS);

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

	col_iter_map.set(col, INVALID_CLASS, VALID_CLASS);

	col_iter_map[col]=i;

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

	col=col_iter_map.push_back(i, VALID_CLASS);

      }

      int_col_map.push_back(col);

      return col;

    }

    /// \e

    Row addRow() {

      int i=_addRow();

      Row row;

      row_iter_map.first(row, INVALID_CLASS);

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

	row_iter_map.set(row, INVALID_CLASS, VALID_CLASS);

	row_iter_map[row]=i;

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

	row=row_iter_map.push_back(i, VALID_CLASS);

      }

      int_row_map.push_back(row);

      return row;

    }

    /// \e

    void eraseCol(const Col& col) {

      col_iter_map.set(col, VALID_CLASS, INVALID_CLASS);

      int cols[2];

      cols[1]=col_iter_map[col];

      _eraseCol(cols[1]);

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

      Col 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];

      }

      int_col_map.erase(int_col_map.begin()+cols[1]);

    }

    /// \e

    void eraseRow(const Row& row) {

      row_iter_map.set(row, VALID_CLASS, INVALID_CLASS);

      int rows[2];

      rows[1]=row_iter_map[row];

      _eraseRow(rows[1]);

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

      Row 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];

      }

      int_row_map.erase(int_row_map.begin()+rows[1]);

    }

    /// \e

    void setColLowerBound(Col col, _Value lo) {

      _setColLowerBound(col_iter_map[col], lo);

    }

    /// \e

    _Value getColLowerBound(Col col) {

      return _getColLowerBound(col_iter_map[col]);

    }

    /// \e

    void setColUpperBound(Col col, _Value up) {

      _setColUpperBound(col_iter_map[col], up);

    }

    /// \e

    _Value getColUpperBound(Col col) {      

      return _getColUpperBound(col_iter_map[col]);

    }

    /// \e

    void setRowLowerBound(Row row, _Value lo) {

      _setRowLowerBound(row_iter_map[row], lo);

    }

    /// \e

    _Value getRowLowerBound(Row row) {

      return _getRowLowerBound(row_iter_map[row]);

    }

    /// \e

    void setRowUpperBound(Row row, _Value up) {

      _setRowUpperBound(row_iter_map[row], up);

    }

    /// \e

    _Value getRowUpperBound(Row row) {      

      return _getRowUpperBound(row_iter_map[row]);

    }

    /// \e

    void setObjCoef(const Col& col, _Value obj_coef) {

      _setObjCoef(col_iter_map[col], obj_coef);

    }

    /// \e

    _Value getObjCoef(const Col& col) {

      return _getObjCoef(col_iter_map[col]);

    }

    //SOLUTION RETRIEVING FUNCTIONS

    /// \e

    _Value getPrimal(const Col& col) {

      return _getPrimal(col_iter_map[col]);

    }    

    //@}

    /*! @name User friend interface

      Problem manipulating functions in the user friend interface

    */

    //@{

    //EXPRESSION TYPES

    /// \e

    typedef Expr<Col, _Value> Expression;

    /// \e

    typedef Expr<Row, _Value> DualExpression;

    /// \e

    typedef Constr<Col, _Value> Constraint;

    //MATRIX MANIPULATING FUNCTIONS

    /// \e

    void setRowCoeffs(Row row, const Expression& expr) {

      std::vector<std::pair<int, _Value> > row_coeffs;

      for(typename Expression::Data::const_iterator i=expr.data.begin(); 

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

	row_coeffs.push_back(std::make_pair

			     (col_iter_map[(*i).first], (*i).second));

      }

      _setRowCoeffs(row_iter_map[row], row_coeffs);

    }

    /// \e 

    void setRow(Row row, const Constraint& constr) {

      setRowCoeffs(row, constr.expr);

      setRowLowerBound(row, constr.lo);

      setRowUpperBound(row, constr.up);

    }

    /// \e 

    Row addRow(const Constraint& constr) {

      Row row=addRow();

      setRowCoeffs(row, constr.expr);

      setRowLowerBound(row, constr.lo);

      setRowUpperBound(row, constr.up);

      return row;

    }

    /// \e

    /// This routine modifies \c expr by only adding to it.

    void getRowCoeffs(Row row, Expression& expr) {

      std::vector<std::pair<int, _Value> > row_coeffs;

      _getRowCoeffs(row_iter_map[row], row_coeffs);

      for(typename std::vector<std::pair<int, _Value> >::const_iterator 

 	    i=row_coeffs.begin(); i!=row_coeffs.end(); ++i) {

 	expr+= (*i).second*int_col_map[(*i).first];

      }

    }

    /// \e

    void setColCoeffs(Col col, const DualExpression& expr) {

      std::vector<std::pair<int, _Value> > col_coeffs;

      for(typename DualExpression::Data::const_iterator i=expr.data.begin(); 

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

	col_coeffs.push_back(std::make_pair

			     (row_iter_map[(*i).first], (*i).second));

      }

      _setColCoeffs(col_iter_map[col], col_coeffs);

    }

    /// \e

    /// This routine modifies \c expr by only adding to it.

    void getColCoeffs(Col col, DualExpression& expr) {

      std::vector<std::pair<int, _Value> > col_coeffs;

      _getColCoeffs(col_iter_map[col], col_coeffs);

      for(typename std::vector<std::pair<int, _Value> >::const_iterator 

 	    i=col_coeffs.begin(); i!=col_coeffs.end(); ++i) {

 	expr+= (*i).second*int_row_map[(*i).first];

      }

    }

    /// \e

    /// \bug ez igy nem jo

    void setObjCoeffs(const Expression& expr) {

      for(typename Expression::Data::const_iterator i=expr.data.begin(); 

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

	setObjCoef((*i).first, (*i).second);

      }

    }

    /// \e

    /// This routine modifies \c expr by only adding to it.

    void getObjCoeffs(Expression& expr) {

      /// FIXME not yet implemented

    }

    //@}

  };

  template <typename _Value>

  const _Value LPSolverBase<_Value>::INF=std::numeric_limits<_Value>::infinity();

  /// \brief Wrapper for GLPK solver

  /// 

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

  class LPGLPK : public LPSolverBase<double> {

  public:

    typedef LPSolverBase<double> Parent;

  public:

    /// \e

    LPX* lp;

  public:

    /// \e

    LPGLPK() : Parent(), 

			lp(lpx_create_prob()) {

      int_row_map.push_back(Row());

      int_col_map.push_back(Col());

      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() { 

      int i=lpx_add_cols(lp, 1);

      _setColLowerBound(i, -INF);

      _setColUpperBound(i, INF);

      return i;

    }

    /// \e

    int _addRow() { 

      int i=lpx_add_rows(lp, 1);

      return i;

    }

    /// \e

    virtual void _setRowCoeffs(int i, 

			       const 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;

      }

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

      delete [] indices;

      delete [] doubles;

    }

    /// \e

    virtual void _getRowCoeffs(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=lpx_get_mat_row(lp, i, indices, doubles);

      for (int i=1; i<=length; ++i) {

	coeffs.push_back(std::make_pair(indices[i], doubles[i]));

      }

      delete [] indices;

      delete [] doubles;

    }

    /// \e

    virtual void _setColCoeffs(int i, 

			       const 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 _getColCoeffs(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=lpx_get_mat_col(lp, i, indices, doubles);

      for (int i=1; i<=length; ++i) {

	coeffs.push_back(std::make_pair(indices[i], doubles[i]));

      }

      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 _setColLowerBound(int i, double lo) {

      if (lo==INF) {

	//FIXME error

      }

      int b=lpx_get_col_type(lp, i);

      double up=lpx_get_col_ub(lp, i);	

      if (lo==-INF) {

	switch (b) {

	case LPX_FR:

	case LPX_LO:

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

	  break;

	case LPX_UP:

	  break;

	case LPX_DB:

	case LPX_FX:

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

	  break;

	default: ;

	  //FIXME error

	}

      } else {

	switch (b) {

	case LPX_FR:

	case LPX_LO:

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

	  break;

	case LPX_UP:	  

	case LPX_DB:

	case LPX_FX:

	  if (lo==up) 

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

	  else 

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

	  break;

	default: ;

	  //FIXME error

	}

      }

    }

    virtual double _getColLowerBound(int i) {

      int b=lpx_get_col_type(lp, i);

      switch (b) {

      case LPX_FR:

	return -INF;

      case LPX_LO:

	return lpx_get_col_lb(lp, i);

      case LPX_UP:

	return -INF;

      case LPX_DB:

      case LPX_FX:

	return lpx_get_col_lb(lp, i);

      default: ;

	//FIXME error

	return 0.0;

      }

    }

    virtual void _setColUpperBound(int i, double up) {

      if (up==-INF) {

	//FIXME error

      }

      int b=lpx_get_col_type(lp, i);

      double lo=lpx_get_col_lb(lp, i);

      if (up==INF) {

	switch (b) {

	case LPX_FR:

	case LPX_LO:

	  break;

	case LPX_UP:

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

	  break;

	case LPX_DB:

	case LPX_FX:

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

	  break;

	default: ;

	  //FIXME error

	}

      } else {

	switch (b) {

	case LPX_FR:

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

	case LPX_LO:

	  if (lo==up) 

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

	  else

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

	  break;

	case LPX_UP:

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

	  break;

	case LPX_DB:

	case LPX_FX:

	  if (lo==up) 

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

	  else 

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

	  break;

	default: ;

	  //FIXME error

	}

      }

    }

    virtual double _getColUpperBound(int i) {

      int b=lpx_get_col_type(lp, i);

      switch (b) {

      case LPX_FR:

      case LPX_LO:

	return INF;

      case LPX_UP:

      case LPX_DB:

      case LPX_FX:

	return lpx_get_col_ub(lp, i);

      default: ;

	//FIXME error

	return 0.0;

      }

    }

    virtual void _setRowLowerBound(int i, double lo) {

      if (lo==INF) {

	//FIXME error

      }

      int b=lpx_get_row_type(lp, i);

      double up=lpx_get_row_ub(lp, i);	

      if (lo==-INF) {

	switch (b) {

	case LPX_FR:

	case LPX_LO:

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

	  break;

	case LPX_UP:

	  break;

	case LPX_DB:

	case LPX_FX:

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

	  break;

	default: ;

	  //FIXME error

	}

      } else {

	switch (b) {

	case LPX_FR:

	case LPX_LO:

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

	  break;

	case LPX_UP:	  

	case LPX_DB:

	case LPX_FX:

	  if (lo==up) 

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

	  else 

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

	  break;

	default: ;

	  //FIXME error

	}

      }

    }

    virtual double _getRowLowerBound(int i) {

      int b=lpx_get_row_type(lp, i);

      switch (b) {

      case LPX_FR:

	return -INF;

      case LPX_LO:

	return lpx_get_row_lb(lp, i);

      case LPX_UP:

	return -INF;

      case LPX_DB:

      case LPX_FX:

	return lpx_get_row_lb(lp, i);

      default: ;

	//FIXME error

	return 0.0;

      }

    }

    virtual void _setRowUpperBound(int i, double up) {

      if (up==-INF) {

	//FIXME error

      }

      int b=lpx_get_row_type(lp, i);

      double lo=lpx_get_row_lb(lp, i);

      if (up==INF) {

	switch (b) {

	case LPX_FR:

	case LPX_LO:

	  break;

	case LPX_UP:

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

	  break;

	case LPX_DB:

	case LPX_FX:

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

	  break;

	default: ;

	  //FIXME error

	}

      } else {

	switch (b) {

	case LPX_FR:

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

	case LPX_LO:

	  if (lo==up) 

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

	  else

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

	  break;

	case LPX_UP:

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

	  break;

	case LPX_DB:

	case LPX_FX:

	  if (lo==up) 

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

	  else 

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

	  break;

	default: ;

	  //FIXME error

	}

      }

    }

    virtual double _getRowUpperBound(int i) {

      int b=lpx_get_row_type(lp, i);

      switch (b) {

      case LPX_FR:

      case LPX_LO:

	return INF;

      case LPX_UP:

      case LPX_DB:

      case LPX_FX:

	return lpx_get_row_ub(lp, i);

      default: ;

	//FIXME error

	return 0.0;

      }

    }

    /// \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.

    int getRowStat(const Row& row) { 

      return lpx_get_row_stat(lp, row_iter_map[row]); 

    }