source: lemon-0.x/src/work/marci/lp/lp_solver_base.h @ 1278:4abea330614d

// -*- c++ -*-
#ifndef LEMON_LP_SOLVER_BASE_H
#define LEMON_LP_SOLVER_BASE_H

///\ingroup misc
///\file

// #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/invalid.h>
#include <expression.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 Class;
    /// 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(Class 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(Class 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.
    Class 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(Class 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[](Class it) { return nodes[it.i].data; }
    /// Returns the data pointed by \c it.
    const T& operator[](Class it) const { return nodes[it.i].data; }
    ///.
    class Class {
      friend class IterablePartition;
    protected:
      int i;
    public:
      /// Default constructor.
      Class() { }
      /// This constructor constructs an iterator which points
      /// to the member of th container indexed by the integer _i.
      Class(const int& _i) : i(_i) { }
      /// Invalid constructor.
      Class(const Invalid&) : i(-1) { }
      friend bool operator<(const Class& x, const Class& y);
      friend std::ostream& operator<<(std::ostream& os, 
                                      const Class& it);
      bool operator==(const Class& node) const {return i == node.i;}
      bool operator!=(const Class& node) const {return i != node.i;}
    };
    friend bool operator<(const Class& x, const Class& y) {
      return (x.i < y.i);
    }
    friend std::ostream& operator<<(std::ostream& os, 
                                    const Class& it) {
      os << it.i;
      return os;
    }
    /// First member of class \c class_id.
    Class& first(Class& it, int class_id) const {
      it.i=tips[class_id].first;
      return it;
    }
    /// Next member.
    Class& next(Class& it) const {
      it.i=nodes[it.i].next;
      return it;
    }
    /// True iff the iterator is valid.
    bool valid(const Class& it) const { return it.i!=-1; }

    class ClassIt : public Class {
      const IterablePartition* iterable_partition;
    public:
      ClassIt() { }
      ClassIt(Invalid i) : Class(i) { }
      ClassIt(const IterablePartition& _iterable_partition, 
              const int& i) : iterable_partition(&_iterable_partition) {
        _iterable_partition.first(*this, i);
      }
      ClassIt(const IterablePartition& _iterable_partition, 
              const Class& _class) : 
        Class(_class), iterable_partition(&_iterable_partition) { }
      ClassIt& operator++() {
        iterable_partition->next(*this);
        return *this;
      }
    };

  };


  /*! \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 IterablePartition<int> Rows;
    /// \e
    typedef IterablePartition<int> Cols;
    /// \e
    typedef _Value Value;
    /// \e
    typedef Rows::Class Row;
    /// \e
    typedef Cols::Class 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;

    //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;
    /// \e
    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;
    /// \e
    virtual void _setCoeff(int col, int row, _Value value) = 0;
    /// \e
    virtual _Value _getCoeff(int col, int row) = 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 _setObjCoeff(int i, _Value obj_coef) = 0;
    /// \e
    virtual _Value _getObjCoeff(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 setCoeff(Col col, Row row, _Value value) {
      _setCoeff(col_iter_map[col], row_iter_map[row], value);
    }
    /// \e
    _Value getCoeff(Col col, Row row) {
      return _getCoeff(col_iter_map[col], row_iter_map[row], value);
    }
    /// \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 setObjCoeff(const Col& col, _Value obj_coef) {
      _setObjCoeff(col_iter_map[col], obj_coef);
    }
    /// \e
    _Value getObjCoeff(const Col& col) {
      return _getObjCoeff(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
    void setObjCoeffs(const Expression& expr) {
      // writing zero everywhere
      for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
        setObjCoeff(it, 0.0);
      // writing the data needed
      for(typename Expression::Data::const_iterator i=expr.data.begin(); 
          i!=expr.data.end(); ++i) {
        setObjCoeff((*i).first, (*i).second);
      }
    }
    /// \e
    /// This routine modifies \c expr by only adding to it.
    void getObjCoeffs(Expression& expr) {
      for(Cols::ClassIt it(col_iter_map, VALID_CLASS); it!=INVALID; ++it)
        expr+=getObjCoeff(it)*it;
    }
    //@}


    /*! @name MIP functions and types (public members)
    */
    //@{
  public:
    /// \e
    virtual void solveBandB() = 0;
    /// \e
    virtual void setLP() = 0;
    /// \e
    virtual void setMIP() = 0;
  protected:
   /// \e
    virtual void _setColCont(int i) = 0;
    /// \e
    virtual void _setColInt(int i) = 0;
    /// \e
    virtual _Value _getMIPPrimal(int i) = 0;
  public:
    /// \e
    void setColCont(Col col) {
      _setColCont(col_iter_map[col]);
    }
    /// \e
    void setColInt(Col col) {
      _setColInt(col_iter_map[col]);
    }
    /// \e
    _Value getMIPPrimal(Col col) {
      return _getMIPPrimal(col_iter_map[col]);
    }
    //@}
  };
  
  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);
    }
    void _setCoeff(int col, int row, double value) {
      /// FIXME not yet implemented
    }
    double _getCoeff(int col, int row) {
      /// FIXME not yet implemented
      return 0.0;
    }
    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);
          break;
        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);
          break;
        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 _getObjCoeff(int i) { 
      return lpx_get_obj_coef(lp, i);
    }
    /// \e
    virtual void _setObjCoeff(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); }
  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]); 
    }
    /// \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.
    int getColStat(const Col& col) { 
      return lpx_get_col_stat(lp, col_iter_map[col]); 
    }
    /// \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;
      }
    }

    // MIP
    /// \e
    void solveBandB() { lpx_integer(lp); }
    /// \e
    void setLP() { lpx_set_class(lp, LPX_LP); }
    /// \e
    void setMIP() { lpx_set_class(lp, LPX_MIP); }
  protected:
    /// \e
    void _setColCont(int i) { lpx_set_col_kind(lp, i, LPX_CV); }
    /// \e
    void _setColInt(int i) { lpx_set_col_kind(lp, i, LPX_IV); }
    /// \e
    double _getMIPPrimal(int i) { return lpx_mip_col_val(lp, i); }
  };
  
  /// @}

} //namespace lemon

#endif //LEMON_LP_SOLVER_BASE_H