RhodeCode

    ///

    /// \pre \ref run() must be called before using this function.

    const PotentialMap& potentialMap() const {

      return *_potential_map;

    }

    /// @}

  private:

    // Initialize internal data structures

    bool init() {

      // Initialize result maps

      if (!_flow_map) {

        _flow_map = new FlowMap(_graph);

        _local_flow = true;

      }

      if (!_potential_map) {

        _potential_map = new PotentialMap(_graph);

        _local_potential = true;

      }

      // Initialize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      int all_node_num = _node_num + 1;

      int all_arc_num = _arc_num + _node_num;

      if (_node_num == 0) return false;

      _arc_ref.resize(_arc_num);

      _source.resize(all_arc_num);

      _target.resize(all_arc_num);

      _cap.resize(all_arc_num);

      _cost.resize(all_arc_num);

      _supply.resize(all_node_num);

      _flow.resize(all_arc_num);

      _pi.resize(all_node_num);

      _parent.resize(all_node_num);

      _pred.resize(all_node_num);

      _forward.resize(all_node_num);

      _thread.resize(all_node_num);

      _rev_thread.resize(all_node_num);

      _succ_num.resize(all_node_num);

      _last_succ.resize(all_node_num);

      _state.resize(all_arc_num);

      // Initialize node related data

      bool valid_supply = true;

      Flow sum_supply = 0;

      if (!_pstsup && !_psupply) {

        _pstsup = true;

        _psource = _ptarget = NodeIt(_graph);

        _pstflow = 0;

      }

      if (_psupply) {

        int i = 0;

        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

          _node_id[n] = i;

          _supply[i] = (*_psupply)[n];

          sum_supply += _supply[i];

        }

        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||

                       (_ptype == LEQ && sum_supply >= 0);

      } else {

        int i = 0;

        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

          _node_id[n] = i;

          _supply[i] = 0;

        }

        _supply[_node_id[_psource]] =  _pstflow;

        _supply[_node_id[_ptarget]] = -_pstflow;

      }

      if (!valid_supply) return false;

      // Infinite capacity value

      Flow inf_cap =

        std::numeric_limits<Flow>::has_infinity ?

        std::numeric_limits<Flow>::infinity() :

        std::numeric_limits<Flow>::max();

      // Initialize artifical cost

      Cost art_cost;

      if (std::numeric_limits<Cost>::is_exact) {

        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;

      } else {

        art_cost = std::numeric_limits<Cost>::min();

        for (int i = 0; i != _arc_num; ++i) {

          if (_cost[i] > art_cost) art_cost = _cost[i];

        }

        art_cost = (art_cost + 1) * _node_num;

      }

      // Run Circulation to check if a feasible solution exists

      typedef ConstMap<Arc, Flow> ConstArcMap;

      FlowNodeMap *csup = NULL;

      bool local_csup = false;

      if (_psupply) {

        csup = _psupply;

      } else {

        csup = new FlowNodeMap(_graph, 0);

        (*csup)[_psource] =  _pstflow;

        (*csup)[_ptarget] = -_pstflow;

        local_csup = true;

      }

      bool circ_result = false;

      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {

        // GEQ problem type

        if (_plower) {

          if (_pupper) {

            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>

              circ(_graph, *_plower, *_pupper, *csup);

            circ_result = circ.run();

          } else {

            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>

            circ_result = circ.run();

          }

        } else {

          if (_pupper) {

            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>

            circ_result = circ.run();

          } else {

            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>

            circ_result = circ.run();

          }

        }

      } else {

        // LEQ problem type

        typedef ReverseDigraph<const GR> RevGraph;

        typedef NegMap<FlowNodeMap> NegNodeMap;

        RevGraph rgraph(_graph);

        NegNodeMap neg_csup(*csup);

        if (_plower) {

          if (_pupper) {

            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>

              circ(rgraph, *_plower, *_pupper, neg_csup);

            circ_result = circ.run();

          } else {

            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>

            circ_result = circ.run();

          }

        } else {

          if (_pupper) {

            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>

            circ_result = circ.run();

          } else {

            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>

            circ_result = circ.run();

          }

        }

      }

      if (local_csup) delete csup;

      if (!circ_result) return false;

      // Set data for the artificial root node

      _root = _node_num;

      _parent[_root] = -1;

      _pred[_root] = -1;

      _thread[_root] = 0;

      _rev_thread[0] = _root;

      _succ_num[_root] = all_node_num;

      _last_succ[_root] = _root - 1;

      _supply[_root] = -sum_supply;

      if (sum_supply < 0) {

        _pi[_root] = -art_cost;

      } else {

        _pi[_root] = art_cost;

      }

      // Store the arcs in a mixed order

      int k = std::max(int(std::sqrt(double(_arc_num))), 10);

      int i = 0;

      for (ArcIt e(_graph); e != INVALID; ++e) {

        _arc_ref[i] = e;

        if ((i += k) >= _arc_num) i = (i % k) + 1;

      }

      // Initialize arc maps

      if (_pupper && _pcost) {

        for (int i = 0; i != _arc_num; ++i) {

          Arc e = _arc_ref[i];

          _source[i] = _node_id[_graph.source(e)];

          _target[i] = _node_id[_graph.target(e)];

          _cap[i] = (*_pupper)[e];

          _cost[i] = (*_pcost)[e];

          _flow[i] = 0;

          _state[i] = STATE_LOWER;

        }

      } else {

        for (int i = 0; i != _arc_num; ++i) {

          Arc e = _arc_ref[i];

          _source[i] = _node_id[_graph.source(e)];

          _target[i] = _node_id[_graph.target(e)];

          _flow[i] = 0;

          _state[i] = STATE_LOWER;

        }

        if (_pupper) {

          for (int i = 0; i != _arc_num; ++i)

            _cap[i] = (*_pupper)[_arc_ref[i]];

        } else {

          for (int i = 0; i != _arc_num; ++i)

            _cap[i] = inf_cap;

        }

        if (_pcost) {

          for (int i = 0; i != _arc_num; ++i)

            _cost[i] = (*_pcost)[_arc_ref[i]];

        } else {

          for (int i = 0; i != _arc_num; ++i)

            _cost[i] = 1;

        }

      }

      // Remove non-zero lower bounds

      if (_plower) {

        for (int i = 0; i != _arc_num; ++i) {

          Flow c = (*_plower)[_arc_ref[i]];

          if (c != 0) {

            _cap[i] -= c;

            _supply[_source[i]] -= c;

            _supply[_target[i]] += c;

          }

        }

      }

      // Add artificial arcs and initialize the spanning tree data structure

      for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

        _thread[u] = u + 1;

        _rev_thread[u + 1] = u;

        _succ_num[u] = 1;

        _last_succ[u] = u;

        _parent[u] = _root;

        _pred[u] = e;

        _cost[e] = art_cost;

        _cap[e] = inf_cap;

        _state[e] = STATE_TREE;

        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {

          _flow[e] = _supply[u];

          _forward[u] = true;

          _pi[u] = -art_cost + _pi[_root];

        } else {

          _flow[e] = -_supply[u];

          _forward[u] = false;

          _pi[u] = art_cost + _pi[_root];

    bool b;

    typename MCF::FlowMap &flow;

    typename MCF::PotentialMap &pot;

  };

};

// Check the feasibility of the given flow (primal soluiton)

template < typename GR, typename LM, typename UM,

           typename SM, typename FM >

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

                const SM& supply, const FM& flow,

                ProblemType type = EQ )

{

  TEMPLATE_DIGRAPH_TYPEDEFS(GR);

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

    if (flow[e] < lower[e] || flow[e] > upper[e]) return false;

  }

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

    typename SM::Value sum = 0;

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

      sum += flow[e];

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

      sum -= flow[e];

    bool b = (type ==  EQ && sum == supply[n]) ||

             (type == GEQ && sum >= supply[n]) ||

             (type == LEQ && sum <= supply[n]);

    if (!b) return false;

  }

  return true;

}

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

// using the "Complementary Slackness" optimality condition

template < typename GR, typename LM, typename UM,

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

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

                     const CM& cost, const SM& supply, const FM& flow, 

                     const PM& pi )

{

  TEMPLATE_DIGRAPH_TYPEDEFS(GR);

  bool opt = true;

  for (ArcIt e(gr); opt && e != INVALID; ++e) {

    typename CM::Value red_cost =

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

    opt = red_cost == 0 ||

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

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

  }

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

    typename SM::Value sum = 0;

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

      sum += flow[e];

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

      sum -= flow[e];

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

  }

  return opt;

}

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

template < typename MCF, typename GR,

           typename LM, typename UM,

           typename CM, typename SM >

void checkMcf( const MCF& mcf, bool mcf_result,

               const GR& gr, const LM& lower, const UM& upper,

               const CM& cost, const SM& supply,

               bool result, typename CM::Value total,

               const std::string &test_id = "",

               ProblemType type = EQ )

{

  check(mcf_result == result, "Wrong result " + test_id);

  if (result) {

    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),

          "The flow is not feasible " + test_id);

    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);

    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),

                         mcf.potentialMap()),

          "Wrong potentials " + test_id);

  }

}

int main()

{

  // Check the interfaces

  {

    typedef int Flow;

    typedef int Cost;

    checkConcept< McfClassConcept<GR, Flow, Cost>,

                  NetworkSimplex<GR, Flow, Cost> >();

  }

  // Run various MCF tests

  typedef ListDigraph Digraph;

  DIGRAPH_TYPEDEFS(ListDigraph);

  // Read the test digraph

  Digraph gr;

  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);

  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);

  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());

  Node v, w;

  std::istringstream input(test_lgf);

  DigraphReader<Digraph>(gr, input)

    .arcMap("cost", c)

    .arcMap("cap", u)

    .arcMap("low1", l1)

    .arcMap("low2", l2)

    .nodeMap("sup1", s1)

    .nodeMap("sup2", s2)

    .nodeMap("sup3", s3)

    .nodeMap("sup4", s4)

    .nodeMap("sup5", s5)

    .node("source", v)

    .node("target", w)

    .run();

  // A. Test NetworkSimplex with the default pivot rule

  {

    NetworkSimplex<Digraph> mcf(gr);

    // Check the equality form

    mcf.upperMap(u).costMap(c);

    checkMcf(mcf, mcf.supplyMap(s1).run(),

             gr, l1, u, c, s1, true,  5240, "#A1");

    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),

             gr, l1, u, c, s2, true,  7620, "#A2");

    mcf.lowerMap(l2);

    checkMcf(mcf, mcf.supplyMap(s1).run(),

             gr, l2, u, c, s1, true,  5970, "#A3");

    checkMcf(mcf, mcf.stSupply(v, w, 27).run(),

             gr, l2, u, c, s2, true,  8010, "#A4");

    mcf.reset();

    checkMcf(mcf, mcf.supplyMap(s1).run(),

             gr, l1, cu, cc, s1, true,  74, "#A5");

    checkMcf(mcf, mcf.lowerMap(l2).stSupply(v, w, 27).run(),

             gr, l2, cu, cc, s2, true,  94, "#A6");

    mcf.reset();

    checkMcf(mcf, mcf.run(),

             gr, l1, cu, cc, s3, true,   0, "#A7");

    checkMcf(mcf, mcf.boundMaps(l2, u).run(),

             gr, l2, u, cc, s3, false,   0, "#A8");

    // Check the GEQ form

    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);

    checkMcf(mcf, mcf.run(),

             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);

    mcf.problemType(mcf.GEQ);

    checkMcf(mcf, mcf.lowerMap(l2).run(),

             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);

    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);

    checkMcf(mcf, mcf.run(),

             gr, l2, u, c, s5, false,    0, "#A11", GEQ);

    // Check the LEQ form

    mcf.reset().problemType(mcf.LEQ);

    mcf.upperMap(u).costMap(c).supplyMap(s5);

    checkMcf(mcf, mcf.run(),

             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);

    checkMcf(mcf, mcf.lowerMap(l2).run(),

             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);

    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);

    checkMcf(mcf, mcf.run(),

             gr, l2, u, c, s4, false,    0, "#A14", LEQ);

  }

  // B. Test NetworkSimplex with each pivot rule

  {

    NetworkSimplex<Digraph> mcf(gr);

    mcf.supplyMap(s1).costMap(c).capacityMap(u).lowerMap(l2);

    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::FIRST_ELIGIBLE),

             gr, l2, u, c, s1, true,  5970, "#B1");

    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BEST_ELIGIBLE),

             gr, l2, u, c, s1, true,  5970, "#B2");

    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::BLOCK_SEARCH),

             gr, l2, u, c, s1, true,  5970, "#B3");

    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::CANDIDATE_LIST),

             gr, l2, u, c, s1, true,  5970, "#B4");

    checkMcf(mcf, mcf.run(NetworkSimplex<Digraph>::ALTERING_LIST),

             gr, l2, u, c, s1, true,  5970, "#B5");

  }

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

 *

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

 *

 * Copyright (C) 2003-2009

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

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

 *

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

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

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

///\ingroup tools

///\file

///\brief DIMACS problem solver.

///

/// This program solves various problems given in DIMACS format.

///

/// See

/// \code

///   dimacs-solver --help

/// \endcode

/// for more info on usage.

#include <iostream>

#include <fstream>

#include <cstring>

#include <lemon/smart_graph.h>

#include <lemon/dimacs.h>

#include <lemon/lgf_writer.h>

#include <lemon/time_measure.h>

#include <lemon/arg_parser.h>

#include <lemon/error.h>

#include <lemon/dijkstra.h>

#include <lemon/preflow.h>

#include <lemon/matching.h>

#include <lemon/network_simplex.h>

using namespace lemon;

typedef SmartDigraph Digraph;

DIGRAPH_TYPEDEFS(Digraph);

typedef SmartGraph Graph;

template<class Value>

void solve_sp(ArgParser &ap, std::istream &is, std::ostream &,

              DimacsDescriptor &desc)

{

  bool report = !ap.given("q");

  Digraph g;

  Node s;

  Digraph::ArcMap<Value> len(g);

  Timer t;

  t.restart();

  readDimacsSp(is, g, len, s, desc);

  if(report) std::cerr << "Read the file: " << t << '\n';

  t.restart();

  Dijkstra<Digraph, Digraph::ArcMap<Value> > dij(g,len);

  if(report) std::cerr << "Setup Dijkstra class: " << t << '\n';

  t.restart();

  dij.run(s);

  if(report) std::cerr << "Run Dijkstra: " << t << '\n';

}

template<class Value>

void solve_max(ArgParser &ap, std::istream &is, std::ostream &,

               Value infty, DimacsDescriptor &desc)

{

  bool report = !ap.given("q");

  Digraph g;

  Node s,t;

  Digraph::ArcMap<Value> cap(g);

  Timer ti;

  ti.restart();

  readDimacsMax(is, g, cap, s, t, infty, desc);

  if(report) std::cerr << "Read the file: " << ti << '\n';

  ti.restart();

  Preflow<Digraph, Digraph::ArcMap<Value> > pre(g,cap,s,t);

  if(report) std::cerr << "Setup Preflow class: " << ti << '\n';

  ti.restart();

  pre.run();

  if(report) std::cerr << "Run Preflow: " << ti << '\n';

  if(report) std::cerr << "\nMax flow value: " << pre.flowValue() << '\n';  

}

template<class Value>

void solve_min(ArgParser &ap, std::istream &is, std::ostream &,

{

  bool report = !ap.given("q");

  Digraph g;

  Digraph::ArcMap<Value> lower(g), cap(g), cost(g);

  Digraph::NodeMap<Value> sup(g);

  Timer ti;

  ti.restart();

  if (report) std::cerr << "Read the file: " << ti << '\n';

  ti.restart();

  NetworkSimplex<Digraph, Value> ns(g);

  ns.lowerMap(lower).capacityMap(cap).costMap(cost).supplyMap(sup);

  if (report) std::cerr << "Setup NetworkSimplex class: " << ti << '\n';

  ti.restart();

}

void solve_mat(ArgParser &ap, std::istream &is, std::ostream &,

              DimacsDescriptor &desc)

{

  bool report = !ap.given("q");

  Graph g;

  Timer ti;

  ti.restart();

  readDimacsMat(is, g, desc);

  if(report) std::cerr << "Read the file: " << ti << '\n';

  ti.restart();

  MaxMatching<Graph> mat(g);

  if(report) std::cerr << "Setup MaxMatching class: " << ti << '\n';

  ti.restart();

  mat.run();

  if(report) std::cerr << "Run MaxMatching: " << ti << '\n';

  if(report) std::cerr << "\nCardinality of max matching: "

                       << mat.matchingSize() << '\n';  

}

template<class Value>

void solve(ArgParser &ap, std::istream &is, std::ostream &os,

           DimacsDescriptor &desc)

{

  std::stringstream iss(static_cast<std::string>(ap["infcap"]));

  Value infty;

  iss >> infty;

  if(iss.fail())

    {

      std::cerr << "Cannot interpret '"

                << static_cast<std::string>(ap["infcap"]) << "' as infinite"

                << std::endl;

      exit(1);

    }

  switch(desc.type)

    {

    case DimacsDescriptor::MIN:

      break;

    case DimacsDescriptor::MAX:

      solve_max<Value>(ap,is,os,infty,desc);

      break;

    case DimacsDescriptor::SP:

      solve_sp<Value>(ap,is,os,desc);

      break;

    case DimacsDescriptor::MAT:

      solve_mat(ap,is,os,desc);

      break;

    default:

      break;

    }

}

int main(int argc, const char *argv[]) {

  typedef SmartDigraph Digraph;

  typedef Digraph::Arc Arc;

  std::string inputName;

  std::string outputName;

  ArgParser ap(argc, argv);

  ap.other("[INFILE [OUTFILE]]",

           "If either the INFILE or OUTFILE file is missing the standard\n"

           "     input/output will be used instead.")

    .boolOption("q", "Do not print any report")

    .boolOption("int","Use 'int' for capacities, costs etc. (default)")

    .optionGroup("datatype","int")

#ifdef HAVE_LONG_LONG

    .boolOption("long","Use 'long long' for capacities, costs etc.")

    .optionGroup("datatype","long")

#endif

    .boolOption("double","Use 'double' for capacities, costs etc.")

    .optionGroup("datatype","double")

    .boolOption("ldouble","Use 'long double' for capacities, costs etc.")

    .optionGroup("datatype","ldouble")

    .onlyOneGroup("datatype")

    .stringOption("infcap","Value used for 'very high' capacities","0")

    .run();

  std::ifstream input;

  std::ofstream output;

  switch(ap.files().size())

    {

    case 2:

      output.open(ap.files()[1].c_str());

      if (!output) {

        throw IoError("Cannot open the file for writing", ap.files()[1]);

      }

    case 1:

      input.open(ap.files()[0].c_str());

      if (!input) {

        throw IoError("File cannot be found", ap.files()[0]);

      }

    case 0:

      break;

    default:

      std::cerr << ap.commandName() << ": too many arguments\n";

      return 1;

    }

  std::istream& is = (ap.files().size()<1 ? std::cin : input);

  std::ostream& os = (ap.files().size()<2 ? std::cout : output);

  DimacsDescriptor desc = dimacsType(is);

  if(!ap.given("q"))

    {

      std::cout << "Problem type: ";

      switch(desc.type)

        {

        case DimacsDescriptor::MIN:

          std::cout << "min";

          break;

        case DimacsDescriptor::MAX:

          std::cout << "max";

          break;

        case DimacsDescriptor::SP:

          std::cout << "sp";

        case DimacsDescriptor::MAT:

          std::cout << "mat";

          break;

        default:

          exit(1);

          break;

        }

      std::cout << "\nNum of nodes: " << desc.nodeNum;

      std::cout << "\nNum of arcs:  " << desc.edgeNum;

      std::cout << "\n\n";

    }

  if(ap.given("double"))

    solve<double>(ap,is,os,desc);

  else if(ap.given("ldouble"))