1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/marci/lp/max_flow_by_lp.cc Tue Aug 17 13:20:46 2004 +0000
1.3 @@ -0,0 +1,233 @@
1.4 +// -*- c++ -*-
1.5 +#include <iostream>
1.6 +#include <fstream>
1.7 +
1.8 +#include <sage_graph.h>
1.9 +#include <hugo/smart_graph.h>
1.10 +#include <hugo/dimacs.h>
1.11 +#include <hugo/time_measure.h>
1.12 +//#include <graph_wrapper.h>
1.13 +#include <hugo/max_flow.h>
1.14 +#include <augmenting_flow.h>
1.15 +//#include <preflow_res.h>
1.16 +#include <for_each_macros.h>
1.17 +#include <lp_solver_wrapper.h>
1.18 +
1.19 +using namespace hugo;
1.20 +
1.21 +// Use a DIMACS max flow file as stdin.
1.22 +// max_flow_demo < dimacs_max_flow_file
1.23 +
1.24 +template<typename Edge, typename EdgeIndexMap>
1.25 +class PrimalMap {
1.26 +protected:
1.27 + LPSolverWrapper* lp;
1.28 + EdgeIndexMap* edge_index_map;
1.29 +public:
1.30 + PrimalMap(LPSolverWrapper& _lp, EdgeIndexMap& _edge_index_map) :
1.31 + lp(&_lp), edge_index_map(&_edge_index_map) { }
1.32 + double operator[](Edge e) const {
1.33 + return lp->getPrimal((*edge_index_map)[e]);
1.34 + }
1.35 +};
1.36 +
1.37 +int main(int, char **) {
1.38 +
1.39 + typedef SageGraph MutableGraph;
1.40 + //typedef SmartGraph Graph;
1.41 + typedef SageGraph Graph;
1.42 + typedef Graph::Node Node;
1.43 + typedef Graph::EdgeIt EdgeIt;
1.44 +
1.45 + Graph g;
1.46 +
1.47 + Node s, t;
1.48 + Graph::EdgeMap<int> cap(g);
1.49 + //readDimacsMaxFlow(std::cin, g, s, t, cap);
1.50 + readDimacs(std::cin, g, cap, s, t);
1.51 + Timer ts;
1.52 + Graph::EdgeMap<int> flow(g); //0 flow
1.53 + MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
1.54 + max_flow_test(g, s, t, cap, flow);
1.55 + AugmentingFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >
1.56 + augmenting_flow_test(g, s, t, cap, flow);
1.57 +
1.58 + Graph::NodeMap<bool> cut(g);
1.59 +
1.60 + {
1.61 + std::cout << "preflow ..." << std::endl;
1.62 + ts.reset();
1.63 + max_flow_test.run();
1.64 + std::cout << "elapsed time: " << ts << std::endl;
1.65 + std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
1.66 + max_flow_test.actMinCut(cut);
1.67 +
1.68 + FOR_EACH_LOC(Graph::EdgeIt, e, g) {
1.69 + if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
1.70 + std::cout << "Slackness does not hold!" << std::endl;
1.71 + if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
1.72 + std::cout << "Slackness does not hold!" << std::endl;
1.73 + }
1.74 + }
1.75 +
1.76 +// {
1.77 +// std::cout << "preflow ..." << std::endl;
1.78 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.79 +// ts.reset();
1.80 +// max_flow_test.preflow(MaxFlow<Graph, int, Graph::EdgeMap<int>, Graph::EdgeMap<int> >::GEN_FLOW);
1.81 +// std::cout << "elapsed time: " << ts << std::endl;
1.82 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
1.83 +
1.84 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
1.85 +// if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
1.86 +// std::cout << "Slackness does not hold!" << std::endl;
1.87 +// if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
1.88 +// std::cout << "Slackness does not hold!" << std::endl;
1.89 +// }
1.90 +// }
1.91 +
1.92 +// {
1.93 +// std::cout << "wrapped preflow ..." << std::endl;
1.94 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.95 +// ts.reset();
1.96 +// pre_flow_res.run();
1.97 +// std::cout << "elapsed time: " << ts << std::endl;
1.98 +// std::cout << "flow value: "<< pre_flow_test.flowValue() << std::endl;
1.99 +// }
1.100 +
1.101 + {
1.102 + std::cout << "physical blocking flow augmentation ..." << std::endl;
1.103 + FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.104 + ts.reset();
1.105 + int i=0;
1.106 + while (augmenting_flow_test.augmentOnBlockingFlow<MutableGraph>()) { ++i; }
1.107 + std::cout << "elapsed time: " << ts << std::endl;
1.108 + std::cout << "number of augmentation phases: " << i << std::endl;
1.109 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
1.110 +
1.111 + FOR_EACH_LOC(Graph::EdgeIt, e, g) {
1.112 + if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
1.113 + std::cout << "Slackness does not hold!" << std::endl;
1.114 + if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
1.115 + std::cout << "Slackness does not hold!" << std::endl;
1.116 + }
1.117 + }
1.118 +
1.119 +// {
1.120 +// std::cout << "faster physical blocking flow augmentation ..." << std::endl;
1.121 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.122 +// ts.reset();
1.123 +// int i=0;
1.124 +// while (max_flow_test.augmentOnBlockingFlow1<MutableGraph>()) { ++i; }
1.125 +// std::cout << "elapsed time: " << ts << std::endl;
1.126 +// std::cout << "number of augmentation phases: " << i << std::endl;
1.127 +// std::cout << "flow value: "<< max_flow_test.flowValue() << std::endl;
1.128 +// }
1.129 +
1.130 + {
1.131 + std::cout << "on-the-fly blocking flow augmentation ..." << std::endl;
1.132 + FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.133 + ts.reset();
1.134 + int i=0;
1.135 + while (augmenting_flow_test.augmentOnBlockingFlow2()) { ++i; }
1.136 + std::cout << "elapsed time: " << ts << std::endl;
1.137 + std::cout << "number of augmentation phases: " << i << std::endl;
1.138 + std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
1.139 +
1.140 + FOR_EACH_LOC(Graph::EdgeIt, e, g) {
1.141 + if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
1.142 + std::cout << "Slackness does not hold!" << std::endl;
1.143 + if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
1.144 + std::cout << "Slackness does not hold!" << std::endl;
1.145 + }
1.146 + }
1.147 +
1.148 +// {
1.149 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
1.150 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.151 +// ts.reset();
1.152 +// int i=0;
1.153 +// while (augmenting_flow_test.augmentOnShortestPath()) { ++i; }
1.154 +// std::cout << "elapsed time: " << ts << std::endl;
1.155 +// std::cout << "number of augmentation phases: " << i << std::endl;
1.156 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
1.157 +
1.158 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
1.159 +// if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
1.160 +// std::cout << "Slackness does not hold!" << std::endl;
1.161 +// if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
1.162 +// std::cout << "Slackness does not hold!" << std::endl;
1.163 +// }
1.164 +// }
1.165 +
1.166 +// {
1.167 +// std::cout << "on-the-fly shortest path augmentation ..." << std::endl;
1.168 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) flow.set(e, 0);
1.169 +// ts.reset();
1.170 +// int i=0;
1.171 +// while (augmenting_flow_test.augmentOnShortestPath2()) { ++i; }
1.172 +// std::cout << "elapsed time: " << ts << std::endl;
1.173 +// std::cout << "number of augmentation phases: " << i << std::endl;
1.174 +// std::cout << "flow value: "<< augmenting_flow_test.flowValue() << std::endl;
1.175 +
1.176 +// FOR_EACH_LOC(Graph::EdgeIt, e, g) {
1.177 +// if (cut[g.tail(e)] && !cut[g.head(e)] && !flow[e]==cap[e])
1.178 +// std::cout << "Slackness does not hold!" << std::endl;
1.179 +// if (!cut[g.tail(e)] && cut[g.head(e)] && flow[e]>0)
1.180 +// std::cout << "Slackness does not hold!" << std::endl;
1.181 +// }
1.182 +// }
1.183 +
1.184 + ts.reset();
1.185 + LPSolverWrapper lp;
1.186 + lp.setMaximize();
1.187 + typedef LPSolverWrapper::ColIt ColIt;
1.188 + typedef LPSolverWrapper::RowIt RowIt;
1.189 + typedef Graph::EdgeMap<ColIt> EdgeIndexMap;
1.190 + EdgeIndexMap edge_index_map(g);
1.191 + PrimalMap<Graph::Edge, EdgeIndexMap> lp_flow(lp, edge_index_map);
1.192 + Graph::EdgeIt e;
1.193 + for (g.first(e); g.valid(e); g.next(e)) {
1.194 + ColIt col_it=lp.addCol();
1.195 + edge_index_map.set(e, col_it);
1.196 + lp.setColBounds(col_it, LPX_DB, 0.0, cap[e]);
1.197 + }
1.198 + Graph::NodeIt n;
1.199 + for (g.first(n); g.valid(n); g.next(n)) {
1.200 + if (n!=s) {
1.201 + //hurokelek miatt
1.202 + Graph::EdgeMap<int> coeffs(g, 0);
1.203 + {
1.204 + Graph::InEdgeIt e;
1.205 + for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]+1);
1.206 + }
1.207 + {
1.208 + Graph::OutEdgeIt e;
1.209 + for (g.first(e, n); g.valid(e); g.next(e)) coeffs.set(e, coeffs[e]-1);
1.210 + }
1.211 + if (n==t) {
1.212 + Graph::EdgeIt e;
1.213 + //std::vector< std::pair<ColIt, double> > row;
1.214 + for (g.first(e); g.valid(e); g.next(e)) {
1.215 + if (coeffs[e]!=0)
1.216 + lp.setObjCoef(edge_index_map[e], coeffs[e]);
1.217 + }
1.218 + } else {
1.219 + RowIt row_it=lp.addRow();
1.220 + Graph::EdgeIt e;
1.221 + std::vector< std::pair<ColIt, double> > row;
1.222 + for (g.first(e); g.valid(e); g.next(e)) {
1.223 + if (coeffs[e]!=0)
1.224 + row.push_back(std::make_pair(edge_index_map[e], coeffs[e]));
1.225 + }
1.226 + lp.setRowCoeffs(row_it, row.begin(), row.end());
1.227 + lp.setRowBounds(row_it, LPX_FX, 0.0, 0.0);
1.228 + }
1.229 + }
1.230 + }
1.231 + lp.solveSimplex();
1.232 + std::cout << "flow value: "<< lp.getObjVal() << std::endl;
1.233 + std::cout << "elapsed time: " << ts << std::endl;
1.234 +
1.235 + return 0;
1.236 +}