source: lemon-0.x/src/work/marci/lp/lp_solver_wrapper_3.h @ 1097:c91e765266d7

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