/* -*- mode: C++; indent-tabs-mode: nil; -*-

 /* -*- mode: C++; indent-tabs-mode: nil; -*-

  *

  *

  * This file is a part of LEMON, a generic C++ optimization library.

  * This file is a part of LEMON, a generic C++ optimization library.

  *

  *

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  * (Egervary Research Group on Combinatorial Optimization, EGRES).

  *

  *

  * Permission to use, modify and distribute this software is granted

  * Permission to use, modify and distribute this software is granted

  * provided that this copyright notice appears in all copies. For

  * provided that this copyright notice appears in all copies. For

   "    8     0    0    0    0    0    3\n"

   "    8     0    0    0    0    0    3\n"

   "    9     3    0    0    0    0    0\n"

   "    9     3    0    0    0    0    0\n"

   "   10    -2    0    0    0   -7   -2\n"

   "   10    -2    0    0    0   -7   -2\n"

   "   11     0    0    0    0  -10    0\n"

   "   11     0    0    0    0  -10    0\n"

   "   12   -20  -27    0  -30  -30  -20\n"

   "   12   -20  -27    0  -30  -30  -20\n"

   "@arcs\n"

   "@arcs\n"

   "       cost  cap low1 low2 low3\n"

   "       cost  cap low1 low2 low3\n"

   " 1  2    70   11    0    8    8\n"

   " 1  2    70   11    0    8    8\n"

   " 1  3   150    3    0    1    0\n"

   " 1  3   150    3    0    1    0\n"

   " 1  4    80   15    0    2    2\n"

   " 1  4    80   15    0    2    2\n"

   template <typename MCF>

   template <typename MCF>

   struct Constraints {

   struct Constraints {

     void constraints() {

     void constraints() {

       checkConcept<concepts::Digraph, GR>();

       checkConcept<concepts::Digraph, GR>();

       const Constraints& me = *this;

       const Constraints& me = *this;

       MCF mcf(me.g);

       MCF mcf(me.g);

       const MCF& const_mcf = mcf;

       const MCF& const_mcf = mcf;

     typedef concepts::ReadMap<Node, Value> NM;

     typedef concepts::ReadMap<Node, Value> NM;

     typedef concepts::ReadMap<Arc, Value> VAM;

     typedef concepts::ReadMap<Arc, Value> VAM;

     typedef concepts::ReadMap<Arc, Cost> CAM;

     typedef concepts::ReadMap<Arc, Cost> CAM;

     typedef concepts::WriteMap<Arc, Value> FlowMap;

     typedef concepts::WriteMap<Arc, Value> FlowMap;

     typedef concepts::WriteMap<Node, Cost> PotMap;

     typedef concepts::WriteMap<Node, Cost> PotMap;

     GR g;

     GR g;

     VAM lower;

     VAM lower;

     VAM upper;

     VAM upper;

     CAM cost;

     CAM cost;

     NM sup;

     NM sup;

 // Check the feasibility of the given potentials (dual soluiton)

 // Check the feasibility of the given potentials (dual soluiton)

 // using the "Complementary Slackness" optimality condition

 // using the "Complementary Slackness" optimality condition

 template < typename GR, typename LM, typename UM,

 template < typename GR, typename LM, typename UM,

            typename CM, typename SM, typename FM, typename PM >

            typename CM, typename SM, typename FM, typename PM >

 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,

 bool checkPotential( const GR& gr, const LM& lower, const UM& upper,

                      const PM& pi, SupplyType type )

                      const PM& pi, SupplyType type )

 {

 {

   TEMPLATE_DIGRAPH_TYPEDEFS(GR);

   TEMPLATE_DIGRAPH_TYPEDEFS(GR);

   bool opt = true;

   bool opt = true;

       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];

       cost[e] + pi[gr.source(e)] - pi[gr.target(e)];

     opt = red_cost == 0 ||

     opt = red_cost == 0 ||

           (red_cost > 0 && flow[e] == lower[e]) ||

           (red_cost > 0 && flow[e] == lower[e]) ||

           (red_cost < 0 && flow[e] == upper[e]);

           (red_cost < 0 && flow[e] == upper[e]);

   }

   }

   for (NodeIt n(gr); opt && n != INVALID; ++n) {

   for (NodeIt n(gr); opt && n != INVALID; ++n) {

     typename SM::Value sum = 0;

     typename SM::Value sum = 0;

     for (OutArcIt e(gr, n); e != INVALID; ++e)

     for (OutArcIt e(gr, n); e != INVALID; ++e)

       sum += flow[e];

       sum += flow[e];

     for (InArcIt e(gr, n); e != INVALID; ++e)

     for (InArcIt e(gr, n); e != INVALID; ++e)

       opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);

       opt = (pi[n] <= 0) && (sum == supply[n] || pi[n] == 0);

     } else {

     } else {

       opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);

       opt = (pi[n] >= 0) && (sum == supply[n] || pi[n] == 0);

     }

     }

   }

   }

   return opt;

   return opt;

 }

 }

 // Check whether the dual cost is equal to the primal cost

 // Check whether the dual cost is equal to the primal cost

 template < typename GR, typename LM, typename UM,

 template < typename GR, typename LM, typename UM,

       dual_cost += lower[a] * cost[a];

       dual_cost += lower[a] * cost[a];

       red_supply[gr.source(a)] -= lower[a];

       red_supply[gr.source(a)] -= lower[a];

       red_supply[gr.target(a)] += lower[a];

       red_supply[gr.target(a)] += lower[a];

     }

     }

   }

   }

   for (NodeIt n(gr); n != INVALID; ++n) {

   for (NodeIt n(gr); n != INVALID; ++n) {

     dual_cost -= red_supply[n] * pi[n];

     dual_cost -= red_supply[n] * pi[n];

   }

   }

   for (ArcIt a(gr); a != INVALID; ++a) {

   for (ArcIt a(gr); a != INVALID; ++a) {

     typename CM::Value red_cost =

     typename CM::Value red_cost =

       cost[a] + pi[gr.source(a)] - pi[gr.target(a)];

       cost[a] + pi[gr.source(a)] - pi[gr.target(a)];

     dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);

     dual_cost -= (upper[a] - lower[a]) * std::max(-red_cost, 0);

   }

   }

   return dual_cost == total;

   return dual_cost == total;

 }

 }

 // Run a minimum cost flow algorithm and check the results

 // Run a minimum cost flow algorithm and check the results

 template < typename MCF, typename GR,

 template < typename MCF, typename GR,

     .nodeMap("sup5", s5)

     .nodeMap("sup5", s5)

     .nodeMap("sup6", s6)

     .nodeMap("sup6", s6)

     .node("source", v)

     .node("source", v)

     .node("target", w)

     .node("target", w)

     .run();

     .run();

   // Build test digraphs with negative costs

   // Build test digraphs with negative costs

   Digraph neg_gr;

   Digraph neg_gr;

   Node n1 = neg_gr.addNode();

   Node n1 = neg_gr.addNode();

   Node n2 = neg_gr.addNode();

   Node n2 = neg_gr.addNode();

   Node n3 = neg_gr.addNode();

   Node n3 = neg_gr.addNode();

   Node n4 = neg_gr.addNode();

   Node n4 = neg_gr.addNode();

   Node n5 = neg_gr.addNode();

   Node n5 = neg_gr.addNode();

   Node n6 = neg_gr.addNode();

   Node n6 = neg_gr.addNode();

   Node n7 = neg_gr.addNode();

   Node n7 = neg_gr.addNode();

   Arc a1 = neg_gr.addArc(n1, n2);

   Arc a1 = neg_gr.addArc(n1, n2);

   Arc a2 = neg_gr.addArc(n1, n3);

   Arc a2 = neg_gr.addArc(n1, n3);

   Arc a3 = neg_gr.addArc(n2, n4);

   Arc a3 = neg_gr.addArc(n2, n4);

   Arc a4 = neg_gr.addArc(n3, n4);

   Arc a4 = neg_gr.addArc(n3, n4);

   Arc a5 = neg_gr.addArc(n3, n2);

   Arc a5 = neg_gr.addArc(n3, n2);

   Arc a6 = neg_gr.addArc(n5, n3);

   Arc a6 = neg_gr.addArc(n5, n3);

   Arc a7 = neg_gr.addArc(n5, n6);

   Arc a7 = neg_gr.addArc(n5, n6);

   Arc a8 = neg_gr.addArc(n6, n7);

   Arc a8 = neg_gr.addArc(n6, n7);

   Arc a9 = neg_gr.addArc(n7, n5);

   Arc a9 = neg_gr.addArc(n7, n5);

   Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);

   Digraph::ArcMap<int> neg_c(neg_gr), neg_l1(neg_gr, 0), neg_l2(neg_gr, 0);

   ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);

   ConstMap<Arc, int> neg_u1(std::numeric_limits<int>::max()), neg_u2(5000);

   Digraph::NodeMap<int> neg_s(neg_gr, 0);

   Digraph::NodeMap<int> neg_s(neg_gr, 0);

   neg_l2[a7] =  1000;

   neg_l2[a7] =  1000;

   neg_l2[a8] = -1000;

   neg_l2[a8] = -1000;

   neg_s[n1] =  100;

   neg_s[n1] =  100;

   neg_s[n4] = -100;

   neg_s[n4] = -100;

   neg_c[a1] =  100;

   neg_c[a1] =  100;

   neg_c[a2] =   30;

   neg_c[a2] =   30;

   neg_c[a3] =   20;

   neg_c[a3] =   20;

   neg_c[a4] =   80;

   neg_c[a4] =   80;

   neg_c[a5] =   50;

   neg_c[a5] =   50;

     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,

     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l1, neg_u2,

       neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");

       neg_c, neg_s, neg_mcf.OPTIMAL, true,  -40000, "#A17");

     neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);

     neg_mcf.reset().lowerMap(neg_l2).costMap(neg_c).supplyMap(neg_s);

     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,

     checkMcf(neg_mcf, neg_mcf.run(), neg_gr, neg_l2, neg_u1,

       neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");

       neg_c, neg_s, neg_mcf.UNBOUNDED, false,    0, "#A18");

     NetworkSimplex<Digraph> negs_mcf(negs_gr);

     NetworkSimplex<Digraph> negs_mcf(negs_gr);

     negs_mcf.costMap(negs_c).supplyMap(negs_s);

     negs_mcf.costMap(negs_c).supplyMap(negs_s);

     checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,

     checkMcf(negs_mcf, negs_mcf.run(), negs_gr, negs_l, negs_u,

       negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);

       negs_c, negs_s, negs_mcf.OPTIMAL, true, -300, "#A19", GEQ);

   }

   }