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(sqrt(_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];