diff -r 86d9dd6da44d lemon/multicommodity_flow.h --- a/lemon/multicommodity_flow.h Tue Jan 03 11:41:43 2012 +0100 +++ b/lemon/multicommodity_flow.h Mon Dec 13 16:31:18 2021 +0000 @@ -2,7 +2,7 @@ * * This file is a part of LEMON, a generic C++ optimization library. * - * Copyright (C) 2003-2009 + * Copyright (C) 2003-2021 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport * (Egervary Research Group on Combinatorial Optimization, EGRES). * @@ -14,510 +14,924 @@ * express or implied, and with no claim as to its suitability for any * purpose. * + * Additional Copyright file multicommodity_flow.h: Lancaster University (2021) */ -#ifndef LEMON_MULTICOMMODITY_FLOW_H -#define LEMON_MULTICOMMODITY_FLOW_H +#ifndef LEMON_MULTICOMMODITY_FLOW_H_ +#define LEMON_MULTICOMMODITY_FLOW_H_ /// \ingroup approx /// /// \file -/// \brief Approximation algorithms for solving multicommodity flow problems. +/// \brief Fleischer's algorithms for finding approximate multicommodity flows. +#include // min +#include // numeric_limits +#include // random_device, mt19937 +#include #include -#include -#include -#include -#include -#include -#include -#include -// -// TODO: -// -// - The algorithms should provide primal and dual solution as well: -// primal objective value, primal solution (flow), primal solution (flow) for each commodity -// dual objective value, dual solution -// Note that these algorithms are approximation methods, so the found primal and dual obj. values -// are not necessarily optimal, ie. the found primal obj. value could be smaller than the found -// dual objective value and their ratio provides a proof for the approximation factor of both -// the primal and the dual solution. -// -// - compute an upper bound for the optimal primal obj. value, e.g. using -// x_i = min ( out_cap(s_i), in_cap(t_i) ) for each commodity i -// the sum of x_i or the minimum of x_i/d_i ratios can be utilized. -// -// - give back as good approximation guarantee r as possible for which: -// dual = r*primal -// -// - revise the implementations and their API -// +#include "lemon/adaptors.h" +#include "lemon/bfs.h" +#include "lemon/concepts/maps.h" +#include "lemon/dijkstra.h" +#include "lemon/maps.h" +#include "lemon/math.h" +#include "lemon/path.h" +#include "lemon/preflow.h" +#include "lemon/static_graph.h" +#include "lemon/tolerance.h" namespace lemon { - /// \addtogroup approx +/// \brief Default traits class of \ref ApproxMCF class. +/// +/// Default traits class of \ref ApproxMCF class. +/// \tparam GR The type of the digraph. +/// \tparam CAPType The type of the capacity values. +/// \tparam CAPMap The type of the capacity map. +/// \tparam V The type for the flows and costs. +/// It must conform to the \ref concepts::ReadMap "ReadMap" concept. +#ifdef DOXYGEN +template +#else +template< + typename GR, + typename CAPMap = typename GR::template ArcMap, + typename V = double> +#endif +struct ApproxMCNFDefaultTraits { + /// The type of the digraph + typedef GR Digraph; + /// The type of the capacity map + typedef CAPMap CapacityMap; + /// The type of the capacity map values. By default is double. + typedef typename CapacityMap::Value CapacityValue; + + /// \brief The type used for flow values and costs. + /// + /// \c CapacityValue must be convertible to \c Value. + /// By default is \c double. + typedef V Value; + typedef typename Digraph::template ArcMap ValueMap; + + /// The tolerance type used for internal computations + typedef lemon::Tolerance Tolerance; + + /// \brief The path type of the decomposed flows + /// + /// It must conform to the \ref lemon::concepts::Path "Path" concept + typedef lemon::Path Path; + typedef typename lemon::Path::ArcIt PathArcIt; + + /// \brief Instantiates a ValueMap. + /// + /// This function instantiates a \ref ValueMap. + /// \param digraph The digraph for which we would like to define + /// the large cost map. This can be used for flows or otherwise. + static ValueMap* createValueMapPtr( + const Digraph& digraph, + const Value& value = 0) { + return new ValueMap(digraph, value); + } + // Same as above just const + static const ValueMap* createValueMapConstPtr( + const Digraph& digraph, + const Value& value = 0) { + return new ValueMap(digraph, value); + } +}; + +/// \addtogroup approx +/// @{ +/// \brief Fleischer's algorithms for the multicommodity flow problem. +/// +/// Default traits class of \ref ApproxMCF class. +/// \tparam GR The type of the digraph. +/// \tparam TR The type of the traits class. +#ifdef DOXYGEN +template +#else +template< + typename GR, + typename CM = typename GR::template ArcMap, + typename TR = ApproxMCNFDefaultTraits> +#endif +class ApproxMCF { + public: + /// The type of the digraph + typedef typename TR::Digraph Digraph; + /// The type of the capacity map + typedef typename TR::CapacityMap CapacityMap; + /// The type of the capacity values + typedef typename TR::CapacityValue CapacityValue; + + /// \brief The large cost type + /// + /// The large cost type used for internal computations. + /// By default is \c double. + typedef typename TR::Value Value; + /// Arc Map with values of type \ref Value + typedef typename TR::ValueMap ValueMap; + + /// The tolerance type + typedef typename TR::Tolerance Tolerance; + + /// \brief The path type of the found cycles + /// + /// The path type of the found cycles. + /// Using the \ref lemon::ApproxMCNFDefaultTraits "default traits class", + /// it is \ref lemon::Path "Path". + typedef typename TR::Path Path; + /// The arc iterator for paths. + typedef typename TR::PathArcIt PathArcIt; + + /// The \ref lemon::ApproxMCNFDefaultTraits "traits class" of the + /// algorithm + typedef TR Traits; + + /// Types of approximation type. + enum ApproxType { + /// Maximum flow using Fleischer's algorithm. + FLEISCHER_MAX, + /// Maximum concurrent flow using Fleischer's algorithm. + FLEISCHER_MAX_CONCURRENT, + /// Minimum cost concurrent flow using binary search with cost-bounded + /// \ref FLEISCHER_MAX_CONCURRENT + BINARY_SEARCH_MIN_COST + }; + + private: + struct PathData { + Path path; + // The amount of flow through the path + Value value; + // The minimum capacity along all edges in the path + CapacityValue min_cap; + // The actual cost of the path + Value cost; + int hash; + + PathData() {} + PathData(const Path& p, Value v, CapacityValue m, Value c, int h = 0) + : path(p), value(v), min_cap(m), cost(c), hash(h) {} + }; + + int hash(const Path& path) const { + int h = path.length(); + for (PathArcIt a(path); a != INVALID; ++a) + h = (h << 4) ^ (h >> 28) ^ _graph.id(a); + return h; + } + + // Matrix indexed by commodity with all the paths found. i.e. [i][p] where i + // is the commodity index and p is the path index. + typedef std::vector> PathMatrix; + // Vector indexed by commodity with \ref ValueMap pointers. + typedef std::vector ValueMapVector; + typedef std::vector ConstValueMapVector; + + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); + + const Digraph& _graph; + const int _num_nodes; + const int _num_arcs; + // Input capacity map + const CapacityMap& _cap; + + // The number of commodities + int _num_commodities; + // The source nodes for each commodity + std::vector _source; + // The target nodes for each commodity + std::vector _target; + // The demands for each commodity (may be 0, or undefined). + std::vector _demand; + // The cost for each commodity (may be 0, or undefined) + // See \ref ValueMapVector. + ConstValueMapVector _cost; + + /** + * Algorithm parameters + */ + // type of approximation algorithm. See \ref ApproxType. + ApproxType _approx_type; + // Stopping criterion + Value _epsilon; + // Whether the internal heuristics will be applied at the end of the + // algorithms + bool _heuristic; + // Whether there is feasible unsplitable flow + bool _feasible_unsplit; + // Objective function of the primal problem + Value _obj; + // Objective function of the dual problem + Value _dual_obj; + bool _local_cost; + // + Value _delta; + // Expected maximum number of iterations + Value _max_iter; + // + Value _limit; + // Number of iterations + int _iter; + // Tolerance function. See \ref lemon::Tolerance. + Tolerance _tol; + + // For the BINARY_SEARCH_MIN_COST case + // Cost of the final solution. + Value _total_cost; + // Cost bounding factor. + Value _phi; + // Cost bound. + Value _cost_bound; + + // The aggregated flow per arc over all commodities + ValueMap _total_flow; + // The value of the length per arc (the dual solution). + ValueMap _len; + // The flow for each commodity. See \ref ValueMapVector. + ValueMapVector _flow; + // The paths found for different commodities. See \ref PathMatrix. + PathMatrix _paths; + + // Solution Attributes + std::vector _final_paths; + std::vector _final_flows; + std::vector _final_costs; + + public: + /// \brief Default constructor + ApproxMCF(const Digraph& gr, const CapacityMap& cp) + : _graph(gr), + _num_nodes(countNodes(_graph)), + _num_arcs(countArcs(_graph)), + _cap(cp), + _num_commodities(0), + _feasible_unsplit(true), + _obj(0), + _dual_obj(0), + _local_cost(true), + _iter(0), + _cost_bound(0), + _total_flow(_graph), + _len(_graph) {} + + ~ApproxMCF() { + for (int i = 0; i < _num_commodities; ++i) { + delete _flow[i]; + _flow[i] = nullptr; + if (_local_cost) { + delete _cost[i]; + _cost[i] = nullptr; + } + } + _flow.clear(); + _cost.clear(); + } + + /// \name Parameters + /// The parameters of the algorithm can be specified using these + /// functions. + /// @{ - /// \brief Implementation of an approximation algorithm for finding - /// a maximum multicommodity flow. + /// \brief Set single commodity parameters. /// - /// \ref MaxMulticommodityFlow implements Lisa Fleischer's algorithm - /// for finding a maximum multicommodity flow. + /// \param s The source node. + /// \param t The target node. + /// \param d The required amount of flow from node \c s to node \c t + /// (i.e. the supply of \c s and the demand of \c t). + /// \param cost_map The cost arc map. /// - /// \tparam Digraph The type of the digraph the algorithm runs on. - /// \tparam CAP The type of the capacity map. - /// The default map type is \ref concepts::Digraph::ArcMap - /// "GR::ArcMap". - /// \tparam V The value type used in the algorithm. - /// The default value type is \c CAP::Value. -#ifdef DOXYGEN - template < typename GR, - typename CAP, - typename V > -#else - template < typename GR, - typename CAP = typename GR::template ArcMap, - typename V = typename CAP::Value > -#endif - class MaxMulticommodityFlow - { - TEMPLATE_DIGRAPH_TYPEDEFS(GR); - - typedef SimplePath SPath; - typedef typename SPath::ArcIt PathArcIt; - - public: - - /// The type of the digraph the algorithm runs on. - typedef GR Digraph; - /// The type of the capacity map. - typedef CAP CapacityMap; - /// The value type of the algorithm. - typedef V Value; - -#ifdef DOXYGEN - /// The type of the flow map (primal solution). - typedef Digraph::ArcMap FlowMap; - /// The type of the length map (dual soluiton). - typedef Digraph::ArcMap LengthMap; -#else - /// The type of the flow map (primal solution). - typedef typename Digraph::template ArcMap FlowMap; - /// The type of the length map (dual soluiton). - typedef typename Digraph::template ArcMap LengthMap; -#endif - - private: - - // The digraph the algorithm runs on - const Digraph &_g; - // The capacity map - const CapacityMap &_cap; - - // The number of commdities - int _k; - // The source nodes for each commodity - std::vector _s; - // The sink nodes for each commodity - std::vector _t; - // The weight of each commodity - // std::vector _w; + /// \return (*this) + ApproxMCF& addCommodity( + const Node& s, + const Node& t, + const CapacityValue& d, + const ValueMap& cost_map) { + _source.push_back(s); + _target.push_back(t); + if (d > 0) + _demand.push_back(d); + _local_cost = false; + _cost.push_back(&cost_map); + ++_num_commodities; + return *this; + } - // The flow map (the primal solution). - FlowMap _flow; - // The length map (the dual solution). - LengthMap _len; - - struct PathData { - SPath path; - Value value; - Value min_cap; - int hash; - - PathData() {} - PathData(const SPath& p, Value v, Value mc, int h = 0) : - path(p), value(v), min_cap(mc), hash(h) {} - }; - - public: - - /// \brief Constructor. - /// - /// Constructor. - /// - /// \param digraph The digraph the algorithm runs on. - /// \param capacity The capacities of the edges. - MaxMulticommodityFlow( const Digraph &digraph, - const CapacityMap &capacity ) : - _g(digraph), _cap(capacity), _k(0), - _flow(_g), _len(_g) - {} + /// @overload no cost map, optional demand + ApproxMCF& + addCommodity(const Node& s, const Node& t, const CapacityValue& d = 0) { + const ValueMap* cost_map_ptr = Traits::createValueMapConstPtr(_graph); + addCommodity(s, t, d, *cost_map_ptr); + _local_cost = true; + return *this; + } - /// \brief Constructor. - /// - /// Constructor. - /// - /// \param gr The digraph the algorithm runs on. - /// \param capacity The capacities of the edges. - /// \param k The number of commodities. - /// \param sources A vector containing the source nodes - /// (its size must be at least \c k). - /// \param sinks A vector containing the sink nodes - /// (its size must be at least \c k). - MaxMulticommodityFlow( const Digraph &gr, - const CapacityMap &capacity, - int k, - const std::vector& sources, - const std::vector& sinks ) : - _g(gr), _cap(capacity), - _k(k), _s(sources), _t(sinks), - _flow(_g), _len(_g) - { - _s.resize(k); - _t.resize(k); - } + /// @} - /// Destructor. - ~MaxMulticommodityFlow() {} + /// \name Execution control - /// \brief Add a new commodity. - /// - /// This function adds a new commodity with the given source and - /// sink nodes. - void addCommodity(const Node& s, const Node& t) { - ++_k; - _s.push_back(s); - _t.push_back(t); - } - - /// \name Execution control - - /// @{ + /// @{ - /// \brief Run the algorithm for finding a maximum multicommodity flow. - /// - /// This function runs the algorithm for finding a maximum multicommodity - /// flow. - /// \param r The approximation factor. It must be greater than 1. - /// \return The effective approximation factor (at least 1), i.e. the - /// ratio of the values of the found dual and primal solutions. - /// Therefore both the primal and dual solutions are r-approximate. - /// - /// \note The smaller parameter \c r is, the slower the algorithm runs, - /// but it might find better solution. According to our tests, using - /// the default value the algorithm performs well and it usually manage - /// to find the exact optimum due to heuristic improvements. - Value run(Value r = 2.0, Value opt = 0) { - - #define MMF_FLEISCHER_IMPR - #define MMF_HEURISTIC_IMPR - //#define MMF_EXP_UPDATE - //#define MMF_DEBUG_OUTPUT - - // Initialize epsilon and delta parameters - Tolerance ftol, ctol; - if (!ftol.less(1.0, r)) return 0; -#ifdef MMF_FLEISCHER_IMPR - #ifdef MMF_EXP_UPDATE - Value epsilon = (r + 1 - 2 * pow(r, 0.5)) / (r - 1); - #else - Value epsilon = (2 * r + 1 - pow(8 * r + 1, 0.5)) / (2 * r); - #endif -#else - #ifdef MMF_EXP_UPDATE - Value epsilon = (2 * r + 1 - pow(8 * r + 1, 0.5)) / (2 * r); - #else - Value epsilon = 1 - pow(r, -0.5); - #endif -#endif - -#ifdef MMF_DEBUG_OUTPUT - std::cout << "r = " << r << "; epsilon = " << epsilon << "\n"; -#endif - - int node_num = countNodes(_g); - Value delta = pow(1 + epsilon, 1 - 1/epsilon) / - pow(node_num, 1/epsilon); - ctol.epsilon(delta * 0.0001); - - // Initialize flow values and lengths - for (ArcIt a(_g); a != INVALID; ++a) { - _flow[a] = 0; - _len[a] = delta; - } - - // Check commodities - for (int i = 0; i < _k; ++i) { - if (_s[i] == _t[i] || !bfs(_g).run(_s[i], _t[i])) { - _s[i] = _s[_k - 1]; - _t[i] = _t[_k - 1]; - --_k; - } - } - if (_k <= 0) return 0; - + /// \brief Run the algorithm for finding a multicommodity flow. + /// + /// This function runs the algorithm for finding a maximum/minimum cost + /// multicommodity flow. + /// + /// \param at The type of approximation algorithm to use. @see \ref + /// ApproxMCF::ApproxType. + /// \param heuristic Whether the internal heuristics will be applied at the + // end of Fleischer's algorithm. + /// \param epsilon The stopping criterion. Must be positive and larger than + /// 0. The smaller the value, the more accurate the result, but the longer the + /// execution time. By default is 0.1. + void run( + const ApproxType& at = FLEISCHER_MAX, + const bool& heuristic = true, + const Value& epsilon = 0.1) { + _approx_type = at; + _epsilon = epsilon; + _heuristic = heuristic; + check(); + switch (_approx_type) { + case FLEISCHER_MAX: + runFleischerMax(); + break; + case BINARY_SEARCH_MIN_COST: + runBinarySearchMinCost(); + break; + case FLEISCHER_MAX_CONCURRENT: + runFleischerMaxConcurrent(); + break; + } + } - // Successive improvement along shortest paths - std::vector > paths(_k); -#ifdef MMF_FLEISCHER_IMPR - Value limit = delta * (1 + epsilon); - int i = 0; -#else - std::vector curr_paths(_k); -#endif - SPath curr_path; - Value curr_dist; - int iter = 0; - while (true) { -#ifdef MMF_FLEISCHER_IMPR - // Find the shortest path for the i-th commodity - // (w.r.t the current length function) - dijkstra(_g, _len).path(curr_path).dist(curr_dist) - .run(_s[i], _t[i]); - - // Check if the path can be used and increase i or the - // length limit if necessary - if (!ctol.less(curr_dist, limit)) { - if (++i == _k) { - if (!ctol.less(limit, 1.0)) break; - i = 0; - limit *= (1 + epsilon); - if (ctol.less(1.0, limit)) limit = 1.0; - } - continue; - } -#else - // Find the shortest path - Value min_dist = 1.0; - int i = -1; - for (int j = 0; j < _k; ++j) { - if ( !dijkstra(_g, _len).path(curr_paths[j]).dist(curr_dist). - run(_s[j], _t[j]) ) { - std::cerr << "\n\nERROR!!! "<::infinity(); - for (PathArcIt a(curr_path); a != INVALID; ++a) { - if (_cap[a] < gamma) gamma = _cap[a]; - } - paths[i].push_back(PathData(curr_path, gamma, gamma, h)); - } else { - gamma = paths[i][j].min_cap; - paths[i][j].value += gamma; - } - - // Modify the arc lengths along the current path -#ifdef MMF_EXP_UPDATE - for (PathArcIt a(curr_path); a != INVALID; ++a) - _len[a] *= exp(epsilon * gamma / _cap[a]); -#else - for (PathArcIt a(curr_path); a != INVALID; ++a) - _len[a] *= 1 + epsilon * gamma / _cap[a]; -#endif - } - - // Setting flow values - for (int i = 0; i < _k; ++i) { - for (int j = 0; j < int(paths[i].size()); ++j) { - Value val = paths[i][j].value; - for (PathArcIt a(paths[i][j].path); a != INVALID; ++a) - _flow[a] += val; + // Allocate/reset flow per commodity and arc + _flow.resize(_num_commodities); + for (int i = 0; i < _num_commodities; ++i) + if (!_flow[i]) { + _flow[i] = Traits::createValueMapPtr(_graph); + } else { + for (ArcIt a(_graph); a != INVALID; ++a) { + (*_flow[i])[a] = 0; } } - // Scaling the solution - Value lambda = 0; - for (ArcIt a(_g); a != INVALID; ++a) { - if (_flow[a] / _cap[a] > lambda) lambda = _flow[a] / _cap[a]; + // Initialise length and flow map + for (ArcIt a(_graph); a != INVALID; ++a) { + _len[a] = _delta; + _total_flow[a] = 0; + } + + // Initialise phi for cost bounding + if (_approx_type == BINARY_SEARCH_MIN_COST && _cost_bound > 0) + _phi = _delta / _cost_bound; + else + _phi = 0; + } + + // Preprocessing checks to ensure that: + // + // 1. Commodities are well-defined + // 2. Input sizes are consistent + // 3. Concurrent case: Feasible split flow exists for each commodity that + // satisfies demand. In the case this doesn't pass, will throw an exception. + // + // TODO(): Group commodities by source, so we only have to solve the shortest + // path once for each group (and length values). + void check() { + LEMON_ASSERT(_tol.positive(_epsilon), "Epsilon must be greater than 0"); + + // Remove commodities with source = sink or unreachable + for (int i = 0; i < _num_commodities; ++i) { + if (_source[i] == _target[i] || + !bfs(_graph).run(_source[i], _target[i])) { + _source[i] = _source[_num_commodities - 1]; + _target[i] = _target[_num_commodities - 1]; + --_num_commodities; + } + } + + const int& s_size = _source.size(); + const int& t_size = _target.size(); + const int& d_size = _demand.size(); + + switch (_approx_type) { + case FLEISCHER_MAX: { + LEMON_ASSERT(_tol.positive(_num_commodities), "No commodities"); + // Check consistency of vector sizes + LEMON_ASSERT( + (_num_commodities == s_size && s_size == t_size && + !_tol.positive(d_size)), + "All vectors should be of size equal to the number of commodities, " + "and have no demand."); + break; } - if (ftol.positive(lambda)) { - lambda = 1 / lambda; - for (ArcIt a(_g); a != INVALID; ++a) _flow[a] *= lambda; - for (int i = 0; i < _k; ++i) { - for (int j = 0; j < int(paths[i].size()); ++j) - paths[i][j].value *= lambda; + case FLEISCHER_MAX_CONCURRENT: + case BINARY_SEARCH_MIN_COST: { + // Check consistency of vector sizes + LEMON_ASSERT( + (_num_commodities == s_size && s_size == t_size && + s_size == d_size), + "All vectors should be of size equal to the number of commodities"); + // Remove commodities with 0 demand. + for (int i = 0; i < _num_commodities; ++i) { + if (!_tol.positive(_demand[i])) { + _source[i] = _source[_num_commodities - 1]; + _target[i] = _target[_num_commodities - 1]; + --_num_commodities; + } + } + LEMON_ASSERT(_tol.positive(_num_commodities), "No commodities"); + + // Check if feasible routing of demand exists for each commodity + for (int i = 0; i < _num_commodities; ++i) { + // Only check commodity if the demand in non-zero + // Check if feasible split flow exists + // Total capacity outgoing at the source + CapacityValue total_cap_source = 0; + CapacityValue max_cap_source = 0; + for (OutArcIt a(_graph, _source[i]); a != INVALID; ++a) { + total_cap_source += _cap[a]; + if (max_cap_source < _cap[a]) + max_cap_source = _cap[a]; + } + + // Total capacity incoming at the target + CapacityValue total_cap_target = 0; + CapacityValue max_cap_target = 0; + for (InArcIt a(_graph, _target[i]); a != INVALID; ++a) { + total_cap_target += _cap[a]; + if (max_cap_target < _cap[a]) + max_cap_target = _cap[a]; + } + + // Check if enough capacity to satisfy the demand + LEMON_ASSERT( + (total_cap_source >= _demand[i] && + total_cap_target >= _demand[i]), + "No feasible (split) flow to satisfy demand for commodity" + + std::to_string(i)); + + // Check if feasible unsplittable flow exists that satisfies demand + if (max_cap_source < _demand[i] || max_cap_target < _demand[i]) + _feasible_unsplit = false; } } + } + } -#ifdef MMF_DEBUG_OUTPUT - Value value1 = 0; - for (int i = 0; i < _k; ++i) { - for (InArcIt a(_g, _t[i]); a != INVALID; ++a) value1 += _flow[a]; - for (OutArcIt a(_g, _t[i]); a != INVALID; ++a) value1 -= _flow[a]; + // Run standard max concurrent algorithm iteratively updating the budget, or + // cost bound, a la binary search to obtain an approximation to the + // multicommodity flow minimum cost + void runBinarySearchMinCost() { + // Run standard max concurrent + runFleischerMaxConcurrent(); + + Value lower_bound = 0, upper_bound = getTotalCost(); + + // Best solution attributes + Value curr_best_cost = std::numeric_limits::infinity(); + ValueMapVector curr_best_flow; + + while ((upper_bound - lower_bound) / upper_bound > _epsilon && + upper_bound >= lower_bound) { + // Run Max Concurrent algorithm with cost bound + _cost_bound = (upper_bound + lower_bound) / 2; + runFleischerMaxConcurrent(); + // Update upper/lower bounds + if (_total_cost < _cost_bound) + upper_bound = _cost_bound; + else + lower_bound = _cost_bound; + // Save current solution if appropriate + if (_total_cost < curr_best_cost) { + curr_best_cost = _total_cost; + curr_best_flow = _flow; } - Value max_opt = value1 / (1 - epsilon); + } + // Update flow + _total_cost = curr_best_cost; + _flow = curr_best_flow; + } - std::cout << "\n"; - std::cout << "Computation finished. Iterations: " << iter << "\n"; - std::cout << "Solution value: " << value1; - if (opt > 0) std::cout << " \t" << value1 / opt * 100 << "%"; - std::cout << std::endl; -#else - opt = opt; // avoid warning -#endif + void runFleischerMaxConcurrent() { + // Successive improvement along shortest paths + init(); + Path curr_path; + Value curr_dist; + while (true) { + // Find the shortest path for the i-th commodity + // (w.r.t the current length function) + for (int i = 0; i < _num_commodities; ++i) { + // Route demand + CapacityValue d_i = _demand[i]; + while (true) { // while (phase) + dijkstra(_graph, _len) + .path(curr_path) + .dist(curr_dist) + .run(_source[i], _target[i]); + // Get + const CapacityValue& routed_flow = updatePaths(curr_path, i, d_i); + // Update demand + d_i -= routed_flow; + updateLen(curr_path, routed_flow, i); + updateDualObjective(); + if (d_i <= 0 || !_tol.less(_dual_obj, 1)) + break; + } // end while (phase) + } + ++_iter; + updateDualObjective(); + if (!_tol.less(_dual_obj, 1)) + break; + } // end while (main) + + // Clean up + scaleFinalFlowsMaxConcurrent(); + + if (_heuristic) + randomisedRounding(); + + setFlowPerCommodity(); + setFinalCost(); + saveFinalPaths(); + } + + void runFleischerMax() { + // Successive improvement along shortest paths + int i = 0; + Path curr_path; + Value curr_dist; + init(); + while (true) { // start while + dijkstra(_graph, _len) + .path(curr_path) + .dist(curr_dist) + .run(_source[i], _target[i]); -#ifdef MMF_DEBUG_OUTPUT - Value tmp = 0; - for (ArcIt a(_g); a != INVALID; ++a) - tmp += _cap[a] * _len[a]; - std::cout << "Final dual: " << tmp << "\n"; -#endif + // Check if the path can be used and increase i or the length limit if + // necessary + if (!_tol.less(curr_dist, _limit)) { + if (++i == _num_commodities) { + if (!_tol.less(_limit, 1.0)) + break; + i = 0; + _limit *= (1 + _epsilon); + if (_tol.less(1.0, _limit)) + _limit = 1.0; + } + continue; + } + const CapacityValue& routed_flow = updatePaths(curr_path, i); + // Update dual values + updateLen(curr_path, routed_flow, i); + ++_iter; + } // end while + + scaleFinalFlowsMaxFlow(); + + if (_heuristic) { + runHeurMax1(); + runHeurMax2(); + } + + setFlowPerCommodity(); + + // Get dual objective + _dual_obj = getDualObjective(); + // Get objective value + _obj = 0; + for (int i = 0; i < _num_commodities; ++i) { + for (InArcIt a(_graph, _target[i]); a != INVALID; ++a) + _obj += (*_flow[i])[a]; + for (OutArcIt a(_graph, _target[i]); a != INVALID; ++a) + _obj -= (*_flow[i])[a]; + } + } -#ifdef MMF_HEURISTIC_IMPR - // Heuristic improvement 1: improve the solution along one of the - // found paths if it is possible - int opi = 0; - for (int i = 0; i < _k; ++i) { - for (int j = 0; j < int(paths[i].size()); ++j) { - Value d = std::numeric_limits::infinity(); - for (PathArcIt a(paths[i][j].path); a != INVALID; ++a) - if (_cap[a] - _flow[a] < d) d = _cap[a] - _flow[a]; - if (ftol.positive(d)) { - paths[i][j].value += d; - for (PathArcIt a(paths[i][j].path); a != INVALID; ++a) - _flow[a] += d; - ++opi; + // Randomised rounding - if flow is split over multiple paths, randomised + // rounding randomly rounds up the largest proportion of flow/cost. + // + // Only runs if flow can be routed without splitting. + // If costs are given, to make expensive paths less likely to occur, we use, + // for each path j a value[j] = flow / cost for the discrete distribution. + // This distribution samples path j with a probability value[j] / sum values + void randomisedRounding() { + // Only works if flow can be routed without splitting. + if (!_feasible_unsplit) + return; + for (int i = 0; i < _num_commodities; ++i) { + const int& p_size = static_cast(_paths[i].size()); + Value total_flow_i = 0; + std::vector values(p_size, 0); + + for (int j = 0; j < p_size; ++j) { + const Value& val = _paths[i][j].value; + const Value& cost = _paths[i][j].cost; + total_flow_i += val; + values[j] = val; + if (_tol.positive(cost)) + values[j] /= cost; + } + // Get random index wrt values[j] / sum of values + int chosen_p_idx; + std::random_device rd; + std::mt19937 gen(rd()); + std::discrete_distribution<> d(values.begin(), values.end()); + chosen_p_idx = d(gen); + // Check routing all demand through chosen_p_idx remains capacity feasible + bool capacity_feasible = true; + // cost is all demand routed through here + for (PathArcIt a(_paths[i][chosen_p_idx].path); a != INVALID; ++a) { + if (_total_flow[a] - _paths[i][chosen_p_idx].value + _demand[i] > + _cap[a]) { + capacity_feasible = false; + break; + } + } + if (capacity_feasible) { + for (int j = 0; j < p_size; ++j) { + // Update total flows and path flow + if (j == chosen_p_idx) { // path j routes demand + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) { + _total_flow[a] = _total_flow[a] - _paths[i][j].value + _demand[i]; + } + _paths[i][j].value = _demand[i]; + } else { // path j flow is 0 + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) { + _total_flow[a] = _total_flow[a] - _paths[i][j].value; + } + _paths[i][j].value = 0; } } } + } + } - // Heuristic improvement 2: improve the solution along new paths - // if it is possible - int npi = 0; - BoolArcMap filter(_g); - for (ArcIt a(_g); a != INVALID; ++a) - filter[a] = ftol.positive(_cap[a] - _flow[a]); - FilterArcs sub_graph(_g, filter); - SPath p; - for (int i = 0; i < _k; ++i) { - while (bfs(sub_graph).path(p).run(_s[i], _t[i])) { - Value d = std::numeric_limits::infinity(); - for (PathArcIt a(p); a != INVALID; ++a) - if (_cap[a] - _flow[a] < d) d = _cap[a] - _flow[a]; - paths[i].push_back(PathData(p, d, d)); - for (PathArcIt a(p); a != INVALID; ++a) { - _flow[a] += d; - filter[a] = ftol.positive(_cap[a] - _flow[a]); - } - ++npi; + // Heuristic improvement 1 (for FLEISCHER_MAX only). + // + // Attempts to improve the solution along one of the found paths if it is + // possible + void runHeurMax1() { + for (int i = 0; i < _num_commodities; ++i) { + for (int j = 0; j < static_cast(_paths[i].size()); ++j) { + CapacityValue d = std::numeric_limits::infinity(); + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) + if (_cap[a] - _total_flow[a] < d) + d = _cap[a] - _total_flow[a]; + if (_tol.positive(d)) { + _paths[i][j].value += d; + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) + _total_flow[a] += d; + } + } + } + } + + // Heuristic improvement 2 (for FLEISCHER_MAX only). + // + // Attempts to improve the solution along new paths + void runHeurMax2() { + BoolArcMap filter(_graph); + for (ArcIt a(_graph); a != INVALID; ++a) + filter[a] = _tol.positive(_cap[a] - _total_flow[a]); + FilterArcs sub_graph(_graph, filter); + Path p; + for (int i = 0; i < _num_commodities; ++i) { + while (bfs(sub_graph).path(p).run(_source[i], _target[i])) { + CapacityValue d = std::numeric_limits::infinity(); + for (PathArcIt a(p); a != INVALID; ++a) + if (_cap[a] - _total_flow[a] < d) + d = _cap[a] - _total_flow[a]; + _paths[i].push_back(PathData(p, d, d, 0, 0)); + for (PathArcIt a(p); a != INVALID; ++a) { + _total_flow[a] += d; + filter[a] = _tol.positive(_cap[a] - _total_flow[a]); } } -#endif - -#ifdef MMF_DEBUG_OUTPUT - Value value2 = 0; - for (int i = 0; i < _k; ++i) { - for (InArcIt a(_g, _t[i]); a != INVALID; ++a) value2 += _flow[a]; - for (OutArcIt a(_g, _t[i]); a != INVALID; ++a) value2 -= _flow[a]; - } - epsilon = 1 - value2 / max_opt; + } + } - #ifdef MMF_HEURISTIC_IMPR - std::cout << "\nHeuristics applied. Old impr: " << opi << " new impr: " << npi << "\n"; - #endif - std::cout << "Solution value: " << value2; - if (opt > 0) std::cout << " \t" << value2 / opt * 100 << "%"; - std::cout << std::endl << std::endl; -#endif + // Update paths for current commodity and return flow routed + Value updatePaths( + const Path& curr_path, + const int& i, + const CapacityValue& d_i = + std::numeric_limits::infinity()) { + int h = hash(curr_path); + int j; // path index + for (j = 0; j < static_cast(_paths[i].size()); ++j) { + bool b = (_paths[i][j].hash == h) && + (_paths[i][j].path.length() == curr_path.length()); + if (b) { + for (PathArcIt a1(_paths[i][j].path), a2(curr_path); + b && a1 != INVALID && a2 != INVALID; + ++a1, ++a2) + b = _graph.id(a1) == _graph.id(a2); + if (b) + break; + } + } + const bool path_seen = !(j == static_cast(_paths[i].size())); - // TODO: better return value - return r; + // Set/modify the flow value for the found path + Value routed_flow = 0; + CapacityValue min_cap; + if (!path_seen) { + // Add new path + Value path_cost = 0; + min_cap = std::numeric_limits::infinity(); + for (PathArcIt a(curr_path); a != INVALID; ++a) { + path_cost += (*_cost[i])[a]; + if (_cap[a] < min_cap) + min_cap = _cap[a]; + } + routed_flow = std::min(min_cap, d_i); + // Save path + _paths[i].push_back( + PathData(curr_path, routed_flow, min_cap, path_cost, h)); + } else { + // Update path flow + min_cap = _paths[i][j].min_cap; + routed_flow = std::min(min_cap, d_i); + _paths[i][j].value += routed_flow; } + return routed_flow; + } - /// @} + // Modify the arc lengths along the current path + void updateLen(const Path& p, const Value& f, const int& i) { + Value path_cost = 0; + for (PathArcIt a(p); a != INVALID; ++a) { + const Value& c_a = (*_cost[i])[a]; + _len[a] *= 1 + _epsilon * f / _cap[a] + c_a * _phi; + path_cost += c_a; + } + // Update phi if needed + if (_tol.positive(_phi)) + _phi *= (1 + (_epsilon * path_cost * f) / _cost_bound); + } - private: + void updateDualObjective() { + _dual_obj = 0; + for (ArcIt a(_graph); a != INVALID; ++a) + _dual_obj += _cap[a] * _len[a]; + } - int hash(const SPath& path) { - int h = path.length(); - for (PathArcIt a(path); a != INVALID; ++a) - h = (h<<4) ^ (h>>28) ^ _g.id(a); - return h; + void setFinalCost() { + _total_cost = 0; + for (int i = 0; i < _num_commodities; ++i) { + for (int j = 0; j < static_cast((_paths[i].size())); ++j) { + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) + _total_cost += (*_cost[i])[a] * (*_flow[i])[a]; + } } + } - public: - - /// \name Query Functions - /// The result of the algorithm can be obtained using these - /// functions.\n - /// \ref run() must be called before using them. - - /// @{ + // Scale flows in the \ref FLEISCHER_MAX_CONCURRENT. + // Flows are scaled by the number of iterations. + void scaleFinalFlowsMaxConcurrent() { + int scaler = _iter; + for (int i = 0; i < _num_commodities; ++i) { + for (int j = 0; j < static_cast(_paths[i].size()); ++j) { + _paths[i][j].value /= scaler; + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) { + _total_flow[a] += _paths[i][j].value; + } + } + } + _obj = (_iter) / scaler; + } - /// \brief Return the flow value on the given arc. - /// - /// This function returns the flow value on the given arc. - /// - /// \pre \ref run() must be called before using this function. - Value flow(const Arc& arc) const { - return _flow[arc]; + // Update total flow map + void updateTotalFlow() { + for (int i = 0; i < _num_commodities; ++i) { + for (int j = 0; j < static_cast(_paths[i].size()); ++j) { + CapacityValue val = _paths[i][j].value; + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) + _total_flow[a] += val; + } } + } - /// \brief Return a const reference to an arc map storing the - /// found flow. - /// - /// This function returns a const reference to an arc map storing - /// the found flow. - /// - /// \pre \ref run() must be called before using this function. - const FlowMap& flowMap() const { - return _flow; + // Scale flows per commodity for the FLEISCHER_MAX case. + // scaler = max {total_flow_a / cap_a} such that all flows are capacity + // feasible + void scaleFinalFlowsMaxFlow() { + updateTotalFlow(); + Value scaler = 0; + for (ArcIt a(_graph); a != INVALID; ++a) { + if (_total_flow[a] / _cap[a] > scaler) + scaler = _total_flow[a] / _cap[a]; + } + if (_tol.positive(scaler)) { + for (ArcIt a(_graph); a != INVALID; ++a) { + _total_flow[a] /= scaler; + } + for (int i = 0; i < _num_commodities; ++i) { + for (int j = 0; j < static_cast(_paths[i].size()); ++j) + _paths[i][j].value /= scaler; + } } + } + + // Set flow decomposed solution for exposure + void setFlowPerCommodity() { + // Setting flow values using saved paths + for (int i = 0; i < _num_commodities; ++i) { + for (int j = 0; j < static_cast(_paths[i].size()); ++j) { + Value val = _paths[i][j].value; + for (PathArcIt a(_paths[i][j].path); a != INVALID; ++a) + (*_flow[i])[a] += val; + } + } + } + + void saveFinalPaths() { + for (int i = 0; i < _num_commodities; ++i) { + const int& p_size = static_cast(_paths[i].size()); + bool saved_path = false; + Value path_flow = 0; + Value path_cost = std::numeric_limits::infinity(); + Path path; - /// \brief Return the length of the given arc (the dual value). - /// - /// This function returns the length of the given arc (the dual value). - /// - /// \pre \ref run() must be called before using this function. - Value length(const Arc& arc) const { - return _len[arc]; + for (int j = 0; j < p_size; ++j) { + const Value& v = _paths[i][j].value; + const Value& c = _paths[i][j].cost; + if (_tol.positive(v)) { + if (saved_path) + std::cout << "[WARNING] multiple paths with positive flow" + << std::endl; + // Keep path with largest flow + if (path_flow < v && c < path_cost) { + path = _paths[i][j].path; + path_flow = v; + path_cost = c; + saved_path = true; + } + } + } + _final_paths[i] = path; + _final_flows[i] = path_flow; + _final_costs[i] = path_cost; } + } + + public: + const Value& flow(const int& k, const Arc& arc) const { + return (*_flow[k])[arc]; + } - /// \brief Return a const reference to an arc map storing the - /// found lengths (the dual solution). - /// - /// This function returns a const reference to an arc map storing - /// the found lengths (the dual solution). - /// - /// \pre \ref run() must be called before using this function. - const LengthMap& lengthMap() const { - return _len; - } + const Value& flow(const int& k, const int& arc_id) const { + return (*_flow[k])[_graph.arcFromId(arc_id)]; + } + + const Value& length(const Arc& arc) const { return _len[arc]; } + + const ValueMap& lengthMap() const { return _len; } + + const Value& getObjective() const { return _obj; } + + const Value& getDualObjective() const { return _dual_obj; } + + const Value& getTotalCost() const { return _total_cost; } - /// @} + const ValueMap& getFlowMap(const int& k) const { return *_flow[k]; } - }; //class MaxMulticommodityFlow + const Path& getPath(const int& k) const { return _final_paths[k]; } + + const Value& getPathCost(const int& k) const { return _final_costs[k]; } - ///@} + const Value& getPathFlow(const int& k) const { return _final_flows[k]; } +}; // namespace lemon +/// @} -} //namespace lemon +} // namespace lemon -#endif //LEMON_MULTICOMMODITY_FLOW_H +#endif // LEMON_MULTICOMMODITY_FLOW_H_ diff -r 86d9dd6da44d test/multicommodity_flow_test.cc --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/multicommodity_flow_test.cc Mon Dec 13 16:31:18 2021 +0000 @@ -0,0 +1,217 @@ +/* -*- mode: C++; indent-tabs-mode: nil; -*- + * + * This file is a part of LEMON, a generic C++ optimization library. + * + * Copyright (C) 2003-2021 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport + * (Egervary Research Group on Combinatorial Optimization, EGRES). + * + * Permission to use, modify and distribute this software is granted + * provided that this copyright notice appears in all copies. For + * precise terms see the accompanying LICENSE file. + * + * This software is provided "AS IS" with no warranty of any kind, + * express or implied, and with no claim as to its suitability for any + * purpose. + * + */ + +#include +#include +#include +#include +#include + +#include +#include + +#include "test/test_tools.h" + +char test_lgf[] = + "@nodes\n" + "label\n" + "0\n" + "1\n" + "2\n" + "3\n" + "@arcs\n" + " label cap cost solution_1 solution_2 solution_1_max " + "solution_2_max\n" + "0 1 0 10.0 1.0 10.0 0.0 10.0 0.0\n" + "0 2 1 10.0 1.0 0.0 10.0 8.3 1.7\n" + "1 3 2 10.0 1.0 10.0 0.0 10.0 0.0\n" + "2 3 3 10.0 1000.0 0.0 0.0 8.3 0.0\n" + "@attributes\n" + "source_1 0\n" + "target_1 3\n" + "source_2 0\n" + "target_2 2\n"; + +using namespace lemon; + +// Returns string for test-naming +template +std::string getEdgeName(const GR& graph, const A& edge) { + return std::to_string(graph.id(graph.source(edge))) + "->" + + std::to_string(graph.id(graph.target(edge))); +} + +template +void checkApproxMCF() { + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); + + Digraph graph; + Node s1, t1, s2, t2; + typedef typename Digraph::template ArcMap DoubleMap; + DoubleMap cap(graph), cost(graph), sol1(graph), sol2(graph), sol1_max(graph), + sol2_max(graph); + Tolerance tol; + tol.epsilon(1e-1); // 1 decimal place + + std::istringstream input(test_lgf); + DigraphReader(graph, input) + .arcMap("cap", cap) + .arcMap("cost", cost) + .arcMap("solution_1", sol1) // flow in the solution for commodity 1 + .arcMap("solution_2", sol2) // flow in the solution for commodity 1 + .arcMap( + "solution_1_max", + sol1_max) // flow in the solution for commodity (max-flow) 1 + .arcMap( + "solution_2_max", + sol2_max) // flow in the solution for commodity (max-flow) 1 + .node("source_1", s1) // source node commodity 1 + .node("target_1", t1) // target node commodity 1 + .node("source_2", s2) // source node commodity 2 + .node("target_2", t2) // target node commodity 2 + .run(); + + // Max + { + ApproxMCF max(graph, cap); + max.addCommodity(s1, t1).addCommodity(s2, t2).run(); + for (ArcIt e(graph); e != INVALID; ++e) { + const std::string test_name = + "[FLEISCHER_MAX] flow on " + getEdgeName(graph, e); + // Check flow for commodity 1 (index 0) to 1 decimal place + check( + !tol.different(max.flow(0, e), sol1_max[e]), + test_name + " commodity 1"); + // Check flow for commodity 2 (index 1) to 1 decimal place + check( + !tol.different(max.flow(1, e), sol2_max[e]), + test_name + " commodity 2"); + } + } + + // Max Concurrent + { + ApproxMCF max(graph, cap); + max.addCommodity(s1, t1, 10.0) + .addCommodity(s2, t2, 10.0) + .run(ApproxMCF::FLEISCHER_MAX_CONCURRENT); + for (ArcIt e(graph); e != INVALID; ++e) { + const std::string test_name = "[FLEISCHER_MAX_CONCURRENT] flow on " + + getEdgeName(graph, e); + // Check commodity 1 (index 0) + check(max.flow(0, e) == sol1[e], test_name + " commodity 1"); + // Check commodity 2 (index 1) + check(max.flow(1, e) == sol2[e], test_name + " commodity 2"); + } + } + + // Min cost + { + ApproxMCF max(graph, cap); + max.addCommodity(s1, t1, 10.0, cost) + .addCommodity(s2, t2, 10.0, cost) + .run(ApproxMCF::BINARY_SEARCH_MIN_COST); + for (ArcIt e(graph); e != INVALID; ++e) { + const std::string test_name = + "[BINARY_SEARCH_MIN_COST] flow on " + getEdgeName(graph, e); + // Check commodity 1 (index 0) + check(max.flow(0, e) == sol1[e], test_name + " commodity 1"); + // Check commodity 2 (index 1) + check(max.flow(1, e) == sol2[e], test_name + " commodity 2"); + } + // Check total cost + check(max.getTotalCost() == 30.0, "min-cost failed"); + } + + // Max Concurrent - 3 commodities + { + ApproxMCF max(graph, cap); + max.addCommodity(s1, t1, 5.0, cost) + .addCommodity(s1, t2, 10.0, cost) + .addCommodity(s1, t1, 5.0, cost) + .run(ApproxMCF::FLEISCHER_MAX_CONCURRENT); + + const DoubleMap& flow = max.getFlowMap(0); + check( + max.getTotalCost() == 30.0, + "3 commodities FLEISCHER_MAX_CONCURRENT failed"); + } +} + +template +void checkApproxMCFDisconnected() { + TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); + + Digraph graph; + typedef typename Digraph::template ArcMap DoubleMap; + DoubleMap cap(graph), cost(graph), sol1(graph), sol2(graph); + + Node n0 = graph.addNode(); + Node n1 = graph.addNode(); + Node n2 = graph.addNode(); + Node n3 = graph.addNode(); + + Arc a; + a = graph.addArc(n0, n1); + cap[a] = 10.0; + cost[a] = 1.0; + sol1[a] = 10.0; + sol2[a] = 0.0; + + a = graph.addArc(n2, n3); + cap[a] = 10.0; + cost[a] = 1.0; + sol1[a] = 0.0; + sol2[a] = 10.0; + + // Max + { + ApproxMCF max(graph, cap); + max.addCommodity(n0, n1).addCommodity(n2, n3).run(); + for (ArcIt e(graph); e != INVALID; ++e) { + const std::string test_name = + "[FLEISCHER_MAX] flow on " + getEdgeName(graph, e); + // Check commodity 1 + check(max.flow(0, e) == sol1[e], test_name + " commodity 1"); + // Check commodity 2 + check(max.flow(1, e) == sol2[e], test_name + " commodity 2"); + } + } + // Max Concurrent + { + ApproxMCF max(graph, cap); + max.addCommodity(n0, n1, 10.0, cost) + .addCommodity(n2, n3, 10.0, cost) + .run(ApproxMCF::FLEISCHER_MAX_CONCURRENT); + for (ArcIt e(graph); e != INVALID; ++e) { + const std::string test_name = "[FLEISCHER_MAX_CONCURRENT] flow on " + + getEdgeName(graph, e); + // Check commodity 1 + check(max.flow(0, e) == sol1[e], test_name + " commodity 1"); + // Check commodity 2 + check(max.flow(1, e) == sol2[e], test_name + " commodity 2"); + } + } +} + +int main() { + checkApproxMCF(); + checkApproxMCF(); + checkApproxMCFDisconnected(); + checkApproxMCFDisconnected(); +}