source: lemon-0.x/src/work/athos/lp/lp_solver_base.h @ 1241:dadc9987c537

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

#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]);
    }
    //@}
  };

} //namespace lemon

#endif //LEMON_LP_SOLVER_BASE_H