# HG changeset patch
# User athos
# Date 1111571381 0
# Node ID 43a3d06e0ee07d8b4cb168924197905121d39469
# Parent  41caee260bd4b6a1c94036f451616f5e43f2486a


diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/expression.h
--- a/src/work/athos/lp/expression.h	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,197 +0,0 @@
-// -*- c++ -*-
-#ifndef LEMON_EXPRESSION_H
-#define LEMON_EXPRESSION_H
-
-#include <iostream>
-#include <map>
-#include <limits>
-
-namespace lemon {
-
-  /*! \brief Linear expression
-
-    \c Expr<_Col,_Value> implements a class of linear expressions with the 
-    operations of addition and multiplication with scalar. 
-
-    \author Marton Makai
-   */
-  template <typename _Col, typename _Value>
-  class Expr;
-
-  template <typename _Col, typename _Value>
-  class Expr {
-//  protected:
-  public:
-    typedef 
-    typename std::map<_Col, _Value> Data; 
-    Data data;
-  public:
-    void simplify() {
-      for (typename Data::iterator i=data.begin(); 
-	   i!=data.end(); ++i) {
-	if ((*i).second==0) data.erase(i);
-      }
-    }
-    Expr() { }
-    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;
-      }
-      simplify();
-      return *this;
-    }
-    Expr& operator+=(const Expr<_Col, _Value>& expr) {
-      for (typename Data::const_iterator j=expr.data.begin(); 
-	   j!=expr.data.end(); ++j) {
-	typename Data::iterator i=data.find((*j).first);
-	if (i==data.end()) {
-	  data.insert(std::make_pair((*j).first, (*j).second));
-	} else {
-	  (*i).second+=(*j).second;
-	}
-      }
-      simplify();
-      return *this;
-    }
-    Expr& operator-=(const Expr<_Col, _Value>& expr) {
-      for (typename Data::const_iterator j=expr.data.begin(); 
-	   j!=expr.data.end(); ++j) {
-	typename Data::iterator i=data.find((*j).first);
-	if (i==data.end()) {
-	  data.insert(std::make_pair((*j).first, -(*j).second));
-	} else {
-	  (*i).second+=-(*j).second;
-	}
-      }
-      simplify();
-      return *this;
-    }
-    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, _Col _col) {
-    Expr<_Col, _Value> tmp(_col);
-    tmp*=_value;
-    tmp.simplify();
-    return tmp;
-  }
-
-  template <typename _Col, typename _Value>
-  Expr<_Col, _Value> operator*(_Value _value, 
-			       const Expr<_Col, _Value>& expr) {
-    Expr<_Col, _Value> tmp(expr);
-    tmp*=_value;
-    tmp.simplify();
-    return tmp;
-  }
-
-  template <typename _Col, typename _Value>
-  Expr<_Col, _Value> operator+(const Expr<_Col, _Value>& expr1, 
-			       const Expr<_Col, _Value>& expr2) {
-    Expr<_Col, _Value> tmp(expr1);
-    tmp+=expr2;
-    tmp.simplify();
-    return tmp;
-  }
-
-  template <typename _Col, typename _Value>
-  Expr<_Col, _Value> operator-(const Expr<_Col, _Value>& expr1, 
-			       const Expr<_Col, _Value>& expr2) {
-    Expr<_Col, _Value> tmp(expr1);
-    tmp-=expr2;
-    tmp.simplify();
-    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;
-  }
-
-  template <typename _Col, typename _Value>
-  class LConstr {
-    //  protected:
-  public:
-    Expr<_Col, _Value> expr;
-    _Value lo;
-  public:
-    LConstr(const Expr<_Col, _Value>& _expr, _Value _lo) : 
-      expr(_expr), lo(_lo) { }
-  };
-  
-  template <typename _Col, typename _Value>
-  LConstr<_Col, _Value> 
-  operator<=(_Value lo, const Expr<_Col, _Value>& expr) {
-    return LConstr<_Col, _Value>(expr, lo);
-  }
-
-  template <typename _Col, typename _Value>
-  class UConstr {
-    //  protected:
-  public:
-    Expr<_Col, _Value> expr;
-    _Value up;
-  public:
-    UConstr(const Expr<_Col, _Value>& _expr, _Value _up) : 
-      expr(_expr), up(_up) { }
-  };
-
-  template <typename _Col, typename _Value>
-  UConstr<_Col, _Value> 
-  operator<=(const Expr<_Col, _Value>& expr, _Value up) {
-    return UConstr<_Col, _Value>(expr, up);
-  }
-
-  template <typename _Col, typename _Value>
-  class Constr {
-    //  protected:
-  public:
-    Expr<_Col, _Value> expr;
-    _Value lo, up;
-  public:
-    Constr(const Expr<_Col, _Value>& _expr, _Value _lo, _Value _up) : 
-      expr(_expr), lo(_lo), up(_up) { }
-    Constr(const LConstr<_Col, _Value>& _lconstr) : 
-      expr(_lconstr.expr), 
-      lo(_lconstr.lo), 
-      up(std::numeric_limits<_Value>::infinity()) { }
-    Constr(const UConstr<_Col, _Value>& _uconstr) : 
-      expr(_uconstr.expr), 
-      lo(-std::numeric_limits<_Value>::infinity()), 
-      up(_uconstr.up) { }
-  };
-
-  template <typename _Col, typename _Value>
-  Constr<_Col, _Value> 
-  operator<=(const LConstr<_Col, _Value>& lconstr, _Value up) {
-    return Constr<_Col, _Value>(lconstr.expr, lconstr.lo, up);
-  }
-
-  template <typename _Col, typename _Value>
-  Constr<_Col, _Value> 
-  operator<=(_Value lo, const UConstr<_Col, _Value>& uconstr) {
-    return Constr<_Col, _Value>(uconstr.expr, lo, uconstr.up);
-  }
-
-  template <typename _Col, typename _Value>
-  Constr<_Col, _Value> 
-  operator==(const Expr<_Col, _Value>& expr, _Value value) {
-    return Constr<_Col, _Value>(expr, value, value);
-  }
-  
-} //namespace lemon
-
-#endif //LEMON_EXPRESSION_H
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/expression_test.cc
--- a/src/work/athos/lp/expression_test.cc	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,23 +0,0 @@
-#include <expression.h>
-#include <iostream>
-#include <string>
-
-using std::cout;
-using std::endl;
-using std::string;
-using namespace lemon;
-
-int main() {
-  Expr<string, double> b;
-  cout << b << endl;
-  Expr<string, double> c("f");
-  cout << c << endl;
-  Expr<string, double> d=8.0*string("g");
-  cout << d << endl;
-  c*=5;
-  cout << c << endl;
-  Expr<string, double> e=c;
-  e+=8.9*9.0*string("l");
-  cout << e << endl;
-  cout << c+d << endl;
-}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/lp_solver_base.h
--- a/src/work/athos/lp/lp_solver_base.h	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,639 +0,0 @@
-// -*- 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
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/lp_solver_glpk.h
--- a/src/work/athos/lp/lp_solver_glpk.h	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,545 +0,0 @@
-// -*- c++ -*-
-#ifndef LEMON_LP_SOLVER_GLPK_H
-#define LEMON_LP_SOLVER_GLPK_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 <lp_solver_base.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 {
-  
-  
-  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 Of course this is not efficient at all, but GLPK knows not more.
-      int change_this;
-      bool get_set_row;
-      //The only thing we can do is optimize on whether working with a row 
-      //or a coloumn
-      int row_num = rowNum();
-      int col_num = colNum();
-      if (col_num<row_num){
-	//There are more rows than coloumns
-	get_set_row=true;
-	int mem_length=1+row_num;
-	int* indices = new int[mem_length];
-	double* doubles = new double[mem_length];
-	int length=lpx_get_mat_col(lp, i, indices, doubles);
-      }else{
-	get_set_row=false;
-	int mem_length=1+col_num;
-	int* indices = new int[mem_length];
-	double* doubles = new double[mem_length];
-	int length=lpx_get_mat_row(lp, i, indices, doubles);
-      }
-      //Itten      
-int* indices = new int[mem_length];
-      double* doubles = new double[mem_length];
-      int length=lpx_get_mat_col(lp, i, indices, doubles);
- 
-      delete [] indices;
-      delete [] doubles;
-
-    }
-    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);
-	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 _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_GLPK_H
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/lp_solver_wrapper.h
--- a/src/work/athos/lp/lp_solver_wrapper.h	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,431 +0,0 @@
-// -*- 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 <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; }
-  };
-  
-  /// \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 LPSolverWrapper {
-  public:
-
-//   class Row {
-//   protected:
-//     int i;
-//   public:
-//     Row() { }
-//     Row(const Invalid&) : i(0) { }
-//     Row(const int& _i) : i(_i) { }
-//     operator int() const { return i; }
-//   };
-//   class RowIt : public Row {
-//   public:
-//     RowIt(const Row& row) : Row(row) { }
-//   };
-
-//   class Col {
-//   protected:
-//     int i;
-//   public:
-//     Col() { }
-//     Col(const Invalid&) : i(0) { }
-//     Col(const int& _i) : i(_i) { }
-//     operator int() const { return i; }
-//   };
-//   class ColIt : public Col {
-//     ColIt(const Col& col) : Col(col) { }
-//   };
-
-  public:
-    ///.
-    LPX* lp;
-    ///.
-    typedef IterablePartition<int>::ClassIt RowIt;
-    ///.
-    IterablePartition<int> row_iter_map;
-    ///.
-    typedef IterablePartition<int>::ClassIt ColIt;
-    ///.
-    IterablePartition<int> col_iter_map;
-    //std::vector<int> row_id_to_lp_row_id;
-    //std::vector<int> col_id_to_lp_col_id;
-    ///.
-    const int VALID_ID;
-    ///.
-    const int INVALID_ID;
-
-  public:
-    ///.
-    LPSolverWrapper() : lp(lpx_create_prob()), 
-			row_iter_map(2), 
-			col_iter_map(2), 
-			//row_id_to_lp_row_id(), col_id_to_lp_col_id(), 
-			VALID_ID(0), INVALID_ID(1) {
-      lpx_set_int_parm(lp, LPX_K_DUAL, 1);
-    }
-    ///.
-    ~LPSolverWrapper() {
-      lpx_delete_prob(lp);
-    }
-    ///.
-    void setMinimize() { 
-      lpx_set_obj_dir(lp, LPX_MIN);
-    }
-    ///.
-    void setMaximize() { 
-      lpx_set_obj_dir(lp, LPX_MAX);
-    }
-    ///.
-    ColIt addCol() {
-      int i=lpx_add_cols(lp, 1);  
-      ColIt col_it;
-      col_iter_map.first(col_it, INVALID_ID);
-      if (col_iter_map.valid(col_it)) { //van hasznalhato hely
-	col_iter_map.set(col_it, INVALID_ID, VALID_ID);
-	col_iter_map[col_it]=i;
-	//col_id_to_lp_col_id[col_iter_map[col_it]]=i;
-      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
-	//col_id_to_lp_col_id.push_back(i);
-	//int j=col_id_to_lp_col_id.size()-1;
-	col_it=col_iter_map.push_back(i, VALID_ID);
-      }
-//    edge_index_map.set(e, i);
-//    lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
-//    lpx_set_obj_coef(lp, i, cost[e]);    
-      return col_it;
-    }
-    ///.
-    RowIt addRow() {
-      int i=lpx_add_rows(lp, 1);  
-      RowIt row_it;
-      row_iter_map.first(row_it, INVALID_ID);
-      if (row_iter_map.valid(row_it)) { //van hasznalhato hely
-	row_iter_map.set(row_it, INVALID_ID, VALID_ID);
-	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_ID);
-      }
-      return row_it;
-    }
-    //pair<RowIt, double>-bol kell megadni egy std range-et
-    ///.
-    template <typename Begin, typename End>
-    void setColCoeffs(const ColIt& col_it, 
-		      Begin begin, End end) {
-      int mem_length=1+lpx_get_num_rows(lp);
-      int* indices = new int[mem_length];
-      double* doubles = new double[mem_length];
-      int length=0;
-      for ( ; begin!=end; ++begin) {
-	++length;
-	indices[length]=row_iter_map[begin->first];
-	doubles[length]=begin->second;
-      }
-      lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
-      delete [] indices;
-      delete [] doubles;
-    }
-    //pair<ColIt, double>-bol kell megadni egy std range-et
-    ///.
-    template <typename Begin, typename End>
-    void setRowCoeffs(const RowIt& row_it, 
-		      Begin begin, End end) {
-      int mem_length=1+lpx_get_num_cols(lp);
-      int* indices = new int[mem_length];
-      double* doubles = new double[mem_length];
-      int length=0;
-      for ( ; begin!=end; ++begin) {
-	++length;
-	indices[length]=col_iter_map[begin->first];
-	doubles[length]=begin->second;
-      }
-      lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
-      delete [] indices;
-      delete [] doubles;
-    }
-    ///.
-    void eraseCol(const ColIt& col_it) {
-      col_iter_map.set(col_it, VALID_ID, INVALID_ID);
-      int cols[2];
-      cols[1]=col_iter_map[col_it];
-      lpx_del_cols(lp, 1, cols);
-      col_iter_map[col_it]=0; //glpk specifikus
-      ColIt it;
-      for (col_iter_map.first(it, VALID_ID); 
-	   col_iter_map.valid(it); col_iter_map.next(it)) {
-	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
-      }
-    }
-    ///.
-    void eraseRow(const RowIt& row_it) {
-      row_iter_map.set(row_it, VALID_ID, INVALID_ID);
-      int rows[2];
-      rows[1]=row_iter_map[row_it];
-      lpx_del_rows(lp, 1, rows);
-      row_iter_map[row_it]=0; //glpk specifikus
-      RowIt it;
-      for (row_iter_map.first(it, VALID_ID); 
-	   row_iter_map.valid(it); row_iter_map.next(it)) {
-	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
-      }
-    }
-    ///.
-    void setColBounds(const ColIt& col_it, int bound_type, 
-		      double lo, double up) {
-      lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
-    }
-    ///.
-    double getObjCoef(const ColIt& col_it) { 
-      return lpx_get_obj_coef(lp, col_iter_map[col_it]);
-    }
-    ///.
-    void setRowBounds(const RowIt& row_it, int bound_type, 
-		      double lo, double up) {
-      lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
-    }
-    ///.
-    void setObjCoef(const ColIt& col_it, double obj_coef) { 
-      lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
-    }
-    ///.
-    void solveSimplex() { lpx_simplex(lp); }
-    ///.
-    void solvePrimalSimplex() { lpx_simplex(lp); }
-    ///.
-    void solveDualSimplex() { lpx_simplex(lp); }
-    ///.
-    double getPrimal(const ColIt& col_it) {
-      return lpx_get_col_prim(lp, col_iter_map[col_it]);
-    }
-    ///.
-    double getObjVal() { return lpx_get_obj_val(lp); }
-    ///.
-    int rowNum() const { return lpx_get_num_rows(lp); }
-    ///.
-    int colNum() const { return lpx_get_num_cols(lp); }
-    ///.
-    int warmUp() { return lpx_warm_up(lp); }
-    ///.
-    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;
-      }
-    }
-    ///.
-    int getPrimalStatus() { return lpx_get_prim_stat(lp); }
-    ///.
-    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;
-      }
-    }
-    ///.
-    int getDualStatus() { return lpx_get_dual_stat(lp); }
-    ///.
-    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]); 
-    }
-    ///.
-    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]); 
-    }
-    ///.
-    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
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/magic_square.cc
--- a/src/work/athos/lp/magic_square.cc	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,99 +0,0 @@
-// -*- c++ -*-
-#include <iostream>
-#include <fstream>
-
-#include <lemon/time_measure.h>
-#include <lp_solver_glpk.h>
-
-using std::cout;
-using std::endl;
-using namespace lemon;
-
-/*
-  On an 1537Mhz PC, the run times with 
-  glpk are the following.
-  for n=3,4, some secondes
-  for n=5, 25 hours
- */
-
-int main(int, char **) {
-  const int n=3;
-  const double row_sum=(1.0+n*n)*n/2;
-  Timer ts;
-  ts.reset();
-  typedef LpGlpk LPSolver;
-  typedef LPSolver::Col Col;
-  LPSolver lp;
-  typedef std::map<std::pair<int, int>, Col> Coords;
-  Coords x;
-  // we create a new variable for each entry 
-  // of the magic square
-  for (int i=1; i<=n; ++i) {
-    for (int j=1; j<=n; ++j) { 
-      Col col=lp.addCol();
-      x[std::make_pair(i,j)]=col;
-      lp.setColLowerBound(col, 1.0);
-      lp.setColUpperBound(col, double(n*n));
-    }
-  }
-  LPSolver::Expression expr3, expr4;
-  for (int i=1; i<=n; ++i) {
-    LPSolver::Expression expr1, expr2;
-    for (int j=1; j<=n; ++j) {
-      expr1+=x[std::make_pair(i, j)];
-      expr2+=x[std::make_pair(j, i)];
-    }
-
-    // sum of rows and columns
-    lp.addRow(expr1==row_sum);
-    lp.addRow(expr2==row_sum);
-      cout <<"Expr1: "<<expr1<<endl;
-      cout <<"Expr2: "<<expr2<<endl;
-
-    expr3+=x[std::make_pair(i, i)];
-    expr4+=x[std::make_pair(i, (n+1)-i)];
-  }
-  cout <<"Expr3: "<<expr3<<endl;
-  cout <<"Expr4: "<<expr4<<endl;
-  // sum of the diagonal entries
-  lp.addRow(expr3==row_sum);
-  lp.addRow(expr4==row_sum);
-  lp.solveSimplex();
-  cout << "elapsed time: " << ts << endl;
-  for (int i=1; i<=n; ++i) {
-    for (int j=1; j<=n; ++j) { 
-      cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)]) 
-	   << endl;
-    }
-  }
-
-
-
-//   // we make new binary variables for each pair of 
-//   // entries of the square to achieve that 
-//   // the values of different entries are different
-//   lp.setMIP();
-//   for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
-//     Coords::const_iterator jt=it; ++jt;
-//     for(; jt!=x.end(); ++jt) {
-//       Col col1=(*it).second;
-//       Col col2=(*jt).second;
-//       Col col=lp.addCol();
-//       lp.setColLowerBound(col, 0.0);
-//       lp.setColUpperBound(col, 1.0);
-//       lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
-//       lp.setColInt(col);
-//     }
-//   }
-//   cout << "elapsed time: " << ts << endl;
-//   lp.solveSimplex();
-//   // let's solve the integer problem
-//   lp.solveBandB();
-//   cout << "elapsed time: " << ts << endl;
-//   for (int i=1; i<=n; ++i) {
-//     for (int j=1; j<=n; ++j) { 
-//       cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)]) 
-// 	   << endl;
-//     }
-//   }
-}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/makefile
--- a/src/work/athos/lp/makefile	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,73 +0,0 @@
-#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
-#GLPKROOT = /usr/local/glpk-4.4
-INCLUDEDIRS ?= -I/home/marci/boost -I/usr/local/cplex/cplex75/include -I../../{marci,alpar,klao,akos,athos} -I. -I../../.. -I../.. -I..# -I$(GLPKROOT)/include
-#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
-CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
-LDFLAGS  =  -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
-
-BINARIES = magic_square max_flow_expression #expression_test max_flow_by_lp# sample sample2 sample11 sample15
-
-#include ../makefile
-
-# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
-# ismert rendszeren :-)  (Misi)
-ifdef GCCVER
-CXX := g++-$(GCCVER)
-else
-CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
-endif
-
-ifdef DEBUG
-CXXFLAGS += -DDEBUG
-endif
-
-CC := $(CXX)
-
-
-all: $(BINARIES)
-
-################
-# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
-# all elo.
-# Kiveve ha siman file.cc -> file  esetrol van szo, amikor is nem kell
-# irni semmit.
-
-#proba: proba.o seged.o
-
-################
-
-
-# .depend dep depend:
-# 	-$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
-
-#makefile: .depend
-#sinclude .depend
-#moert nem megy az eredeti /usr/bin/ld-vel?
-
-# %: %.o
-# 	$(CXX) -o $@ $< $(LDFLAGS)
-
-# %.o: %.cc
-# 	$(CXX) $(CXXFLAGS) -c $<
-
-%: %.cc
-	$(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
-
-sample11prof: sample11prof.o
-	 $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
-sample11prof.o: sample11.cc
-	$(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
-
-# sample.o: sample.cc
-# 	$(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
-
-# sample2: sample2.o
-# 	$(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
-# sample2.o: sample2.cc
-# 	$(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
-
-
-clean:
-	$(RM) *.o $(BINARIES) .depend
-
-.PHONY: all clean dep depend
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/max_flow_by_lp.cc
--- a/src/work/athos/lp/max_flow_by_lp.cc	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,186 +0,0 @@
-// -*- c++ -*-
-#include <iostream>
-#include <fstream>
-
-#include <lemon/smart_graph.h>
-#include <lemon/list_graph.h>
-#include <lemon/dimacs.h>
-#include <lemon/time_measure.h>
-//#include <graph_wrapper.h>
-#include <lemon/preflow.h>
-#include <augmenting_flow.h>
-//#include <preflow_res.h>
-//#include <lp_solver_wrapper_2.h>
-#include <min_cost_gen_flow.h>
-
-// Use a DIMACS max flow file as stdin.
-// max_flow_demo < dimacs_max_flow_file
-
-using namespace lemon;
-
-int main(int, char **) {
-
-  typedef ListGraph MutableGraph;
-  typedef ListGraph Graph;
-  typedef Graph::Node Node;
-  typedef Graph::Edge Edge;
-  typedef Graph::EdgeIt EdgeIt;
-
-  Graph g;
-
-  Node s, t;
-  Graph::EdgeMap<int> cap(g);
-  //readDimacsMaxFlow(std::cin, g, s, t, cap);
-  readDimacs(std::cin, g, cap, s, t);
-  Timer ts;
-  Graph::EdgeMap<int> flow(g); //0 flow
-  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
-    max_flow_test(g, s, t, cap, flow);
-  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
-    augmenting_flow_test(g, s, t, cap, flow);
-  
-  Graph::NodeMap<bool> cut(g);
-
-  {
-    std::cout << "preflow ..." << std::endl;
-    ts.reset();
-    max_flow_test.run();
-    std::cout << "elapsed time: " << ts << std::endl;
-    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
-    max_flow_test.minCut(cut);
-
-    for (EdgeIt e(g); e!=INVALID; ++e) {
-      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
-	std::cout << "Slackness does not hold!" << std::endl;
-      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
-	std::cout << "Slackness does not hold!" << std::endl;
-    }
-  }
-
-//   {
-//     std::cout << "preflow ..." << std::endl;
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
-//     ts.reset();
-//     max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
-//     std::cout << "elapsed time: " << ts << std::endl;
-//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
-
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
-//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
-// 	std::cout << "Slackness does not hold!" << std::endl;
-//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
-// 	std::cout << "Slackness does not hold!" << std::endl;
-//     }
-//   }
-
-//   {
-//     std::cout << "wrapped preflow ..." << std::endl;
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
-//     ts.reset();
-//     pre_flow_res.run();
-//     std::cout << "elapsed time: " << ts << std::endl;
-//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
-//   }
-
-  {
-    std::cout << "physical blocking flow augmentation ..." << std::endl;
-    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
-    ts.reset();
-    int i=0;
-    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
-    std::cout << "elapsed time: " << ts << std::endl;
-    std::cout << "number of augmentation phases: " << i << std::endl; 
-    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
-
-    for (EdgeIt e(g); e!=INVALID; ++e) {
-      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
-	std::cout << "Slackness does not hold!" << std::endl;
-      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
-	std::cout << "Slackness does not hold!" << std::endl;
-    }
-  }
-
-//   {
-//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
-//     ts.reset();
-//     int i=0;
-//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
-//     std::cout << "elapsed time: " << ts << std::endl;
-//     std::cout << "number of augmentation phases: " << i << std::endl; 
-//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
-//   }
-
-  {
-    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
-    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
-    ts.reset();
-    int i=0;
-    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
-    std::cout << "elapsed time: " << ts << std::endl;
-    std::cout << "number of augmentation phases: " << i << std::endl; 
-    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
-
-    for (EdgeIt e(g); e!=INVALID; ++e) {
-      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
-	std::cout << "Slackness does not hold!" << std::endl;
-      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
-	std::cout << "Slackness does not hold!" << std::endl;
-    }
-  }
-
-//   {
-//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
-//     ts.reset();
-//     int i=0;
-//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
-//     std::cout << "elapsed time: " << ts << std::endl;
-//     std::cout << "number of augmentation phases: " << i << std::endl; 
-//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
-
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
-//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
-// 	std::cout << "Slackness does not hold!" << std::endl;
-//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
-// 	std::cout << "Slackness does not hold!" << std::endl;
-//     }
-//   }
-
-//   {
-//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
-//     ts.reset();
-//     int i=0;
-//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
-//     std::cout << "elapsed time: " << ts << std::endl;
-//     std::cout << "number of augmentation phases: " << i << std::endl; 
-//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
-
-//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
-//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
-// 	std::cout << "Slackness does not hold!" << std::endl;
-//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
-// 	std::cout << "Slackness does not hold!" << std::endl;
-//     }
-//   }
-
-  ts.reset();
-
-  Edge e=g.addEdge(t, s);
-  Graph::EdgeMap<int> cost(g, 0);
-  cost.set(e, -1);
-  cap.set(e, 10000);
-  typedef ConstMap<Node, int> Excess;
-  Excess excess(0);
-  typedef ConstMap<Edge, int> LCap;
-  LCap lcap(0);
-
-  MinCostGenFlow<Graph, int, Excess, LCap> 
-    min_cost(g, excess, lcap, cap, flow, cost); 
-  min_cost.feasible();
-  min_cost.runByLP();
-
-  std::cout << "elapsed time: " << ts << std::endl;
-  std::cout << "flow value: "<< flow[e] << std::endl;
-}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/max_flow_expression.cc
--- a/src/work/athos/lp/max_flow_expression.cc	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,109 +0,0 @@
-// -*- c++ -*-
-#include <iostream>
-#include <fstream>
-
-#include <lemon/graph_utils.h>
-#include <lemon/smart_graph.h>
-#include <lemon/list_graph.h>
-#include <lemon/dimacs.h>
-#include <lemon/time_measure.h>
-#include <lp_solver_glpk.h>
-
-using std::cout;
-using std::endl;
-using namespace lemon;
-
-template<typename Edge, typename EdgeIndexMap> 
-class PrimalMap {
-protected:
-  LpGlpk* lp;
-  EdgeIndexMap* edge_index_map;
-public:
-  PrimalMap(LpGlpk& _lp, EdgeIndexMap& _edge_index_map) : 
-    lp(&_lp), edge_index_map(&_edge_index_map) { }
-  double operator[](Edge e) const { 
-    return lp->getPrimal((*edge_index_map)[e]);
-  }
-};
-
-// Use a DIMACS max flow file as stdin.
-// max_flow_expresion < dimacs_max_flow_file
-
-int main(int, char **) {
-
-  typedef ListGraph Graph;
-  typedef Graph::Node Node;
-  typedef Graph::Edge Edge;
-  typedef Graph::EdgeIt EdgeIt;
-
-  Graph g;
-  
-  Node s, t;
-  Graph::EdgeMap<int> cap(g);
-  readDimacs(std::cin, g, cap, s, t);
-  Timer ts;
-  
-  typedef LpGlpk LPSolver;
-  LPSolver lp;
-  lp.setMaximize();
-  typedef LPSolver::Col Col;
-  typedef LPSolver::Row Row;
-  typedef Graph::EdgeMap<Col> EdgeIndexMap;
-  typedef Graph::NodeMap<Row> NodeIndexMap;
-  EdgeIndexMap edge_index_map(g);
-  NodeIndexMap node_index_map(g);
-  PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
-
-  // nonnegativity of flow and capacity function
-  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
-    Col col=lp.addCol();
-    edge_index_map.set(e, col);
-    // interesting property in GLPK:
-    // if you change the order of the following two lines, the 
-    // two runs of GLPK are extremely different
-      lp.setColLowerBound(col, 0);
-      lp.setColUpperBound(col, cap[e]);
-  }
-  
-  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
-    LPSolver::Expression expr;
-    for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
-      expr+=edge_index_map[e];
-    for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
-      expr-=edge_index_map[e];
-    // cost function
-    if (n==s) {
-      lp.setObjCoeffs(expr);      
-    }
-    // flow conservation constraints
-    if ((n!=s) && (n!=t)) {
-      node_index_map.set(n, lp.addRow(expr == 0.0));
-    }
-  }
-  lp.solveSimplex();
-  cout << "elapsed time: " << ts << endl;
-//   cout << "rows:" << endl;
-//   for (
-//        LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
-//        i!=INVALID;
-//        ++i) { 
-//     cout << i << " ";
-//   }
-//   cout << endl;
-//   cout << "cols:" << endl;
-//   for (
-//        LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
-//        i!=INVALID;
-//        ++i) { 
-//     cout << i << " ";
-//   }
-//   cout << endl;
-  lp.setMIP();
-  cout << "elapsed time: " << ts << endl;
-  for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
-    lp.setColInt(it);
-  }
-  cout << "elapsed time: " << ts << endl;
-  lp.solveBandB();
-  cout << "elapsed time: " << ts << endl;
-}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp/min_cost_gen_flow.h
--- a/src/work/athos/lp/min_cost_gen_flow.h	Tue Mar 22 16:49:30 2005 +0000
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,268 +0,0 @@
-// -*- c++ -*-
-#ifndef LEMON_MIN_COST_GEN_FLOW_H
-#define LEMON_MIN_COST_GEN_FLOW_H
-#include <iostream>
-//#include <fstream>
-
-#include <lemon/smart_graph.h>
-#include <lemon/list_graph.h>
-//#include <lemon/dimacs.h>
-//#include <lemon/time_measure.h>
-//#include <graph_wrapper.h>
-#include <lemon/preflow.h>
-#include <lemon/min_cost_flow.h>
-//#include <augmenting_flow.h>
-//#include <preflow_res.h>
-#include <work/marci/merge_node_graph_wrapper.h>
-#include <work/marci/lp/lp_solver_wrapper_3.h>
-
-namespace lemon {
-
-  template<typename Edge, typename EdgeIndexMap> 
-  class PrimalMap {
-  protected:
-    LPGLPK* lp;
-    EdgeIndexMap* edge_index_map;
-  public:
-    PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
-      lp(&_lp), edge_index_map(&_edge_index_map) { }
-    double operator[](Edge e) const { 
-      return lp->getPrimal((*edge_index_map)[e]);
-    }
-  };
-
-  // excess: rho-delta egyelore csak =0-ra.
-  template <typename Graph, typename Num,
-	    typename Excess=typename Graph::template NodeMap<Num>, 
-	    typename LCapMap=typename Graph::template EdgeMap<Num>,
-	    typename CapMap=typename Graph::template EdgeMap<Num>,
-            typename FlowMap=typename Graph::template EdgeMap<Num>,
-	    typename CostMap=typename Graph::template EdgeMap<Num> >
-  class MinCostGenFlow {
-  protected:
-    const Graph& g;
-    const Excess& excess;
-    const LCapMap& lcapacity;
-    const CapMap& capacity;
-    FlowMap& flow;
-    const CostMap& cost;
-  public:
-    MinCostGenFlow(const Graph& _g, const Excess& _excess, 
-		   const LCapMap& _lcapacity, const CapMap& _capacity, 
-		   FlowMap& _flow, 
-		   const CostMap& _cost) :
-      g(_g), excess(_excess), lcapacity(_lcapacity),
-      capacity(_capacity), flow(_flow), cost(_cost) { }
-    bool feasible() {
-      //      std::cout << "making new vertices..." << std::endl; 
-      typedef ListGraph Graph2;
-      Graph2 g2;
-      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
-      //      std::cout << "merging..." << std::endl; 
-      GW gw(g, g2);
-      typename GW::Node s(INVALID, g2.addNode(), true);
-      typename GW::Node t(INVALID, g2.addNode(), true);
-      typedef SmartGraph Graph3;
-      //      std::cout << "making extender graph..." << std::endl; 
-      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
-//       {
-// 	checkConcept<StaticGraph, GWW>();   
-//       }
-      GWW gww(gw);
-      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
-      GWWW gwww(gw, gww);
-
-      //      std::cout << "making new edges..." << std::endl; 
-      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
-
-      for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
-	translated_cap.set(typename GWWW::Edge(e,INVALID,false), 
-			   capacity[e]-lcapacity[e]);
-	//	cout << "t_cap " << gw.id(e) << " " 
-	//	     << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
-      }
-
-      Num expected=0;
-
-      //      std::cout << "making new edges 2..." << std::endl; 
-      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
-	Num a=0;
-	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
-	  a+=lcapacity[e];
-	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
-	  a-=lcapacity[e];
-	if (excess[n]>a) {
-	  typename GWW::Edge e=
-	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
-	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
-			     excess[n]-a);
-	  //	  std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
-	}
-	if (excess[n]<a) {
-	  typename GWW::Edge e=
-	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
-	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
-			     a-excess[n]);
-	  expected+=a-excess[n];
-	  //	  std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
-	}
-      }
-
-      //      std::cout << "preflow..." << std::endl; 
-      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
-      Preflow<GWWW, Num> preflow(gwww, s, t, 
-				 translated_cap, translated_flow);
-      preflow.run();
-      //      std::cout << "fv: " << preflow.flowValue() << std::endl; 
-      //      std::cout << "expected: " << expected << std::endl; 
-
-      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
-	typename GW::Edge ew(e, INVALID, false);
-	typename GWWW::Edge ewww(ew, INVALID, false);
-	flow.set(e, translated_flow[ewww]+lcapacity[e]);
-      }
-      return (preflow.flowValue()>=expected);
-    }
-    // for nonnegative costs
-    bool run() {
-      //      std::cout << "making new vertices..." << std::endl; 
-      typedef ListGraph Graph2;
-      Graph2 g2;
-      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
-      //      std::cout << "merging..." << std::endl; 
-      GW gw(g, g2);
-      typename GW::Node s(INVALID, g2.addNode(), true);
-      typename GW::Node t(INVALID, g2.addNode(), true);
-      typedef SmartGraph Graph3;
-      //      std::cout << "making extender graph..." << std::endl; 
-      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
-//       {
-// 	checkConcept<StaticGraph, GWW>();   
-//       }
-      GWW gww(gw);
-      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
-      GWWW gwww(gw, gww);
-
-      //      std::cout << "making new edges..." << std::endl; 
-      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
-
-      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
-	typename GW::Edge ew(e, INVALID, false);
-	typename GWWW::Edge ewww(ew, INVALID, false);
-	translated_cap.set(ewww, capacity[e]-lcapacity[e]);
-	//	cout << "t_cap " << g.id(e) << " " 
-	//	     << translated_cap[ewww] << endl;
-      }
-
-      Num expected=0;
-
-      //      std::cout << "making new edges 2..." << std::endl; 
-      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
-	//	std::cout << "node: " << g.id(n) << std::endl;
-	Num a=0;
-	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
-	  a+=lcapacity[e];
-	  //	  std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
-	}
-	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
-	  a-=lcapacity[e];
-	  //	  std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
-	}
-	//	std::cout << "excess " << g.id(n) << ": " << a << std::endl;
-	if (0>a) {
-	  typename GWW::Edge e=
-	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
-	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
-			     -a);
-	  //	  std::cout << g.id(n) << "->t " << -a << std::endl;
-	}
-	if (0<a) {
-	  typename GWW::Edge e=
-	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
-	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
-			     a);
-	  expected+=a;
-	  //	  std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
-	}
-      }
-
-      //      std::cout << "preflow..." << std::endl; 
-      typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
-      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
-	translated_cost.set(typename GWWW::Edge(
-        typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
-      }
-      //      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
-      MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>, 
-      typename GWWW::template EdgeMap<Num> > 
-      min_cost_flow(gwww, translated_cost, translated_cap, 
-		    s, t);
-      while (min_cost_flow.augment()) { }
-      std::cout << "fv: " << min_cost_flow.flowValue() << std::endl; 
-      std::cout << "expected: " << expected << std::endl; 
-
-      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
-	typename GW::Edge ew(e, INVALID, false);
-	typename GWWW::Edge ewww(ew, INVALID, false);
-	//	std::cout << g.id(e) << " " << flow[e] << std::endl;
-	flow.set(e, lcapacity[e]+
-		 min_cost_flow.getFlow()[ewww]);
-      }
-      return (min_cost_flow.flowValue()>=expected);
-    }
-    void runByLP() {
-      typedef LPGLPK LPSolver;
-      LPSolver lp;
-      lp.setMinimize();
-      typedef LPSolver::ColIt ColIt;
-      typedef LPSolver::RowIt RowIt;
-      typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
-      EdgeIndexMap edge_index_map(g);
-      PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
-      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
-	ColIt col_it=lp.addCol();
-	edge_index_map.set(e, col_it);
-	if (lcapacity[e]==capacity[e])
-	  lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
-	else 
-	  lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
-	lp.setObjCoef(col_it, cost[e]);
-      }
-      LPSolver::ColIt col_it;
-      for (lp.col_iter_map.first(col_it, lp.VALID_CLASS); 
-	   lp.col_iter_map.valid(col_it); 
-	   lp.col_iter_map.next(col_it)) {
-//	std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
-      }
-      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
-	typename Graph::template EdgeMap<Num> coeffs(g, 0);
-	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
-	coeffs.set(e, coeffs[e]+1);
-	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
-	coeffs.set(e, coeffs[e]-1);
-	RowIt row_it=lp.addRow();
-	typename std::vector< std::pair<ColIt, double> > row;
-	//std::cout << "node:" <<g.id(n)<<std::endl;
-	for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
-	  if (coeffs[e]!=0) {
-	    //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
-	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
-	  }
-	}
-	//std::cout << std::endl;
-	//std::cout << " " << g.id(n) << " " << row.size() << std::endl;
-	lp.setRowCoeffs(row_it, row.begin(), row.end());
-	lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
-      }
-      lp.solveSimplex();
-      //std::cout << lp.colNum() << std::endl;
-      //std::cout << lp.rowNum() << std::endl;
-      //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
-      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) 
-      flow.set(e, lp_flow[e]);
-    }
-  };
-
-} // namespace lemon
-
-#endif //LEMON_MIN_COST_GEN_FLOW_H
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/expression.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/expression.h	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,197 @@
+// -*- c++ -*-
+#ifndef LEMON_EXPRESSION_H
+#define LEMON_EXPRESSION_H
+
+#include <iostream>
+#include <map>
+#include <limits>
+
+namespace lemon {
+
+  /*! \brief Linear expression
+
+    \c Expr<_Col,_Value> implements a class of linear expressions with the 
+    operations of addition and multiplication with scalar. 
+
+    \author Marton Makai
+   */
+  template <typename _Col, typename _Value>
+  class Expr;
+
+  template <typename _Col, typename _Value>
+  class Expr {
+//  protected:
+  public:
+    typedef 
+    typename std::map<_Col, _Value> Data; 
+    Data data;
+  public:
+    void simplify() {
+      for (typename Data::iterator i=data.begin(); 
+	   i!=data.end(); ++i) {
+	if ((*i).second==0) data.erase(i);
+      }
+    }
+    Expr() { }
+    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;
+      }
+      simplify();
+      return *this;
+    }
+    Expr& operator+=(const Expr<_Col, _Value>& expr) {
+      for (typename Data::const_iterator j=expr.data.begin(); 
+	   j!=expr.data.end(); ++j) {
+	typename Data::iterator i=data.find((*j).first);
+	if (i==data.end()) {
+	  data.insert(std::make_pair((*j).first, (*j).second));
+	} else {
+	  (*i).second+=(*j).second;
+	}
+      }
+      simplify();
+      return *this;
+    }
+    Expr& operator-=(const Expr<_Col, _Value>& expr) {
+      for (typename Data::const_iterator j=expr.data.begin(); 
+	   j!=expr.data.end(); ++j) {
+	typename Data::iterator i=data.find((*j).first);
+	if (i==data.end()) {
+	  data.insert(std::make_pair((*j).first, -(*j).second));
+	} else {
+	  (*i).second+=-(*j).second;
+	}
+      }
+      simplify();
+      return *this;
+    }
+    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, _Col _col) {
+    Expr<_Col, _Value> tmp(_col);
+    tmp*=_value;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator*(_Value _value, 
+			       const Expr<_Col, _Value>& expr) {
+    Expr<_Col, _Value> tmp(expr);
+    tmp*=_value;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator+(const Expr<_Col, _Value>& expr1, 
+			       const Expr<_Col, _Value>& expr2) {
+    Expr<_Col, _Value> tmp(expr1);
+    tmp+=expr2;
+    tmp.simplify();
+    return tmp;
+  }
+
+  template <typename _Col, typename _Value>
+  Expr<_Col, _Value> operator-(const Expr<_Col, _Value>& expr1, 
+			       const Expr<_Col, _Value>& expr2) {
+    Expr<_Col, _Value> tmp(expr1);
+    tmp-=expr2;
+    tmp.simplify();
+    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;
+  }
+
+  template <typename _Col, typename _Value>
+  class LConstr {
+    //  protected:
+  public:
+    Expr<_Col, _Value> expr;
+    _Value lo;
+  public:
+    LConstr(const Expr<_Col, _Value>& _expr, _Value _lo) : 
+      expr(_expr), lo(_lo) { }
+  };
+  
+  template <typename _Col, typename _Value>
+  LConstr<_Col, _Value> 
+  operator<=(_Value lo, const Expr<_Col, _Value>& expr) {
+    return LConstr<_Col, _Value>(expr, lo);
+  }
+
+  template <typename _Col, typename _Value>
+  class UConstr {
+    //  protected:
+  public:
+    Expr<_Col, _Value> expr;
+    _Value up;
+  public:
+    UConstr(const Expr<_Col, _Value>& _expr, _Value _up) : 
+      expr(_expr), up(_up) { }
+  };
+
+  template <typename _Col, typename _Value>
+  UConstr<_Col, _Value> 
+  operator<=(const Expr<_Col, _Value>& expr, _Value up) {
+    return UConstr<_Col, _Value>(expr, up);
+  }
+
+  template <typename _Col, typename _Value>
+  class Constr {
+    //  protected:
+  public:
+    Expr<_Col, _Value> expr;
+    _Value lo, up;
+  public:
+    Constr(const Expr<_Col, _Value>& _expr, _Value _lo, _Value _up) : 
+      expr(_expr), lo(_lo), up(_up) { }
+    Constr(const LConstr<_Col, _Value>& _lconstr) : 
+      expr(_lconstr.expr), 
+      lo(_lconstr.lo), 
+      up(std::numeric_limits<_Value>::infinity()) { }
+    Constr(const UConstr<_Col, _Value>& _uconstr) : 
+      expr(_uconstr.expr), 
+      lo(-std::numeric_limits<_Value>::infinity()), 
+      up(_uconstr.up) { }
+  };
+
+  template <typename _Col, typename _Value>
+  Constr<_Col, _Value> 
+  operator<=(const LConstr<_Col, _Value>& lconstr, _Value up) {
+    return Constr<_Col, _Value>(lconstr.expr, lconstr.lo, up);
+  }
+
+  template <typename _Col, typename _Value>
+  Constr<_Col, _Value> 
+  operator<=(_Value lo, const UConstr<_Col, _Value>& uconstr) {
+    return Constr<_Col, _Value>(uconstr.expr, lo, uconstr.up);
+  }
+
+  template <typename _Col, typename _Value>
+  Constr<_Col, _Value> 
+  operator==(const Expr<_Col, _Value>& expr, _Value value) {
+    return Constr<_Col, _Value>(expr, value, value);
+  }
+  
+} //namespace lemon
+
+#endif //LEMON_EXPRESSION_H
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/expression_test.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/expression_test.cc	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,23 @@
+#include <expression.h>
+#include <iostream>
+#include <string>
+
+using std::cout;
+using std::endl;
+using std::string;
+using namespace lemon;
+
+int main() {
+  Expr<string, double> b;
+  cout << b << endl;
+  Expr<string, double> c("f");
+  cout << c << endl;
+  Expr<string, double> d=8.0*string("g");
+  cout << d << endl;
+  c*=5;
+  cout << c << endl;
+  Expr<string, double> e=c;
+  e+=8.9*9.0*string("l");
+  cout << e << endl;
+  cout << c+d << endl;
+}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/lp_solver_base.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/lp_solver_base.h	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,639 @@
+// -*- 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
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/lp_solver_glpk.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/lp_solver_glpk.h	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,545 @@
+// -*- c++ -*-
+#ifndef LEMON_LP_SOLVER_GLPK_H
+#define LEMON_LP_SOLVER_GLPK_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 <lp_solver_base.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 {
+  
+  
+  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 Of course this is not efficient at all, but GLPK knows not more.
+      int change_this;
+      bool get_set_row;
+      //The only thing we can do is optimize on whether working with a row 
+      //or a coloumn
+      int row_num = rowNum();
+      int col_num = colNum();
+      if (col_num<row_num){
+	//There are more rows than coloumns
+	get_set_row=true;
+	int mem_length=1+row_num;
+	int* indices = new int[mem_length];
+	double* doubles = new double[mem_length];
+	int length=lpx_get_mat_col(lp, i, indices, doubles);
+      }else{
+	get_set_row=false;
+	int mem_length=1+col_num;
+	int* indices = new int[mem_length];
+	double* doubles = new double[mem_length];
+	int length=lpx_get_mat_row(lp, i, indices, doubles);
+      }
+      //Itten      
+int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=lpx_get_mat_col(lp, i, indices, doubles);
+ 
+      delete [] indices;
+      delete [] doubles;
+
+    }
+    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);
+	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 _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_GLPK_H
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/lp_solver_wrapper.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/lp_solver_wrapper.h	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,431 @@
+// -*- 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 <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; }
+  };
+  
+  /// \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 LPSolverWrapper {
+  public:
+
+//   class Row {
+//   protected:
+//     int i;
+//   public:
+//     Row() { }
+//     Row(const Invalid&) : i(0) { }
+//     Row(const int& _i) : i(_i) { }
+//     operator int() const { return i; }
+//   };
+//   class RowIt : public Row {
+//   public:
+//     RowIt(const Row& row) : Row(row) { }
+//   };
+
+//   class Col {
+//   protected:
+//     int i;
+//   public:
+//     Col() { }
+//     Col(const Invalid&) : i(0) { }
+//     Col(const int& _i) : i(_i) { }
+//     operator int() const { return i; }
+//   };
+//   class ColIt : public Col {
+//     ColIt(const Col& col) : Col(col) { }
+//   };
+
+  public:
+    ///.
+    LPX* lp;
+    ///.
+    typedef IterablePartition<int>::ClassIt RowIt;
+    ///.
+    IterablePartition<int> row_iter_map;
+    ///.
+    typedef IterablePartition<int>::ClassIt ColIt;
+    ///.
+    IterablePartition<int> col_iter_map;
+    //std::vector<int> row_id_to_lp_row_id;
+    //std::vector<int> col_id_to_lp_col_id;
+    ///.
+    const int VALID_ID;
+    ///.
+    const int INVALID_ID;
+
+  public:
+    ///.
+    LPSolverWrapper() : lp(lpx_create_prob()), 
+			row_iter_map(2), 
+			col_iter_map(2), 
+			//row_id_to_lp_row_id(), col_id_to_lp_col_id(), 
+			VALID_ID(0), INVALID_ID(1) {
+      lpx_set_int_parm(lp, LPX_K_DUAL, 1);
+    }
+    ///.
+    ~LPSolverWrapper() {
+      lpx_delete_prob(lp);
+    }
+    ///.
+    void setMinimize() { 
+      lpx_set_obj_dir(lp, LPX_MIN);
+    }
+    ///.
+    void setMaximize() { 
+      lpx_set_obj_dir(lp, LPX_MAX);
+    }
+    ///.
+    ColIt addCol() {
+      int i=lpx_add_cols(lp, 1);  
+      ColIt col_it;
+      col_iter_map.first(col_it, INVALID_ID);
+      if (col_iter_map.valid(col_it)) { //van hasznalhato hely
+	col_iter_map.set(col_it, INVALID_ID, VALID_ID);
+	col_iter_map[col_it]=i;
+	//col_id_to_lp_col_id[col_iter_map[col_it]]=i;
+      } else { //a cucc vegere kell inzertalni mert nincs szabad hely
+	//col_id_to_lp_col_id.push_back(i);
+	//int j=col_id_to_lp_col_id.size()-1;
+	col_it=col_iter_map.push_back(i, VALID_ID);
+      }
+//    edge_index_map.set(e, i);
+//    lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
+//    lpx_set_obj_coef(lp, i, cost[e]);    
+      return col_it;
+    }
+    ///.
+    RowIt addRow() {
+      int i=lpx_add_rows(lp, 1);  
+      RowIt row_it;
+      row_iter_map.first(row_it, INVALID_ID);
+      if (row_iter_map.valid(row_it)) { //van hasznalhato hely
+	row_iter_map.set(row_it, INVALID_ID, VALID_ID);
+	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_ID);
+      }
+      return row_it;
+    }
+    //pair<RowIt, double>-bol kell megadni egy std range-et
+    ///.
+    template <typename Begin, typename End>
+    void setColCoeffs(const ColIt& col_it, 
+		      Begin begin, End end) {
+      int mem_length=1+lpx_get_num_rows(lp);
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=0;
+      for ( ; begin!=end; ++begin) {
+	++length;
+	indices[length]=row_iter_map[begin->first];
+	doubles[length]=begin->second;
+      }
+      lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
+      delete [] indices;
+      delete [] doubles;
+    }
+    //pair<ColIt, double>-bol kell megadni egy std range-et
+    ///.
+    template <typename Begin, typename End>
+    void setRowCoeffs(const RowIt& row_it, 
+		      Begin begin, End end) {
+      int mem_length=1+lpx_get_num_cols(lp);
+      int* indices = new int[mem_length];
+      double* doubles = new double[mem_length];
+      int length=0;
+      for ( ; begin!=end; ++begin) {
+	++length;
+	indices[length]=col_iter_map[begin->first];
+	doubles[length]=begin->second;
+      }
+      lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
+      delete [] indices;
+      delete [] doubles;
+    }
+    ///.
+    void eraseCol(const ColIt& col_it) {
+      col_iter_map.set(col_it, VALID_ID, INVALID_ID);
+      int cols[2];
+      cols[1]=col_iter_map[col_it];
+      lpx_del_cols(lp, 1, cols);
+      col_iter_map[col_it]=0; //glpk specifikus
+      ColIt it;
+      for (col_iter_map.first(it, VALID_ID); 
+	   col_iter_map.valid(it); col_iter_map.next(it)) {
+	if (col_iter_map[it]>cols[1]) --col_iter_map[it];
+      }
+    }
+    ///.
+    void eraseRow(const RowIt& row_it) {
+      row_iter_map.set(row_it, VALID_ID, INVALID_ID);
+      int rows[2];
+      rows[1]=row_iter_map[row_it];
+      lpx_del_rows(lp, 1, rows);
+      row_iter_map[row_it]=0; //glpk specifikus
+      RowIt it;
+      for (row_iter_map.first(it, VALID_ID); 
+	   row_iter_map.valid(it); row_iter_map.next(it)) {
+	if (row_iter_map[it]>rows[1]) --row_iter_map[it];
+      }
+    }
+    ///.
+    void setColBounds(const ColIt& col_it, int bound_type, 
+		      double lo, double up) {
+      lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
+    }
+    ///.
+    double getObjCoef(const ColIt& col_it) { 
+      return lpx_get_obj_coef(lp, col_iter_map[col_it]);
+    }
+    ///.
+    void setRowBounds(const RowIt& row_it, int bound_type, 
+		      double lo, double up) {
+      lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
+    }
+    ///.
+    void setObjCoef(const ColIt& col_it, double obj_coef) { 
+      lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
+    }
+    ///.
+    void solveSimplex() { lpx_simplex(lp); }
+    ///.
+    void solvePrimalSimplex() { lpx_simplex(lp); }
+    ///.
+    void solveDualSimplex() { lpx_simplex(lp); }
+    ///.
+    double getPrimal(const ColIt& col_it) {
+      return lpx_get_col_prim(lp, col_iter_map[col_it]);
+    }
+    ///.
+    double getObjVal() { return lpx_get_obj_val(lp); }
+    ///.
+    int rowNum() const { return lpx_get_num_rows(lp); }
+    ///.
+    int colNum() const { return lpx_get_num_cols(lp); }
+    ///.
+    int warmUp() { return lpx_warm_up(lp); }
+    ///.
+    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;
+      }
+    }
+    ///.
+    int getPrimalStatus() { return lpx_get_prim_stat(lp); }
+    ///.
+    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;
+      }
+    }
+    ///.
+    int getDualStatus() { return lpx_get_dual_stat(lp); }
+    ///.
+    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]); 
+    }
+    ///.
+    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]); 
+    }
+    ///.
+    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
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/magic_square.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/magic_square.cc	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,99 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/time_measure.h>
+#include <lp_solver_glpk.h>
+
+using std::cout;
+using std::endl;
+using namespace lemon;
+
+/*
+  On an 1537Mhz PC, the run times with 
+  glpk are the following.
+  for n=3,4, some secondes
+  for n=5, 25 hours
+ */
+
+int main(int, char **) {
+  const int n=3;
+  const double row_sum=(1.0+n*n)*n/2;
+  Timer ts;
+  ts.reset();
+  typedef LpGlpk LPSolver;
+  typedef LPSolver::Col Col;
+  LPSolver lp;
+  typedef std::map<std::pair<int, int>, Col> Coords;
+  Coords x;
+  // we create a new variable for each entry 
+  // of the magic square
+  for (int i=1; i<=n; ++i) {
+    for (int j=1; j<=n; ++j) { 
+      Col col=lp.addCol();
+      x[std::make_pair(i,j)]=col;
+      lp.setColLowerBound(col, 1.0);
+      lp.setColUpperBound(col, double(n*n));
+    }
+  }
+  LPSolver::Expression expr3, expr4;
+  for (int i=1; i<=n; ++i) {
+    LPSolver::Expression expr1, expr2;
+    for (int j=1; j<=n; ++j) {
+      expr1+=x[std::make_pair(i, j)];
+      expr2+=x[std::make_pair(j, i)];
+    }
+
+    // sum of rows and columns
+    lp.addRow(expr1==row_sum);
+    lp.addRow(expr2==row_sum);
+      cout <<"Expr1: "<<expr1<<endl;
+      cout <<"Expr2: "<<expr2<<endl;
+
+    expr3+=x[std::make_pair(i, i)];
+    expr4+=x[std::make_pair(i, (n+1)-i)];
+  }
+  cout <<"Expr3: "<<expr3<<endl;
+  cout <<"Expr4: "<<expr4<<endl;
+  // sum of the diagonal entries
+  lp.addRow(expr3==row_sum);
+  lp.addRow(expr4==row_sum);
+  lp.solveSimplex();
+  cout << "elapsed time: " << ts << endl;
+  for (int i=1; i<=n; ++i) {
+    for (int j=1; j<=n; ++j) { 
+      cout << "x("<<i<<","<<j<<")="<<lp.getPrimal(x[std::make_pair(i,j)]) 
+	   << endl;
+    }
+  }
+
+
+
+//   // we make new binary variables for each pair of 
+//   // entries of the square to achieve that 
+//   // the values of different entries are different
+//   lp.setMIP();
+//   for (Coords::const_iterator it=x.begin(); it!=x.end(); ++it) {
+//     Coords::const_iterator jt=it; ++jt;
+//     for(; jt!=x.end(); ++jt) {
+//       Col col1=(*it).second;
+//       Col col2=(*jt).second;
+//       Col col=lp.addCol();
+//       lp.setColLowerBound(col, 0.0);
+//       lp.setColUpperBound(col, 1.0);
+//       lp.addRow(double(-n*n+1.0)<=1.0*col2-1.0*col1-double(n*n)*col<=-1.0);
+//       lp.setColInt(col);
+//     }
+//   }
+//   cout << "elapsed time: " << ts << endl;
+//   lp.solveSimplex();
+//   // let's solve the integer problem
+//   lp.solveBandB();
+//   cout << "elapsed time: " << ts << endl;
+//   for (int i=1; i<=n; ++i) {
+//     for (int j=1; j<=n; ++j) { 
+//       cout << "x("<<i<<","<<j<<")="<<lp.getMIPPrimal(x[std::make_pair(i,j)]) 
+// 	   << endl;
+//     }
+//   }
+}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/makefile
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/makefile	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,73 @@
+#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
+#GLPKROOT = /usr/local/glpk-4.4
+INCLUDEDIRS ?= -I/home/marci/boost -I/usr/local/cplex/cplex75/include -I../../{marci,alpar,klao,akos,athos} -I. -I../../.. -I../.. -I..# -I$(GLPKROOT)/include
+#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
+CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
+LDFLAGS  =  -lglpk#-lcplex -lm -lpthread -lilocplex -L/usr/local/cplex/cplex75/lib/i86_linux2_glibc2.2_gcc3.0/static_mt# -L$(GLPKROOT)/lib
+
+BINARIES = magic_square max_flow_expression #expression_test max_flow_by_lp# sample sample2 sample11 sample15
+
+#include ../makefile
+
+# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
+# ismert rendszeren :-)  (Misi)
+ifdef GCCVER
+CXX := g++-$(GCCVER)
+else
+CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
+endif
+
+ifdef DEBUG
+CXXFLAGS += -DDEBUG
+endif
+
+CC := $(CXX)
+
+
+all: $(BINARIES)
+
+################
+# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
+# all elo.
+# Kiveve ha siman file.cc -> file  esetrol van szo, amikor is nem kell
+# irni semmit.
+
+#proba: proba.o seged.o
+
+################
+
+
+# .depend dep depend:
+# 	-$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
+
+#makefile: .depend
+#sinclude .depend
+#moert nem megy az eredeti /usr/bin/ld-vel?
+
+# %: %.o
+# 	$(CXX) -o $@ $< $(LDFLAGS)
+
+# %.o: %.cc
+# 	$(CXX) $(CXXFLAGS) -c $<
+
+%: %.cc
+	$(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
+
+sample11prof: sample11prof.o
+	 $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
+sample11prof.o: sample11.cc
+	$(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
+
+# sample.o: sample.cc
+# 	$(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
+
+# sample2: sample2.o
+# 	$(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
+# sample2.o: sample2.cc
+# 	$(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
+
+
+clean:
+	$(RM) *.o $(BINARIES) .depend
+
+.PHONY: all clean dep depend
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/max_flow_by_lp.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/max_flow_by_lp.cc	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,186 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+//#include <graph_wrapper.h>
+#include <lemon/preflow.h>
+#include <augmenting_flow.h>
+//#include <preflow_res.h>
+//#include <lp_solver_wrapper_2.h>
+#include <min_cost_gen_flow.h>
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_demo < dimacs_max_flow_file
+
+using namespace lemon;
+
+int main(int, char **) {
+
+  typedef ListGraph MutableGraph;
+  typedef ListGraph Graph;
+  typedef Graph::Node Node;
+  typedef Graph::Edge Edge;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph g;
+
+  Node s, t;
+  Graph::EdgeMap<int> cap(g);
+  //readDimacsMaxFlow(std::cin, g, s, t, cap);
+  readDimacs(std::cin, g, cap, s, t);
+  Timer ts;
+  Graph::EdgeMap<int> flow(g); //0 flow
+  Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
+    max_flow_test(g, s, t, cap, flow);
+  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
+    augmenting_flow_test(g, s, t, cap, flow);
+  
+  Graph::NodeMap<bool> cut(g);
+
+  {
+    std::cout << "preflow ..." << std::endl;
+    ts.reset();
+    max_flow_test.run();
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+    max_flow_test.minCut(cut);
+
+    for (EdgeIt e(g); e!=INVALID; ++e) {
+      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+	std::cout << "Slackness does not hold!" << std::endl;
+    }
+  }
+
+//   {
+//     std::cout << "preflow ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     max_flow_test.preflow(Preflow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+//   {
+//     std::cout << "wrapped preflow ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     pre_flow_res.run();
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
+//   }
+
+  {
+    std::cout << "physical blocking flow augmentation ..." << std::endl;
+    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
+    ts.reset();
+    int i=0;
+    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "number of augmentation phases: " << i << std::endl; 
+    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+    for (EdgeIt e(g); e!=INVALID; ++e) {
+      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+	std::cout << "Slackness does not hold!" << std::endl;
+    }
+  }
+
+//   {
+//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "number of augmentation phases: " << i << std::endl; 
+//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
+//   }
+
+  {
+    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
+    for (EdgeIt e(g); e!=INVALID; ++e) flow.set(e, 0);
+    ts.reset();
+    int i=0;
+    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
+    std::cout << "elapsed time: " << ts << std::endl;
+    std::cout << "number of augmentation phases: " << i << std::endl; 
+    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+    for (EdgeIt e(g); e!=INVALID; ++e) {
+      if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+	std::cout << "Slackness does not hold!" << std::endl;
+    }
+  }
+
+//   {
+//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "number of augmentation phases: " << i << std::endl; 
+//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+//   {
+//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
+//     ts.reset();
+//     int i=0;
+//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
+//     std::cout << "elapsed time: " << ts << std::endl;
+//     std::cout << "number of augmentation phases: " << i << std::endl; 
+//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
+
+//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+//       if (cut[g.source(e)] && !cut[g.target(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.source(e)] && cut[g.target(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+  ts.reset();
+
+  Edge e=g.addEdge(t, s);
+  Graph::EdgeMap<int> cost(g, 0);
+  cost.set(e, -1);
+  cap.set(e, 10000);
+  typedef ConstMap<Node, int> Excess;
+  Excess excess(0);
+  typedef ConstMap<Edge, int> LCap;
+  LCap lcap(0);
+
+  MinCostGenFlow<Graph, int, Excess, LCap> 
+    min_cost(g, excess, lcap, cap, flow, cost); 
+  min_cost.feasible();
+  min_cost.runByLP();
+
+  std::cout << "elapsed time: " << ts << std::endl;
+  std::cout << "flow value: "<< flow[e] << std::endl;
+}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/max_flow_expression.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/max_flow_expression.cc	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,109 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <lemon/graph_utils.h>
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+#include <lemon/dimacs.h>
+#include <lemon/time_measure.h>
+#include <lp_solver_glpk.h>
+
+using std::cout;
+using std::endl;
+using namespace lemon;
+
+template<typename Edge, typename EdgeIndexMap> 
+class PrimalMap {
+protected:
+  LpGlpk* lp;
+  EdgeIndexMap* edge_index_map;
+public:
+  PrimalMap(LpGlpk& _lp, EdgeIndexMap& _edge_index_map) : 
+    lp(&_lp), edge_index_map(&_edge_index_map) { }
+  double operator[](Edge e) const { 
+    return lp->getPrimal((*edge_index_map)[e]);
+  }
+};
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_expresion < dimacs_max_flow_file
+
+int main(int, char **) {
+
+  typedef ListGraph Graph;
+  typedef Graph::Node Node;
+  typedef Graph::Edge Edge;
+  typedef Graph::EdgeIt EdgeIt;
+
+  Graph g;
+  
+  Node s, t;
+  Graph::EdgeMap<int> cap(g);
+  readDimacs(std::cin, g, cap, s, t);
+  Timer ts;
+  
+  typedef LpGlpk LPSolver;
+  LPSolver lp;
+  lp.setMaximize();
+  typedef LPSolver::Col Col;
+  typedef LPSolver::Row Row;
+  typedef Graph::EdgeMap<Col> EdgeIndexMap;
+  typedef Graph::NodeMap<Row> NodeIndexMap;
+  EdgeIndexMap edge_index_map(g);
+  NodeIndexMap node_index_map(g);
+  PrimalMap<Edge, EdgeIndexMap> flow(lp, edge_index_map);
+
+  // nonnegativity of flow and capacity function
+  for (Graph::EdgeIt e(g); e!=INVALID; ++e) {
+    Col col=lp.addCol();
+    edge_index_map.set(e, col);
+    // interesting property in GLPK:
+    // if you change the order of the following two lines, the 
+    // two runs of GLPK are extremely different
+      lp.setColLowerBound(col, 0);
+      lp.setColUpperBound(col, cap[e]);
+  }
+  
+  for (Graph::NodeIt n(g); n!=INVALID; ++n) {
+    LPSolver::Expression expr;
+    for (Graph::OutEdgeIt e(g, n); e!=INVALID; ++e)
+      expr+=edge_index_map[e];
+    for (Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+      expr-=edge_index_map[e];
+    // cost function
+    if (n==s) {
+      lp.setObjCoeffs(expr);      
+    }
+    // flow conservation constraints
+    if ((n!=s) && (n!=t)) {
+      node_index_map.set(n, lp.addRow(expr == 0.0));
+    }
+  }
+  lp.solveSimplex();
+  cout << "elapsed time: " << ts << endl;
+//   cout << "rows:" << endl;
+//   for (
+//        LPSolver::Rows::ClassIt i(lp.row_iter_map, 0);
+//        i!=INVALID;
+//        ++i) { 
+//     cout << i << " ";
+//   }
+//   cout << endl;
+//   cout << "cols:" << endl;
+//   for (
+//        LPSolver::Cols::ClassIt i(lp.col_iter_map, 0);
+//        i!=INVALID;
+//        ++i) { 
+//     cout << i << " ";
+//   }
+//   cout << endl;
+  lp.setMIP();
+  cout << "elapsed time: " << ts << endl;
+  for (LPSolver::Cols::ClassIt it(lp.col_iter_map ,1); it!=INVALID; ++it) {
+    lp.setColInt(it);
+  }
+  cout << "elapsed time: " << ts << endl;
+  lp.solveBandB();
+  cout << "elapsed time: " << ts << endl;
+}
diff -r 41caee260bd4 -r 43a3d06e0ee0 src/work/athos/lp_old/min_cost_gen_flow.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/athos/lp_old/min_cost_gen_flow.h	Wed Mar 23 09:49:41 2005 +0000
@@ -0,0 +1,268 @@
+// -*- c++ -*-
+#ifndef LEMON_MIN_COST_GEN_FLOW_H
+#define LEMON_MIN_COST_GEN_FLOW_H
+#include <iostream>
+//#include <fstream>
+
+#include <lemon/smart_graph.h>
+#include <lemon/list_graph.h>
+//#include <lemon/dimacs.h>
+//#include <lemon/time_measure.h>
+//#include <graph_wrapper.h>
+#include <lemon/preflow.h>
+#include <lemon/min_cost_flow.h>
+//#include <augmenting_flow.h>
+//#include <preflow_res.h>
+#include <work/marci/merge_node_graph_wrapper.h>
+#include <work/marci/lp/lp_solver_wrapper_3.h>
+
+namespace lemon {
+
+  template<typename Edge, typename EdgeIndexMap> 
+  class PrimalMap {
+  protected:
+    LPGLPK* lp;
+    EdgeIndexMap* edge_index_map;
+  public:
+    PrimalMap(LPGLPK& _lp, EdgeIndexMap& _edge_index_map) : 
+      lp(&_lp), edge_index_map(&_edge_index_map) { }
+    double operator[](Edge e) const { 
+      return lp->getPrimal((*edge_index_map)[e]);
+    }
+  };
+
+  // excess: rho-delta egyelore csak =0-ra.
+  template <typename Graph, typename Num,
+	    typename Excess=typename Graph::template NodeMap<Num>, 
+	    typename LCapMap=typename Graph::template EdgeMap<Num>,
+	    typename CapMap=typename Graph::template EdgeMap<Num>,
+            typename FlowMap=typename Graph::template EdgeMap<Num>,
+	    typename CostMap=typename Graph::template EdgeMap<Num> >
+  class MinCostGenFlow {
+  protected:
+    const Graph& g;
+    const Excess& excess;
+    const LCapMap& lcapacity;
+    const CapMap& capacity;
+    FlowMap& flow;
+    const CostMap& cost;
+  public:
+    MinCostGenFlow(const Graph& _g, const Excess& _excess, 
+		   const LCapMap& _lcapacity, const CapMap& _capacity, 
+		   FlowMap& _flow, 
+		   const CostMap& _cost) :
+      g(_g), excess(_excess), lcapacity(_lcapacity),
+      capacity(_capacity), flow(_flow), cost(_cost) { }
+    bool feasible() {
+      //      std::cout << "making new vertices..." << std::endl; 
+      typedef ListGraph Graph2;
+      Graph2 g2;
+      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
+      //      std::cout << "merging..." << std::endl; 
+      GW gw(g, g2);
+      typename GW::Node s(INVALID, g2.addNode(), true);
+      typename GW::Node t(INVALID, g2.addNode(), true);
+      typedef SmartGraph Graph3;
+      //      std::cout << "making extender graph..." << std::endl; 
+      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
+//       {
+// 	checkConcept<StaticGraph, GWW>();   
+//       }
+      GWW gww(gw);
+      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
+      GWWW gwww(gw, gww);
+
+      //      std::cout << "making new edges..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
+
+      for (typename GW::EdgeIt e(gw); e!=INVALID; ++e) {
+	translated_cap.set(typename GWWW::Edge(e,INVALID,false), 
+			   capacity[e]-lcapacity[e]);
+	//	cout << "t_cap " << gw.id(e) << " " 
+	//	     << translated_cap[typename GWWW::Edge(e,INVALID,false)] << endl;
+      }
+
+      Num expected=0;
+
+      //      std::cout << "making new edges 2..." << std::endl; 
+      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
+	Num a=0;
+	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+	  a+=lcapacity[e];
+	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
+	  a-=lcapacity[e];
+	if (excess[n]>a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     excess[n]-a);
+	  //	  std::cout << g.id(n) << "->t " << excess[n]-a << std::endl;
+	}
+	if (excess[n]<a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     a-excess[n]);
+	  expected+=a-excess[n];
+	  //	  std::cout << "s->" << g.id(n) << " "<< a-excess[n] <<std:: endl;
+	}
+      }
+
+      //      std::cout << "preflow..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
+      Preflow<GWWW, Num> preflow(gwww, s, t, 
+				 translated_cap, translated_flow);
+      preflow.run();
+      //      std::cout << "fv: " << preflow.flowValue() << std::endl; 
+      //      std::cout << "expected: " << expected << std::endl; 
+
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	typename GW::Edge ew(e, INVALID, false);
+	typename GWWW::Edge ewww(ew, INVALID, false);
+	flow.set(e, translated_flow[ewww]+lcapacity[e]);
+      }
+      return (preflow.flowValue()>=expected);
+    }
+    // for nonnegative costs
+    bool run() {
+      //      std::cout << "making new vertices..." << std::endl; 
+      typedef ListGraph Graph2;
+      Graph2 g2;
+      typedef MergeEdgeGraphWrapper<const Graph, Graph2> GW;
+      //      std::cout << "merging..." << std::endl; 
+      GW gw(g, g2);
+      typename GW::Node s(INVALID, g2.addNode(), true);
+      typename GW::Node t(INVALID, g2.addNode(), true);
+      typedef SmartGraph Graph3;
+      //      std::cout << "making extender graph..." << std::endl; 
+      typedef NewEdgeSetGraphWrapper2<GW, Graph3> GWW;
+//       {
+// 	checkConcept<StaticGraph, GWW>();   
+//       }
+      GWW gww(gw);
+      typedef AugmentingGraphWrapper<GW, GWW> GWWW;
+      GWWW gwww(gw, gww);
+
+      //      std::cout << "making new edges..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_cap(gwww);
+
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	typename GW::Edge ew(e, INVALID, false);
+	typename GWWW::Edge ewww(ew, INVALID, false);
+	translated_cap.set(ewww, capacity[e]-lcapacity[e]);
+	//	cout << "t_cap " << g.id(e) << " " 
+	//	     << translated_cap[ewww] << endl;
+      }
+
+      Num expected=0;
+
+      //      std::cout << "making new edges 2..." << std::endl; 
+      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
+	//	std::cout << "node: " << g.id(n) << std::endl;
+	Num a=0;
+	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e) {
+	  a+=lcapacity[e];
+	  //	  std::cout << "bee: " << g.id(e) << " " << lcapacity[e] << std::endl;
+	}
+	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) {
+	  a-=lcapacity[e];
+	  //	  std::cout << "kie: " << g.id(e) << " " << lcapacity[e] << std::endl;
+	}
+	//	std::cout << "excess " << g.id(n) << ": " << a << std::endl;
+	if (0>a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(typename GW::Node(n,INVALID,false), t);
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     -a);
+	  //	  std::cout << g.id(n) << "->t " << -a << std::endl;
+	}
+	if (0<a) {
+	  typename GWW::Edge e=
+	    gww.addEdge(s, typename GW::Node(n,INVALID,false));
+	  translated_cap.set(typename GWWW::Edge(INVALID, e, true), 
+			     a);
+	  expected+=a;
+	  //	  std::cout << "s->" << g.id(n) << " "<< a <<std:: endl;
+	}
+      }
+
+      //      std::cout << "preflow..." << std::endl; 
+      typename GWWW::template EdgeMap<Num> translated_cost(gwww, 0);
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	translated_cost.set(typename GWWW::Edge(
+        typename GW::Edge(e, INVALID, false), INVALID, false), cost[e]);
+      }
+      //      typename GWWW::template EdgeMap<Num> translated_flow(gwww, 0);
+      MinCostFlow<GWWW, typename GWWW::template EdgeMap<Num>, 
+      typename GWWW::template EdgeMap<Num> > 
+      min_cost_flow(gwww, translated_cost, translated_cap, 
+		    s, t);
+      while (min_cost_flow.augment()) { }
+      std::cout << "fv: " << min_cost_flow.flowValue() << std::endl; 
+      std::cout << "expected: " << expected << std::endl; 
+
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	typename GW::Edge ew(e, INVALID, false);
+	typename GWWW::Edge ewww(ew, INVALID, false);
+	//	std::cout << g.id(e) << " " << flow[e] << std::endl;
+	flow.set(e, lcapacity[e]+
+		 min_cost_flow.getFlow()[ewww]);
+      }
+      return (min_cost_flow.flowValue()>=expected);
+    }
+    void runByLP() {
+      typedef LPGLPK LPSolver;
+      LPSolver lp;
+      lp.setMinimize();
+      typedef LPSolver::ColIt ColIt;
+      typedef LPSolver::RowIt RowIt;
+      typedef typename Graph::template EdgeMap<ColIt> EdgeIndexMap;
+      EdgeIndexMap edge_index_map(g);
+      PrimalMap<typename Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	ColIt col_it=lp.addCol();
+	edge_index_map.set(e, col_it);
+	if (lcapacity[e]==capacity[e])
+	  lp.setColBounds(col_it, LPSolver::FIXED, lcapacity[e], capacity[e]);
+	else 
+	  lp.setColBounds(col_it, LPSolver::DOUBLE, lcapacity[e], capacity[e]);
+	lp.setObjCoef(col_it, cost[e]);
+      }
+      LPSolver::ColIt col_it;
+      for (lp.col_iter_map.first(col_it, lp.VALID_CLASS); 
+	   lp.col_iter_map.valid(col_it); 
+	   lp.col_iter_map.next(col_it)) {
+//	std::cout << "ize " << lp.col_iter_map[col_it] << std::endl;
+      }
+      for (typename Graph::NodeIt n(g); n!=INVALID; ++n) {
+	typename Graph::template EdgeMap<Num> coeffs(g, 0);
+	for (typename Graph::InEdgeIt e(g, n); e!=INVALID; ++e)
+	coeffs.set(e, coeffs[e]+1);
+	for (typename Graph::OutEdgeIt e(g, n); e!=INVALID; ++e) 
+	coeffs.set(e, coeffs[e]-1);
+	RowIt row_it=lp.addRow();
+	typename std::vector< std::pair<ColIt, double> > row;
+	//std::cout << "node:" <<g.id(n)<<std::endl;
+	for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) {
+	  if (coeffs[e]!=0) {
+	    //std::cout << " edge:" <<g.id(e)<<" "<<coeffs[e];
+	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
+	  }
+	}
+	//std::cout << std::endl;
+	//std::cout << " " << g.id(n) << " " << row.size() << std::endl;
+	lp.setRowCoeffs(row_it, row.begin(), row.end());
+	lp.setRowBounds(row_it, LPSolver::FIXED, 0.0, 0.0);
+      }
+      lp.solveSimplex();
+      //std::cout << lp.colNum() << std::endl;
+      //std::cout << lp.rowNum() << std::endl;
+      //std::cout << "flow value: "<< lp.getObjVal() << std::endl;
+      for (typename Graph::EdgeIt e(g); e!=INVALID; ++e) 
+      flow.set(e, lp_flow[e]);
+    }
+  };
+
+} // namespace lemon
+
+#endif //LEMON_MIN_COST_GEN_FLOW_H