# HG changeset patch
# User marci
# Date 1092748846 0
# Node ID 615aca7091d264b0ad0e959e2584a1e439854ff2
# Parent  151b5754c7c6e60e234d7f1a7c1399da51b80448
An experimental LPSolverWrapper class which uses glpk. For a short
demo, max flow problems are solved with it. This demo does not
demonstrates, but the main aims of this class are row and column
generation capabilities, i.e. to be a core for easily
implementable branch-and-cut a column generetion algorithms.

diff -r 151b5754c7c6 -r 615aca7091d2 src/work/marci/lp/lp_solver_wrapper.h
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/marci/lp/lp_solver_wrapper.h	Tue Aug 17 13:20:46 2004 +0000
@@ -0,0 +1,382 @@
+// -*- c++ -*-
+#ifndef HUGO_LP_SOLVER_WRAPPER_H
+#define HUGO_LP_SOLVER_WRAPPER
+
+// #include <stdio.h>
+#include <stdlib.h>
+// #include <stdio>
+//#include <stdlib>
+#include "glpk.h"
+
+#include <iostream>
+#include <vector>
+#include <string>
+#include <list>
+#include <memory>
+#include <utility>
+
+//#include <sage_graph.h>
+//#include <hugo/list_graph.h>
+//#include <hugo/graph_wrapper.h>
+#include <hugo/invalid.h>
+//#include <bfs_dfs.h>
+//#include <stp.h>
+//#include <hugo/max_flow.h>
+//#include <augmenting_flow.h>
+//#include <iter_map.h>
+
+using std::cout;
+using std::cin;
+using std::endl;
+
+namespace hugo {
+
+  /// \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 hugo 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 hugo 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 hugo wrapper for glpk.
+  /// Later other LP-solvers will be wrapped into hugo.
+  /// 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);
+    }
+    void setObjCoef(const ColIt& col_it, double obj_coef) { 
+      lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
+    }
+    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 RowIt& row_it, double obj_coef) { 
+//     lpx_set_obj_coef(lp, row_iter_map[row_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 hugo
+
+#endif //HUGO_LP_SOLVER_WRAPPER_H
diff -r 151b5754c7c6 -r 615aca7091d2 src/work/marci/lp/makefile
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/marci/lp/makefile	Tue Aug 17 13:20:46 2004 +0000
@@ -0,0 +1,73 @@
+#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
+GLPKROOT = /usr/local/glpk-4.4
+INCLUDEDIRS ?= -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  = -L$(GLPKROOT)/lib -lglpk
+
+BINARIES = 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 151b5754c7c6 -r 615aca7091d2 src/work/marci/lp/max_flow_by_lp.cc
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/src/work/marci/lp/max_flow_by_lp.cc	Tue Aug 17 13:20:46 2004 +0000
@@ -0,0 +1,233 @@
+// -*- c++ -*-
+#include <iostream>
+#include <fstream>
+
+#include <sage_graph.h>
+#include <hugo/smart_graph.h>
+#include <hugo/dimacs.h>
+#include <hugo/time_measure.h>
+//#include <graph_wrapper.h>
+#include <hugo/max_flow.h>
+#include <augmenting_flow.h>
+//#include <preflow_res.h>
+#include <for_each_macros.h>
+#include <lp_solver_wrapper.h>
+
+using namespace hugo;
+
+// Use a DIMACS max flow file as stdin.
+// max_flow_demo < dimacs_max_flow_file
+
+template<typename Edge, typename EdgeIndexMap> 
+class PrimalMap {
+protected:
+  LPSolverWrapper* lp;
+  EdgeIndexMap* edge_index_map;
+public:
+  PrimalMap(LPSolverWrapper& _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]);
+  }
+};
+
+int main(int, char **) {
+
+  typedef SageGraph MutableGraph;
+  //typedef SmartGraph Graph;
+  typedef SageGraph Graph;
+  typedef Graph::Node Node;
+  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
+  MaxFlow<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.actMinCut(cut);
+
+    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
+      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.tail(e)] && cut[g.head(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(MaxFlow<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.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.tail(e)] && cut[g.head(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_EACH_LOC(Graph::EdgeIt, e, g) 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_EACH_LOC(Graph::EdgeIt, e, g) {
+      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.tail(e)] && cut[g.head(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_EACH_LOC(Graph::EdgeIt, e, g) 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_EACH_LOC(Graph::EdgeIt, e, g) {
+      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
+	std::cout << "Slackness does not hold!" << std::endl;
+      if (!cut[g.tail(e)] && cut[g.head(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.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.tail(e)] && cut[g.head(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.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
+// 	std::cout << "Slackness does not hold!" << std::endl;
+//     }
+//   }
+
+  ts.reset();
+  LPSolverWrapper lp;
+  lp.setMaximize();
+  typedef LPSolverWrapper::ColIt ColIt;
+  typedef LPSolverWrapper::RowIt RowIt;
+  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
+  EdgeIndexMap edge_index_map(g);
+  PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
+  Graph::EdgeIt e;
+  for (g.first(e); g.valid(e); g.next(e)) {
+    ColIt col_it=lp.addCol();
+    edge_index_map.set(e, col_it);
+    lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
+  }
+  Graph::NodeIt n;
+  for (g.first(n); g.valid(n); g.next(n)) {
+    if (n!=s) {
+      //hurokelek miatt
+      Graph::EdgeMap<int> coeffs(g, 0);
+      {
+	Graph::InEdgeIt e;
+	for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]+1);
+      }
+      {
+	Graph::OutEdgeIt e;
+	for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]-1);
+      }
+      if (n==t) {
+	Graph::EdgeIt e;
+	//std::vector< std::pair<ColIt, double> > row;
+	for (g.first(e); g.valid(e); g.next(e)) {
+	  if (coeffs[e]!=0) 
+	    lp.setObjCoef(edge_index_map[e], coeffs[e]);
+	}
+      } else  {
+	RowIt row_it=lp.addRow();
+	Graph::EdgeIt e;
+	std::vector< std::pair<ColIt, double> > row;
+	for (g.first(e); g.valid(e); g.next(e)) {
+	  if (coeffs[e]!=0) 
+	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
+	}	
+	lp.setRowCoeffs(row_it, row.begin(), row.end());
+	lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
+      }
+    }
+  }
+  lp.solveSimplex();
+  std::cout << "flow value: "<< lp.getObjVal() << std::endl;
+  std::cout << "elapsed time: " << ts << std::endl;
+
+  return 0;
+}