1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/test/min_cost_flow_test.cc Tue Apr 21 15:18:54 2009 +0100
1.3 @@ -0,0 +1,339 @@
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 +
1.25 +#include <lemon/list_graph.h>
1.26 +#include <lemon/lgf_reader.h>
1.27 +
1.28 +#include <lemon/network_simplex.h>
1.29 +
1.30 +#include <lemon/concepts/digraph.h>
1.31 +#include <lemon/concept_check.h>
1.32 +
1.33 +#include "test_tools.h"
1.34 +
1.35 +using namespace lemon;
1.36 +
1.37 +char test_lgf[] =
1.38 + "@nodes\n"
1.39 + "label sup1 sup2 sup3 sup4 sup5\n"
1.40 + " 1 20 27 0 20 30\n"
1.41 + " 2 -4 0 0 -8 -3\n"
1.42 + " 3 0 0 0 0 0\n"
1.43 + " 4 0 0 0 0 0\n"
1.44 + " 5 9 0 0 6 11\n"
1.45 + " 6 -6 0 0 -5 -6\n"
1.46 + " 7 0 0 0 0 0\n"
1.47 + " 8 0 0 0 0 3\n"
1.48 + " 9 3 0 0 0 0\n"
1.49 + " 10 -2 0 0 -7 -2\n"
1.50 + " 11 0 0 0 -10 0\n"
1.51 + " 12 -20 -27 0 -30 -20\n"
1.52 + "\n"
1.53 + "@arcs\n"
1.54 + " cost cap low1 low2\n"
1.55 + " 1 2 70 11 0 8\n"
1.56 + " 1 3 150 3 0 1\n"
1.57 + " 1 4 80 15 0 2\n"
1.58 + " 2 8 80 12 0 0\n"
1.59 + " 3 5 140 5 0 3\n"
1.60 + " 4 6 60 10 0 1\n"
1.61 + " 4 7 80 2 0 0\n"
1.62 + " 4 8 110 3 0 0\n"
1.63 + " 5 7 60 14 0 0\n"
1.64 + " 5 11 120 12 0 0\n"
1.65 + " 6 3 0 3 0 0\n"
1.66 + " 6 9 140 4 0 0\n"
1.67 + " 6 10 90 8 0 0\n"
1.68 + " 7 1 30 5 0 0\n"
1.69 + " 8 12 60 16 0 4\n"
1.70 + " 9 12 50 6 0 0\n"
1.71 + "10 12 70 13 0 5\n"
1.72 + "10 2 100 7 0 0\n"
1.73 + "10 7 60 10 0 0\n"
1.74 + "11 10 20 14 0 6\n"
1.75 + "12 11 30 10 0 0\n"
1.76 + "\n"
1.77 + "@attributes\n"
1.78 + "source 1\n"
1.79 + "target 12\n";
1.80 +
1.81 +
1.82 +enum ProblemType {
1.83 + EQ,
1.84 + GEQ,
1.85 + LEQ
1.86 +};
1.87 +
1.88 +// Check the interface of an MCF algorithm
1.89 +template <typename GR, typename Flow, typename Cost>
1.90 +class McfClassConcept
1.91 +{
1.92 +public:
1.93 +
1.94 + template <typename MCF>
1.95 + struct Constraints {
1.96 + void constraints() {
1.97 + checkConcept<concepts::Digraph, GR>();
1.98 +
1.99 + MCF mcf(g);
1.100 +
1.101 + b = mcf.reset()
1.102 + .lowerMap(lower)
1.103 + .upperMap(upper)
1.104 + .capacityMap(upper)
1.105 + .boundMaps(lower, upper)
1.106 + .costMap(cost)
1.107 + .supplyMap(sup)
1.108 + .stSupply(n, n, k)
1.109 + .flowMap(flow)
1.110 + .potentialMap(pot)
1.111 + .run();
1.112 +
1.113 + const MCF& const_mcf = mcf;
1.114 +
1.115 + const typename MCF::FlowMap &fm = const_mcf.flowMap();
1.116 + const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
1.117 +
1.118 + v = const_mcf.totalCost();
1.119 + double x = const_mcf.template totalCost<double>();
1.120 + v = const_mcf.flow(a);
1.121 + v = const_mcf.potential(n);
1.122 +
1.123 + ignore_unused_variable_warning(fm);
1.124 + ignore_unused_variable_warning(pm);
1.125 + ignore_unused_variable_warning(x);
1.126 + }
1.127 +
1.128 + typedef typename GR::Node Node;
1.129 + typedef typename GR::Arc Arc;
1.130 + typedef concepts::ReadMap<Node, Flow> NM;
1.131 + typedef concepts::ReadMap<Arc, Flow> FAM;
1.132 + typedef concepts::ReadMap<Arc, Cost> CAM;
1.133 +
1.134 + const GR &g;
1.135 + const FAM &lower;
1.136 + const FAM &upper;
1.137 + const CAM &cost;
1.138 + const NM ⊃
1.139 + const Node &n;
1.140 + const Arc &a;
1.141 + const Flow &k;
1.142 + Flow v;
1.143 + bool b;
1.144 +
1.145 + typename MCF::FlowMap &flow;
1.146 + typename MCF::PotentialMap &pot;
1.147 + };
1.148 +
1.149 +};
1.150 +
1.151 +
1.152 +// Check the feasibility of the given flow (primal soluiton)
1.153 +template < typename GR, typename LM, typename UM,
1.154 + typename SM, typename FM >
1.155 +bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
1.156 + const SM& supply, const FM& flow,
1.157 + ProblemType type = EQ )
1.158 +{
1.159 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.160 +
1.161 + for (ArcIt e(gr); e != INVALID; ++e) {
1.162 + if (flow[e] < lower[e] || flow[e] > upper[e]) return false;
1.163 + }
1.164 +
1.165 + for (NodeIt n(gr); n != INVALID; ++n) {
1.166 + typename SM::Value sum = 0;
1.167 + for (OutArcIt e(gr, n); e != INVALID; ++e)
1.168 + sum += flow[e];
1.169 + for (InArcIt e(gr, n); e != INVALID; ++e)
1.170 + sum -= flow[e];
1.171 + bool b = (type == EQ && sum == supply[n]) ||
1.172 + (type == GEQ && sum >= supply[n]) ||
1.173 + (type == LEQ && sum <= supply[n]);
1.174 + if (!b) return false;
1.175 + }
1.176 +
1.177 + return true;
1.178 +}
1.179 +
1.180 +// Check the feasibility of the given potentials (dual soluiton)
1.181 +// using the "Complementary Slackness" optimality condition
1.182 +template < typename GR, typename LM, typename UM,
1.183 + typename CM, typename SM, typename FM, typename PM >
1.184 +bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
1.185 + const CM& cost, const SM& supply, const FM& flow,
1.186 + const PM& pi )
1.187 +{
1.188 + TEMPLATE_DIGRAPH_TYPEDEFS(GR);
1.189 +
1.190 + bool opt = true;
1.191 + for (ArcIt e(gr); opt && e != INVALID; ++e) {
1.192 + typename CM::Value red_cost =
1.193 + cost[e] + pi[gr.source(e)] - pi[gr.target(e)];
1.194 + opt = red_cost == 0 ||
1.195 + (red_cost > 0 && flow[e] == lower[e]) ||
1.196 + (red_cost < 0 && flow[e] == upper[e]);
1.197 + }
1.198 +
1.199 + for (NodeIt n(gr); opt && n != INVALID; ++n) {
1.200 + typename SM::Value sum = 0;
1.201 + for (OutArcIt e(gr, n); e != INVALID; ++e)
1.202 + sum += flow[e];
1.203 + for (InArcIt e(gr, n); e != INVALID; ++e)
1.204 + sum -= flow[e];
1.205 + opt = (sum == supply[n]) || (pi[n] == 0);
1.206 + }
1.207 +
1.208 + return opt;
1.209 +}
1.210 +
1.211 +// Run a minimum cost flow algorithm and check the results
1.212 +template < typename MCF, typename GR,
1.213 + typename LM, typename UM,
1.214 + typename CM, typename SM >
1.215 +void checkMcf( const MCF& mcf, bool mcf_result,
1.216 + const GR& gr, const LM& lower, const UM& upper,
1.217 + const CM& cost, const SM& supply,
1.218 + bool result, typename CM::Value total,
1.219 + const std::string &test_id = "",
1.220 + ProblemType type = EQ )
1.221 +{
1.222 + check(mcf_result == result, "Wrong result " + test_id);
1.223 + if (result) {
1.224 + check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
1.225 + "The flow is not feasible " + test_id);
1.226 + check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
1.227 + check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
1.228 + mcf.potentialMap()),
1.229 + "Wrong potentials " + test_id);
1.230 + }
1.231 +}
1.232 +
1.233 +int main()
1.234 +{
1.235 + // Check the interfaces
1.236 + {
1.237 + typedef int Flow;
1.238 + typedef int Cost;
1.239 + // TODO: This typedef should be enabled if the standard maps are
1.240 + // reference maps in the graph concepts (See #190).
1.241 +/**/
1.242 + //typedef concepts::Digraph GR;
1.243 + typedef ListDigraph GR;
1.244 +/**/
1.245 + checkConcept< McfClassConcept<GR, Flow, Cost>,
1.246 + NetworkSimplex<GR, Flow, Cost> >();
1.247 + }
1.248 +
1.249 + // Run various MCF tests
1.250 + typedef ListDigraph Digraph;
1.251 + DIGRAPH_TYPEDEFS(ListDigraph);
1.252 +
1.253 + // Read the test digraph
1.254 + Digraph gr;
1.255 + Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
1.256 + Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
1.257 + ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
1.258 + Node v, w;
1.259 +
1.260 + std::istringstream input(test_lgf);
1.261 + DigraphReader<Digraph>(gr, input)
1.262 + .arcMap("cost", c)
1.263 + .arcMap("cap", u)
1.264 + .arcMap("low1", l1)
1.265 + .arcMap("low2", l2)
1.266 + .nodeMap("sup1", s1)
1.267 + .nodeMap("sup2", s2)
1.268 + .nodeMap("sup3", s3)
1.269 + .nodeMap("sup4", s4)
1.270 + .nodeMap("sup5", s5)
1.271 + .node("source", v)
1.272 + .node("target", w)
1.273 + .run();
1.274 +
1.275 + // A. Test NetworkSimplex with the default pivot rule
1.276 + {
1.277 + NetworkSimplex<Digraph> mcf(gr);
1.278 +
1.279 + // Check the equality form
1.280 + mcf.upperMap(u).costMap(c);
1.281 + checkMcf(mcf, mcf.supplyMap(s1).run(),
1.282 + gr, l1, u, c, s1, true, 5240, "#A1");
1.283 + checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
1.284 + gr, l1, u, c, s2, true, 7620, "#A2");
1.285 + mcf.lowerMap(l2);
1.286 + checkMcf(mcf, mcf.supplyMap(s1).run(),
1.287 + gr, l2, u, c, s1, true, 5970, "#A3");
1.288 + checkMcf(mcf, mcf.stSupply(v, w, 27).run(),
1.289 + gr, l2, u, c, s2, true, 8010, "#A4");
1.290 + mcf.reset();
1.291 + checkMcf(mcf, mcf.supplyMap(s1).run(),
1.292 + gr, l1, cu, cc, s1, true, 74, "#A5");
1.293 + checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),
1.294 + gr, l2, cu, cc, s2, true, 94, "#A6");
1.295 + mcf.reset();
1.296 + checkMcf(mcf, mcf.run(),
1.297 + gr, l1, cu, cc, s3, true, 0, "#A7");
1.298 + checkMcf(mcf, mcf.boundMaps(l2, u).run(),
1.299 + gr, l2, u, cc, s3, false, 0, "#A8");
1.300 +
1.301 + // Check the GEQ form
1.302 + mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
1.303 + checkMcf(mcf, mcf.run(),
1.304 + gr, l1, u, c, s4, true, 3530, "#A9", GEQ);
1.305 + mcf.problemType(mcf.GEQ);
1.306 + checkMcf(mcf, mcf.lowerMap(l2).run(),
1.307 + gr, l2, u, c, s4, true, 4540, "#A10", GEQ);
1.308 + mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
1.309 + checkMcf(mcf, mcf.run(),
1.310 + gr, l2, u, c, s5, false, 0, "#A11", GEQ);
1.311 +
1.312 + // Check the LEQ form
1.313 + mcf.reset().problemType(mcf.LEQ);
1.314 + mcf.upperMap(u).costMap(c).supplyMap(s5);
1.315 + checkMcf(mcf, mcf.run(),
1.316 + gr, l1, u, c, s5, true, 5080, "#A12", LEQ);
1.317 + checkMcf(mcf, mcf.lowerMap(l2).run(),
1.318 + gr, l2, u, c, s5, true, 5930, "#A13", LEQ);
1.319 + mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
1.320 + checkMcf(mcf, mcf.run(),
1.321 + gr, l2, u, c, s4, false, 0, "#A14", LEQ);
1.322 + }
1.323 +
1.324 + // B. Test NetworkSimplex with each pivot rule
1.325 + {
1.326 + NetworkSimplex<Digraph> mcf(gr);
1.327 + mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);
1.328 +
1.329 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),
1.330 + gr, l2, u, c, s1, true, 5970, "#B1");
1.331 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),
1.332 + gr, l2, u, c, s1, true, 5970, "#B2");
1.333 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),
1.334 + gr, l2, u, c, s1, true, 5970, "#B3");
1.335 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),
1.336 + gr, l2, u, c, s1, true, 5970, "#B4");
1.337 + checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),
1.338 + gr, l2, u, c, s1, true, 5970, "#B5");
1.339 + }
1.340 +
1.341 + return 0;
1.342 +}