1.1 --- a/demo/Makefile.am Fri Nov 21 10:49:39 2008 +0000
1.2 +++ b/demo/Makefile.am Fri Nov 21 14:42:47 2008 +0000
1.3 @@ -1,16 +1,19 @@
1.4 EXTRA_DIST += \
1.5 demo/CMakeLists.txt \
1.6 + demo/circulation-input.lgf \
1.7 demo/digraph.lgf
1.8
1.9 if WANT_DEMO
1.10
1.11 noinst_PROGRAMS += \
1.12 demo/arg_parser_demo \
1.13 + demo/circulation_demo \
1.14 demo/graph_to_eps_demo \
1.15 demo/lgf_demo
1.16
1.17 endif WANT_DEMO
1.18
1.19 demo_arg_parser_demo_SOURCES = demo/arg_parser_demo.cc
1.20 +demo_circulation_demo_SOURCES = demo/circulation_demo.cc
1.21 demo_graph_to_eps_demo_SOURCES = demo/graph_to_eps_demo.cc
1.22 demo_lgf_demo_SOURCES = demo/lgf_demo.cc
2.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
2.2 +++ b/demo/circulation-input.lgf Fri Nov 21 14:42:47 2008 +0000
2.3 @@ -0,0 +1,32 @@
2.4 +@nodes
2.5 +coordinates_x coordinates_y delta label
2.6 +-396.638 -311.798 0 0
2.7 +154.409 -214.714 13 1
2.8 +-378.119 -135.808 0 2
2.9 +-138.182 -58.0452 0 3
2.10 +55 -76.1018 0 4
2.11 +-167.302 169.88 0 5
2.12 +71.6876 38.7452 0 6
2.13 +-328.784 257.777 0 7
2.14 +354.242 67.9628 -13 8
2.15 +147 266 0 9
2.16 +@edges
2.17 + label lo_cap up_cap
2.18 +0 1 0 0 20
2.19 +0 2 1 0 0
2.20 +1 1 2 0 3
2.21 +1 2 3 0 8
2.22 +1 3 4 0 8
2.23 +2 5 5 0 5
2.24 +3 2 6 0 5
2.25 +3 5 7 0 5
2.26 +3 6 8 0 5
2.27 +4 3 9 0 3
2.28 +5 7 10 0 3
2.29 +5 6 11 0 10
2.30 +5 8 12 0 10
2.31 +6 8 13 0 8
2.32 +8 9 14 0 20
2.33 +8 1 15 0 5
2.34 +9 5 16 0 5
2.35 +@end
3.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
3.2 +++ b/demo/circulation_demo.cc Fri Nov 21 14:42:47 2008 +0000
3.3 @@ -0,0 +1,107 @@
3.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
3.5 + *
3.6 + * This file is a part of LEMON, a generic C++ optimization library.
3.7 + *
3.8 + * Copyright (C) 2003-2008
3.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
3.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
3.11 + *
3.12 + * Permission to use, modify and distribute this software is granted
3.13 + * provided that this copyright notice appears in all copies. For
3.14 + * precise terms see the accompanying LICENSE file.
3.15 + *
3.16 + * This software is provided "AS IS" with no warranty of any kind,
3.17 + * express or implied, and with no claim as to its suitability for any
3.18 + * purpose.
3.19 + *
3.20 + */
3.21 +
3.22 +///\ingroup demos
3.23 +///\file
3.24 +///\brief Demonstrating the usage of LEMON's General Flow algorithm
3.25 +///
3.26 +/// This demo program reads a general network circulation problem from the
3.27 +/// file 'circulation-input.lgf', runs the preflow based algorithm for
3.28 +/// finding a feasible solution and writes the output
3.29 +/// to 'circulation-input.lgf'
3.30 +///
3.31 +/// \include circulation_demo.cc
3.32 +
3.33 +#include <iostream>
3.34 +
3.35 +#include <lemon/list_graph.h>
3.36 +#include <lemon/circulation.h>
3.37 +#include <lemon/lgf_reader.h>
3.38 +#include <lemon/lgf_writer.h>
3.39 +
3.40 +using namespace lemon;
3.41 +
3.42 +
3.43 +int main (int, char*[])
3.44 +{
3.45 +
3.46 + typedef ListDigraph Digraph;
3.47 + typedef Digraph::Node Node;
3.48 + typedef Digraph::NodeIt NodeIt;
3.49 + typedef Digraph::Arc Arc;
3.50 + typedef Digraph::ArcIt ArcIt;
3.51 + typedef Digraph::ArcMap<int> ArcMap;
3.52 + typedef Digraph::NodeMap<int> NodeMap;
3.53 + typedef Digraph::NodeMap<double> DNodeMap;
3.54 +
3.55 + Digraph g;
3.56 + ArcMap lo(g);
3.57 + ArcMap up(g);
3.58 + NodeMap delta(g);
3.59 + NodeMap nid(g);
3.60 + ArcMap eid(g);
3.61 + DNodeMap cx(g);
3.62 + DNodeMap cy(g);
3.63 +
3.64 + DigraphReader<Digraph>(g,"circulation-input.lgf").
3.65 + arcMap("lo_cap", lo).
3.66 + arcMap("up_cap", up).
3.67 + nodeMap("delta", delta).
3.68 + arcMap("label", eid).
3.69 + nodeMap("label", nid).
3.70 + nodeMap("coordinates_x", cx).
3.71 + nodeMap("coordinates_y", cy).
3.72 + run();
3.73 +
3.74 + Circulation<Digraph> gen(g,lo,up,delta);
3.75 + bool ret=gen.run();
3.76 + if(ret)
3.77 + {
3.78 + std::cout << "\nA feasible flow has been found.\n";
3.79 + if(!gen.checkFlow()) std::cerr << "Oops!!!\n";
3.80 + DigraphWriter<Digraph>(g, "circulation-output.lgf").
3.81 + arcMap("lo_cap", lo).
3.82 + arcMap("up_cap", up).
3.83 + arcMap("flow", gen.flowMap()).
3.84 + nodeMap("delta", delta).
3.85 + arcMap("label", eid).
3.86 + nodeMap("label", nid).
3.87 + nodeMap("coordinates_x", cx).
3.88 + nodeMap("coordinates_y", cy).
3.89 + run();
3.90 + }
3.91 + else {
3.92 + std::cout << "\nThere is no such a flow\n";
3.93 + Digraph::NodeMap<int> bar(g);
3.94 + gen.barrierMap(bar);
3.95 + if(!gen.checkBarrier()) std::cerr << "Dual Oops!!!\n";
3.96 +
3.97 + DigraphWriter<Digraph>(g, "circulation-output.lgf").
3.98 + arcMap("lo_cap", lo).
3.99 + arcMap("up_cap", up).
3.100 + nodeMap("barrier", bar).
3.101 + nodeMap("delta", delta).
3.102 + arcMap("label", eid).
3.103 + nodeMap("label", nid).
3.104 + nodeMap("coordinates_x", cx).
3.105 + nodeMap("coordinates_y", cy).
3.106 + run();
3.107 + }
3.108 + std::cout << "The output is written to file 'circulation-output.lgf'\n\n";
3.109 +
3.110 +}
4.1 --- a/lemon/Makefile.am Fri Nov 21 10:49:39 2008 +0000
4.2 +++ b/lemon/Makefile.am Fri Nov 21 14:42:47 2008 +0000
4.3 @@ -20,6 +20,7 @@
4.4 lemon/assert.h \
4.5 lemon/bfs.h \
4.6 lemon/bin_heap.h \
4.7 + lemon/circulation.h \
4.8 lemon/color.h \
4.9 lemon/concept_check.h \
4.10 lemon/counter.h \
5.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
5.2 +++ b/lemon/circulation.h Fri Nov 21 14:42:47 2008 +0000
5.3 @@ -0,0 +1,656 @@
5.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
5.5 + *
5.6 + * This file is a part of LEMON, a generic C++ optimization library.
5.7 + *
5.8 + * Copyright (C) 2003-2008
5.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
5.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
5.11 + *
5.12 + * Permission to use, modify and distribute this software is granted
5.13 + * provided that this copyright notice appears in all copies. For
5.14 + * precise terms see the accompanying LICENSE file.
5.15 + *
5.16 + * This software is provided "AS IS" with no warranty of any kind,
5.17 + * express or implied, and with no claim as to its suitability for any
5.18 + * purpose.
5.19 + *
5.20 + */
5.21 +
5.22 +#ifndef LEMON_CIRCULATION_H
5.23 +#define LEMON_CIRCULATION_H
5.24 +
5.25 +#include <iostream>
5.26 +#include <queue>
5.27 +#include <lemon/tolerance.h>
5.28 +#include <lemon/elevator.h>
5.29 +
5.30 +///\ingroup max_flow
5.31 +///\file
5.32 +///\brief Push-prelabel algorithm for finding a feasible circulation.
5.33 +///
5.34 +namespace lemon {
5.35 +
5.36 + /// \brief Default traits class of Circulation class.
5.37 + ///
5.38 + /// Default traits class of Circulation class.
5.39 + /// \param _Graph Digraph type.
5.40 + /// \param _CapacityMap Type of capacity map.
5.41 + template <typename _Graph, typename _LCapMap,
5.42 + typename _UCapMap, typename _DeltaMap>
5.43 + struct CirculationDefaultTraits {
5.44 +
5.45 + /// \brief The digraph type the algorithm runs on.
5.46 + typedef _Graph Digraph;
5.47 +
5.48 + /// \brief The type of the map that stores the circulation lower
5.49 + /// bound.
5.50 + ///
5.51 + /// The type of the map that stores the circulation lower bound.
5.52 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
5.53 + typedef _LCapMap LCapMap;
5.54 +
5.55 + /// \brief The type of the map that stores the circulation upper
5.56 + /// bound.
5.57 + ///
5.58 + /// The type of the map that stores the circulation upper bound.
5.59 + /// It must meet the \ref concepts::ReadMap "ReadMap" concept.
5.60 + typedef _UCapMap UCapMap;
5.61 +
5.62 + /// \brief The type of the map that stores the upper bound of
5.63 + /// node excess.
5.64 + ///
5.65 + /// The type of the map that stores the lower bound of node
5.66 + /// excess. It must meet the \ref concepts::ReadMap "ReadMap"
5.67 + /// concept.
5.68 + typedef _DeltaMap DeltaMap;
5.69 +
5.70 + /// \brief The type of the length of the arcs.
5.71 + typedef typename DeltaMap::Value Value;
5.72 +
5.73 + /// \brief The map type that stores the flow values.
5.74 + ///
5.75 + /// The map type that stores the flow values.
5.76 + /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.
5.77 + typedef typename Digraph::template ArcMap<Value> FlowMap;
5.78 +
5.79 + /// \brief Instantiates a FlowMap.
5.80 + ///
5.81 + /// This function instantiates a \ref FlowMap.
5.82 + /// \param digraph The digraph, to which we would like to define
5.83 + /// the flow map.
5.84 + static FlowMap* createFlowMap(const Digraph& digraph) {
5.85 + return new FlowMap(digraph);
5.86 + }
5.87 +
5.88 + /// \brief The eleavator type used by Circulation algorithm.
5.89 + ///
5.90 + /// The elevator type used by Circulation algorithm.
5.91 + ///
5.92 + /// \sa Elevator
5.93 + /// \sa LinkedElevator
5.94 + typedef lemon::Elevator<Digraph, typename Digraph::Node> Elevator;
5.95 +
5.96 + /// \brief Instantiates an Elevator.
5.97 + ///
5.98 + /// This function instantiates a \ref Elevator.
5.99 + /// \param digraph The digraph, to which we would like to define
5.100 + /// the elevator.
5.101 + /// \param max_level The maximum level of the elevator.
5.102 + static Elevator* createElevator(const Digraph& digraph, int max_level) {
5.103 + return new Elevator(digraph, max_level);
5.104 + }
5.105 +
5.106 + /// \brief The tolerance used by the algorithm
5.107 + ///
5.108 + /// The tolerance used by the algorithm to handle inexact computation.
5.109 + typedef lemon::Tolerance<Value> Tolerance;
5.110 +
5.111 + };
5.112 +
5.113 + ///Push-relabel algorithm for the Network Circulation Problem.
5.114 +
5.115 + /**
5.116 + \ingroup max_flow
5.117 + This class implements a push-relabel algorithm
5.118 + or the Network Circulation Problem.
5.119 + The exact formulation of this problem is the following.
5.120 + \f[\sum_{e\in\rho(v)}x(e)-\sum_{e\in\delta(v)}x(e)\leq
5.121 + -delta(v)\quad \forall v\in V \f]
5.122 + \f[ lo(e)\leq x(e) \leq up(e) \quad \forall e\in E \f]
5.123 + */
5.124 + template<class _Graph,
5.125 + class _LCapMap=typename _Graph::template ArcMap<int>,
5.126 + class _UCapMap=_LCapMap,
5.127 + class _DeltaMap=typename _Graph::template NodeMap<
5.128 + typename _UCapMap::Value>,
5.129 + class _Traits=CirculationDefaultTraits<_Graph, _LCapMap,
5.130 + _UCapMap, _DeltaMap> >
5.131 + class Circulation {
5.132 +
5.133 + typedef _Traits Traits;
5.134 + typedef typename Traits::Digraph Digraph;
5.135 + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
5.136 +
5.137 + typedef typename Traits::Value Value;
5.138 +
5.139 + typedef typename Traits::LCapMap LCapMap;
5.140 + typedef typename Traits::UCapMap UCapMap;
5.141 + typedef typename Traits::DeltaMap DeltaMap;
5.142 + typedef typename Traits::FlowMap FlowMap;
5.143 + typedef typename Traits::Elevator Elevator;
5.144 + typedef typename Traits::Tolerance Tolerance;
5.145 +
5.146 + typedef typename Digraph::template NodeMap<Value> ExcessMap;
5.147 +
5.148 + const Digraph &_g;
5.149 + int _node_num;
5.150 +
5.151 + const LCapMap *_lo;
5.152 + const UCapMap *_up;
5.153 + const DeltaMap *_delta;
5.154 +
5.155 + FlowMap *_flow;
5.156 + bool _local_flow;
5.157 +
5.158 + Elevator* _level;
5.159 + bool _local_level;
5.160 +
5.161 + ExcessMap* _excess;
5.162 +
5.163 + Tolerance _tol;
5.164 + int _el;
5.165 +
5.166 + public:
5.167 +
5.168 + typedef Circulation Create;
5.169 +
5.170 + ///\name Named template parameters
5.171 +
5.172 + ///@{
5.173 +
5.174 + template <typename _FlowMap>
5.175 + struct DefFlowMapTraits : public Traits {
5.176 + typedef _FlowMap FlowMap;
5.177 + static FlowMap *createFlowMap(const Digraph&) {
5.178 + LEMON_ASSERT(false, "FlowMap is not initialized");
5.179 + return 0; // ignore warnings
5.180 + }
5.181 + };
5.182 +
5.183 + /// \brief \ref named-templ-param "Named parameter" for setting
5.184 + /// FlowMap type
5.185 + ///
5.186 + /// \ref named-templ-param "Named parameter" for setting FlowMap
5.187 + /// type
5.188 + template <typename _FlowMap>
5.189 + struct DefFlowMap
5.190 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
5.191 + DefFlowMapTraits<_FlowMap> > {
5.192 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
5.193 + DefFlowMapTraits<_FlowMap> > Create;
5.194 + };
5.195 +
5.196 + template <typename _Elevator>
5.197 + struct DefElevatorTraits : public Traits {
5.198 + typedef _Elevator Elevator;
5.199 + static Elevator *createElevator(const Digraph&, int) {
5.200 + LEMON_ASSERT(false, "Elevator is not initialized");
5.201 + return 0; // ignore warnings
5.202 + }
5.203 + };
5.204 +
5.205 + /// \brief \ref named-templ-param "Named parameter" for setting
5.206 + /// Elevator type
5.207 + ///
5.208 + /// \ref named-templ-param "Named parameter" for setting Elevator
5.209 + /// type
5.210 + template <typename _Elevator>
5.211 + struct DefElevator
5.212 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
5.213 + DefElevatorTraits<_Elevator> > {
5.214 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
5.215 + DefElevatorTraits<_Elevator> > Create;
5.216 + };
5.217 +
5.218 + template <typename _Elevator>
5.219 + struct DefStandardElevatorTraits : public Traits {
5.220 + typedef _Elevator Elevator;
5.221 + static Elevator *createElevator(const Digraph& digraph, int max_level) {
5.222 + return new Elevator(digraph, max_level);
5.223 + }
5.224 + };
5.225 +
5.226 + /// \brief \ref named-templ-param "Named parameter" for setting
5.227 + /// Elevator type
5.228 + ///
5.229 + /// \ref named-templ-param "Named parameter" for setting Elevator
5.230 + /// type. The Elevator should be standard constructor interface, ie.
5.231 + /// the digraph and the maximum level should be passed to it.
5.232 + template <typename _Elevator>
5.233 + struct DefStandardElevator
5.234 + : public Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
5.235 + DefStandardElevatorTraits<_Elevator> > {
5.236 + typedef Circulation<Digraph, LCapMap, UCapMap, DeltaMap,
5.237 + DefStandardElevatorTraits<_Elevator> > Create;
5.238 + };
5.239 +
5.240 + /// @}
5.241 +
5.242 + protected:
5.243 +
5.244 + Circulation() {}
5.245 +
5.246 + public:
5.247 +
5.248 + /// The constructor of the class.
5.249 +
5.250 + /// The constructor of the class.
5.251 + /// \param g The digraph the algorithm runs on.
5.252 + /// \param lo The lower bound capacity of the arcs.
5.253 + /// \param up The upper bound capacity of the arcs.
5.254 + /// \param delta The lower bound on node excess.
5.255 + Circulation(const Digraph &g,const LCapMap &lo,
5.256 + const UCapMap &up,const DeltaMap &delta)
5.257 + : _g(g), _node_num(),
5.258 + _lo(&lo),_up(&up),_delta(&delta),_flow(0),_local_flow(false),
5.259 + _level(0), _local_level(false), _excess(0), _el() {}
5.260 +
5.261 + /// Destrcutor.
5.262 + ~Circulation() {
5.263 + destroyStructures();
5.264 + }
5.265 +
5.266 + private:
5.267 +
5.268 + void createStructures() {
5.269 + _node_num = _el = countNodes(_g);
5.270 +
5.271 + if (!_flow) {
5.272 + _flow = Traits::createFlowMap(_g);
5.273 + _local_flow = true;
5.274 + }
5.275 + if (!_level) {
5.276 + _level = Traits::createElevator(_g, _node_num);
5.277 + _local_level = true;
5.278 + }
5.279 + if (!_excess) {
5.280 + _excess = new ExcessMap(_g);
5.281 + }
5.282 + }
5.283 +
5.284 + void destroyStructures() {
5.285 + if (_local_flow) {
5.286 + delete _flow;
5.287 + }
5.288 + if (_local_level) {
5.289 + delete _level;
5.290 + }
5.291 + if (_excess) {
5.292 + delete _excess;
5.293 + }
5.294 + }
5.295 +
5.296 + public:
5.297 +
5.298 + /// Sets the lower bound capacity map.
5.299 +
5.300 + /// Sets the lower bound capacity map.
5.301 + /// \return \c (*this)
5.302 + Circulation& lowerCapMap(const LCapMap& map) {
5.303 + _lo = ↦
5.304 + return *this;
5.305 + }
5.306 +
5.307 + /// Sets the upper bound capacity map.
5.308 +
5.309 + /// Sets the upper bound capacity map.
5.310 + /// \return \c (*this)
5.311 + Circulation& upperCapMap(const LCapMap& map) {
5.312 + _up = ↦
5.313 + return *this;
5.314 + }
5.315 +
5.316 + /// Sets the lower bound map on excess.
5.317 +
5.318 + /// Sets the lower bound map on excess.
5.319 + /// \return \c (*this)
5.320 + Circulation& deltaMap(const DeltaMap& map) {
5.321 + _delta = ↦
5.322 + return *this;
5.323 + }
5.324 +
5.325 + /// Sets the flow map.
5.326 +
5.327 + /// Sets the flow map.
5.328 + /// \return \c (*this)
5.329 + Circulation& flowMap(FlowMap& map) {
5.330 + if (_local_flow) {
5.331 + delete _flow;
5.332 + _local_flow = false;
5.333 + }
5.334 + _flow = ↦
5.335 + return *this;
5.336 + }
5.337 +
5.338 + /// Returns the flow map.
5.339 +
5.340 + /// \return The flow map.
5.341 + ///
5.342 + const FlowMap& flowMap() {
5.343 + return *_flow;
5.344 + }
5.345 +
5.346 + /// Sets the elevator.
5.347 +
5.348 + /// Sets the elevator.
5.349 + /// \return \c (*this)
5.350 + Circulation& elevator(Elevator& elevator) {
5.351 + if (_local_level) {
5.352 + delete _level;
5.353 + _local_level = false;
5.354 + }
5.355 + _level = &elevator;
5.356 + return *this;
5.357 + }
5.358 +
5.359 + /// Returns the elevator.
5.360 +
5.361 + /// \return The elevator.
5.362 + ///
5.363 + const Elevator& elevator() {
5.364 + return *_level;
5.365 + }
5.366 +
5.367 + /// Sets the tolerance used by algorithm.
5.368 +
5.369 + /// Sets the tolerance used by algorithm.
5.370 + ///
5.371 + Circulation& tolerance(const Tolerance& tolerance) const {
5.372 + _tol = tolerance;
5.373 + return *this;
5.374 + }
5.375 +
5.376 + /// Returns the tolerance used by algorithm.
5.377 +
5.378 + /// Returns the tolerance used by algorithm.
5.379 + ///
5.380 + const Tolerance& tolerance() const {
5.381 + return tolerance;
5.382 + }
5.383 +
5.384 + /// \name Execution control
5.385 + /// The simplest way to execute the algorithm is to use one of the
5.386 + /// member functions called \c run().
5.387 + /// \n
5.388 + /// If you need more control on initial solution or execution then
5.389 + /// you have to call one \ref init() function and then the start()
5.390 + /// function.
5.391 +
5.392 + ///@{
5.393 +
5.394 + /// Initializes the internal data structures.
5.395 +
5.396 + /// Initializes the internal data structures. This function sets
5.397 + /// all flow values to the lower bound.
5.398 + /// \return This function returns false if the initialization
5.399 + /// process found a barrier.
5.400 + void init()
5.401 + {
5.402 + createStructures();
5.403 +
5.404 + for(NodeIt n(_g);n!=INVALID;++n) {
5.405 + _excess->set(n, (*_delta)[n]);
5.406 + }
5.407 +
5.408 + for (ArcIt e(_g);e!=INVALID;++e) {
5.409 + _flow->set(e, (*_lo)[e]);
5.410 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_flow)[e]);
5.411 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_flow)[e]);
5.412 + }
5.413 +
5.414 + // global relabeling tested, but in general case it provides
5.415 + // worse performance for random digraphs
5.416 + _level->initStart();
5.417 + for(NodeIt n(_g);n!=INVALID;++n)
5.418 + _level->initAddItem(n);
5.419 + _level->initFinish();
5.420 + for(NodeIt n(_g);n!=INVALID;++n)
5.421 + if(_tol.positive((*_excess)[n]))
5.422 + _level->activate(n);
5.423 + }
5.424 +
5.425 + /// Initializes the internal data structures.
5.426 +
5.427 + /// Initializes the internal data structures. This functions uses
5.428 + /// greedy approach to construct the initial solution.
5.429 + void greedyInit()
5.430 + {
5.431 + createStructures();
5.432 +
5.433 + for(NodeIt n(_g);n!=INVALID;++n) {
5.434 + _excess->set(n, (*_delta)[n]);
5.435 + }
5.436 +
5.437 + for (ArcIt e(_g);e!=INVALID;++e) {
5.438 + if (!_tol.positive((*_excess)[_g.target(e)] + (*_up)[e])) {
5.439 + _flow->set(e, (*_up)[e]);
5.440 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_up)[e]);
5.441 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_up)[e]);
5.442 + } else if (_tol.positive((*_excess)[_g.target(e)] + (*_lo)[e])) {
5.443 + _flow->set(e, (*_lo)[e]);
5.444 + _excess->set(_g.target(e), (*_excess)[_g.target(e)] + (*_lo)[e]);
5.445 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - (*_lo)[e]);
5.446 + } else {
5.447 + Value fc = -(*_excess)[_g.target(e)];
5.448 + _flow->set(e, fc);
5.449 + _excess->set(_g.target(e), 0);
5.450 + _excess->set(_g.source(e), (*_excess)[_g.source(e)] - fc);
5.451 + }
5.452 + }
5.453 +
5.454 + _level->initStart();
5.455 + for(NodeIt n(_g);n!=INVALID;++n)
5.456 + _level->initAddItem(n);
5.457 + _level->initFinish();
5.458 + for(NodeIt n(_g);n!=INVALID;++n)
5.459 + if(_tol.positive((*_excess)[n]))
5.460 + _level->activate(n);
5.461 + }
5.462 +
5.463 + ///Starts the algorithm
5.464 +
5.465 + ///This function starts the algorithm.
5.466 + ///\return This function returns true if it found a feasible circulation.
5.467 + ///
5.468 + ///\sa barrier()
5.469 + bool start()
5.470 + {
5.471 +
5.472 + Node act;
5.473 + Node bact=INVALID;
5.474 + Node last_activated=INVALID;
5.475 + while((act=_level->highestActive())!=INVALID) {
5.476 + int actlevel=(*_level)[act];
5.477 + int mlevel=_node_num;
5.478 + Value exc=(*_excess)[act];
5.479 +
5.480 + for(OutArcIt e(_g,act);e!=INVALID; ++e) {
5.481 + Node v = _g.target(e);
5.482 + Value fc=(*_up)[e]-(*_flow)[e];
5.483 + if(!_tol.positive(fc)) continue;
5.484 + if((*_level)[v]<actlevel) {
5.485 + if(!_tol.less(fc, exc)) {
5.486 + _flow->set(e, (*_flow)[e] + exc);
5.487 + _excess->set(v, (*_excess)[v] + exc);
5.488 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
5.489 + _level->activate(v);
5.490 + _excess->set(act,0);
5.491 + _level->deactivate(act);
5.492 + goto next_l;
5.493 + }
5.494 + else {
5.495 + _flow->set(e, (*_up)[e]);
5.496 + _excess->set(v, (*_excess)[v] + fc);
5.497 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
5.498 + _level->activate(v);
5.499 + exc-=fc;
5.500 + }
5.501 + }
5.502 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
5.503 + }
5.504 + for(InArcIt e(_g,act);e!=INVALID; ++e) {
5.505 + Node v = _g.source(e);
5.506 + Value fc=(*_flow)[e]-(*_lo)[e];
5.507 + if(!_tol.positive(fc)) continue;
5.508 + if((*_level)[v]<actlevel) {
5.509 + if(!_tol.less(fc, exc)) {
5.510 + _flow->set(e, (*_flow)[e] - exc);
5.511 + _excess->set(v, (*_excess)[v] + exc);
5.512 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
5.513 + _level->activate(v);
5.514 + _excess->set(act,0);
5.515 + _level->deactivate(act);
5.516 + goto next_l;
5.517 + }
5.518 + else {
5.519 + _flow->set(e, (*_lo)[e]);
5.520 + _excess->set(v, (*_excess)[v] + fc);
5.521 + if(!_level->active(v) && _tol.positive((*_excess)[v]))
5.522 + _level->activate(v);
5.523 + exc-=fc;
5.524 + }
5.525 + }
5.526 + else if((*_level)[v]<mlevel) mlevel=(*_level)[v];
5.527 + }
5.528 +
5.529 + _excess->set(act, exc);
5.530 + if(!_tol.positive(exc)) _level->deactivate(act);
5.531 + else if(mlevel==_node_num) {
5.532 + _level->liftHighestActiveToTop();
5.533 + _el = _node_num;
5.534 + return false;
5.535 + }
5.536 + else {
5.537 + _level->liftHighestActive(mlevel+1);
5.538 + if(_level->onLevel(actlevel)==0) {
5.539 + _el = actlevel;
5.540 + return false;
5.541 + }
5.542 + }
5.543 + next_l:
5.544 + ;
5.545 + }
5.546 + return true;
5.547 + }
5.548 +
5.549 + /// Runs the circulation algorithm.
5.550 +
5.551 + /// Runs the circulation algorithm.
5.552 + /// \note fc.run() is just a shortcut of the following code.
5.553 + /// \code
5.554 + /// fc.greedyInit();
5.555 + /// return fc.start();
5.556 + /// \endcode
5.557 + bool run() {
5.558 + greedyInit();
5.559 + return start();
5.560 + }
5.561 +
5.562 + /// @}
5.563 +
5.564 + /// \name Query Functions
5.565 + /// The result of the %Circulation algorithm can be obtained using
5.566 + /// these functions.
5.567 + /// \n
5.568 + /// Before the use of these functions,
5.569 + /// either run() or start() must be called.
5.570 +
5.571 + ///@{
5.572 +
5.573 + /**
5.574 + \brief Returns a barrier
5.575 +
5.576 + Barrier is a set \e B of nodes for which
5.577 + \f[ \sum_{v\in B}-delta(v)<
5.578 + \sum_{e\in\rho(B)}lo(e)-\sum_{e\in\delta(B)}up(e) \f]
5.579 + holds. The existence of a set with this property prooves that a feasible
5.580 + flow cannot exists.
5.581 + \sa checkBarrier()
5.582 + \sa run()
5.583 + */
5.584 + template<class GT>
5.585 + void barrierMap(GT &bar)
5.586 + {
5.587 + for(NodeIt n(_g);n!=INVALID;++n)
5.588 + bar.set(n, (*_level)[n] >= _el);
5.589 + }
5.590 +
5.591 + ///Returns true if the node is in the barrier
5.592 +
5.593 + ///Returns true if the node is in the barrier
5.594 + ///\sa barrierMap()
5.595 + bool barrier(const Node& node)
5.596 + {
5.597 + return (*_level)[node] >= _el;
5.598 + }
5.599 +
5.600 + /// \brief Returns the flow on the arc.
5.601 + ///
5.602 + /// Sets the \c flowMap to the flow on the arcs. This method can
5.603 + /// be called after the second phase of algorithm.
5.604 + Value flow(const Arc& arc) const {
5.605 + return (*_flow)[arc];
5.606 + }
5.607 +
5.608 + /// @}
5.609 +
5.610 + /// \name Checker Functions
5.611 + /// The feasibility of the results can be checked using
5.612 + /// these functions.
5.613 + /// \n
5.614 + /// Before the use of these functions,
5.615 + /// either run() or start() must be called.
5.616 +
5.617 + ///@{
5.618 +
5.619 + ///Check if the \c flow is a feasible circulation
5.620 + bool checkFlow() {
5.621 + for(ArcIt e(_g);e!=INVALID;++e)
5.622 + if((*_flow)[e]<(*_lo)[e]||(*_flow)[e]>(*_up)[e]) return false;
5.623 + for(NodeIt n(_g);n!=INVALID;++n)
5.624 + {
5.625 + Value dif=-(*_delta)[n];
5.626 + for(InArcIt e(_g,n);e!=INVALID;++e) dif-=(*_flow)[e];
5.627 + for(OutArcIt e(_g,n);e!=INVALID;++e) dif+=(*_flow)[e];
5.628 + if(_tol.negative(dif)) return false;
5.629 + }
5.630 + return true;
5.631 + }
5.632 +
5.633 + ///Check whether or not the last execution provides a barrier
5.634 +
5.635 + ///Check whether or not the last execution provides a barrier
5.636 + ///\sa barrier()
5.637 + bool checkBarrier()
5.638 + {
5.639 + Value delta=0;
5.640 + for(NodeIt n(_g);n!=INVALID;++n)
5.641 + if(barrier(n))
5.642 + delta-=(*_delta)[n];
5.643 + for(ArcIt e(_g);e!=INVALID;++e)
5.644 + {
5.645 + Node s=_g.source(e);
5.646 + Node t=_g.target(e);
5.647 + if(barrier(s)&&!barrier(t)) delta+=(*_up)[e];
5.648 + else if(barrier(t)&&!barrier(s)) delta-=(*_lo)[e];
5.649 + }
5.650 + return _tol.negative(delta);
5.651 + }
5.652 +
5.653 + /// @}
5.654 +
5.655 + };
5.656 +
5.657 +}
5.658 +
5.659 +#endif