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.
1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/marci/lp/lp_solver_wrapper.h Tue Aug 17 13:20:46 2004 +0000
1.3 @@ -0,0 +1,382 @@
1.4 +// -*- c++ -*-
1.5 +#ifndef HUGO_LP_SOLVER_WRAPPER_H
1.6 +#define HUGO_LP_SOLVER_WRAPPER
1.7 +
1.8 +// #include <stdio.h>
1.9 +#include <stdlib.h>
1.10 +// #include <stdio>
1.11 +//#include <stdlib>
1.12 +#include "glpk.h"
1.13 +
1.14 +#include <iostream>
1.15 +#include <vector>
1.16 +#include <string>
1.17 +#include <list>
1.18 +#include <memory>
1.19 +#include <utility>
1.20 +
1.21 +//#include <sage_graph.h>
1.22 +//#include <hugo/list_graph.h>
1.23 +//#include <hugo/graph_wrapper.h>
1.24 +#include <hugo/invalid.h>
1.25 +//#include <bfs_dfs.h>
1.26 +//#include <stp.h>
1.27 +//#include <hugo/max_flow.h>
1.28 +//#include <augmenting_flow.h>
1.29 +//#include <iter_map.h>
1.30 +
1.31 +using std::cout;
1.32 +using std::cin;
1.33 +using std::endl;
1.34 +
1.35 +namespace hugo {
1.36 +
1.37 + /// \brief A partitioned vector with iterable classes.
1.38 + ///
1.39 + /// This class implements a container in which the data is stored in an
1.40 + /// stl vector, the range is partitioned into sets and each set is
1.41 + /// doubly linked in a list.
1.42 + /// That is, each class is iterable by hugo iterators, and any member of
1.43 + /// the vector can bo moved to an other class.
1.44 + template <typename T>
1.45 + class IterablePartition {
1.46 + protected:
1.47 + struct Node {
1.48 + T data;
1.49 + int prev; //invalid az -1
1.50 + int next;
1.51 + };
1.52 + std::vector<Node> nodes;
1.53 + struct Tip {
1.54 + int first;
1.55 + int last;
1.56 + };
1.57 + std::vector<Tip> tips;
1.58 + public:
1.59 + /// The classes are indexed by integers from \c 0 to \c classNum()-1.
1.60 + int classNum() const { return tips.size(); }
1.61 + /// This hugo style iterator iterates through a class.
1.62 + class ClassIt;
1.63 + /// Constructor. The number of classes is to be given which is fixed
1.64 + /// over the life of the container.
1.65 + /// The partition classes are indexed from 0 to class_num-1.
1.66 + IterablePartition(int class_num) {
1.67 + for (int i=0; i<class_num; ++i) {
1.68 + Tip t;
1.69 + t.first=t.last=-1;
1.70 + tips.push_back(t);
1.71 + }
1.72 + }
1.73 + protected:
1.74 + void befuz(ClassIt it, int class_id) {
1.75 + if (tips[class_id].first==-1) {
1.76 + if (tips[class_id].last==-1) {
1.77 + nodes[it.i].prev=nodes[it.i].next=-1;
1.78 + tips[class_id].first=tips[class_id].last=it.i;
1.79 + }
1.80 + } else {
1.81 + nodes[it.i].prev=tips[class_id].last;
1.82 + nodes[it.i].next=-1;
1.83 + nodes[tips[class_id].last].next=it.i;
1.84 + tips[class_id].last=it.i;
1.85 + }
1.86 + }
1.87 + void kifuz(ClassIt it, int class_id) {
1.88 + if (tips[class_id].first==it.i) {
1.89 + if (tips[class_id].last==it.i) {
1.90 + tips[class_id].first=tips[class_id].last=-1;
1.91 + } else {
1.92 + tips[class_id].first=nodes[it.i].next;
1.93 + nodes[nodes[it.i].next].prev=-1;
1.94 + }
1.95 + } else {
1.96 + if (tips[class_id].last==it.i) {
1.97 + tips[class_id].last=nodes[it.i].prev;
1.98 + nodes[nodes[it.i].prev].next=-1;
1.99 + } else {
1.100 + nodes[nodes[it.i].next].prev=nodes[it.i].prev;
1.101 + nodes[nodes[it.i].prev].next=nodes[it.i].next;
1.102 + }
1.103 + }
1.104 + }
1.105 + public:
1.106 + /// A new element with data \c t is pushed into the vector and into class
1.107 + /// \c class_id.
1.108 + ClassIt push_back(const T& t, int class_id) {
1.109 + Node n;
1.110 + n.data=t;
1.111 + nodes.push_back(n);
1.112 + int i=nodes.size()-1;
1.113 + befuz(i, class_id);
1.114 + return i;
1.115 + }
1.116 + /// A member is moved to an other class.
1.117 + void set(ClassIt it, int old_class_id, int new_class_id) {
1.118 + kifuz(it.i, old_class_id);
1.119 + befuz(it.i, new_class_id);
1.120 + }
1.121 + /// Returns the data pointed by \c it.
1.122 + T& operator[](ClassIt it) { return nodes[it.i].data; }
1.123 + /// Returns the data pointed by \c it.
1.124 + const T& operator[](ClassIt it) const { return nodes[it.i].data; }
1.125 + class ClassIt {
1.126 + friend class IterablePartition;
1.127 + protected:
1.128 + int i;
1.129 + public:
1.130 + /// Default constructor.
1.131 + ClassIt() { }
1.132 + /// This constructor constructs an iterator which points
1.133 + /// to the member of th container indexed by the integer _i.
1.134 + ClassIt(const int& _i) : i(_i) { }
1.135 + /// Invalid constructor.
1.136 + ClassIt(const Invalid&) : i(-1) { }
1.137 + };
1.138 + /// First member of class \c class_id.
1.139 + ClassIt& first(ClassIt& it, int class_id) const {
1.140 + it.i=tips[class_id].first;
1.141 + return it;
1.142 + }
1.143 + /// Next member.
1.144 + ClassIt& next(ClassIt& it) const {
1.145 + it.i=nodes[it.i].next;
1.146 + return it;
1.147 + }
1.148 + /// True iff the iterator is valid.
1.149 + bool valid(const ClassIt& it) const { return it.i!=-1; }
1.150 + };
1.151 +
1.152 + /// \brief Wrappers for LP solvers
1.153 + ///
1.154 + /// This class implements a hugo wrapper for glpk.
1.155 + /// Later other LP-solvers will be wrapped into hugo.
1.156 + /// The aim of this class is to give a general surface to different
1.157 + /// solvers, i.e. it makes possible to write algorithms using LP's,
1.158 + /// in which the solver can be changed to an other one easily.
1.159 + class LPSolverWrapper {
1.160 + public:
1.161 +
1.162 +// class Row {
1.163 +// protected:
1.164 +// int i;
1.165 +// public:
1.166 +// Row() { }
1.167 +// Row(const Invalid&) : i(0) { }
1.168 +// Row(const int& _i) : i(_i) { }
1.169 +// operator int() const { return i; }
1.170 +// };
1.171 +// class RowIt : public Row {
1.172 +// public:
1.173 +// RowIt(const Row& row) : Row(row) { }
1.174 +// };
1.175 +
1.176 +// class Col {
1.177 +// protected:
1.178 +// int i;
1.179 +// public:
1.180 +// Col() { }
1.181 +// Col(const Invalid&) : i(0) { }
1.182 +// Col(const int& _i) : i(_i) { }
1.183 +// operator int() const { return i; }
1.184 +// };
1.185 +// class ColIt : public Col {
1.186 +// ColIt(const Col& col) : Col(col) { }
1.187 +// };
1.188 +
1.189 + public:
1.190 + LPX* lp;
1.191 + typedef IterablePartition<int>::ClassIt RowIt;
1.192 + IterablePartition<int> row_iter_map;
1.193 + typedef IterablePartition<int>::ClassIt ColIt;
1.194 + IterablePartition<int> col_iter_map;
1.195 + //std::vector<int> row_id_to_lp_row_id;
1.196 + //std::vector<int> col_id_to_lp_col_id;
1.197 + const int VALID_ID;
1.198 + const int INVALID_ID;
1.199 +
1.200 + public:
1.201 + LPSolverWrapper() : lp(lpx_create_prob()),
1.202 + row_iter_map(2),
1.203 + col_iter_map(2),
1.204 + //row_id_to_lp_row_id(), col_id_to_lp_col_id(),
1.205 + VALID_ID(0), INVALID_ID(1) {
1.206 + lpx_set_int_parm(lp, LPX_K_DUAL, 1);
1.207 + }
1.208 + ~LPSolverWrapper() {
1.209 + lpx_delete_prob(lp);
1.210 + }
1.211 + void setMinimize() {
1.212 + lpx_set_obj_dir(lp, LPX_MIN);
1.213 + }
1.214 + void setMaximize() {
1.215 + lpx_set_obj_dir(lp, LPX_MAX);
1.216 + }
1.217 + ColIt addCol() {
1.218 + int i=lpx_add_cols(lp, 1);
1.219 + ColIt col_it;
1.220 + col_iter_map.first(col_it, INVALID_ID);
1.221 + if (col_iter_map.valid(col_it)) { //van hasznalhato hely
1.222 + col_iter_map.set(col_it, INVALID_ID, VALID_ID);
1.223 + col_iter_map[col_it]=i;
1.224 + //col_id_to_lp_col_id[col_iter_map[col_it]]=i;
1.225 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
1.226 + //col_id_to_lp_col_id.push_back(i);
1.227 + //int j=col_id_to_lp_col_id.size()-1;
1.228 + col_it=col_iter_map.push_back(i, VALID_ID);
1.229 + }
1.230 +// edge_index_map.set(e, i);
1.231 +// lpx_set_col_bnds(lp, i, LPX_DB, 0.0, 1.0);
1.232 +// lpx_set_obj_coef(lp, i, cost[e]);
1.233 + return col_it;
1.234 + }
1.235 + RowIt addRow() {
1.236 + int i=lpx_add_rows(lp, 1);
1.237 + RowIt row_it;
1.238 + row_iter_map.first(row_it, INVALID_ID);
1.239 + if (row_iter_map.valid(row_it)) { //van hasznalhato hely
1.240 + row_iter_map.set(row_it, INVALID_ID, VALID_ID);
1.241 + row_iter_map[row_it]=i;
1.242 + } else { //a cucc vegere kell inzertalni mert nincs szabad hely
1.243 + row_it=row_iter_map.push_back(i, VALID_ID);
1.244 + }
1.245 + return row_it;
1.246 + }
1.247 + //pair<RowIt, double>-bol kell megadni egy std range-et
1.248 + template <typename Begin, typename End>
1.249 + void setColCoeffs(const ColIt& col_it,
1.250 + Begin begin, End end) {
1.251 + int mem_length=1+lpx_get_num_rows(lp);
1.252 + int* indices = new int[mem_length];
1.253 + double* doubles = new double[mem_length];
1.254 + int length=0;
1.255 + for ( ; begin!=end; ++begin) {
1.256 + ++length;
1.257 + indices[length]=row_iter_map[begin->first];
1.258 + doubles[length]=begin->second;
1.259 + }
1.260 + lpx_set_mat_col(lp, col_iter_map[col_it], length, indices, doubles);
1.261 + delete [] indices;
1.262 + delete [] doubles;
1.263 + }
1.264 + //pair<ColIt, double>-bol kell megadni egy std range-et
1.265 + template <typename Begin, typename End>
1.266 + void setRowCoeffs(const RowIt& row_it,
1.267 + Begin begin, End end) {
1.268 + int mem_length=1+lpx_get_num_cols(lp);
1.269 + int* indices = new int[mem_length];
1.270 + double* doubles = new double[mem_length];
1.271 + int length=0;
1.272 + for ( ; begin!=end; ++begin) {
1.273 + ++length;
1.274 + indices[length]=col_iter_map[begin->first];
1.275 + doubles[length]=begin->second;
1.276 + }
1.277 + lpx_set_mat_row(lp, row_iter_map[row_it], length, indices, doubles);
1.278 + delete [] indices;
1.279 + delete [] doubles;
1.280 + }
1.281 + void eraseCol(const ColIt& col_it) {
1.282 + col_iter_map.set(col_it, VALID_ID, INVALID_ID);
1.283 + int cols[2];
1.284 + cols[1]=col_iter_map[col_it];
1.285 + lpx_del_cols(lp, 1, cols);
1.286 + col_iter_map[col_it]=0; //glpk specifikus
1.287 + ColIt it;
1.288 + for (col_iter_map.first(it, VALID_ID);
1.289 + col_iter_map.valid(it); col_iter_map.next(it)) {
1.290 + if (col_iter_map[it]>cols[1]) --col_iter_map[it];
1.291 + }
1.292 + }
1.293 + void eraseRow(const RowIt& row_it) {
1.294 + row_iter_map.set(row_it, VALID_ID, INVALID_ID);
1.295 + int rows[2];
1.296 + rows[1]=row_iter_map[row_it];
1.297 + lpx_del_rows(lp, 1, rows);
1.298 + row_iter_map[row_it]=0; //glpk specifikus
1.299 + RowIt it;
1.300 + for (row_iter_map.first(it, VALID_ID);
1.301 + row_iter_map.valid(it); row_iter_map.next(it)) {
1.302 + if (row_iter_map[it]>rows[1]) --row_iter_map[it];
1.303 + }
1.304 + }
1.305 + void setColBounds(const ColIt& col_it, int bound_type,
1.306 + double lo, double up) {
1.307 + lpx_set_col_bnds(lp, col_iter_map[col_it], bound_type, lo, up);
1.308 + }
1.309 + void setObjCoef(const ColIt& col_it, double obj_coef) {
1.310 + lpx_set_obj_coef(lp, col_iter_map[col_it], obj_coef);
1.311 + }
1.312 + void setRowBounds(const RowIt& row_it, int bound_type,
1.313 + double lo, double up) {
1.314 + lpx_set_row_bnds(lp, row_iter_map[row_it], bound_type, lo, up);
1.315 + }
1.316 +// void setObjCoef(const RowIt& row_it, double obj_coef) {
1.317 +// lpx_set_obj_coef(lp, row_iter_map[row_it], obj_coef);
1.318 +// }
1.319 + void solveSimplex() { lpx_simplex(lp); }
1.320 + void solvePrimalSimplex() { lpx_simplex(lp); }
1.321 + void solveDualSimplex() { lpx_simplex(lp); }
1.322 + double getPrimal(const ColIt& col_it) {
1.323 + return lpx_get_col_prim(lp, col_iter_map[col_it]);
1.324 + }
1.325 + double getObjVal() { return lpx_get_obj_val(lp); }
1.326 + int rowNum() const { return lpx_get_num_rows(lp); }
1.327 + int colNum() const { return lpx_get_num_cols(lp); }
1.328 + int warmUp() { return lpx_warm_up(lp); }
1.329 + void printWarmUpStatus(int i) {
1.330 + switch (i) {
1.331 + case LPX_E_OK: cout << "LPX_E_OK" << endl; break;
1.332 + case LPX_E_EMPTY: cout << "LPX_E_EMPTY" << endl; break;
1.333 + case LPX_E_BADB: cout << "LPX_E_BADB" << endl; break;
1.334 + case LPX_E_SING: cout << "LPX_E_SING" << endl; break;
1.335 + }
1.336 + }
1.337 + int getPrimalStatus() { return lpx_get_prim_stat(lp); }
1.338 + void printPrimalStatus(int i) {
1.339 + switch (i) {
1.340 + case LPX_P_UNDEF: cout << "LPX_P_UNDEF" << endl; break;
1.341 + case LPX_P_FEAS: cout << "LPX_P_FEAS" << endl; break;
1.342 + case LPX_P_INFEAS: cout << "LPX_P_INFEAS" << endl; break;
1.343 + case LPX_P_NOFEAS: cout << "LPX_P_NOFEAS" << endl; break;
1.344 + }
1.345 + }
1.346 + int getDualStatus() { return lpx_get_dual_stat(lp); }
1.347 + void printDualStatus(int i) {
1.348 + switch (i) {
1.349 + case LPX_D_UNDEF: cout << "LPX_D_UNDEF" << endl; break;
1.350 + case LPX_D_FEAS: cout << "LPX_D_FEAS" << endl; break;
1.351 + case LPX_D_INFEAS: cout << "LPX_D_INFEAS" << endl; break;
1.352 + case LPX_D_NOFEAS: cout << "LPX_D_NOFEAS" << endl; break;
1.353 + }
1.354 + }
1.355 + /// Returns the status of the slack variable assigned to row \c row_it.
1.356 + int getRowStat(const RowIt& row_it) {
1.357 + return lpx_get_row_stat(lp, row_iter_map[row_it]);
1.358 + }
1.359 + void printRowStatus(int i) {
1.360 + switch (i) {
1.361 + case LPX_BS: cout << "LPX_BS" << endl; break;
1.362 + case LPX_NL: cout << "LPX_NL" << endl; break;
1.363 + case LPX_NU: cout << "LPX_NU" << endl; break;
1.364 + case LPX_NF: cout << "LPX_NF" << endl; break;
1.365 + case LPX_NS: cout << "LPX_NS" << endl; break;
1.366 + }
1.367 + }
1.368 + /// Returns the status of the variable assigned to column \c col_it.
1.369 + int getColStat(const ColIt& col_it) {
1.370 + return lpx_get_col_stat(lp, col_iter_map[col_it]);
1.371 + }
1.372 + void printColStatus(int i) {
1.373 + switch (i) {
1.374 + case LPX_BS: cout << "LPX_BS" << endl; break;
1.375 + case LPX_NL: cout << "LPX_NL" << endl; break;
1.376 + case LPX_NU: cout << "LPX_NU" << endl; break;
1.377 + case LPX_NF: cout << "LPX_NF" << endl; break;
1.378 + case LPX_NS: cout << "LPX_NS" << endl; break;
1.379 + }
1.380 + }
1.381 + };
1.382 +
1.383 +} //namespace hugo
1.384 +
1.385 +#endif //HUGO_LP_SOLVER_WRAPPER_H
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/work/marci/lp/makefile Tue Aug 17 13:20:46 2004 +0000
2.3 @@ -0,0 +1,73 @@
2.4 +#INCLUDEDIRS ?= -I.. -I. -I./{marci,jacint,alpar,klao,akos}
2.5 +GLPKROOT = /usr/local/glpk-4.4
2.6 +INCLUDEDIRS ?= -I../../{marci,alpar,klao,akos,athos} -I. -I../../.. -I../.. -I.. -I$(GLPKROOT)/include
2.7 +#INCLUDEDIRS ?= -I../.. -I../.. -I../../{marci,jacint,alpar,klao,akos} -I/usr/local/glpk-4.4/include
2.8 +CXXFLAGS = -g -O2 -W -Wall $(INCLUDEDIRS) -ansi -pedantic
2.9 +LDFLAGS = -L$(GLPKROOT)/lib -lglpk
2.10 +
2.11 +BINARIES = max_flow_by_lp sample sample2 sample11 sample15
2.12 +
2.13 +#include ../makefile
2.14 +
2.15 +# Hat, ez elismerem, hogy nagyon ronda, de mukodik minden altalam
2.16 +# ismert rendszeren :-) (Misi)
2.17 +ifdef GCCVER
2.18 +CXX := g++-$(GCCVER)
2.19 +else
2.20 +CXX := $(shell type -p g++-3.3 || type -p g++-3.2 || type -p g++-3.0 || type -p g++-3 || echo g++)
2.21 +endif
2.22 +
2.23 +ifdef DEBUG
2.24 +CXXFLAGS += -DDEBUG
2.25 +endif
2.26 +
2.27 +CC := $(CXX)
2.28 +
2.29 +
2.30 +all: $(BINARIES)
2.31 +
2.32 +################
2.33 +# Minden binarishoz egy sor, felsorolva, hogy mely object file-okbol
2.34 +# all elo.
2.35 +# Kiveve ha siman file.cc -> file esetrol van szo, amikor is nem kell
2.36 +# irni semmit.
2.37 +
2.38 +#proba: proba.o seged.o
2.39 +
2.40 +################
2.41 +
2.42 +
2.43 +# .depend dep depend:
2.44 +# -$(CXX) $(CXXFLAGS) -M $(BINARIES:=.cc) > .depend
2.45 +
2.46 +#makefile: .depend
2.47 +#sinclude .depend
2.48 +#moert nem megy az eredeti /usr/bin/ld-vel?
2.49 +
2.50 +# %: %.o
2.51 +# $(CXX) -o $@ $< $(LDFLAGS)
2.52 +
2.53 +# %.o: %.cc
2.54 +# $(CXX) $(CXXFLAGS) -c $<
2.55 +
2.56 +%: %.cc
2.57 + $(CXX) $(CXXFLAGS) -o $@ $< $(LDFLAGS)
2.58 +
2.59 +sample11prof: sample11prof.o
2.60 + $(CXX) -pg -o sample11prof sample11prof.o -L$(GLPKROOT)/lib -lglpk
2.61 +sample11prof.o: sample11.cc
2.62 + $(CXX) -pg $(CXXFLAGS) -c -o sample11prof.o sample11.cc
2.63 +
2.64 +# sample.o: sample.cc
2.65 +# $(CXX) $(CXXFLAGS) -c -o sample.o sample.cc
2.66 +
2.67 +# sample2: sample2.o
2.68 +# $(CXX) -o sample2 sample2.o -L/usr/local/glpk-4.4/lib -lglpk
2.69 +# sample2.o: sample2.cc
2.70 +# $(CXX) $(CXXFLAGS) -c -o sample2.o sample2.cc
2.71 +
2.72 +
2.73 +clean:
2.74 + $(RM) *.o $(BINARIES) .depend
2.75 +
2.76 +.PHONY: all clean dep depend
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/src/work/marci/lp/max_flow_by_lp.cc Tue Aug 17 13:20:46 2004 +0000
3.3 @@ -0,0 +1,233 @@
3.4 +// -*- c++ -*-
3.5 +#include <iostream>
3.6 +#include <fstream>
3.7 +
3.8 +#include <sage_graph.h>
3.9 +#include <hugo/smart_graph.h>
3.10 +#include <hugo/dimacs.h>
3.11 +#include <hugo/time_measure.h>
3.12 +//#include <graph_wrapper.h>
3.13 +#include <hugo/max_flow.h>
3.14 +#include <augmenting_flow.h>
3.15 +//#include <preflow_res.h>
3.16 +#include <for_each_macros.h>
3.17 +#include <lp_solver_wrapper.h>
3.18 +
3.19 +using namespace hugo;
3.20 +
3.21 +// Use a DIMACS max flow file as stdin.
3.22 +// max_flow_demo < dimacs_max_flow_file
3.23 +
3.24 +template<typename Edge, typename EdgeIndexMap>
3.25 +class PrimalMap {
3.26 +protected:
3.27 + LPSolverWrapper* lp;
3.28 + EdgeIndexMap* edge_index_map;
3.29 +public:
3.30 + PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) :
3.31 + lp(&_lp), edge_index_map(&_edge_index_map) { }
3.32 + double operator[](Edge e) const {
3.33 + return lp->getPrimal((*edge_index_map)[e]);
3.34 + }
3.35 +};
3.36 +
3.37 +int main(int, char **) {
3.38 +
3.39 + typedef SageGraph MutableGraph;
3.40 + //typedef SmartGraph Graph;
3.41 + typedef SageGraph Graph;
3.42 + typedef Graph::Node Node;
3.43 + typedef Graph::EdgeIt EdgeIt;
3.44 +
3.45 + Graph g;
3.46 +
3.47 + Node s, t;
3.48 + Graph::EdgeMap<int> cap(g);
3.49 + //readDimacsMaxFlow(std::cin, g, s, t, cap);
3.50 + readDimacs(std::cin, g, cap, s, t);
3.51 + Timer ts;
3.52 + Graph::EdgeMap<int> flow(g); //0 flow
3.53 + MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
3.54 + max_flow_test(g, s, t, cap, flow);
3.55 + AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
3.56 + augmenting_flow_test(g, s, t, cap, flow);
3.57 +
3.58 + Graph::NodeMap<bool> cut(g);
3.59 +
3.60 + {
3.61 + std::cout << "preflow ..." << std::endl;
3.62 + ts.reset();
3.63 + max_flow_test.run();
3.64 + std::cout << "elapsed time: " << ts << std::endl;
3.65 + std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
3.66 + max_flow_test.actMinCut(cut);
3.67 +
3.68 + FOR_EACH_LOC(Graph::EdgeIt, e, g) {
3.69 + if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
3.70 + std::cout << "Slackness does not hold!" << std::endl;
3.71 + if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
3.72 + std::cout << "Slackness does not hold!" << std::endl;
3.73 + }
3.74 + }
3.75 +
3.76 +// {
3.77 +// std::cout << "preflow ..." << std::endl;
3.78 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.79 +// ts.reset();
3.80 +// max_flow_test.preflow(MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
3.81 +// std::cout << "elapsed time: " << ts << std::endl;
3.82 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
3.83 +
3.84 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
3.85 +// if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
3.86 +// std::cout << "Slackness does not hold!" << std::endl;
3.87 +// if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
3.88 +// std::cout << "Slackness does not hold!" << std::endl;
3.89 +// }
3.90 +// }
3.91 +
3.92 +// {
3.93 +// std::cout << "wrapped preflow ..." << std::endl;
3.94 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.95 +// ts.reset();
3.96 +// pre_flow_res.run();
3.97 +// std::cout << "elapsed time: " << ts << std::endl;
3.98 +// std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
3.99 +// }
3.100 +
3.101 + {
3.102 + std::cout << "physical blocking flow augmentation ..." << std::endl;
3.103 + FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.104 + ts.reset();
3.105 + int i=0;
3.106 + while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
3.107 + std::cout << "elapsed time: " << ts << std::endl;
3.108 + std::cout << "number of augmentation phases: " << i << std::endl;
3.109 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
3.110 +
3.111 + FOR_EACH_LOC(Graph::EdgeIt, e, g) {
3.112 + if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
3.113 + std::cout << "Slackness does not hold!" << std::endl;
3.114 + if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
3.115 + std::cout << "Slackness does not hold!" << std::endl;
3.116 + }
3.117 + }
3.118 +
3.119 +// {
3.120 +// std::cout << "faster physical blocking flow augmentation ..." << std::endl;
3.121 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.122 +// ts.reset();
3.123 +// int i=0;
3.124 +// while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
3.125 +// std::cout << "elapsed time: " << ts << std::endl;
3.126 +// std::cout << "number of augmentation phases: " << i << std::endl;
3.127 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
3.128 +// }
3.129 +
3.130 + {
3.131 + std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
3.132 + FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.133 + ts.reset();
3.134 + int i=0;
3.135 + while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
3.136 + std::cout << "elapsed time: " << ts << std::endl;
3.137 + std::cout << "number of augmentation phases: " << i << std::endl;
3.138 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
3.139 +
3.140 + FOR_EACH_LOC(Graph::EdgeIt, e, g) {
3.141 + if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
3.142 + std::cout << "Slackness does not hold!" << std::endl;
3.143 + if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
3.144 + std::cout << "Slackness does not hold!" << std::endl;
3.145 + }
3.146 + }
3.147 +
3.148 +// {
3.149 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
3.150 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.151 +// ts.reset();
3.152 +// int i=0;
3.153 +// while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
3.154 +// std::cout << "elapsed time: " << ts << std::endl;
3.155 +// std::cout << "number of augmentation phases: " << i << std::endl;
3.156 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
3.157 +
3.158 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
3.159 +// if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
3.160 +// std::cout << "Slackness does not hold!" << std::endl;
3.161 +// if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
3.162 +// std::cout << "Slackness does not hold!" << std::endl;
3.163 +// }
3.164 +// }
3.165 +
3.166 +// {
3.167 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
3.168 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
3.169 +// ts.reset();
3.170 +// int i=0;
3.171 +// while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
3.172 +// std::cout << "elapsed time: " << ts << std::endl;
3.173 +// std::cout << "number of augmentation phases: " << i << std::endl;
3.174 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
3.175 +
3.176 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
3.177 +// if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
3.178 +// std::cout << "Slackness does not hold!" << std::endl;
3.179 +// if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
3.180 +// std::cout << "Slackness does not hold!" << std::endl;
3.181 +// }
3.182 +// }
3.183 +
3.184 + ts.reset();
3.185 + LPSolverWrapper lp;
3.186 + lp.setMaximize();
3.187 + typedef LPSolverWrapper::ColIt ColIt;
3.188 + typedef LPSolverWrapper::RowIt RowIt;
3.189 + typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
3.190 + EdgeIndexMap edge_index_map(g);
3.191 + PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
3.192 + Graph::EdgeIt e;
3.193 + for (g.first(e); g.valid(e); g.next(e)) {
3.194 + ColIt col_it=lp.addCol();
3.195 + edge_index_map.set(e, col_it);
3.196 + lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
3.197 + }
3.198 + Graph::NodeIt n;
3.199 + for (g.first(n); g.valid(n); g.next(n)) {
3.200 + if (n!=s) {
3.201 + //hurokelek miatt
3.202 + Graph::EdgeMap<int> coeffs(g, 0);
3.203 + {
3.204 + Graph::InEdgeIt e;
3.205 + for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]+1);
3.206 + }
3.207 + {
3.208 + Graph::OutEdgeIt e;
3.209 + for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]-1);
3.210 + }
3.211 + if (n==t) {
3.212 + Graph::EdgeIt e;
3.213 + //std::vector< std::pair<ColIt, double> > row;
3.214 + for (g.first(e); g.valid(e); g.next(e)) {
3.215 + if (coeffs[e]!=0)
3.216 + lp.setObjCoef(edge_index_map[e], coeffs[e]);
3.217 + }
3.218 + } else {
3.219 + RowIt row_it=lp.addRow();
3.220 + Graph::EdgeIt e;
3.221 + std::vector< std::pair<ColIt, double> > row;
3.222 + for (g.first(e); g.valid(e); g.next(e)) {
3.223 + if (coeffs[e]!=0)
3.224 + row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
3.225 + }
3.226 + lp.setRowCoeffs(row_it, row.begin(), row.end());
3.227 + lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
3.228 + }
3.229 + }
3.230 + }
3.231 + lp.solveSimplex();
3.232 + std::cout << "flow value: "<< lp.getObjVal() << std::endl;
3.233 + std::cout << "elapsed time: " << ts << std::endl;
3.234 +
3.235 + return 0;
3.236 +}