Bug fix in CapacityScaling.
3 * This file is a part of LEMON, a generic C++ optimization library
5 * Copyright (C) 2003-2008
6 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 * (Egervary Research Group on Combinatorial Optimization, EGRES).
9 * Permission to use, modify and distribute this software is granted
10 * provided that this copyright notice appears in all copies. For
11 * precise terms see the accompanying LICENSE file.
13 * This software is provided "AS IS" with no warranty of any kind,
14 * express or implied, and with no claim as to its suitability for any
22 #include <lemon/list_graph.h>
23 #include <lemon/smart_graph.h>
24 #include <lemon/graph_reader.h>
25 #include <lemon/dimacs.h>
26 #include <lemon/time_measure.h>
28 #include <lemon/cycle_canceling.h>
29 #include <lemon/capacity_scaling.h>
30 #include <lemon/cost_scaling.h>
31 #include <lemon/network_simplex.h>
33 #include <lemon/min_cost_flow.h>
34 #include <lemon/min_cost_max_flow.h>
36 #include "test_tools.h"
38 using namespace lemon;
40 // Checks the feasibility of a flow
41 template < typename Graph, typename LowerMap, typename CapacityMap,
42 typename SupplyMap, typename FlowMap >
43 bool checkFlow( const Graph& gr,
44 const LowerMap& lower, const CapacityMap& upper,
45 const SupplyMap& supply, const FlowMap& flow )
47 GRAPH_TYPEDEFS(typename Graph);
48 for (EdgeIt e(gr); e != INVALID; ++e)
49 if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
51 for (NodeIt n(gr); n != INVALID; ++n) {
52 typename SupplyMap::Value sum = 0;
53 for (OutEdgeIt e(gr, n); e != INVALID; ++e)
55 for (InEdgeIt e(gr, n); e != INVALID; ++e)
57 if (sum != supply[n]) return false;
63 // Checks the optimalitiy of a flow
64 template < typename Graph, typename LowerMap, typename CapacityMap,
65 typename CostMap, typename FlowMap, typename PotentialMap >
66 bool checkOptimality( const Graph& gr, const LowerMap& lower,
67 const CapacityMap& upper, const CostMap& cost,
68 const FlowMap& flow, const PotentialMap& pi )
70 GRAPH_TYPEDEFS(typename Graph);
71 // Checking the Complementary Slackness optimality condition
73 for (EdgeIt e(gr); e != INVALID; ++e) {
74 typename CostMap::Value red_cost =
75 cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
76 opt = red_cost == 0 ||
77 (red_cost > 0 && flow[e] == lower[e]) ||
78 (red_cost < 0 && flow[e] == upper[e]);
84 // Runs a minimum cost flow algorithm and checks the results
85 template < typename MinCostFlowImpl, typename Graph,
86 typename LowerMap, typename CapacityMap,
87 typename CostMap, typename SupplyMap >
88 void checkMcf( std::string test_id,
89 const MinCostFlowImpl& mcf, const Graph& gr,
90 const LowerMap& lower, const CapacityMap& upper,
91 const CostMap& cost, const SupplyMap& supply,
92 bool mcf_result, bool result,
93 typename CostMap::Value total )
95 check(mcf_result == result, "Wrong result " + test_id);
97 check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
98 "The flow is not feasible " + test_id);
99 check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
100 check(checkOptimality(gr, lower, upper, cost, mcf.flowMap(), mcf.potentialMap()),
101 "Wrong potentials " + test_id);
107 // Various tests on a small graph
109 typedef ListGraph Graph;
110 GRAPH_TYPEDEFS(ListGraph);
112 // Reading the test graph
114 Graph::EdgeMap<int> c(gr), l1(gr), l2(gr), u(gr);
115 Graph::NodeMap<int> s1(gr), s2(gr), s3(gr);
120 fname = std::string(getenv("srcdir"));
122 fname += "/test/min_cost_flow_test.lgf";
124 std::ifstream input(fname.c_str());
125 check(input, "Input file '" << fname << "' not found");
126 GraphReader<Graph>(input, gr).
127 readEdgeMap("cost", c).
128 readEdgeMap("capacity", u).
129 readEdgeMap("lower1", l1).
130 readEdgeMap("lower2", l2).
131 readNodeMap("supply1", s1).
132 readNodeMap("supply2", s2).
133 readNodeMap("supply3", s3).
134 readNode("source", v).
135 readNode("target", w).
139 // Testing CapacityScaling (scaling enabled)
141 CapacityScaling<Graph> mcf1(gr,u,c,s1); checkMcf("#A1",mcf1,gr,l1,u,c,s1,mcf1.run(),true, 0);
142 CapacityScaling<Graph> mcf2(gr,u,c,s2); checkMcf("#A2",mcf2,gr,l1,u,c,s2,mcf2.run(),true, 5240);
143 CapacityScaling<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#A3",mcf3,gr,l1,u,c,s3,mcf3.run(),true, 7620);
144 CapacityScaling<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#A4",mcf4,gr,l2,u,c,s1,mcf4.run(),false, 0);
145 CapacityScaling<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#A5",mcf5,gr,l2,u,c,s2,mcf5.run(),true, 5970);
146 CapacityScaling<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#A6",mcf6,gr,l2,u,c,s3,mcf6.run(),true, 8010);
148 // Testing CapacityScaling (scaling disabled)
150 CapacityScaling<Graph> mcf1(gr,u,c,s1); checkMcf("#B1",mcf1,gr,l1,u,c,s1,mcf1.run(false),true, 0);
151 CapacityScaling<Graph> mcf2(gr,u,c,s2); checkMcf("#B2",mcf2,gr,l1,u,c,s2,mcf2.run(false),true, 5240);
152 CapacityScaling<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#B3",mcf3,gr,l1,u,c,s3,mcf3.run(false),true, 7620);
153 CapacityScaling<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#B4",mcf4,gr,l2,u,c,s1,mcf4.run(false),false, 0);
154 CapacityScaling<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#B5",mcf5,gr,l2,u,c,s2,mcf5.run(false),true, 5970);
155 CapacityScaling<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#B6",mcf6,gr,l2,u,c,s3,mcf6.run(false),true, 8010);
158 // Testing CostScaling
160 CostScaling<Graph> mcf1(gr,u,c,s1); checkMcf("#C1",mcf1,gr,l1,u,c,s1,mcf1.run(),true, 0);
161 CostScaling<Graph> mcf2(gr,u,c,s2); checkMcf("#C2",mcf2,gr,l1,u,c,s2,mcf2.run(),true, 5240);
162 CostScaling<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#C3",mcf3,gr,l1,u,c,s3,mcf3.run(),true, 7620);
163 CostScaling<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#C4",mcf4,gr,l2,u,c,s1,mcf4.run(),false, 0);
164 CostScaling<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#C5",mcf5,gr,l2,u,c,s2,mcf5.run(),true, 5970);
165 CostScaling<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#C6",mcf6,gr,l2,u,c,s3,mcf6.run(),true, 8010);
168 // Testing NetworkSimplex (with the default pivot rule)
170 NetworkSimplex<Graph> mcf1(gr,u,c,s1); checkMcf("#D1",mcf1,gr,l1,u,c,s1,mcf1.run(),true, 0);
171 NetworkSimplex<Graph> mcf2(gr,u,c,s2); checkMcf("#D2",mcf2,gr,l1,u,c,s2,mcf2.run(),true, 5240);
172 NetworkSimplex<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#D3",mcf3,gr,l1,u,c,s3,mcf3.run(),true, 7620);
173 NetworkSimplex<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#D4",mcf4,gr,l2,u,c,s1,mcf4.run(),false, 0);
174 NetworkSimplex<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#D5",mcf5,gr,l2,u,c,s2,mcf5.run(),true, 5970);
175 NetworkSimplex<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#D6",mcf6,gr,l2,u,c,s3,mcf6.run(),true, 8010);
177 // Testing NetworkSimplex (with FIRST_ELIGIBLE_PIVOT)
179 NetworkSimplex<Graph>::PivotRuleEnum pr = NetworkSimplex<Graph>::FIRST_ELIGIBLE_PIVOT;
180 NetworkSimplex<Graph> mcf1(gr,u,c,s1); checkMcf("#E1",mcf1,gr,l1,u,c,s1,mcf1.run(pr),true, 0);
181 NetworkSimplex<Graph> mcf2(gr,u,c,s2); checkMcf("#E2",mcf2,gr,l1,u,c,s2,mcf2.run(pr),true, 5240);
182 NetworkSimplex<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#E3",mcf3,gr,l1,u,c,s3,mcf3.run(pr),true, 7620);
183 NetworkSimplex<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#E4",mcf4,gr,l2,u,c,s1,mcf4.run(pr),false, 0);
184 NetworkSimplex<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#E5",mcf5,gr,l2,u,c,s2,mcf5.run(pr),true, 5970);
185 NetworkSimplex<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#E6",mcf6,gr,l2,u,c,s3,mcf6.run(pr),true, 8010);
187 // Testing NetworkSimplex (with BEST_ELIGIBLE_PIVOT)
189 NetworkSimplex<Graph>::PivotRuleEnum pr = NetworkSimplex<Graph>::BEST_ELIGIBLE_PIVOT;
190 NetworkSimplex<Graph> mcf1(gr,u,c,s1); checkMcf("#F1",mcf1,gr,l1,u,c,s1,mcf1.run(pr),true, 0);
191 NetworkSimplex<Graph> mcf2(gr,u,c,s2); checkMcf("#F2",mcf2,gr,l1,u,c,s2,mcf2.run(pr),true, 5240);
192 NetworkSimplex<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#F3",mcf3,gr,l1,u,c,s3,mcf3.run(pr),true, 7620);
193 NetworkSimplex<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#F4",mcf4,gr,l2,u,c,s1,mcf4.run(pr),false, 0);
194 NetworkSimplex<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#F5",mcf5,gr,l2,u,c,s2,mcf5.run(pr),true, 5970);
195 NetworkSimplex<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#F6",mcf6,gr,l2,u,c,s3,mcf6.run(pr),true, 8010);
197 // Testing NetworkSimplex (with BLOCK_SEARCH_PIVOT)
199 NetworkSimplex<Graph>::PivotRuleEnum pr = NetworkSimplex<Graph>::BLOCK_SEARCH_PIVOT;
200 NetworkSimplex<Graph> mcf1(gr,u,c,s1); checkMcf("#G1",mcf1,gr,l1,u,c,s1,mcf1.run(pr),true, 0);
201 NetworkSimplex<Graph> mcf2(gr,u,c,s2); checkMcf("#G2",mcf2,gr,l1,u,c,s2,mcf2.run(pr),true, 5240);
202 NetworkSimplex<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#G3",mcf3,gr,l1,u,c,s3,mcf3.run(pr),true, 7620);
203 NetworkSimplex<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#G4",mcf4,gr,l2,u,c,s1,mcf4.run(pr),false, 0);
204 NetworkSimplex<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#G5",mcf5,gr,l2,u,c,s2,mcf5.run(pr),true, 5970);
205 NetworkSimplex<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#G6",mcf6,gr,l2,u,c,s3,mcf6.run(pr),true, 8010);
207 // Testing NetworkSimplex (with LIMITED_SEARCH_PIVOT)
209 NetworkSimplex<Graph>::PivotRuleEnum pr = NetworkSimplex<Graph>::LIMITED_SEARCH_PIVOT;
210 NetworkSimplex<Graph> mcf1(gr,u,c,s1); checkMcf("#H1",mcf1,gr,l1,u,c,s1,mcf1.run(pr),true, 0);
211 NetworkSimplex<Graph> mcf2(gr,u,c,s2); checkMcf("#H2",mcf2,gr,l1,u,c,s2,mcf2.run(pr),true, 5240);
212 NetworkSimplex<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#H3",mcf3,gr,l1,u,c,s3,mcf3.run(pr),true, 7620);
213 NetworkSimplex<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#H4",mcf4,gr,l2,u,c,s1,mcf4.run(pr),false, 0);
214 NetworkSimplex<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#H5",mcf5,gr,l2,u,c,s2,mcf5.run(pr),true, 5970);
215 NetworkSimplex<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#H6",mcf6,gr,l2,u,c,s3,mcf6.run(pr),true, 8010);
217 // Testing NetworkSimplex (with CANDIDATE_LIST_PIVOT)
219 NetworkSimplex<Graph>::PivotRuleEnum pr = NetworkSimplex<Graph>::CANDIDATE_LIST_PIVOT;
220 NetworkSimplex<Graph> mcf1(gr,u,c,s1); checkMcf("#I1",mcf1,gr,l1,u,c,s1,mcf1.run(pr),true, 0);
221 NetworkSimplex<Graph> mcf2(gr,u,c,s2); checkMcf("#I2",mcf2,gr,l1,u,c,s2,mcf2.run(pr),true, 5240);
222 NetworkSimplex<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#I3",mcf3,gr,l1,u,c,s3,mcf3.run(pr),true, 7620);
223 NetworkSimplex<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#I4",mcf4,gr,l2,u,c,s1,mcf4.run(pr),false, 0);
224 NetworkSimplex<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#I5",mcf5,gr,l2,u,c,s2,mcf5.run(pr),true, 5970);
225 NetworkSimplex<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#I6",mcf6,gr,l2,u,c,s3,mcf6.run(pr),true, 8010);
228 // Testing CycleCanceling (with BellmanFord)
230 CycleCanceling<Graph> mcf1(gr,u,c,s1); checkMcf("#J1",mcf1,gr,l1,u,c,s1,mcf1.run(),true, 0);
231 CycleCanceling<Graph> mcf2(gr,u,c,s2); checkMcf("#J2",mcf2,gr,l1,u,c,s2,mcf2.run(),true, 5240);
232 CycleCanceling<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#J3",mcf3,gr,l1,u,c,s3,mcf3.run(),true, 7620);
233 CycleCanceling<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#J4",mcf4,gr,l2,u,c,s1,mcf4.run(),false, 0);
234 CycleCanceling<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#J5",mcf5,gr,l2,u,c,s2,mcf5.run(),true, 5970);
235 CycleCanceling<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#J6",mcf6,gr,l2,u,c,s3,mcf6.run(),true, 8010);
237 // Testing CycleCanceling (with MinMeanCycle)
239 CycleCanceling<Graph> mcf1(gr,u,c,s1); checkMcf("#K1",mcf1,gr,l1,u,c,s1,mcf1.run(true),true, 0);
240 CycleCanceling<Graph> mcf2(gr,u,c,s2); checkMcf("#K2",mcf2,gr,l1,u,c,s2,mcf2.run(true),true, 5240);
241 CycleCanceling<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#K3",mcf3,gr,l1,u,c,s3,mcf3.run(true),true, 7620);
242 CycleCanceling<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#K4",mcf4,gr,l2,u,c,s1,mcf4.run(true),false, 0);
243 CycleCanceling<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#K5",mcf5,gr,l2,u,c,s2,mcf5.run(true),true, 5970);
244 CycleCanceling<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#K6",mcf6,gr,l2,u,c,s3,mcf6.run(true),true, 8010);
247 // Testing MinCostFlow
249 MinCostFlow<Graph> mcf1(gr,u,c,s1); checkMcf("#L1",mcf1,gr,l1,u,c,s1,mcf1.run(),true, 0);
250 MinCostFlow<Graph> mcf2(gr,u,c,s2); checkMcf("#L2",mcf2,gr,l1,u,c,s2,mcf2.run(),true, 5240);
251 MinCostFlow<Graph> mcf3(gr,u,c,v,w,27); checkMcf("#L3",mcf3,gr,l1,u,c,s3,mcf3.run(),true, 7620);
252 MinCostFlow<Graph> mcf4(gr,l2,u,c,s1); checkMcf("#L4",mcf4,gr,l2,u,c,s1,mcf4.run(),false, 0);
253 MinCostFlow<Graph> mcf5(gr,l2,u,c,s2); checkMcf("#L5",mcf5,gr,l2,u,c,s2,mcf5.run(),true, 5970);
254 MinCostFlow<Graph> mcf6(gr,l2,u,c,v,w,27); checkMcf("#L6",mcf6,gr,l2,u,c,s3,mcf6.run(),true, 8010);
257 // Testing MinCostMaxFlow
259 MinCostMaxFlow<Graph> mcmf(gr,u,c,v,w);
261 checkMcf("#M1",mcmf,gr,l1,u,c,s3,true,true,7620);
265 // Benchmark test on a DIMACS network
267 typedef SmartGraph Graph;
268 GRAPH_TYPEDEFS(SmartGraph);
270 // Reading the test graph
272 Graph::EdgeMap<int> lower(graph), capacity(graph), cost(graph);
273 Graph::NodeMap<int> supply(graph);
277 fname = std::string(getenv("srcdir"));
279 fname += "/test/min_cost_flow_test.net";
281 std::ifstream input(fname.c_str());
282 check(input, "Input file '" << fname << "' not found");
283 readDimacs(input, graph, lower, capacity, cost, supply);
289 NetworkSimplex<Graph> mcf(graph, lower, capacity, cost, supply);
290 bool res = mcf.run();
292 checkMcf("#T3", mcf, graph, lower, capacity, cost, supply, res, true, 196587626);
293 std::cout << "NetworkSimplex";
294 std::cout << std::endl << t << std::endl;
299 CapacityScaling<Graph> mcf(graph, lower, capacity, cost, supply);
300 bool res = mcf.run();
302 checkMcf("#T1", mcf, graph, lower, capacity, cost, supply, res, true, 196587626);
303 std::cout << "CapacityScaling";
304 std::cout << std::endl << t << std::endl;
309 CostScaling<Graph> mcf(graph, lower, capacity, cost, supply);
310 bool res = mcf.run();
312 checkMcf("#T2", mcf, graph, lower, capacity, cost, supply, res, true, 196587626);
313 std::cout << "CostScaling";
314 std::cout << std::endl << t << std::endl;
319 CycleCanceling<Graph> mcf(graph, lower, capacity, cost, supply);
320 bool res = mcf.run();
322 checkMcf("#T4", mcf, graph, lower, capacity, cost, supply, res, true, 196587626);
323 std::cout << "CycleCanceling";
324 std::cout << std::endl << t << std::endl;