src/work/marci/lp/max_flow_by_lp.cc
author deba
Mon, 04 Oct 2004 17:13:21 +0000
changeset 937 d4e911acef3d
parent 764 615aca7091d2
child 986 e997802b855c
permissions -rw-r--r--
Revert backport changes -r1230.
marci@764
     1
// -*- c++ -*-
marci@764
     2
#include <iostream>
marci@764
     3
#include <fstream>
marci@764
     4
marci@764
     5
#include <sage_graph.h>
alpar@921
     6
#include <lemon/smart_graph.h>
alpar@921
     7
#include <lemon/dimacs.h>
alpar@921
     8
#include <lemon/time_measure.h>
marci@764
     9
//#include <graph_wrapper.h>
alpar@921
    10
#include <lemon/max_flow.h>
marci@764
    11
#include <augmenting_flow.h>
marci@764
    12
//#include <preflow_res.h>
marci@764
    13
#include <for_each_macros.h>
marci@764
    14
#include <lp_solver_wrapper.h>
marci@764
    15
alpar@921
    16
using namespace lemon;
marci@764
    17
marci@764
    18
// Use a DIMACS max flow file as stdin.
marci@764
    19
// max_flow_demo < dimacs_max_flow_file
marci@764
    20
marci@764
    21
template<typename Edge, typename EdgeIndexMap> 
marci@764
    22
class PrimalMap {
marci@764
    23
protected:
marci@764
    24
  LPSolverWrapper* lp;
marci@764
    25
  EdgeIndexMap* edge_index_map;
marci@764
    26
public:
marci@764
    27
  PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) : 
marci@764
    28
    lp(&_lp), edge_index_map(&_edge_index_map) { }
marci@764
    29
  double operator[](Edge e) const { 
marci@764
    30
    return lp->getPrimal((*edge_index_map)[e]);
marci@764
    31
  }
marci@764
    32
};
marci@764
    33
marci@764
    34
int main(int, char **) {
marci@764
    35
marci@764
    36
  typedef SageGraph MutableGraph;
marci@764
    37
  //typedef SmartGraph Graph;
marci@764
    38
  typedef SageGraph Graph;
marci@764
    39
  typedef Graph::Node Node;
marci@764
    40
  typedef Graph::EdgeIt EdgeIt;
marci@764
    41
marci@764
    42
  Graph g;
marci@764
    43
marci@764
    44
  Node s, t;
marci@764
    45
  Graph::EdgeMap<int> cap(g);
marci@764
    46
  //readDimacsMaxFlow(std::cin, g, s, t, cap);
marci@764
    47
  readDimacs(std::cin, g, cap, s, t);
marci@764
    48
  Timer ts;
marci@764
    49
  Graph::EdgeMap<int> flow(g); //0 flow
marci@764
    50
  MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
marci@764
    51
    max_flow_test(g, s, t, cap, flow);
marci@764
    52
  AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> > 
marci@764
    53
    augmenting_flow_test(g, s, t, cap, flow);
marci@764
    54
  
marci@764
    55
  Graph::NodeMap<bool> cut(g);
marci@764
    56
marci@764
    57
  {
marci@764
    58
    std::cout << "preflow ..." << std::endl;
marci@764
    59
    ts.reset();
marci@764
    60
    max_flow_test.run();
marci@764
    61
    std::cout << "elapsed time: " << ts << std::endl;
marci@764
    62
    std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
marci@764
    63
    max_flow_test.actMinCut(cut);
marci@764
    64
marci@764
    65
    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
marci@764
    66
      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
marci@764
    67
	std::cout << "Slackness does not hold!" << std::endl;
marci@764
    68
      if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
marci@764
    69
	std::cout << "Slackness does not hold!" << std::endl;
marci@764
    70
    }
marci@764
    71
  }
marci@764
    72
marci@764
    73
//   {
marci@764
    74
//     std::cout << "preflow ..." << std::endl;
marci@764
    75
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
    76
//     ts.reset();
marci@764
    77
//     max_flow_test.preflow(MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
marci@764
    78
//     std::cout << "elapsed time: " << ts << std::endl;
marci@764
    79
//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
marci@764
    80
marci@764
    81
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
marci@764
    82
//       if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
marci@764
    83
// 	std::cout << "Slackness does not hold!" << std::endl;
marci@764
    84
//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
marci@764
    85
// 	std::cout << "Slackness does not hold!" << std::endl;
marci@764
    86
//     }
marci@764
    87
//   }
marci@764
    88
marci@764
    89
//   {
marci@764
    90
//     std::cout << "wrapped preflow ..." << std::endl;
marci@764
    91
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
    92
//     ts.reset();
marci@764
    93
//     pre_flow_res.run();
marci@764
    94
//     std::cout << "elapsed time: " << ts << std::endl;
marci@764
    95
//     std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
marci@764
    96
//   }
marci@764
    97
marci@764
    98
  {
marci@764
    99
    std::cout << "physical blocking flow augmentation ..." << std::endl;
marci@764
   100
    FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
   101
    ts.reset();
marci@764
   102
    int i=0;
marci@764
   103
    while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
marci@764
   104
    std::cout << "elapsed time: " << ts << std::endl;
marci@764
   105
    std::cout << "number of augmentation phases: " << i << std::endl; 
marci@764
   106
    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
marci@764
   107
marci@764
   108
    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
marci@764
   109
      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
marci@764
   110
	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   111
      if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
marci@764
   112
	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   113
    }
marci@764
   114
  }
marci@764
   115
marci@764
   116
//   {
marci@764
   117
//     std::cout << "faster physical blocking flow augmentation ..." << std::endl;
marci@764
   118
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
   119
//     ts.reset();
marci@764
   120
//     int i=0;
marci@764
   121
//     while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
marci@764
   122
//     std::cout << "elapsed time: " << ts << std::endl;
marci@764
   123
//     std::cout << "number of augmentation phases: " << i << std::endl; 
marci@764
   124
//     std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
marci@764
   125
//   }
marci@764
   126
marci@764
   127
  {
marci@764
   128
    std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
marci@764
   129
    FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
   130
    ts.reset();
marci@764
   131
    int i=0;
marci@764
   132
    while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
marci@764
   133
    std::cout << "elapsed time: " << ts << std::endl;
marci@764
   134
    std::cout << "number of augmentation phases: " << i << std::endl; 
marci@764
   135
    std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
marci@764
   136
marci@764
   137
    FOR_EACH_LOC(Graph::EdgeIt, e, g) {
marci@764
   138
      if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
marci@764
   139
	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   140
      if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
marci@764
   141
	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   142
    }
marci@764
   143
  }
marci@764
   144
marci@764
   145
//   {
marci@764
   146
//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
marci@764
   147
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
   148
//     ts.reset();
marci@764
   149
//     int i=0;
marci@764
   150
//     while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
marci@764
   151
//     std::cout << "elapsed time: " << ts << std::endl;
marci@764
   152
//     std::cout << "number of augmentation phases: " << i << std::endl; 
marci@764
   153
//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
marci@764
   154
marci@764
   155
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
marci@764
   156
//       if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
marci@764
   157
// 	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   158
//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
marci@764
   159
// 	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   160
//     }
marci@764
   161
//   }
marci@764
   162
marci@764
   163
//   {
marci@764
   164
//     std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
marci@764
   165
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
marci@764
   166
//     ts.reset();
marci@764
   167
//     int i=0;
marci@764
   168
//     while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
marci@764
   169
//     std::cout << "elapsed time: " << ts << std::endl;
marci@764
   170
//     std::cout << "number of augmentation phases: " << i << std::endl; 
marci@764
   171
//     std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
marci@764
   172
marci@764
   173
//     FOR_EACH_LOC(Graph::EdgeIt, e, g) {
marci@764
   174
//       if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e]) 
marci@764
   175
// 	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   176
//       if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0) 
marci@764
   177
// 	std::cout << "Slackness does not hold!" << std::endl;
marci@764
   178
//     }
marci@764
   179
//   }
marci@764
   180
marci@764
   181
  ts.reset();
marci@764
   182
  LPSolverWrapper lp;
marci@764
   183
  lp.setMaximize();
marci@764
   184
  typedef LPSolverWrapper::ColIt ColIt;
marci@764
   185
  typedef LPSolverWrapper::RowIt RowIt;
marci@764
   186
  typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
marci@764
   187
  EdgeIndexMap edge_index_map(g);
marci@764
   188
  PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
marci@764
   189
  Graph::EdgeIt e;
marci@764
   190
  for (g.first(e); g.valid(e); g.next(e)) {
marci@764
   191
    ColIt col_it=lp.addCol();
marci@764
   192
    edge_index_map.set(e, col_it);
marci@764
   193
    lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
marci@764
   194
  }
marci@764
   195
  Graph::NodeIt n;
marci@764
   196
  for (g.first(n); g.valid(n); g.next(n)) {
marci@764
   197
    if (n!=s) {
marci@764
   198
      //hurokelek miatt
marci@764
   199
      Graph::EdgeMap<int> coeffs(g, 0);
marci@764
   200
      {
marci@764
   201
	Graph::InEdgeIt e;
marci@764
   202
	for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]+1);
marci@764
   203
      }
marci@764
   204
      {
marci@764
   205
	Graph::OutEdgeIt e;
marci@764
   206
	for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]-1);
marci@764
   207
      }
marci@764
   208
      if (n==t) {
marci@764
   209
	Graph::EdgeIt e;
marci@764
   210
	//std::vector< std::pair<ColIt, double> > row;
marci@764
   211
	for (g.first(e); g.valid(e); g.next(e)) {
marci@764
   212
	  if (coeffs[e]!=0) 
marci@764
   213
	    lp.setObjCoef(edge_index_map[e], coeffs[e]);
marci@764
   214
	}
marci@764
   215
      } else  {
marci@764
   216
	RowIt row_it=lp.addRow();
marci@764
   217
	Graph::EdgeIt e;
marci@764
   218
	std::vector< std::pair<ColIt, double> > row;
marci@764
   219
	for (g.first(e); g.valid(e); g.next(e)) {
marci@764
   220
	  if (coeffs[e]!=0) 
marci@764
   221
	    row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
marci@764
   222
	}	
marci@764
   223
	lp.setRowCoeffs(row_it, row.begin(), row.end());
marci@764
   224
	lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
marci@764
   225
      }
marci@764
   226
    }
marci@764
   227
  }
marci@764
   228
  lp.solveSimplex();
marci@764
   229
  std::cout << "flow value: "<< lp.getObjVal() << std::endl;
marci@764
   230
  std::cout << "elapsed time: " << ts << std::endl;
marci@764
   231
marci@764
   232
  return 0;
marci@764
   233
}