1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/work/jacint/dimacs_jgraph.hh Thu Feb 26 11:38:51 2004 +0000
1.3 @@ -0,0 +1,61 @@
1.4 +#ifndef DIMACS_HH
1.5 +#define DIMACS_HH
1.6 +
1.7 +#include <iostream>
1.8 +#include <string>
1.9 +#include <vector>
1.10 +
1.11 +namespace hugo {
1.12 +
1.13 + template<typename Graph, typename CapacityMap>
1.14 + void readDimacsMaxFlow(std::istream& is, Graph &G, typename Graph::TrivNodeIt &s, typename Graph::TrivNodeIt &t, CapacityMap& capacity) {
1.15 + G.clear();
1.16 + int cap;
1.17 + char d;
1.18 + std::string problem;
1.19 + char c;
1.20 + int i, j;
1.21 + std::string str;
1.22 + int n, m;
1.23 + std::vector<typename Graph::TrivNodeIt> nodes;
1.24 + while (is>>c) {
1.25 + switch (c) {
1.26 + case 'c': //comment
1.27 + getline(is, str);
1.28 + break;
1.29 + case 't': //type
1.30 + getline(is, str);
1.31 + break;
1.32 + case 'p': //problem definition
1.33 + is >> problem >> n >> m;
1.34 + getline(is, str);
1.35 + nodes.resize(n+1);
1.36 + for (int k=1; k<=n; ++k) nodes[k]=G.addNode();
1.37 + break;
1.38 + case 'n': //node definition
1.39 + if (problem=="sp") { //shortest path problem
1.40 + is >> i;
1.41 + getline(is, str);
1.42 + s=nodes[i];
1.43 + }
1.44 + if (problem=="max") { //max flow problem
1.45 + is >> i >> d;
1.46 + getline(is, str);
1.47 + if (d=='s') s=nodes[i];
1.48 + if (d=='t') t=nodes[i];
1.49 + }
1.50 + break;
1.51 + case 'a':
1.52 + is >> i >> j >> cap;
1.53 + getline(is, str);
1.54 + typename Graph::TrivEdgeIt e=G.addEdge(nodes[i], nodes[j]);
1.55 + capacity.update();
1.56 + capacity.set(e, cap);
1.57 + break;
1.58 + }
1.59 + }
1.60 + }
1.61 +
1.62 +} //namespace marci
1.63 +
1.64 +#endif //DIMACS_HH
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/src/work/jacint/j_graph.h Thu Feb 26 11:38:51 2004 +0000
2.3 @@ -0,0 +1,363 @@
2.4 +// -*- mode:C++ -*-
2.5 +
2.6 +#ifndef J_GRAPH_H
2.7 +#define J_GRAPH_H
2.8 +
2.9 +#include <iostream>
2.10 +#include <vector>
2.11 +
2.12 +namespace hugo {
2.13 +
2.14 + class JGraph {
2.15 +
2.16 + struct NodeT
2.17 + {
2.18 + int in, out;
2.19 + NodeT() : in(), out() {}
2.20 + };
2.21 +
2.22 + struct EdgeT
2.23 + {
2.24 + int head, tail, in, out;
2.25 + };
2.26 +
2.27 +
2.28 + std::vector<NodeT> nodes;
2.29 + std::vector<EdgeT> edges;
2.30 +
2.31 +
2.32 + /* template <typename Key> class DynMapBase
2.33 + {
2.34 + protected:
2.35 + SmartGraph* G;
2.36 + public:
2.37 + virtual void add(const Key k) = NULL;
2.38 + virtual void erase(const Key k) = NULL;
2.39 + DynMapBase(SmartGraph &_G) : G(&_G) {}
2.40 + virtual ~DynMapBase() {}
2.41 + friend class SmartGraph;
2.42 + };
2.43 +
2.44 + template <typename T> class DynEdgeMap;
2.45 + template <typename T> class DynEdgeMap;
2.46 + */
2.47 +
2.48 + public:
2.49 +
2.50 + class TrivNodeIt;
2.51 + class TrivEdgeIt;
2.52 + class NodeIt;
2.53 + class EdgeIt;
2.54 + class OutEdgeIt;
2.55 + class InEdgeIt;
2.56 +
2.57 + /* protected:
2.58 + std::vector<DynMapBase<NodeIt> * > dyn_node_maps;
2.59 + std::vector<DynMapBase<EdgeIt> * > dyn_edge_maps;
2.60 + */
2.61 +
2.62 + template <typename T> class NodeMap;
2.63 + template <typename T> class EdgeMap;
2.64 +
2.65 + public:
2.66 +
2.67 + JGraph() : nodes(1), edges(1) {
2.68 + edges[0].head=0;
2.69 + edges[0].tail=0;
2.70 + edges[0].in=0; //numEdges (numNodes is maintained in nodes[0].in)
2.71 + edges[0].out=0;
2.72 + }
2.73 +
2.74 +
2.75 + /* ~SmartGraph()
2.76 + {
2.77 + for(std::vector<DynMapBase<NodeIt> * >::iterator i=dyn_node_maps.begin();
2.78 + i!=dyn_node_maps.end(); ++i) (**i).G=NULL;
2.79 + for(std::vector<DynMapBase<EdgeIt> * >::iterator i=dyn_edge_maps.begin();
2.80 + i!=dyn_edge_maps.end(); ++i) (**i).G=NULL;
2.81 + }*/
2.82 +
2.83 + int numNodes() const { return nodes[0].in; }
2.84 + int numEdges() const { return edges[0].in; }
2.85 +
2.86 + TrivNodeIt tail(TrivEdgeIt e) const { return edges[e.n].tail; }
2.87 + TrivNodeIt head(TrivEdgeIt e) const { return edges[e.n].head; }
2.88 +
2.89 + /* EachNodeIt& getFirst(EachNodeIt& v) const {
2.90 + v=EachNodeIt(*this); return v; }
2.91 + EachEdgeIt& getFirst(EachEdgeIt& e) const {
2.92 + e=EachEdgeIt(*this); return e; }
2.93 + OutEdgeIt& getFirst(OutEdgeIt& e, const NodeIt v) const {
2.94 + e=OutEdgeIt(*this,v); return e; }
2.95 + InEdgeIt& getFirst(InEdgeIt& e, const NodeIt v) const {
2.96 + e=InEdgeIt(*this,v); return e; }
2.97 + */
2.98 +
2.99 + NodeIt firstNode() const {
2.100 + return NodeIt( numNodes() );
2.101 + }
2.102 + EdgeIt firstEdge() const {
2.103 + return EdgeIt( numEdges() );
2.104 + }
2.105 + InEdgeIt firstIn(TrivNodeIt v) const {
2.106 + return InEdgeIt(nodes[v.n].in);
2.107 + }
2.108 + OutEdgeIt firstOut(TrivNodeIt v) const {
2.109 + return OutEdgeIt(nodes[v.n].out);
2.110 + }
2.111 +
2.112 +
2.113 +
2.114 + void next(NodeIt& it) const {--it.n;}
2.115 + void next(OutEdgeIt& it) const { it.n=edges[it.n].out; }
2.116 + void next(InEdgeIt& it) const { it.n=edges[it.n].in; }
2.117 + void next(EdgeIt& it) const {--it.n;}
2.118 +
2.119 + int id(TrivNodeIt v) const { return v.n; }
2.120 + int id(TrivEdgeIt e) const { return e.n; }
2.121 +
2.122 + TrivNodeIt addNode() {
2.123 + TrivNodeIt n;
2.124 + n.n=++nodes[0].in;
2.125 + nodes.push_back(NodeT()); //FIXME
2.126 +
2.127 + /* for(std::vector<DynMapBase<NodeIt> * >::iterator i=dyn_node_maps.begin();
2.128 + i!=dyn_node_maps.end(); ++i) (**i).add(n.n);*/
2.129 +
2.130 + return n;
2.131 + }
2.132 +
2.133 +
2.134 + TrivEdgeIt addEdge(TrivNodeIt u, TrivNodeIt v) {
2.135 + if ( u && v ) {
2.136 + TrivEdgeIt e;
2.137 + e.n=++edges[0].in;
2.138 + edges.push_back(EdgeT()); //FIXME: Hmmm...
2.139 + edges[e.n].tail=u.n; edges[e.n].head=v.n;
2.140 + edges[e.n].out=nodes[u.n].out;
2.141 + edges[e.n].in=nodes[v.n].in;
2.142 + nodes[u.n].out=nodes[v.n].in=e.n;
2.143 + return e;
2.144 + }
2.145 + else return TrivEdgeIt();
2.146 +
2.147 + /* for(std::vector<DynMapBase<EdgeIt> * >::iterator i=dyn_edge_maps.begin();
2.148 + i!=dyn_edge_maps.end(); ++i) (**i).add(e.n);*/
2.149 +
2.150 +
2.151 + }
2.152 +
2.153 +
2.154 + void clear() {
2.155 + nodes.resize(1);
2.156 + edges.resize(1); edges[0].in=0;
2.157 + }
2.158 +
2.159 +
2.160 + class TrivNodeIt {
2.161 + friend class JGraph;
2.162 + template <typename T> friend class NodeMap;
2.163 +
2.164 + friend class TrivEdgeIt;
2.165 + friend class OutEdgeIt;
2.166 + friend class InEdgeIt;
2.167 +
2.168 + protected:
2.169 + int n;
2.170 + public:
2.171 + TrivNodeIt() : n() {}
2.172 + TrivNodeIt(int number) {n=number;}
2.173 + bool operator==(const TrivNodeIt i) const {return n==i.n;}
2.174 + bool operator!=(const TrivNodeIt i) const {return n!=i.n;}
2.175 +
2.176 + operator bool() const { return n!=0; }
2.177 +
2.178 + };
2.179 +
2.180 +
2.181 + class NodeIt : public TrivNodeIt {
2.182 + friend class JGraph;
2.183 + public:
2.184 + NodeIt() : TrivNodeIt() { }
2.185 + NodeIt(int i) : TrivNodeIt(i) { }
2.186 + operator bool() const { return n!=0; }
2.187 + };
2.188 +
2.189 +
2.190 + class TrivEdgeIt {
2.191 + friend class JGraph;
2.192 + template <typename T> friend class EdgeMap;
2.193 +
2.194 + friend class TrivNodeIt;
2.195 + friend class NodeIt;
2.196 + protected:
2.197 + int n;
2.198 + public:
2.199 + TrivEdgeIt() : n() { }
2.200 + TrivEdgeIt(int number) {n=number;}
2.201 + bool operator==(const TrivEdgeIt i) const {return n==i.n;}
2.202 + bool operator!=(const TrivEdgeIt i) const {return n!=i.n;}
2.203 + operator bool() const { return n!=0; }
2.204 +
2.205 + };
2.206 +
2.207 +
2.208 + class EdgeIt : public TrivEdgeIt {
2.209 + friend class JGraph;
2.210 + public:
2.211 + EdgeIt() : TrivEdgeIt() { }
2.212 + EdgeIt(int i) : TrivEdgeIt(i) { }
2.213 + operator bool() const { return n!=0; }
2.214 + };
2.215 +
2.216 +
2.217 + class OutEdgeIt : public TrivEdgeIt {
2.218 + friend class JGraph;
2.219 + public:
2.220 + OutEdgeIt() : TrivEdgeIt() { }
2.221 + OutEdgeIt(int i) : TrivEdgeIt(i) { }
2.222 + };
2.223 +
2.224 +
2.225 + class InEdgeIt : public TrivEdgeIt {
2.226 + friend class JGraph;
2.227 + public:
2.228 + InEdgeIt() : TrivEdgeIt() { }
2.229 + InEdgeIt(int i) : TrivEdgeIt(i) { }
2.230 + };
2.231 +
2.232 +
2.233 + // Map types
2.234 +
2.235 + template <typename T>
2.236 + class NodeMap {
2.237 + const JGraph& G;
2.238 + std::vector<T> container;
2.239 + public:
2.240 + NodeMap(const JGraph& _G) : G(_G), container( G.numNodes()+1 ) { }
2.241 + NodeMap(const JGraph& _G, T a) :
2.242 + G(_G), container(G.numNodes()+1, a) { }
2.243 + void set(TrivNodeIt n, T a) { container[n.n]=a; }
2.244 + T get(TrivNodeIt n) const { return container[n.n]; }
2.245 + /*T& operator[](NodeIt n) { return container[n.n]; }
2.246 + const T& operator[](NodeIt n) const { return container[n.n]; }*/
2.247 + void update() { container.resize(G.numNodes()+1); }
2.248 + void update(T a) { container.resize(G.numNodes()+1, a); }
2.249 + };
2.250 +
2.251 + template <typename T>
2.252 + class EdgeMap {
2.253 + const JGraph& G;
2.254 + std::vector<T> container;
2.255 + public:
2.256 + EdgeMap(const JGraph& _G) : G(_G), container(G.numEdges()+1) { }
2.257 + EdgeMap(const JGraph& _G, T a) :
2.258 + G(_G), container(G.numEdges()+1, a) { }
2.259 + void set(TrivEdgeIt e, T a) { container[e.n]=a; }
2.260 + T get(TrivEdgeIt e) const { return container[e.n]; }
2.261 + /*T& operator[](EdgeIt e) { return container[e.n]; }
2.262 + const T& operator[](EdgeIt e) const { return container[e.n]; } */
2.263 + void update() { container.resize(G.numEdges()+1); }
2.264 + void update(T a) { container.resize(G.numEdges()+1, a); }
2.265 + };
2.266 +
2.267 + /*template <typename T> class DynNodeMap : public DynMapBase<NodeIt>
2.268 + {
2.269 + std::vector<T> container;
2.270 +
2.271 + public:
2.272 + typedef T ValueType;
2.273 + typedef NodeIt KeyType;
2.274 +
2.275 + DynNodeMap(JGraph &_G) :
2.276 + DynMapBase<NodeIt>(_G), container(_G.maxNodeId())
2.277 + {
2.278 + //FIXME: What if there are empty Id's?
2.279 + //FIXME: Can I use 'this' in a constructor?
2.280 + G->dyn_node_maps.push_back(this);
2.281 + }
2.282 + ~DynNodeMap()
2.283 + {
2.284 + if(G) {
2.285 + std::vector<DynMapBase<NodeIt>* >::iterator i;
2.286 + for(i=G->dyn_node_maps.begin();
2.287 + i!=G->dyn_node_maps.end() && *i!=this; ++i) ;
2.288 + //if(*i==this) G->dyn_node_maps.erase(i); //FIXME: Way too slow...
2.289 + //A better way to do that: (Is this really important?)
2.290 + if(*i==this) {
2.291 + *i=G->dyn_node_maps.back();
2.292 + G->dyn_node_maps.pop_back();
2.293 + }
2.294 + }
2.295 + }
2.296 +
2.297 + void add(const NodeIt k)
2.298 + {
2.299 + if(k.n>=container.size()) container.resize(k.n+1);
2.300 + }
2.301 + void erase(const NodeIt k)
2.302 + {
2.303 + //FIXME: Please implement me.
2.304 + }
2.305 +
2.306 + void set(NodeIt n, T a) { container[n.n]=a; }
2.307 + T get(NodeIt n) const { return container[n.n]; }
2.308 + T& operator[](NodeIt n) { return container[n.n]; }
2.309 + const T& operator[](NodeIt n) const { return container[n.n]; }
2.310 +
2.311 + void update() {} //Useless for DynMaps
2.312 + void update(T a) {} //Useless for DynMaps
2.313 + };
2.314 +
2.315 + template <typename T> class DynEdgeMap : public DynMapBase<EdgeIt>
2.316 + {
2.317 + std::vector<T> container;
2.318 +
2.319 + public:
2.320 + typedef T ValueType;
2.321 + typedef EdgeIt KeyType;
2.322 +
2.323 + DynEdgeMap(JGraph &_G) :
2.324 + DynMapBase<EdgeIt>(_G), container(_G.maxEdgeId())
2.325 + {
2.326 + //FIXME: What if there are empty Id's?
2.327 + //FIXME: Can I use 'this' in a constructor?
2.328 + G->dyn_edge_maps.push_back(this);
2.329 + }
2.330 + ~DynEdgeMap()
2.331 + {
2.332 + if(G) {
2.333 + std::vector<DynMapBase<EdgeIt>* >::iterator i;
2.334 + for(i=G->dyn_edge_maps.begin();
2.335 + i!=G->dyn_edge_maps.end() && *i!=this; ++i) ;
2.336 + //if(*i==this) G->dyn_edge_maps.erase(i); //Way too slow...
2.337 + //A better way to do that: (Is this really important?)
2.338 + if(*i==this) {
2.339 + *i=G->dyn_edge_maps.back();
2.340 + G->dyn_edge_maps.pop_back();
2.341 + }
2.342 + }
2.343 + }
2.344 +
2.345 + void add(const EdgeIt k)
2.346 + {
2.347 + if(k.n>=int(container.size())) container.resize(k.n+1);
2.348 + }
2.349 + void erase(const EdgeIt k)
2.350 + {
2.351 + //FIXME: Please implement me.
2.352 + }
2.353 +
2.354 + void set(EdgeIt n, T a) { container[n.n]=a; }
2.355 + T get(EdgeIt n) const { return container[n.n]; }
2.356 + T& operator[](EdgeIt n) { return container[n.n]; }
2.357 + const T& operator[](EdgeIt n) const { return container[n.n]; }
2.358 +
2.359 + void update() {} //Useless for DynMaps
2.360 + void update(T a) {} //Useless for DynMaps
2.361 + };*/
2.362 +
2.363 + };
2.364 +} //namespace hugo
2.365 +
2.366 +#endif //J_GRAPH_H
3.1 --- a/src/work/jacint/makefile Wed Feb 25 15:27:17 2004 +0000
3.2 +++ b/src/work/jacint/makefile Thu Feb 26 11:38:51 2004 +0000
3.3 @@ -3,36 +3,44 @@
3.4 CXXFLAGS = -W -Wall -ansi -pedantic
3.5 LEDAROOT = /ledasrc/LEDA-4.1
3.6
3.7 -BINARIES = preflow preflow_param preflow_hl0 preflow_hl1 preflow_hl2 preflow_hl3 preflow_hl4 preflow_max_flow
3.8 +BINARIES = preflow_jgraph
3.9
3.10 all: $(BINARIES)
3.11
3.12 makefile: .depend
3.13 sinclude .depend
3.14
3.15 +
3.16 +
3.17 +preflow_jgraph:
3.18 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_jgraph preflow_jgraph.cc
3.19 +
3.20 +preflowalpar:
3.21 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../alpar -I../jacint -o preflowalpar preflowalpar.cc
3.22 +
3.23 preflow_max_flow:
3.24 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_max_flow preflow_max_flow.cc
3.25 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_max_flow preflow_max_flow.cc
3.26
3.27 preflow:
3.28 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow preflow.cc
3.29 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow preflow.cc
3.30
3.31 preflow_hl0:
3.32 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_hl0 preflow_hl0.cc
3.33 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_hl0 preflow_hl0.cc
3.34
3.35 preflow_hl1:
3.36 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_hl1 preflow_hl1.cc
3.37 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_hl1 preflow_hl1.cc
3.38
3.39 preflow_hl2:
3.40 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_hl2 preflow_hl2.cc
3.41 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_hl2 preflow_hl2.cc
3.42
3.43 preflow_hl3:
3.44 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_hl3 preflow_hl3.cc
3.45 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_hl3 preflow_hl3.cc
3.46
3.47 preflow_hl4:
3.48 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_hl4 preflow_hl4.cc
3.49 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_hl4 preflow_hl4.cc
3.50
3.51 preflow_param:
3.52 - $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../marci -I../jacint -o preflow_param preflow_param.cc
3.53 + $(CXX3) $(CXXFLAGS) -O3 -I. -I.. -I../jacint -o preflow_param preflow_param.cc
3.54
3.55 clean:
3.56 $(RM) *.o $(BINARIES) .depend
4.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
4.2 +++ b/src/work/jacint/preflow_jgraph.cc Thu Feb 26 11:38:51 2004 +0000
4.3 @@ -0,0 +1,72 @@
4.4 +#include <iostream>
4.5 +#include <fstream>
4.6 +
4.7 +#include <j_graph.h>
4.8 +#include <dimacs_jgraph.hh>
4.9 +#include <preflow_jgraph.h>
4.10 +#include <time_measure.h>
4.11 +
4.12 +using namespace hugo;
4.13 +
4.14 +// Use a DIMACS max flow file as stdin.
4.15 +// read_dimacs_demo < dimacs_max_flow_file
4.16 +int main(int, char **) {
4.17 + typedef JGraph::TrivNodeIt TrivNodeIt;
4.18 + typedef JGraph::EdgeIt EdgeIt;
4.19 +
4.20 + JGraph G;
4.21 + TrivNodeIt s, t;
4.22 + JGraph::EdgeMap<int> cap(G);
4.23 + readDimacsMaxFlow(std::cin, G, s, t, cap);
4.24 +
4.25 + std::cout << "preflow demo jacint ..." << std::endl;
4.26 +
4.27 + double mintime=1000000;
4.28 +
4.29 + for ( int i=1; i!=11; ++i ) {
4.30 + double pre_time=currTime();
4.31 + preflow<JGraph, int> max_flow_test(G, s, t, cap);
4.32 + double post_time=currTime();
4.33 + if ( mintime > post_time-pre_time ) mintime = post_time-pre_time;
4.34 + }
4.35 +
4.36 + double pre_time=currTime();
4.37 + preflow<JGraph, int> max_flow_test(G, s, t, cap);
4.38 + double post_time=currTime();
4.39 +
4.40 + JGraph::NodeMap<bool> cut(G);
4.41 + max_flow_test.minCut(cut);
4.42 + int min_cut_value=0;
4.43 + for(EdgeIt e=G.firstEdge(); e; G.next(e)) {
4.44 + if (cut.get(G.tail(e)) && !cut.get(G.head(e))) min_cut_value+=cap.get(e);
4.45 + }
4.46 +
4.47 + JGraph::NodeMap<bool> cut1(G);
4.48 + max_flow_test.minMinCut(cut1);
4.49 + int min_min_cut_value=0;
4.50 + for(EdgeIt e=G.firstEdge(); e; G.next(e)) {
4.51 + if (cut.get(G.tail(e)) && !cut.get(G.head(e)))
4.52 + min_min_cut_value+=cap.get(e);
4.53 + }
4.54 +
4.55 + JGraph::NodeMap<bool> cut2(G);
4.56 + max_flow_test.maxMinCut(cut2);
4.57 + int max_min_cut_value=0;
4.58 + for(EdgeIt e=G.firstEdge(); e; G.next(e)) {
4.59 + if (cut2.get(G.tail(e)) && !cut2.get(G.head(e)))
4.60 + max_min_cut_value+=cap.get(e);
4.61 + }
4.62 +
4.63 + std::cout << "min time of 10 runs: " << mintime << " sec"<< std::endl;
4.64 + std::cout << "phase 0: " << max_flow_test.time-pre_time
4.65 + << " sec"<< std::endl;
4.66 + std::cout << "phase 1: " << post_time-max_flow_test.time
4.67 + << " sec"<< std::endl;
4.68 + std::cout << "flow value: "<< max_flow_test.maxFlow() << std::endl;
4.69 + std::cout << "min cut value: "<< min_cut_value << std::endl;
4.70 + std::cout << "min min cut value: "<< min_min_cut_value << std::endl;
4.71 + std::cout << "max min cut value: "<< max_min_cut_value <<
4.72 + std::endl<< std::endl;
4.73 +
4.74 + return 0;
4.75 +}
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/src/work/jacint/preflow_jgraph.h Thu Feb 26 11:38:51 2004 +0000
5.3 @@ -0,0 +1,536 @@
5.4 +// -*- C++ -*-
5.5 +/*
5.6 +preflow.h with 'j_graph interface'
5.7 +by jacint.
5.8 +Heuristics:
5.9 + 2 phase
5.10 + gap
5.11 + list 'level_list' on the nodes on level i implemented by hand
5.12 + stack 'active' on the active nodes on level i implemented by hand
5.13 + runs heuristic 'highest label' for H1*n relabels
5.14 + runs heuristic 'bound decrease' for H0*n relabels, starts with 'highest label'
5.15 +
5.16 +Parameters H0 and H1 are initialized to 20 and 10.
5.17 +
5.18 +The best preflow I could ever write.
5.19 +
5.20 +The constructor runs the algorithm.
5.21 +
5.22 +Members:
5.23 +
5.24 +T maxFlow() : returns the value of a maximum flow
5.25 +
5.26 +T flowOnEdge(EdgeIt e) : for a fixed maximum flow x it returns x(e)
5.27 +
5.28 +FlowMap Flow() : returns the fixed maximum flow x
5.29 +
5.30 +void minMinCut(CutMap& M) : sets M to the characteristic vector of the
5.31 + minimum min cut. M should be a map of bools initialized to false.
5.32 +
5.33 +void maxMinCut(CutMap& M) : sets M to the characteristic vector of the
5.34 + maximum min cut. M should be a map of bools initialized to false.
5.35 +
5.36 +void minCut(CutMap& M) : sets M to the characteristic vector of
5.37 + a min cut. M should be a map of bools initialized to false.
5.38 +
5.39 +*/
5.40 +
5.41 +#ifndef PREFLOW_H
5.42 +#define PREFLOW_H
5.43 +
5.44 +#define H0 20
5.45 +#define H1 1
5.46 +
5.47 +#include <vector>
5.48 +#include <queue>
5.49 +
5.50 +#include<iostream>
5.51 +
5.52 +#include <time_measure.h>
5.53 +
5.54 +namespace hugo {
5.55 +
5.56 + template <typename Graph, typename T,
5.57 + typename FlowMap=typename Graph::EdgeMap<T>,
5.58 + typename CapMap=typename Graph::EdgeMap<T> >
5.59 + class preflow {
5.60 +
5.61 + typedef typename Graph::TrivNodeIt NodeIt;
5.62 + typedef typename Graph::TrivEdgeIt EdgeIt;
5.63 + typedef typename Graph::NodeIt EachNodeIt;
5.64 + typedef typename Graph::OutEdgeIt OutEdgeIt;
5.65 + typedef typename Graph::InEdgeIt InEdgeIt;
5.66 +
5.67 + Graph& G;
5.68 + NodeIt s;
5.69 + NodeIt t;
5.70 + FlowMap flow;
5.71 + CapMap& capacity;
5.72 + T value;
5.73 +
5.74 + public:
5.75 + double time;
5.76 + preflow(Graph& _G, NodeIt _s, NodeIt _t, CapMap& _capacity ) :
5.77 + G(_G), s(_s), t(_t), flow(_G, 0), capacity(_capacity)
5.78 + {
5.79 +
5.80 + bool phase=0; //phase 0 is the 1st phase, phase 1 is the 2nd
5.81 + int n=G.numNodes();
5.82 + int heur0=(int)(H0*n); //time while running 'bound decrease'
5.83 + int heur1=(int)(H1*n); //time while running 'highest label'
5.84 + int heur=heur1; //starting time interval (#of relabels)
5.85 + bool what_heur=1;
5.86 + /*
5.87 + what_heur is 0 in case 'bound decrease'
5.88 + and 1 in case 'highest label'
5.89 + */
5.90 + bool end=false;
5.91 + /*
5.92 + Needed for 'bound decrease', 'true'
5.93 + means no active nodes are above bound b.
5.94 + */
5.95 + int relabel=0;
5.96 + int k=n-2; //bound on the highest level under n containing a node
5.97 + int b=k; //bound on the highest level under n of an active node
5.98 +
5.99 + typename Graph::NodeMap<int> level(G,n);
5.100 + typename Graph::NodeMap<T> excess(G);
5.101 +
5.102 + std::vector<NodeIt> active(n);
5.103 + typename Graph::NodeMap<NodeIt> next(G);
5.104 + //Stack of the active nodes in level i < n.
5.105 + //We use it in both phases.
5.106 +
5.107 + typename Graph::NodeMap<NodeIt> left(G);
5.108 + typename Graph::NodeMap<NodeIt> right(G);
5.109 + std::vector<NodeIt> level_list(n);
5.110 + /*
5.111 + List of the nodes in level i<n.
5.112 + */
5.113 +
5.114 + /*Reverse_bfs from t, to find the starting level.*/
5.115 + level.set(t,0);
5.116 + std::queue<NodeIt> bfs_queue;
5.117 + bfs_queue.push(t);
5.118 +
5.119 + while (!bfs_queue.empty()) {
5.120 +
5.121 + NodeIt v=bfs_queue.front();
5.122 + bfs_queue.pop();
5.123 + int l=level.get(v)+1;
5.124 +
5.125 + for(InEdgeIt e=G.firstIn(v); e; G.next(e)) {
5.126 + NodeIt w=G.tail(e);
5.127 + if ( level.get(w) == n && w != s ) {
5.128 + bfs_queue.push(w);
5.129 + NodeIt first=level_list[l];
5.130 + if ( first ) left.set(first,w);
5.131 + right.set(w,first);
5.132 + level_list[l]=w;
5.133 + level.set(w, l);
5.134 + }
5.135 + }
5.136 + }
5.137 +
5.138 + level.set(s,n);
5.139 +
5.140 +
5.141 + /* Starting flow. It is everywhere 0 at the moment. */
5.142 + for(OutEdgeIt e=G.firstOut(s); e; G.next(e))
5.143 + {
5.144 + T c=capacity.get(e);
5.145 + if ( c == 0 ) continue;
5.146 + NodeIt w=G.head(e);
5.147 + if ( level.get(w) < n ) {
5.148 + if ( excess.get(w) == 0 && w!=t ) {
5.149 + next.set(w,active[level.get(w)]);
5.150 + active[level.get(w)]=w;
5.151 + }
5.152 + flow.set(e, c);
5.153 + excess.set(w, excess.get(w)+c);
5.154 + }
5.155 + }
5.156 +
5.157 + /*
5.158 + End of preprocessing
5.159 + */
5.160 +
5.161 +
5.162 +
5.163 + /*
5.164 + Push/relabel on the highest level active nodes.
5.165 + */
5.166 + while ( true ) {
5.167 +
5.168 + if ( b == 0 ) {
5.169 + if ( phase ) break;
5.170 +
5.171 + if ( !what_heur && !end && k > 0 ) {
5.172 + b=k;
5.173 + end=true;
5.174 + } else {
5.175 + phase=1;
5.176 + time=currTime();
5.177 + level.set(s,0);
5.178 + std::queue<NodeIt> bfs_queue;
5.179 + bfs_queue.push(s);
5.180 +
5.181 + while (!bfs_queue.empty()) {
5.182 +
5.183 + NodeIt v=bfs_queue.front();
5.184 + bfs_queue.pop();
5.185 + int l=level.get(v)+1;
5.186 +
5.187 + for(InEdgeIt e=G.firstIn(v); e; G.next(e)) {
5.188 + if ( capacity.get(e) == flow.get(e) ) continue;
5.189 + NodeIt u=G.tail(e);
5.190 + if ( level.get(u) >= n ) {
5.191 + bfs_queue.push(u);
5.192 + level.set(u, l);
5.193 + if ( excess.get(u) > 0 ) {
5.194 + next.set(u,active[l]);
5.195 + active[l]=u;
5.196 + }
5.197 + }
5.198 + }
5.199 +
5.200 + for(OutEdgeIt e=G.firstOut(v); e; G.next(e)) {
5.201 + if ( 0 == flow.get(e) ) continue;
5.202 + NodeIt u=G.head(e);
5.203 + if ( level.get(u) >= n ) {
5.204 + bfs_queue.push(u);
5.205 + level.set(u, l);
5.206 + if ( excess.get(u) > 0 ) {
5.207 + next.set(u,active[l]);
5.208 + active[l]=u;
5.209 + }
5.210 + }
5.211 + }
5.212 + }
5.213 + b=n-2;
5.214 + }
5.215 +
5.216 + }
5.217 +
5.218 +
5.219 + if ( !active[b] ) --b;
5.220 + else {
5.221 + end=false;
5.222 +
5.223 + NodeIt w=active[b];
5.224 + active[b]=next.get(w);
5.225 + int lev=level.get(w);
5.226 + T exc=excess.get(w);
5.227 + int newlevel=n; //bound on the next level of w
5.228 +
5.229 + for(OutEdgeIt e=G.firstOut(w); e; G.next(e)) {
5.230 +
5.231 + if ( flow.get(e) == capacity.get(e) ) continue;
5.232 + NodeIt v=G.head(e);
5.233 + //e=wv
5.234 +
5.235 + if( lev > level.get(v) ) {
5.236 + /*Push is allowed now*/
5.237 +
5.238 + if ( excess.get(v)==0 && v!=t && v!=s ) {
5.239 + int lev_v=level.get(v);
5.240 + next.set(v,active[lev_v]);
5.241 + active[lev_v]=v;
5.242 + }
5.243 +
5.244 + T cap=capacity.get(e);
5.245 + T flo=flow.get(e);
5.246 + T remcap=cap-flo;
5.247 +
5.248 + if ( remcap >= exc ) {
5.249 + /*A nonsaturating push.*/
5.250 +
5.251 + flow.set(e, flo+exc);
5.252 + excess.set(v, excess.get(v)+exc);
5.253 + exc=0;
5.254 + break;
5.255 +
5.256 + } else {
5.257 + /*A saturating push.*/
5.258 +
5.259 + flow.set(e, cap);
5.260 + excess.set(v, excess.get(v)+remcap);
5.261 + exc-=remcap;
5.262 + }
5.263 + } else if ( newlevel > level.get(v) ){
5.264 + newlevel = level.get(v);
5.265 + }
5.266 +
5.267 + } //for out edges wv
5.268 +
5.269 +
5.270 + if ( exc > 0 ) {
5.271 + for( InEdgeIt e=G.firstIn(w); e; G.next(e)) {
5.272 +
5.273 + if( flow.get(e) == 0 ) continue;
5.274 + NodeIt v=G.tail(e);
5.275 + //e=vw
5.276 +
5.277 + if( lev > level.get(v) ) {
5.278 + /*Push is allowed now*/
5.279 +
5.280 + if ( excess.get(v)==0 && v!=t && v!=s ) {
5.281 + int lev_v=level.get(v);
5.282 + next.set(v,active[lev_v]);
5.283 + active[lev_v]=v;
5.284 + }
5.285 +
5.286 + T flo=flow.get(e);
5.287 +
5.288 + if ( flo >= exc ) {
5.289 + /*A nonsaturating push.*/
5.290 +
5.291 + flow.set(e, flo-exc);
5.292 + excess.set(v, excess.get(v)+exc);
5.293 + exc=0;
5.294 + break;
5.295 + } else {
5.296 +/*A saturating push.*/
5.297 +
5.298 + excess.set(v, excess.get(v)+flo);
5.299 + exc-=flo;
5.300 + flow.set(e,0);
5.301 + }
5.302 + } else if ( newlevel > level.get(v) ) {
5.303 + newlevel = level.get(v);
5.304 + }
5.305 + } //for in edges vw
5.306 +
5.307 + } // if w still has excess after the out edge for cycle
5.308 +
5.309 + excess.set(w, exc);
5.310 +
5.311 + /*
5.312 + Relabel
5.313 + */
5.314 +
5.315 +
5.316 + if ( exc > 0 ) {
5.317 + //now 'lev' is the old level of w
5.318 +
5.319 + if ( phase ) {
5.320 + level.set(w,++newlevel);
5.321 + next.set(w,active[newlevel]);
5.322 + active[newlevel]=w;
5.323 + b=newlevel;
5.324 + } else {
5.325 + //unlacing starts
5.326 + NodeIt right_n=right.get(w);
5.327 + NodeIt left_n=left.get(w);
5.328 +
5.329 + if ( right_n ) {
5.330 + if ( left_n ) {
5.331 + right.set(left_n, right_n);
5.332 + left.set(right_n, left_n);
5.333 + } else {
5.334 + level_list[lev]=right_n;
5.335 + left.set(right_n, NodeIt());
5.336 + }
5.337 + } else {
5.338 + if ( left_n ) {
5.339 + right.set(left_n, NodeIt());
5.340 + } else {
5.341 + level_list[lev]=NodeIt();
5.342 +
5.343 + }
5.344 + }
5.345 + //unlacing ends
5.346 +
5.347 + //gapping starts
5.348 + if ( !level_list[lev] ) {
5.349 +
5.350 + for (int i=lev; i!=k ; ) {
5.351 + NodeIt v=level_list[++i];
5.352 + while ( v ) {
5.353 + level.set(v,n);
5.354 + v=right.get(v);
5.355 + }
5.356 + level_list[i]=NodeIt();
5.357 + if ( !what_heur ) active[i]=NodeIt();
5.358 + }
5.359 +
5.360 + level.set(w,n);
5.361 + b=lev-1;
5.362 + k=b;
5.363 + //gapping ends
5.364 + } else {
5.365 +
5.366 + if ( newlevel == n ) level.set(w,n);
5.367 + else {
5.368 + level.set(w,++newlevel);
5.369 + next.set(w,active[newlevel]);
5.370 + active[newlevel]=w;
5.371 + if ( what_heur ) b=newlevel;
5.372 + if ( k < newlevel ) ++k;
5.373 + NodeIt first=level_list[newlevel];
5.374 + if ( first ) left.set(first,w);
5.375 + right.set(w,first);
5.376 + left.set(w,NodeIt());
5.377 + level_list[newlevel]=w;
5.378 + }
5.379 + }
5.380 +
5.381 +
5.382 + ++relabel;
5.383 + if ( relabel >= heur ) {
5.384 + relabel=0;
5.385 + if ( what_heur ) {
5.386 + what_heur=0;
5.387 + heur=heur0;
5.388 + end=false;
5.389 + } else {
5.390 + what_heur=1;
5.391 + heur=heur1;
5.392 + b=k;
5.393 + }
5.394 + }
5.395 + } //phase 0
5.396 +
5.397 +
5.398 + } // if ( exc > 0 )
5.399 +
5.400 +
5.401 + } // if stack[b] is nonempty
5.402 +
5.403 + } // while(true)
5.404 +
5.405 +
5.406 + value = excess.get(t);
5.407 + /*Max flow value.*/
5.408 +
5.409 + } //void run()
5.410 +
5.411 +
5.412 +
5.413 +
5.414 +
5.415 + /*
5.416 + Returns the maximum value of a flow.
5.417 + */
5.418 +
5.419 + T maxFlow() {
5.420 + return value;
5.421 + }
5.422 +
5.423 +
5.424 +
5.425 + /*
5.426 + For the maximum flow x found by the algorithm,
5.427 + it returns the flow value on edge e, i.e. x(e).
5.428 + */
5.429 +
5.430 + T flowOnEdge(EdgeIt e) {
5.431 + return flow.get(e);
5.432 + }
5.433 +
5.434 +
5.435 +
5.436 + FlowMap Flow() {
5.437 + return flow;
5.438 + }
5.439 +
5.440 +
5.441 +
5.442 + void Flow(FlowMap& _flow ) {
5.443 + for(EachNodeIt v=G.firstNode() ; v; G.next(v))
5.444 + _flow.set(v,flow.get(v));
5.445 + }
5.446 +
5.447 +
5.448 +
5.449 + /*
5.450 + Returns the minimum min cut, by a bfs from s in the residual graph.
5.451 + */
5.452 +
5.453 + template<typename _CutMap>
5.454 + void minMinCut(_CutMap& M) {
5.455 +
5.456 + std::queue<NodeIt> queue;
5.457 +
5.458 + M.set(s,true);
5.459 + queue.push(s);
5.460 +
5.461 + while (!queue.empty()) {
5.462 + NodeIt w=queue.front();
5.463 + queue.pop();
5.464 +
5.465 + for(OutEdgeIt e=G.firstOut(w) ; e; G.next(e)) {
5.466 + NodeIt v=G.head(e);
5.467 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
5.468 + queue.push(v);
5.469 + M.set(v, true);
5.470 + }
5.471 + }
5.472 +
5.473 + for(InEdgeIt e=G.firstIn(w) ; e; G.next(e)) {
5.474 + NodeIt v=G.tail(e);
5.475 + if (!M.get(v) && flow.get(e) > 0 ) {
5.476 + queue.push(v);
5.477 + M.set(v, true);
5.478 + }
5.479 + }
5.480 + }
5.481 + }
5.482 +
5.483 +
5.484 +
5.485 + /*
5.486 + Returns the maximum min cut, by a reverse bfs
5.487 + from t in the residual graph.
5.488 + */
5.489 +
5.490 + template<typename _CutMap>
5.491 + void maxMinCut(_CutMap& M) {
5.492 +
5.493 + std::queue<NodeIt> queue;
5.494 +
5.495 + M.set(t,true);
5.496 + queue.push(t);
5.497 +
5.498 + while (!queue.empty()) {
5.499 + NodeIt w=queue.front();
5.500 + queue.pop();
5.501 +
5.502 + for(InEdgeIt e=G.firstIn(w) ; e; G.next(e)) {
5.503 + NodeIt v=G.tail(e);
5.504 + if (!M.get(v) && flow.get(e) < capacity.get(e) ) {
5.505 + queue.push(v);
5.506 + M.set(v, true);
5.507 + }
5.508 + }
5.509 +
5.510 + for(OutEdgeIt e=G.firstOut(w) ; e; G.next(e)) {
5.511 + NodeIt v=G.head(e);
5.512 + if (!M.get(v) && flow.get(e) > 0 ) {
5.513 + queue.push(v);
5.514 + M.set(v, true);
5.515 + }
5.516 + }
5.517 + }
5.518 +
5.519 + for(EachNodeIt v=G.firstNode() ; v; G.next(v)) {
5.520 + M.set(v, !M.get(v));
5.521 + }
5.522 +
5.523 + }
5.524 +
5.525 +
5.526 +
5.527 + template<typename CutMap>
5.528 + void minCut(CutMap& M) {
5.529 + minMinCut(M);
5.530 + }
5.531 +
5.532 +
5.533 + };
5.534 +}//namespace marci
5.535 +#endif
5.536 +
5.537 +
5.538 +
5.539 +