test/min_cost_flow_test.cc
changeset 711 cc61d09f053b
parent 689 111698359429
child 716 4faca85d40e6
     1.1 --- a/test/min_cost_flow_test.cc	Tue May 12 12:06:40 2009 +0200
     1.2 +++ b/test/min_cost_flow_test.cc	Tue May 12 12:08:06 2009 +0200
     1.3 @@ -174,7 +174,7 @@
     1.4             typename CM, typename SM, typename FM, typename PM >
     1.5  bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
     1.6                       const CM& cost, const SM& supply, const FM& flow, 
     1.7 -                     const PM& pi )
     1.8 +                     const PM& pi, SupplyType type )
     1.9  {
    1.10    TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    1.11  
    1.12 @@ -193,12 +193,50 @@
    1.13        sum += flow[e];
    1.14      for (InArcIt e(gr, n); e != INVALID; ++e)
    1.15        sum -= flow[e];
    1.16 -    opt = (sum == supply[n]) || (pi[n] == 0);
    1.17 +    if (type != LEQ) {
    1.18 +      opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);
    1.19 +    } else {
    1.20 +      opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);
    1.21 +    }
    1.22    }
    1.23    
    1.24    return opt;
    1.25  }
    1.26  
    1.27 +// Check whether the dual cost is equal to the primal cost
    1.28 +template < typename GR, typename LM, typename UM,
    1.29 +           typename CM, typename SM, typename PM >
    1.30 +bool checkDualCost( const GR& gr, const LM& lower, const UM& upper,
    1.31 +                    const CM& cost, const SM& supply, const PM& pi,
    1.32 +                    typename CM::Value total )
    1.33 +{
    1.34 +  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    1.35 +
    1.36 +  typename CM::Value dual_cost = 0;
    1.37 +  SM red_supply(gr);
    1.38 +  for (NodeIt n(gr); n != INVALID; ++n) {
    1.39 +    red_supply[n] = supply[n];
    1.40 +  }
    1.41 +  for (ArcIt a(gr); a != INVALID; ++a) {
    1.42 +    if (lower[a] != 0) {
    1.43 +      dual_cost += lower[a] * cost[a];
    1.44 +      red_supply[gr.source(a)] -= lower[a];
    1.45 +      red_supply[gr.target(a)] += lower[a];
    1.46 +    }
    1.47 +  }
    1.48 +  
    1.49 +  for (NodeIt n(gr); n != INVALID; ++n) {
    1.50 +    dual_cost -= red_supply[n] * pi[n];
    1.51 +  }
    1.52 +  for (ArcIt a(gr); a != INVALID; ++a) {
    1.53 +    typename CM::Value red_cost =
    1.54 +      cost[a] + pi[gr.source(a)] - pi[gr.target(a)];
    1.55 +    dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);
    1.56 +  }
    1.57 +  
    1.58 +  return dual_cost == total;
    1.59 +}
    1.60 +
    1.61  // Run a minimum cost flow algorithm and check the results
    1.62  template < typename MCF, typename GR,
    1.63             typename LM, typename UM,
    1.64 @@ -220,8 +258,10 @@
    1.65      check(checkFlow(gr, lower, upper, supply, flow, type),
    1.66            "The flow is not feasible " + test_id);
    1.67      check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
    1.68 -    check(checkPotential(gr, lower, upper, cost, supply, flow, pi),
    1.69 +    check(checkPotential(gr, lower, upper, cost, supply, flow, pi, type),
    1.70            "Wrong potentials " + test_id);
    1.71 +    check(checkDualCost(gr, lower, upper, cost, supply, pi, total),
    1.72 +          "Wrong dual cost " + test_id);
    1.73    }
    1.74  }
    1.75  
    1.76 @@ -266,45 +306,56 @@
    1.77      .node("target", w)
    1.78      .run();
    1.79    
    1.80 -  // Build a test digraph for testing negative costs
    1.81 -  Digraph ngr;
    1.82 -  Node n1 = ngr.addNode();
    1.83 -  Node n2 = ngr.addNode();
    1.84 -  Node n3 = ngr.addNode();
    1.85 -  Node n4 = ngr.addNode();
    1.86 -  Node n5 = ngr.addNode();
    1.87 -  Node n6 = ngr.addNode();
    1.88 -  Node n7 = ngr.addNode();
    1.89 +  // Build test digraphs with negative costs
    1.90 +  Digraph neg_gr;
    1.91 +  Node n1 = neg_gr.addNode();
    1.92 +  Node n2 = neg_gr.addNode();
    1.93 +  Node n3 = neg_gr.addNode();
    1.94 +  Node n4 = neg_gr.addNode();
    1.95 +  Node n5 = neg_gr.addNode();
    1.96 +  Node n6 = neg_gr.addNode();
    1.97 +  Node n7 = neg_gr.addNode();
    1.98    
    1.99 -  Arc a1 = ngr.addArc(n1, n2);
   1.100 -  Arc a2 = ngr.addArc(n1, n3);
   1.101 -  Arc a3 = ngr.addArc(n2, n4);
   1.102 -  Arc a4 = ngr.addArc(n3, n4);
   1.103 -  Arc a5 = ngr.addArc(n3, n2);
   1.104 -  Arc a6 = ngr.addArc(n5, n3);
   1.105 -  Arc a7 = ngr.addArc(n5, n6);
   1.106 -  Arc a8 = ngr.addArc(n6, n7);
   1.107 -  Arc a9 = ngr.addArc(n7, n5);
   1.108 +  Arc a1 = neg_gr.addArc(n1, n2);
   1.109 +  Arc a2 = neg_gr.addArc(n1, n3);
   1.110 +  Arc a3 = neg_gr.addArc(n2, n4);
   1.111 +  Arc a4 = neg_gr.addArc(n3, n4);
   1.112 +  Arc a5 = neg_gr.addArc(n3, n2);
   1.113 +  Arc a6 = neg_gr.addArc(n5, n3);
   1.114 +  Arc a7 = neg_gr.addArc(n5, n6);
   1.115 +  Arc a8 = neg_gr.addArc(n6, n7);
   1.116 +  Arc a9 = neg_gr.addArc(n7, n5);
   1.117    
   1.118 -  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);
   1.119 -  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);
   1.120 -  Digraph::NodeMap<int> ns(ngr, 0);
   1.121 +  Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);
   1.122 +  ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);
   1.123 +  Digraph::NodeMap<int> neg_s(neg_gr, 0);
   1.124    
   1.125 -  nl2[a7] =  1000;
   1.126 -  nl2[a8] = -1000;
   1.127 +  neg_l2[a7] =  1000;
   1.128 +  neg_l2[a8] = -1000;
   1.129    
   1.130 -  ns[n1] =  100;
   1.131 -  ns[n4] = -100;
   1.132 +  neg_s[n1] =  100;
   1.133 +  neg_s[n4] = -100;
   1.134    
   1.135 -  nc[a1] =  100;
   1.136 -  nc[a2] =   30;
   1.137 -  nc[a3] =   20;
   1.138 -  nc[a4] =   80;
   1.139 -  nc[a5] =   50;
   1.140 -  nc[a6] =   10;
   1.141 -  nc[a7] =   80;
   1.142 -  nc[a8] =   30;
   1.143 -  nc[a9] = -120;
   1.144 +  neg_c[a1] =  100;
   1.145 +  neg_c[a2] =   30;
   1.146 +  neg_c[a3] =   20;
   1.147 +  neg_c[a4] =   80;
   1.148 +  neg_c[a5] =   50;
   1.149 +  neg_c[a6] =   10;
   1.150 +  neg_c[a7] =   80;
   1.151 +  neg_c[a8] =   30;
   1.152 +  neg_c[a9] = -120;
   1.153 +
   1.154 +  Digraph negs_gr;
   1.155 +  Digraph::NodeMap<int> negs_s(negs_gr);
   1.156 +  Digraph::ArcMap<int> negs_c(negs_gr);
   1.157 +  ConstMap<Arc, int> negs_l(0), negs_u(1000);
   1.158 +  n1 = negs_gr.addNode();
   1.159 +  n2 = negs_gr.addNode();
   1.160 +  negs_s[n1] = 100;
   1.161 +  negs_s[n2] = -300;
   1.162 +  negs_c[negs_gr.addArc(n1, n2)] = -1;
   1.163 +
   1.164  
   1.165    // A. Test NetworkSimplex with the default pivot rule
   1.166    {
   1.167 @@ -342,7 +393,7 @@
   1.168      mcf.supplyType(mcf.GEQ);
   1.169      checkMcf(mcf, mcf.lowerMap(l2).run(),
   1.170               gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);
   1.171 -    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);
   1.172 +    mcf.supplyMap(s6);
   1.173      checkMcf(mcf, mcf.run(),
   1.174               gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);
   1.175  
   1.176 @@ -353,20 +404,26 @@
   1.177               gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);
   1.178      checkMcf(mcf, mcf.lowerMap(l2).run(),
   1.179               gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);
   1.180 -    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);
   1.181 +    mcf.supplyMap(s5);
   1.182      checkMcf(mcf, mcf.run(),
   1.183               gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);
   1.184  
   1.185      // Check negative costs
   1.186 -    NetworkSimplex<Digraph> nmcf(ngr);
   1.187 -    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);
   1.188 -    checkMcf(nmcf, nmcf.run(),
   1.189 -      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");
   1.190 -    checkMcf(nmcf, nmcf.upperMap(nu2).run(),
   1.191 -      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");
   1.192 -    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);
   1.193 -    checkMcf(nmcf, nmcf.run(),
   1.194 -      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");
   1.195 +    NetworkSimplex<Digraph> neg_mcf(neg_gr);
   1.196 +    neg_mcf.lowerMap(neg_l1).costMap(neg_c).supplyMap(neg_s);
   1.197 +    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u1,
   1.198 +      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A16");
   1.199 +    neg_mcf.upperMap(neg_u2);
   1.200 +    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,
   1.201 +      neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");
   1.202 +    neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);
   1.203 +    checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,
   1.204 +      neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");
   1.205 +      
   1.206 +    NetworkSimplex<Digraph> negs_mcf(negs_gr);
   1.207 +    negs_mcf.costMap(negs_c).supplyMap(negs_s);
   1.208 +    checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,
   1.209 +      negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);
   1.210    }
   1.211  
   1.212    // B. Test NetworkSimplex with each pivot rule