1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/test/min_cost_flow_test.cc Thu Dec 10 17:05:35 2009 +0100
1.3 @@ -0,0 +1,450 @@
1.4 +/* -*- mode: C++; indent-tabs-mode: nil; -*-
1.5 + *
1.6 + * This file is a part of LEMON, a generic C++ optimization library.
1.7 + *
1.8 + * Copyright (C) 2003-2009
1.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
1.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
1.11 + *
1.12 + * Permission to use, modify and distribute this software is granted
1.13 + * provided that this copyright notice appears in all copies. For
1.14 + * precise terms see the accompanying LICENSE file.
1.15 + *
1.16 + * This software is provided "AS IS" with no warranty of any kind,
1.17 + * express or implied, and with no claim as to its suitability for any
1.18 + * purpose.
1.19 + *
1.20 + */
1.21 +
1.22 +#include <iostream>
1.23 +#include <fstream>
1.24 +#include <limits>
1.25 +
1.26 +#include <lemon/list_graph.h>
1.27 +#include <lemon/lgf_reader.h>
1.28 +
1.29 +#include <lemon/network_simplex.h>
1.30 +
1.31 +#include <lemon/concepts/digraph.h>
1.32 +#include <lemon/concept_check.h>
1.33 +
1.34 +#include "test_tools.h"
1.35 +
1.36 +using namespace lemon;
1.37 +
1.38 +char test_lgf[] =
1.39 + "@nodes\n"
1.40 + "label sup1 sup2 sup3 sup4 sup5 sup6\n"
1.41 + " 1 20 27 0 30 20 30\n"
1.42 + " 2 -4 0 0 0 -8 -3\n"
1.43 + " 3 0 0 0 0 0 0\n"
1.44 + " 4 0 0 0 0 0 0\n"
1.45 + " 5 9 0 0 0 6 11\n"
1.46 + " 6 -6 0 0 0 -5 -6\n"
1.47 + " 7 0 0 0 0 0 0\n"
1.48 + " 8 0 0 0 0 0 3\n"
1.49 + " 9 3 0 0 0 0 0\n"
1.50 + " 10 -2 0 0 0 -7 -2\n"
1.51 + " 11 0 0 0 0 -10 0\n"
1.52 + " 12 -20 -27 0 -30 -30 -20\n"
1.53 + "\n"
1.54 + "@arcs\n"
1.55 + " cost cap low1 low2 low3\n"
1.56 + " 1 2 70 11 0 8 8\n"
1.57 + " 1 3 150 3 0 1 0\n"
1.58 + " 1 4 80 15 0 2 2\n"
1.59 + " 2 8 80 12 0 0 0\n"
1.60 + " 3 5 140 5 0 3 1\n"
1.61 + " 4 6 60 10 0 1 0\n"
1.62 + " 4 7 80 2 0 0 0\n"
1.63 + " 4 8 110 3 0 0 0\n"
1.64 + " 5 7 60 14 0 0 0\n"
1.65 + " 5 11 120 12 0 0 0\n"
1.66 + " 6 3 0 3 0 0 0\n"
1.67 + " 6 9 140 4 0 0 0\n"
1.68 + " 6 10 90 8 0 0 0\n"
1.69 + " 7 1 30 5 0 0 -5\n"
1.70 + " 8 12 60 16 0 4 3\n"
1.71 + " 9 12 50 6 0 0 0\n"
1.72 + "10 12 70 13 0 5 2\n"
1.73 + "10 2 100 7 0 0 0\n"
1.74 + "10 7 60 10 0 0 -3\n"
1.75 + "11 10 20 14 0 6 -20\n"
1.76 + "12 11 30 10 0 0 -10\n"
1.77 + "\n"
1.78 + "@attributes\n"
1.79 + "source 1\n"
1.80 + "target 12\n";
1.81 +
1.82 +
1.83 +enum SupplyType {
1.84 + EQ,
1.85 + GEQ,
1.86 + LEQ
1.87 +};
1.88 +
1.89 +// Check the interface of an MCF algorithm
1.90 +template <typename GR, typename Value, typename Cost>
1.91 +class McfClassConcept
1.92 +{
1.93 +public:
1.94 +
1.95 + template <typename MCF>
1.96 + struct Constraints {
1.97 + void constraints() {
1.98 + checkConcept<concepts::Digraph, GR>();
1.99 +
1.100 + const Constraints& me = *this;
1.101 +
1.102 + MCF mcf(me.g);
1.103 + const MCF& const_mcf = mcf;
1.104 +
1.105 + b = mcf.reset()
1.106 + .lowerMap(me.lower)
1.107 + .upperMap(me.upper)
1.108 + .costMap(me.cost)
1.109 + .supplyMap(me.sup)
1.110 + .stSupply(me.n, me.n, me.k)
1.111 + .run();
1.112 +
1.113 + c = const_mcf.totalCost();
1.114 + x = const_mcf.template totalCost<double>();
1.115 + v = const_mcf.flow(me.a);
1.116 + c = const_mcf.potential(me.n);
1.117 + const_mcf.flowMap(fm);
1.118 + const_mcf.potentialMap(pm);
1.119 + }
1.120 +
1.121 + typedef typename GR::Node Node;
1.122 + typedef typename GR::Arc Arc;
1.123 + typedef concepts::ReadMap<Node, Value> NM;
1.124 + typedef concepts::ReadMap<Arc, Value> VAM;
1.125 + typedef concepts::ReadMap<Arc, Cost> CAM;
1.126 + typedef concepts::WriteMap<Arc, Value> FlowMap;
1.127 + typedef concepts::WriteMap<Node, Cost> PotMap;
1.128 +
1.129 + GR g;
1.130 + VAM lower;
1.131 + VAM upper;
1.132 + CAM cost;
1.133 + NM sup;
1.134 + Node n;
1.135 + Arc a;
1.136 + Value k;
1.137 +
1.138 + FlowMap fm;
1.139 + PotMap pm;
1.140 + bool b;
1.141 + double x;
1.142 + typename MCF::Value v;
1.143 + typename MCF::Cost c;
1.144 + };
1.145 +
1.146 +};
1.147 +
1.148 +
1.149 +// Check the feasibility of the given flow (primal soluiton)
1.150 +template < typename GR, typename LM, typename UM,
1.151 + typename SM, typename FM >
1.152 +bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
1.153 + const SM& supply, const FM& flow,
1.154 + SupplyType type = EQ )
1.155 +{
1.156 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.157 +
1.158 + for (ArcIt e(gr); e != INVALID; ++e) {
1.159 + if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
1.160 + }
1.161 +
1.162 + for (NodeIt n(gr); n != INVALID; ++n) {
1.163 + typename SM::Value sum = 0;
1.164 + for (OutArcIt e(gr, n); e != INVALID; ++e)
1.165 + sum += flow[e];
1.166 + for (InArcIt e(gr, n); e != INVALID; ++e)
1.167 + sum -= flow[e];
1.168 + bool b = (type == EQ && sum == supply[n]) ||
1.169 + (type == GEQ && sum >= supply[n]) ||
1.170 + (type == LEQ && sum <= supply[n]);
1.171 + if (!b) return false;
1.172 + }
1.173 +
1.174 + return true;
1.175 +}
1.176 +
1.177 +// Check the feasibility of the given potentials (dual soluiton)
1.178 +// using the "Complementary Slackness" optimality condition
1.179 +template < typename GR, typename LM, typename UM,
1.180 + typename CM, typename SM, typename FM, typename PM >
1.181 +bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
1.182 + const CM& cost, const SM& supply, const FM& flow,
1.183 + const PM& pi, SupplyType type )
1.184 +{
1.185 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.186 +
1.187 + bool opt = true;
1.188 + for (ArcIt e(gr); opt && e != INVALID; ++e) {
1.189 + typename CM::Value red_cost =
1.190 + cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
1.191 + opt = red_cost == 0 ||
1.192 + (red_cost > 0 && flow[e] == lower[e]) ||
1.193 + (red_cost < 0 && flow[e] == upper[e]);
1.194 + }
1.195 +
1.196 + for (NodeIt n(gr); opt && n != INVALID; ++n) {
1.197 + typename SM::Value sum = 0;
1.198 + for (OutArcIt e(gr, n); e != INVALID; ++e)
1.199 + sum += flow[e];
1.200 + for (InArcIt e(gr, n); e != INVALID; ++e)
1.201 + sum -= flow[e];
1.202 + if (type != LEQ) {
1.203 + opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
1.204 + } else {
1.205 + opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
1.206 + }
1.207 + }
1.208 +
1.209 + return opt;
1.210 +}
1.211 +
1.212 +// Check whether the dual cost is equal to the primal cost
1.213 +template < typename GR, typename LM, typename UM,
1.214 + typename CM, typename SM, typename PM >
1.215 +bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
1.216 + const CM& cost, const SM& supply, const PM& pi,
1.217 + typename CM::Value total )
1.218 +{
1.219 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.220 +
1.221 + typename CM::Value dual_cost = 0;
1.222 + SM red_supply(gr);
1.223 + for (NodeIt n(gr); n != INVALID; ++n) {
1.224 + red_supply[n] = supply[n];
1.225 + }
1.226 + for (ArcIt a(gr); a != INVALID; ++a) {
1.227 + if (lower[a] != 0) {
1.228 + dual_cost += lower[a] * cost[a];
1.229 + red_supply[gr.source(a)] -= lower[a];
1.230 + red_supply[gr.target(a)] += lower[a];
1.231 + }
1.232 + }
1.233 +
1.234 + for (NodeIt n(gr); n != INVALID; ++n) {
1.235 + dual_cost -= red_supply[n] * pi[n];
1.236 + }
1.237 + for (ArcIt a(gr); a != INVALID; ++a) {
1.238 + typename CM::Value red_cost =
1.239 + cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
1.240 + dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
1.241 + }
1.242 +
1.243 + return dual_cost == total;
1.244 +}
1.245 +
1.246 +// Run a minimum cost flow algorithm and check the results
1.247 +template < typename MCF, typename GR,
1.248 + typename LM, typename UM,
1.249 + typename CM, typename SM,
1.250 + typename PT >
1.251 +void checkMcf( const MCF& mcf, PT mcf_result,
1.252 + const GR& gr, const LM& lower, const UM& upper,
1.253 + const CM& cost, const SM& supply,
1.254 + PT result, bool optimal, typename CM::Value total,
1.255 + const std::string &test_id = "",
1.256 + SupplyType type = EQ )
1.257 +{
1.258 + check(mcf_result == result, "Wrong result " + test_id);
1.259 + if (optimal) {
1.260 + typename GR::template ArcMap<typename SM::Value> flow(gr);
1.261 + typename GR::template NodeMap<typename CM::Value> pi(gr);
1.262 + mcf.flowMap(flow);
1.263 + mcf.potentialMap(pi);
1.264 + check(checkFlow(gr, lower, upper, supply, flow, type),
1.265 + "The flow is not feasible " + test_id);
1.266 + check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
1.267 + check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
1.268 + "Wrong potentials " + test_id);
1.269 + check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
1.270 + "Wrong dual cost " + test_id);
1.271 + }
1.272 +}
1.273 +
1.274 +int main()
1.275 +{
1.276 + // Check the interfaces
1.277 + {
1.278 + typedef concepts::Digraph GR;
1.279 + checkConcept< McfClassConcept<GR, int, int>,
1.280 + NetworkSimplex<GR> >();
1.281 + checkConcept< McfClassConcept<GR, double, double>,
1.282 + NetworkSimplex<GR, double> >();
1.283 + checkConcept< McfClassConcept<GR, int, double>,
1.284 + NetworkSimplex<GR, int, double> >();
1.285 + }
1.286 +
1.287 + // Run various MCF tests
1.288 + typedef ListDigraph Digraph;
1.289 + DIGRAPH_TYPEDEFS(ListDigraph);
1.290 +
1.291 + // Read the test digraph
1.292 + Digraph gr;
1.293 + Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);
1.294 + Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr), s6(gr);
1.295 + ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
1.296 + Node v, w;
1.297 +
1.298 + std::istringstream input(test_lgf);
1.299 + DigraphReader<Digraph>(gr, input)
1.300 + .arcMap("cost", c)
1.301 + .arcMap("cap", u)
1.302 + .arcMap("low1", l1)
1.303 + .arcMap("low2", l2)
1.304 + .arcMap("low3", l3)
1.305 + .nodeMap("sup1", s1)
1.306 + .nodeMap("sup2", s2)
1.307 + .nodeMap("sup3", s3)
1.308 + .nodeMap("sup4", s4)
1.309 + .nodeMap("sup5", s5)
1.310 + .nodeMap("sup6", s6)
1.311 + .node("source", v)
1.312 + .node("target", w)
1.313 + .run();
1.314 +
1.315 + // Build test digraphs with negative costs
1.316 + Digraph neg_gr;
1.317 + Node n1 = neg_gr.addNode();
1.318 + Node n2 = neg_gr.addNode();
1.319 + Node n3 = neg_gr.addNode();
1.320 + Node n4 = neg_gr.addNode();
1.321 + Node n5 = neg_gr.addNode();
1.322 + Node n6 = neg_gr.addNode();
1.323 + Node n7 = neg_gr.addNode();
1.324 +
1.325 + Arc a1 = neg_gr.addArc(n1, n2);
1.326 + Arc a2 = neg_gr.addArc(n1, n3);
1.327 + Arc a3 = neg_gr.addArc(n2, n4);
1.328 + Arc a4 = neg_gr.addArc(n3, n4);
1.329 + Arc a5 = neg_gr.addArc(n3, n2);
1.330 + Arc a6 = neg_gr.addArc(n5, n3);
1.331 + Arc a7 = neg_gr.addArc(n5, n6);
1.332 + Arc a8 = neg_gr.addArc(n6, n7);
1.333 + Arc a9 = neg_gr.addArc(n7, n5);
1.334 +
1.335 + Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
1.336 + ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
1.337 + Digraph::NodeMap<int> neg_s(neg_gr, 0);
1.338 +
1.339 + neg_l2[a7] = 1000;
1.340 + neg_l2[a8] = -1000;
1.341 +
1.342 + neg_s[n1] = 100;
1.343 + neg_s[n4] = -100;
1.344 +
1.345 + neg_c[a1] = 100;
1.346 + neg_c[a2] = 30;
1.347 + neg_c[a3] = 20;
1.348 + neg_c[a4] = 80;
1.349 + neg_c[a5] = 50;
1.350 + neg_c[a6] = 10;
1.351 + neg_c[a7] = 80;
1.352 + neg_c[a8] = 30;
1.353 + neg_c[a9] = -120;
1.354 +
1.355 + Digraph negs_gr;
1.356 + Digraph::NodeMap<int> negs_s(negs_gr);
1.357 + Digraph::ArcMap<int> negs_c(negs_gr);
1.358 + ConstMap<Arc, int> negs_l(0), negs_u(1000);
1.359 + n1 = negs_gr.addNode();
1.360 + n2 = negs_gr.addNode();
1.361 + negs_s[n1] = 100;
1.362 + negs_s[n2] = -300;
1.363 + negs_c[negs_gr.addArc(n1, n2)] = -1;
1.364 +
1.365 +
1.366 + // A. Test NetworkSimplex with the default pivot rule
1.367 + {
1.368 + NetworkSimplex<Digraph> mcf(gr);
1.369 +
1.370 + // Check the equality form
1.371 + mcf.upperMap(u).costMap(c);
1.372 + checkMcf(mcf, mcf.supplyMap(s1).run(),
1.373 + gr, l1, u, c, s1, mcf.OPTIMAL, true, 5240, "#A1");
1.374 + checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
1.375 + gr, l1, u, c, s2, mcf.OPTIMAL, true, 7620, "#A2");
1.376 + mcf.lowerMap(l2);
1.377 + checkMcf(mcf, mcf.supplyMap(s1).run(),
1.378 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#A3");
1.379 + checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
1.380 + gr, l2, u, c, s2, mcf.OPTIMAL, true, 8010, "#A4");
1.381 + mcf.reset();
1.382 + checkMcf(mcf, mcf.supplyMap(s1).run(),
1.383 + gr, l1, cu, cc, s1, mcf.OPTIMAL, true, 74, "#A5");
1.384 + checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
1.385 + gr, l2, cu, cc, s2, mcf.OPTIMAL, true, 94, "#A6");
1.386 + mcf.reset();
1.387 + checkMcf(mcf, mcf.run(),
1.388 + gr, l1, cu, cc, s3, mcf.OPTIMAL, true, 0, "#A7");
1.389 + checkMcf(mcf, mcf.lowerMap(l2).upperMap(u).run(),
1.390 + gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");
1.391 + mcf.reset().lowerMap(l3).upperMap(u).costMap(c).supplyMap(s4);
1.392 + checkMcf(mcf, mcf.run(),
1.393 + gr, l3, u, c, s4, mcf.OPTIMAL, true, 6360, "#A9");
1.394 +
1.395 + // Check the GEQ form
1.396 + mcf.reset().upperMap(u).costMap(c).supplyMap(s5);
1.397 + checkMcf(mcf, mcf.run(),
1.398 + gr, l1, u, c, s5, mcf.OPTIMAL, true, 3530, "#A10", GEQ);
1.399 + mcf.supplyType(mcf.GEQ);
1.400 + checkMcf(mcf, mcf.lowerMap(l2).run(),
1.401 + gr, l2, u, c, s5, mcf.OPTIMAL, true, 4540, "#A11", GEQ);
1.402 + mcf.supplyMap(s6);
1.403 + checkMcf(mcf, mcf.run(),
1.404 + gr, l2, u, c, s6, mcf.INFEASIBLE, false, 0, "#A12", GEQ);
1.405 +
1.406 + // Check the LEQ form
1.407 + mcf.reset().supplyType(mcf.LEQ);
1.408 + mcf.upperMap(u).costMap(c).supplyMap(s6);
1.409 + checkMcf(mcf, mcf.run(),
1.410 + gr, l1, u, c, s6, mcf.OPTIMAL, true, 5080, "#A13", LEQ);
1.411 + checkMcf(mcf, mcf.lowerMap(l2).run(),
1.412 + gr, l2, u, c, s6, mcf.OPTIMAL, true, 5930, "#A14", LEQ);
1.413 + mcf.supplyMap(s5);
1.414 + checkMcf(mcf, mcf.run(),
1.415 + gr, l2, u, c, s5, mcf.INFEASIBLE, false, 0, "#A15", LEQ);
1.416 +
1.417 + // Check negative costs
1.418 + NetworkSimplex<Digraph> neg_mcf(neg_gr);
1.419 + neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
1.420 + checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
1.421 + neg_c, neg_s, neg_mcf.UNBOUNDED, false, 0, "#A16");
1.422 + neg_mcf.upperMap(neg_u2);
1.423 + checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
1.424 + neg_c, neg_s, neg_mcf.OPTIMAL, true, -40000, "#A17");
1.425 + neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
1.426 + checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
1.427 + neg_c, neg_s, neg_mcf.UNBOUNDED, false, 0, "#A18");
1.428 +
1.429 + NetworkSimplex<Digraph> negs_mcf(negs_gr);
1.430 + negs_mcf.costMap(negs_c).supplyMap(negs_s);
1.431 + checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
1.432 + negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
1.433 + }
1.434 +
1.435 + // B. Test NetworkSimplex with each pivot rule
1.436 + {
1.437 + NetworkSimplex<Digraph> mcf(gr);
1.438 + mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);
1.439 +
1.440 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
1.441 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B1");
1.442 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
1.443 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B2");
1.444 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
1.445 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B3");
1.446 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
1.447 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B4");
1.448 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
1.449 + gr, l2, u, c, s1, mcf.OPTIMAL, true, 5970, "#B5");
1.450 + }
1.451 +
1.452 + return 0;
1.453 +}