tpl_netgraph.H

# ifndef TPL_NETGRAPH_H
# define TPL_NETGRAPH_H

# include <limits>
# include <tpl_dynDlist.H>
# include <tpl_dynListStack.H>
# include <tpl_dynBinHeap.H>
# include <tpl_dynSetTree.H>
# include <tpl_random_queue.H>
# include <Set.H>
# include <tpl_graph_utils.H>
# include <tpl_find_path.H>

using namespace Aleph;

namespace Aleph {

template <class N> class Res_F;


  template <typename Arc_Info, typename F_Type = double>
class Net_Arc : public Graph_Aarc<Arc_Info>
{
public:

  typedef F_Type Flow_Type;

  Flow_Type cap; 

  Flow_Type flow; 

  bool check_arc() const { return flow >= 0 and flow <= cap; }

  Net_Arc * img_arc;          

  bool is_residual;      

  Net_Arc() : img_arc(NULL), is_residual(false) { /* empty */ } 

  Net_Arc(Net_Arc * net_arc) 
    : Graph_Aarc<Arc_Info>(*net_arc), cap(net_arc->cap), flow(net_arc->flow), 
      img_arc(NULL), is_residual(false) 
  {
    /* empty */ 
  }

  Net_Arc(const Arc_Info & info)
    :  Graph_Aarc<Arc_Info>(info), img_arc(NULL), is_residual(false) 
  { /* empty */ } 

  Net_Arc & operator = (Net_Arc & arc)
  {
    if (this == &arc) 
      return *this;

    this->get_info() = arc.get_info();
    cap             = arc.cap;
    flow            = arc.flow;
    img_arc         = NULL;
    is_residual     = false;

    return *this;
  }
};

    template <typename Node_Info, typename F_Type = double>
class Net_Node : public Graph_Anode<Node_Info>
{
public:

  typedef F_Type Flow_Type;

  size_t in_degree;   

  Flow_Type out_cap;  

  Flow_Type in_cap;   

  Flow_Type out_flow; 

  Flow_Type in_flow;  

  Net_Node() : in_degree(0), out_cap(0), in_cap(0), out_flow(0), in_flow(0) 
  {
    // empty
  }

  Net_Node(const Node_Info & node_info) 
    : Graph_Anode<Node_Info>(node_info), in_degree(0),
      out_cap(0), in_cap(0), out_flow(0), in_flow(0) 
  {
    /* empty */ 
  }

  Net_Node(Net_Node * net_node) 
    : Graph_Anode<Node_Info>(*net_node), in_degree(net_node->in_degree),
      out_cap(net_node->out_cap), in_cap(net_node->in_cap), 
      out_flow(net_node->out_flow), in_flow(net_node->in_flow) 
  {
    /* empty */ 
  }
};

  template <class NodeT = Net_Node<Empty_Class, double>, 
            class ArcT  = Net_Arc<Empty_Class, double> >
class Net_Graph : public Array_Digraph<NodeT, ArcT>
{
public: 

  typedef Array_Digraph<NodeT, ArcT> Digraph;

  typedef ArcT Arc;

  typedef NodeT Node;

  typedef typename Arc::Flow_Type Flow_Type;

  typedef typename Node::Node_Type Node_Type;

  typedef typename Arc::Arc_Type Arc_Type;

public:

  mutable Flow_Type Infinity;

  Flow_Type get_in_cap(Node * node) const { return node->in_cap; }

  Flow_Type get_out_cap(Node * node) const { return node->out_cap; }

  size_t get_in_degree(Node * node) const { return node->in_degree; }

  size_t get_out_degree(Node * node) const { return get_num_arcs(node); }

  Flow_Type get_out_flow(Node * node) const { return node->out_flow; }

  Flow_Type get_in_flow(Node * node) const { return node->in_flow; }

  bool is_source(Node * node) { return src_nodes.count(node) > 0; }

  bool is_sink(Node * node)  { return sink_nodes.count(node) > 0; }

  bool is_connected(Node * node) const 
  {
    return node->in_degree != 0 or node->num_arcs != 0;
  }

  bool check_node(Node * node) 
  { 
    if (not is_connected(node)) 
      return false;
    
    if (is_sink(node)) 
      return node->out_flow == 0 and node->in_flow >= 0;

    if (is_source(node))
      return node->in_flow == 0 and node->out_flow >= 0;

    Flow_Type sum_out_cap = 0;
    Flow_Type sum_out_flow = 0;
    for (Node_Arc_Iterator<Net_Graph> it(node); it.has_curr(); it.next())
      {
        Arc * arc = it.get_current();
        if (arc->is_residual)
          continue;

        if (not arc->check_arc())
          return false;

        sum_out_cap += arc->cap;
        sum_out_flow += arc->flow;
      }

    if (sum_out_cap != node->out_cap or sum_out_flow != node->out_flow)
      return false;

    if (is_source(node))
      return true;

    return node->out_flow == node->in_flow;
  }

private:

  Aleph::set<Node*> src_nodes;
  Aleph::set<Node*> sink_nodes;

public:

  Aleph::set<Node*> & get_src_nodes() { return src_nodes; }

  Aleph::set<Node*> & get_sink_nodes() { return sink_nodes; }

public:

  void make_super_source()
  {
    if (src_nodes.size() == 1) 
      return;

    if (src_nodes.size() == 0)
      throw std::domain_error("network has not source node (it has cicles)");

    DynDlist<Node*> src_list;
    for (typename Aleph::set<Node*>::iterator it = src_nodes.begin(); 
         it != src_nodes.end(); ++it)
      src_list.append(*it);

    Node * super_source = insert_node(); 
    while (not src_list.is_empty())
      {
        Node * p = src_list.remove_first();
        insert_arc(super_source, p, get_out_cap(p));
      }

    with_super_source = true;
  }

  void unmake_super_source()
  {
    if (not with_super_source) 
      return;

    I(src_nodes.size() == 1);

    remove_node(*src_nodes.begin());
    with_super_source = false;
  }

  void make_super_sink()
  {
    if (sink_nodes.size() == 1) 
      return;

    if (sink_nodes.size() == 0)
      throw std::domain_error("network has not sink node (it has cicles)");

    DynDlist<Node*> sink_list;
    for (typename Aleph::set<Node*>::iterator it = sink_nodes.begin(); 
         it != sink_nodes.end(); ++it)
      sink_list.append(*it);

    Node * super_sink = insert_node(); 
    while (not sink_list.is_empty())
      {
        Node * p = sink_list.remove_first();
        insert_arc(p, super_sink, get_in_cap(p));
      }

    with_super_sink = true;
  }

  void unmake_super_sink()
  {
    if (not with_super_sink) 
      return;

    I(sink_nodes.size() == 1);

    remove_node(*sink_nodes.begin());
    with_super_sink = false;
  }

  void make_super_nodes()
  {
    make_super_source();
    try
      {
        make_super_sink();
      }
    catch (std::bad_alloc)
      {
        unmake_super_source();
      }
  }

  void unmake_super_nodes()
  {
    unmake_super_source();
    unmake_super_sink();
  }

  Node * get_source() { return *get_src_nodes().begin(); }

  Node * get_sink() { return *get_sink_nodes().begin(); }

  Node * insert_node(const Node_Type & node_info)
  {
    Node * p = Digraph::insert_node(node_info); 
    try
      {
        src_nodes.insert(p);
        try
          {
            sink_nodes.insert(p);
            return p;
          }
        catch (bad_alloc)
          {
            src_nodes.erase(p);
            Digraph::remove_node(p);
            throw;
          }
      }
    catch (bad_alloc)
      {
        Digraph::remove_node(p);
        throw; // propagar excepción 
      }
  }

  Node * insert_node() { return insert_node(Node_Type()); }

  Node * insert_node(Node * p)
  {
    Digraph::insert_node(p); 
    try
      {
        src_nodes.insert(p);
        try
          {
            sink_nodes.insert(p);
            return p;
          }
        catch (bad_alloc)
          {
            src_nodes.erase(p);
            Digraph::remove_node(p);
            throw;
          }
      }
    catch (bad_alloc)
      {
        Digraph::remove_node(p);
        throw; // propagar excepción 
      }
  }

  virtual Arc * insert_arc(Node * src_node, Node * tgt_node, 
                           const typename Arc::Arc_Type & arc_info,
                           const Flow_Type & cap, const Flow_Type & flow)
  {     // inserción en clase base
    Arc * arc = Digraph::insert_arc(src_node, tgt_node, arc_info); 

    src_nodes.erase(tgt_node);  // actualización de source/sink
    sink_nodes.erase(src_node); 

    arc->cap  = cap;            // ajuste capacidad y flujo de arco
    arc->flow = flow;

    if (not arc->check_arc())
      throw std::overflow_error("flow is greater than capacity");

    tgt_node->in_degree++;   // actualizar info de control de nodo
    src_node->out_cap  += cap;
    tgt_node->in_cap   += cap;
    src_node->out_flow += flow;   
    tgt_node->in_flow  += flow;

    return arc;
  }

  Arc * connect_arc(Arc * arc)
  {     
    Digraph::connect_arc(arc); 

    Node * src = this->get_src_node(arc);
    Node * tgt = this->get_tgt_node(arc);

    src_nodes.erase(tgt);        // elimina destino de src_nodes
    sink_nodes.erase(src);       // elimina fuente de sink_nodes

    tgt->in_degree++;            // actualiza información de control de nodo
    src->out_cap  += arc->cap;
    tgt->in_cap   += arc->cap;
    src->out_flow += arc->flow;   
    tgt->in_flow  += arc->flow;
    
    return arc;
  }

       virtual Arc * 
  insert_arc(Node * src_node, Node * tgt_node, const Flow_Type & cap)
  {
    return insert_arc(src_node, tgt_node, Arc_Type(), cap, 0);
  }

  virtual Arc * insert_arc(Node * src_node, Node * tgt_node, 
                           const typename Arc::Arc_Type & arc_info)
  {
    return insert_arc(src_node, tgt_node, arc_info, 0, 0);
  }

  virtual Arc * insert_arc(Node * src_node, Node * tgt_node)
  {
    return insert_arc(src_node, tgt_node, Arc_Type(), 0, 0);
  }

  virtual void remove_arc(Arc * arc)
  {
    Node * src = this->get_src_node(arc);
    Node * tgt = this->get_tgt_node(arc);

    if (--(tgt->in_degree) == 0)
      src_nodes.insert(tgt);  // tgt deviene un nodo fuente 

    src->out_cap  -= arc->cap; // actualizar caps y flujos en src y tgt
    src->out_flow -= arc->flow;
    tgt->in_cap   -= arc->cap;
    tgt->in_flow  -= arc->flow;

    Digraph::remove_arc(arc); // eliminar en clase base

    if (this->get_num_arcs(src) == 0)
      sink_nodes.insert(src); // src deviene un nodo sumidero
  }

  void disconnect_arc(Arc * arc)
  {
    Node * src = this->get_src_node(arc);
    Node * tgt = this->get_tgt_node(arc);
    if (--(tgt->in_degree) == 0)
      src_nodes.insert(tgt); // tgt deviene un nodo fuente 

    src->out_cap  -= arc->cap;   // actualizar caps y flujos en src y tgt
    src->out_flow -= arc->flow;
    tgt->in_cap   -= arc->cap;
    tgt->in_flow  -= arc->flow;

    Digraph::disconnect_arc(arc); // desconeción en clase base

    if (this->get_num_arcs(src) == 0)
      sink_nodes.insert(src); // src deviene un nodo sumidero
  }

  virtual void remove_node(Node * p)          
  {
    Digraph::remove_node(p);  // eliminación en clase base
    src_nodes.erase(p);
    sink_nodes.erase(p);
  }

  Net_Graph(Digraph & digraph)
  : with_super_source(false), with_super_sink(false), residual_net(false)
  {
    Net_Graph::Net_Graph();     // inicializa atributos 
    copy_graph(*this, digraph, true);  // copia mapeada
  }

  Net_Graph(Net_Graph & net) 
    : Array_Digraph<NodeT, ArcT>::Array_Digraph(),
      Infinity(numeric_limits<typename Arc::Flow_Type>::max()),
      with_super_source(net.with_super_source),
      with_super_sink(net.with_super_sink), residual_net(net.residual_net)
  {
    copy_graph(*this, net, false);            // copia sin mapear
  }

  ~Net_Graph()
  {
    if (residual_net)
      unmake_residual_net();
  }

  void set_cap(Arc * arc, const Flow_Type & cap) 
  {
    if (cap < arc->flow)
      throw std::out_of_range("capacity value is smaller than flow");

    const Flow_Type old_cap = arc->cap;
    arc->cap = cap;

    Node * src_node = get_src_node(arc);
    src_node->out_cap -= old_cap;
    src_node->out_cap += cap;

    Node * tgt_node = get_tgt_node(arc);
    tgt_node->in_cap -= old_cap;
    tgt_node->in_cap += cap;
  }

  void set_flow(Arc * arc, const Flow_Type & flow)
  {
    if (flow > arc->cap)
      throw std::out_of_range("flow value is greater than capacity");

    const Flow_Type old_flow = arc->flow;
    arc->flow = flow;

    Node * src_node = get_src_node(arc);
    src_node->out_flow -= old_flow;
    src_node->out_flow += flow;

    Node * tgt_node = get_tgt_node(arc);
    tgt_node->in_flow -= old_flow;
    tgt_node->in_flow += flow;
  }

  const Flow_Type & get_flow(Arc * arc) const { return arc->flow; }

  const Flow_Type & get_cap(Arc * arc) const { return arc->cap; }

  void reset()
  {
    for (Arc_Iterator<Net_Graph> it(*this); it.has_current(); it.next())
      it.get_current()->flow = 0;

    for (Node_Iterator<Net_Graph> it(*this); it.has_current(); it.next())
      {
        Node * p = it.get_current();
        p->in_flow = p->out_flow = 0;
      }
  }

  bool check_network() 
  {
    for (Node_Iterator<Net_Graph> it(*this); it.has_current(); it.next())
      if (not check_node(it.get_current()))
        return false;

    return true; // todos los nodos son válidos
  }

  Flow_Type flow_value() { return get_source()->out_flow; }

  bool with_super_source; 

  bool with_super_sink;

  bool residual_net;

  Net_Graph() 
    : Infinity(numeric_limits<typename Arc::Flow_Type>::max()),
      with_super_source(false), with_super_sink(false), residual_net(false)
  { /* empty */ }

  void insert_residual_arc(Arc * arc)
  {
    Arc * res_arc = Digraph::insert_arc(this->get_tgt_node(arc), 
                                        this->get_src_node(arc)); 

    res_arc->is_residual = true;
    res_arc->img_arc     = arc;
    res_arc->cap         = arc->cap;
    res_arc->flow        = arc->cap - arc->flow;

    arc->img_arc = res_arc;
  }

private:

  void remove_residual_arc(Arc * arc)
  {    
    I(arc->is_residual);
    
    Digraph::remove_arc(arc);
  }

public:

  void make_residual_net()
  {
    try
      {
        if (residual_net)
          throw std::domain_error("Residual net is currently computed");

        size_t n = this->get_num_arcs(); // num arcos residuales a insertar
        for (Arc_Iterator<Net_Graph> it(*this); it.has_curr() and n > 0; 
             it.next())
          {
            Arc * arc = it.get_current_arc();
            if (not arc->is_residual)
              {
                insert_residual_arc(arc);
                --n; 
              }
          }

        residual_net = true;
      }
    catch (...)
      {
        unmake_residual_net();
      }
  }

  void unmake_residual_net()
  {
    for (Arc_Iterator<Net_Graph> it(*this); it.has_current(); )
      {
        Arc * arc = it.get_current_arc();
        if (arc->is_residual)
          {
            it.next();
            remove_residual_arc(arc);
          }
        else
          it.next();
      }
    residual_net = false;
  }

public:

  bool is_there_residual_net() const { return residual_net; }

  Arc * get_residual_arc(Arc * a) 
  {
    if (not residual_net)
      throw std::domain_error("Residual net is not computed");

    return a->img_arc; 
  }

  void increase_out_flow(Node * p, const Flow_Type & flow)
  {
    if (is_sink(p))
      return;
    
    p->out_flow += flow;
  }

  void decrease_out_flow(Node * p, const Flow_Type & flow)
  {
    if (is_sink(p))
      return;

    p->out_flow -= flow;
  }

  void increase_in_flow(Node * p, const Flow_Type & flow)
  {
    if (is_source(p))
      return;

    p->in_flow += flow;
  }

  void decrease_in_flow(Node * p, const Flow_Type & flow)
  {
    if (is_source(p))
      return;

    p->in_flow -= flow;
  }
};

    template <class N> 
class Res_F
{
public:

  bool operator () (typename N::Arc * a) const
  {
    return (a->cap - a->flow) != 0;
  }
};

    template <class Net>
typename Net::Flow_Type increase_flow(Net & net, Path<Net> & path)
{
  if (not net.residual_net)
    throw std::domain_error("Residual net is not created");

  typename Net::Flow_Type slack = net.Infinity; // eslabón mínimo

      // calcular el eslabón mínimo del camino de aumento
  for (typename Path<Net>::Iterator it(path); it.has_current_arc(); it.next())
    {
      typename Net::Arc * arc = it.get_current_arc();
      const typename Net::Flow_Type w = arc->cap - arc->flow;

      I(w >= 0);

      if (w < slack)
        slack = w;
    }

  if (slack == 0)
    return slack;

      // aumentar el flujo de la red por el camino de aumento
  for (typename Path<Net>::Iterator it(path); it.has_current_arc(); it.next())
    {
      typename Net::Arc * arc = it.get_current_arc();
      typename Net::Arc * img = (typename Net::Arc*) arc->img_arc;
      arc->flow += slack; 
      img->flow -= slack;

      if (not arc->is_residual)
        {
          net.increase_out_flow(net.get_src_node(arc), slack);
          net.increase_in_flow(net.get_tgt_node(arc), slack);
        }
      else
        {
          net.decrease_in_flow(net.get_src_node(arc), slack);
          net.decrease_out_flow(net.get_tgt_node(arc), slack);
        }

      I(arc->check_arc());
    }

  return slack;
}

    template <class Net>
bool has_aumenting_path(Net & net) 
{
  if (not net.residual_net)
    throw std::domain_error("Residual net has not been generated");

  return test_for_path<Net, Res_F<Net> >(net, net.get_source(), net.get_sink());
}

    template <class Net, template <class, class> class Find_Path>
  typename Net::Flow_Type 
aumenting_path_maximum_flow(Net & net, const bool & leave_residual = false)
{
  net.make_super_nodes();
  try
    {
      net.make_residual_net();
    }
  catch (std::bad_alloc)
    {
      net.unmake_super_nodes();
      throw;
    }

  typename Net::Node * source = net.get_source();
  typename Net::Node * sink   = net.get_sink();
  try
    {
      while (true) // mientras exista un camino de aumento
        {
          Path<Net> path(net);
          if (not Find_Path<Net, Res_F<Net> > () (net, source, sink, path))
            break;

          increase_flow <Net> (net, path);
        }
    }
  catch (std::bad_alloc)
    {
      net.unmake_residual_net();
      net.unmake_super_nodes();
      throw;
    }

  const typename Net::Flow_Type ret_val = source->out_flow;
  if (not leave_residual)
    { 
      net.unmake_residual_net();
      net.unmake_super_nodes();
    }

  return ret_val;
}

      template <class Net> typename Net::Flow_Type 
ford_fulkerson_maximum_flow(Net & net, const bool & leave_residual = false)
{
  return aumenting_path_maximum_flow<Net, Find_Path_Depth_First>
    (net, leave_residual);
} 

    template <class Net> class 
Ford_Fulkerson_Maximum_Flow 
{
public:

      typename Net::Flow_Type 
  operator () (Net & net, const bool & leave_residual = false) const
  {
    return ford_fulkerson_maximum_flow(net, leave_residual);
  }
};

      template <class Net> typename Net::Flow_Type
edmonds_karp_maximum_flow(Net & net, const bool & leave_residual = false)
{
  return aumenting_path_maximum_flow<Net, Find_Path_Breadth_First>
    (net, leave_residual);
}

    template <class Net> 
class Edmonds_Karp_Maximum_Flow
{
public:

      typename Net::Flow_Type 
  operator () (Net & net, 
               const bool & leave_residual = false) const
  {
    return edmonds_karp_maximum_flow(net, leave_residual);
  }
};

    template <class Net> static
bool is_node_active(typename Net::Node * p) 
{
  return p->in_flow > p->out_flow; 
}

    template <class Net> static
long & node_height(typename Net::Node * p) { return NODE_COUNTER(p); }

    template <class N> 
class Res_Arc
{
public:

  bool operator () (typename N::Arc * a) const
  {
    return a->is_residual;
  }
};

    template <class Net> static
bool initial_height(Net & net, typename Net::Node * p, typename Net::Arc * a)
{
  node_height<Net>(p) = node_height<Net>(net.get_src_node(a)) + 1;
  return false;
}

    template <class Net> static 
void init_height_in_nodes(Net & net)
{
  breadth_first_traversal <Net, Res_Arc<Net>> 
    (net, net.get_sink(), &initial_height); 
}

    template <class Q_Type> static
void put_in_active_queue(Q_Type & q, typename Q_Type::Item_Type & p)
{
  if (NODE_BITS(p).get_bit(Aleph::Maximum_Flow))
    return;

  NODE_BITS(p).set_bit(Aleph::Maximum_Flow, true);
  q.put(p);
}

    template <class Q_Type> static 
typename Q_Type::Item_Type get_from_active_queue(Q_Type & q)
{
  typename Q_Type::Item_Type p = q.get();

  I(NODE_BITS(p).get_bit(Aleph::Maximum_Flow));

  NODE_BITS(p).set_bit(Aleph::Maximum_Flow, false);

  return p;
}


   template <class Q_Type> static
void remove_from_active_queue(Q_Type & q, typename Q_Type::Item_Type & p)
{
  NODE_BITS(p).set_bit(Aleph::Maximum_Flow, false);
  q.remove(p);
}

    template <class Net, class Q_Type> typename Net::Flow_Type
generic_preflow_vertex_push_maximum_flow
    (Net & net, const bool & leave_residual = false)
{
  net.make_super_nodes();
  try
    {
      net.make_residual_net();
      net.reset_counter_nodes();
    }
  catch (std::bad_alloc)
    {
      net.unmake_super_nodes();
      throw;
    }

  init_height_in_nodes(net); 
  typename Net::Node * source = net.get_source();
  typename Net::Node * sink   = net.get_sink();

  Q_Type q;  // instancia el conjunto de nodos activos
  try
    {
      for (Node_Arc_Iterator<Net, Res_F<Net> > it(source); it.has_curr(); 
           it.next())
        {
          typename Net::Arc * arc  = it.get_current_arc();  
          typename Net::Node * tgt = net.get_tgt_node(arc); 
          arc->flow = tgt->in_flow = arc->cap - arc->flow; // inunde arco
          arc->img_arc->flow = 0; 
          if (tgt != sink)
            put_in_active_queue(q, tgt); // tgt deviene activo
        }

      source->out_flow = source->out_cap;

      while (not q.is_empty()) // mientras haya nodos activos
        {
          typename Net::Node * src = get_from_active_queue(q);
          typename Net::Flow_Type excess = src->in_flow - src->out_flow;

          for (Node_Arc_Iterator<Net, Res_F<Net> > it(src); 
               it.has_curr() and excess > 0; it.next()) 
            {
              typename Net::Arc * arc  = it.get_current_arc();
              typename Net::Node * tgt = net.get_tgt_node(arc);
              if (node_height<Net>(src) != node_height<Net>(tgt) + 1)
                continue; // el nodo no es elegible

              const typename Net::Flow_Type flow_avail_in_arc = 
                arc->cap - arc->flow;   
              typename Net::Flow_Type flow_to_push = 
                std::min(flow_avail_in_arc, excess);

              arc->flow          += flow_to_push;
              arc->img_arc->flow -= flow_to_push;
              if (arc->is_residual)
                {
                  net.decrease_out_flow(tgt, flow_to_push);
                  net.decrease_in_flow(src, flow_to_push);
                }
              else
                {
                  net.increase_out_flow(src, flow_to_push);
                  net.increase_in_flow(tgt, flow_to_push);
                }

              excess -= flow_to_push;

              if (is_node_active<Net>(tgt) and tgt != sink and tgt != source) 
                put_in_active_queue(q, tgt);
        }

          if (excess > 0) // ¿src aún sigue activo?
            { // sí ==> incremente h(src) y re-inserte en q
              node_height<Net>(src)++;
              put_in_active_queue(q, src);
            }
        }

      const typename Net::Flow_Type ret_val = source->out_flow; 

      if (not leave_residual)
        {
          net.unmake_residual_net();
          net.unmake_super_nodes();
        }

      return ret_val;
    } 
  catch (...)
    {
      net.unmake_residual_net();
      net.unmake_super_nodes();
      throw;
    }
}

    template <class Net> typename Net::Flow_Type
fifo_preflow_maximum_flow
    (Net & net, const bool & leave_residual = false)
{
  return generic_preflow_vertex_push_maximum_flow
    <Net,DynListQueue<typename Net::Node*> >(net, leave_residual);
}

    template <class Net> 
class Fifo_Preflow_Maximum_Flow
{
public:
       typename Net::Flow_Type 
   operator () (Net & net, 
             const bool & leave_residual = false) const
  {
    return fifo_preflow_maximum_flow(net, leave_residual);
  }
};

    template <class Net> 
struct Compare_Height
{
  bool operator () (typename Net::Node * n1, typename Net::Node * n2) const
  {
    return node_height<Net>(n1) > node_height<Net>(n2);
  }
};

    template <class Net> typename Net::Flow_Type
heap_preflow_maximum_flow(Net & net, const bool & leave_residual = false)
{
  return generic_preflow_vertex_push_maximum_flow 
    <Net, DynBinHeap<typename Net::Node*, Compare_Height<Net> > > 
     (net, leave_residual);
}

   template <class Net> 
class Heap_Preflow_Maximum_Flow
{
public:

      typename Net::Flow_Type
  operator () (Net & net, 
               const bool & leave_residual = false) const
  {
    return heap_preflow_maximum_flow(net, leave_residual);
  }
};

    template <class Net> typename Net::Flow_Type
random_preflow_maximum_flow(Net & net, const bool & leave_residual = false)
{
  return generic_preflow_vertex_push_maximum_flow
    <Net, Random_Set<typename Net::Node*> > (net, leave_residual);
}

template <class Net> class Random_Preflow_Maximum_Flow
{
public:
       typename Net::Flow_Type 
   operator () (Net & net,
                const bool & leave_residual = false) const
  {
    return random_preflow_maximum_flow(net, leave_residual);
  }
};


    template <class Net> static
bool has_arcs_in_active_queue(typename Net::Node * p)
{
  for (Node_Arc_Iterator<Net, Res_F<Net> > it(p); it.has_current(); it.next()) 
    if (IS_ARC_VISITED(it.get_current(), Aleph::Maximum_Flow))
      return true;

  return false;
}

    template <class Net, class QN_Type, class QA_Type> 
    typename Net::Flow_Type
generic_preflow_edge_maximum_flow(Net & net, 
                                  const bool & leave_residual = false)
{
  net.make_super_nodes();
  try
    {
      net.make_residual_net();
      net.reset_counter_nodes();
    }
  catch (std::bad_alloc)
    {
       net.unmake_super_nodes();
       throw;
    }

  QA_Type q;  // conjunto de arcos elegibles
  QN_Type p;  // conjunto de nodos activos
  typename Net::Node * source = net.get_source();
  typename Net::Node * sink   = net.get_sink();
  init_height_in_nodes(net); 

  for (Node_Arc_Iterator<Net, Res_F<Net> > i(source); i.has_curr(); i.next())
    {
      typename Net::Arc * arc  = i.get_current_arc();  
      typename Net::Node * tgt = net.get_tgt_node(arc); 
      arc->flow = tgt->in_flow = arc->cap - arc->flow;  // inundar arco
      arc->img_arc->flow       = 0;     

          // meter arcos elegibles de tgt
      for (Node_Arc_Iterator<Net, Res_F<Net> > j(tgt); j.has_curr(); j.next())
        {
          typename Net::Arc * a  = j.get_current_arc(); 
          if (node_height<Net>(tgt) == 
              node_height<Net>(net.get_tgt_node(a)) + 1)
            put_in_active_queue(q, a); 
        }
      put_in_active_queue(p, tgt); // mete en cola de nodos activos
    }
  source->out_flow = source->out_cap;

  while (true)
    {
      while (not q.is_empty()) // mientras haya arcos elegibles
        {
          typename Net::Arc * arc  = get_from_active_queue(q);
          typename Net::Node * src = net.get_src_node(arc);
          typename Net::Node * tgt = net.get_tgt_node(arc);
          if (node_height<Net>(src) != node_height<Net>(tgt) + 1) 
            continue; // No ==> ignore arco

          typename Net::Flow_Type excess = src->in_flow - src->out_flow;
          if (excess == 0) // ¿fuente descargado en iteración anterior?
            {     // sí ==> ignore arco
              remove_from_active_queue(p, src); 
              continue; // avance al siguiente arco elegible
            }
          
          const typename Net::Flow_Type push = 
            std::min(excess, arc->cap - arc->flow);
          arc->flow          += push;
          arc->img_arc->flow -= push;
          if (arc->is_residual)
            {
              net.decrease_out_flow(net.get_tgt_node(arc), push); 
              net.decrease_in_flow(net.get_src_node(arc), push);
            }
          else
            {
              net.increase_in_flow(net.get_tgt_node(arc), push);
              net.increase_out_flow(net.get_src_node(arc), push); 
            }
          excess -= push;

          if (is_node_active<Net>(tgt) and tgt != source and tgt != sink) 
            {
              for (Node_Arc_Iterator<Net, Res_F<Net> > it(tgt); 
                   it.has_curr(); it.next()) 
                {
                  typename Net::Arc * a  = it.get_current_arc();
                  if (node_height<Net>(tgt) == 
                      node_height<Net>(net.get_tgt_node(a)) + 1)
                    put_in_active_queue(q, a);
                }

              put_in_active_queue(p, tgt);
            } 

          if (excess == 0) // ¿se descargó el fuente?
            {     // sí ==> fuente deviene inactivo
              remove_from_active_queue(p, src);
              continue;
            }

          if (src != source and src != sink and   // ¿empuje no saturante?
              not has_arcs_in_active_queue<Net>(src)) 
            { // sí ==> incremente h(active) y re-inserte arc en q
              remove_from_active_queue(p, src);
              node_height<Net>(src)++;
              put_in_active_queue(p, src);

              for (Node_Arc_Iterator<Net> it(src); it.has_current(); it.next()) 
                {
                  typename Net::Arc * a = it.get_current_arc(); // sale de src
                  if ((node_height<Net>(src) == 
                       node_height<Net>(net.get_tgt_node(a)) + 1) 
                      and (a->cap - a->flow > 0)) 
                    put_in_active_queue(q, a);

                  typename Net::Arc * i = a->img_arc; // arco entrante a src
                  if ((i->cap - i->flow > 0) and // ¿i es elegible?
                      (node_height<Net>(net.get_src_node(i)) == 
                       node_height<Net>(src) + 1))
                    put_in_active_queue(q, i);
                }
            }
        }

      if (p.is_empty()) // ¿quedan nodos activos?
        break; // ya no hay nodos activos

      typename Net::Node * src = get_from_active_queue(p);
      node_height<Net>(src)++;
      put_in_active_queue(p, src);

      for (Node_Arc_Iterator<Net> it(src); it.has_current(); it.next()) 
        {
          typename Net::Arc * a = it.get_current_arc(); // sale de src
          if ((node_height<Net>(src) == 
               node_height<Net>(net.get_tgt_node(a)) + 1) 
              and (a->cap - a->flow > 0)) 
            put_in_active_queue(q, a);

          typename Net::Arc * i = a->img_arc; // arco entrante a src
          if ((i->cap - i->flow > 0) and // ¿i es elegible?
              (node_height<Net>(net.get_src_node(i)) == 
               node_height<Net>(src) + 1))
            put_in_active_queue(q, i);
        }
    } // fin  while (true)

  if (not leave_residual)
    {
      net.unmake_residual_net();
      I(source->out_flow == sink->in_flow);
      net.unmake_super_nodes();
    }

  return source->out_flow;
}
    template <class Net> typename Net::Flow_Type
depth_first_preflow_edge_maximum_flow
  (Net & net, const bool & leave_residual = false)
{
  typedef DynSetTreap<typename Net::Node*> Active;
  typedef DynListStack<typename Net::Arc*> Stack;

  return generic_preflow_edge_maximum_flow <Net,Active,Stack>
    (net, leave_residual);
}

template <class Net> class Depth_First_Preflow_Maximum_Flow
{
public:

       typename Net::Flow_Type 
   operator () (Net & net, 
                const bool & leave_residual = false) const
  {
    return depth_first_preflow_edge_maximum_flow(net, leave_residual);
  }
};

    template <class Net> typename Net::Flow_Type
breadth_first_preflow_edge_maximum_flow
   (Net & net, const bool leave_residual = false)
{
  typedef DynSetTreap<typename Net::Node*> Active;
  typedef DynListQueue<typename Net::Arc*> Queue;
  return generic_preflow_edge_maximum_flow<Net, Active, Queue>
    (net, leave_residual);
}

template <class Net> class Breadth_First_Preflow_Maximum_Flow
{
public:
       typename Net::Flow_Type 
   operator () (Net & net, 
                const bool & leave_residual = false) const
  {
    return breadth_first_preflow_edge_maximum_flow(net, leave_residual);
  }
};

template <class Net> struct Compare_Arc
{
  bool operator () (typename Net::Arc * a1, typename Net::Arc * a2) const
  {
    typename Net::Node * src1 = (typename Net::Node *) a1->src_node;
    typename Net::Node * src2 = (typename Net::Node *) a2->src_node;
    if (src1->counter == src2->counter)
      return a1->cap - a1->flow > a2->cap - a2->flow;

    return src1->counter > src2->counter;
  }
};
    template <class Net> typename Net::Flow_Type
priority_first_preflow_edge_maximum_flow
  (Net & net, const bool & leave_residual = false)
{
  typedef DynSetTreap<typename Net::Node*> Active;
  typedef DynBinHeap<typename Net::Arc*, Compare_Arc<Net> > Prio_Queue;
  return generic_preflow_edge_maximum_flow <Net, Active, Prio_Queue> 
    (net, leave_residual);
}

template <class Net> class Priority_First_Preflow_Maximum_Flow
{
public:
       typename Net::Flow_Type
   operator () (Net & net, 
                const bool & leave_residual = false) const
  {
    return priority_first_preflow_edge_maximum_flow(net, leave_residual);
  }
};

    template <class Net> typename Net::Flow_Type
random_first_preflow_edge_maximum_flow
  (Net & net, const bool & leave_residual = false)
{
  typedef DynSetTreap<typename Net::Node*> Active;
  typedef Random_Set <typename Net::Arc*> Random_Queue;
  return generic_preflow_edge_maximum_flow <Net, Active, Random_Queue> 
    (net, leave_residual);
}

template <class Net> class Random_First_Preflow_Maximum_Flow
{
public:
       typename Net::Flow_Type 
   operator () (Net & net,
                const bool & leave_residual = false) const
  {
    return random_first_preflow_edge_maximum_flow(net, leave_residual);
  }
};


static void * vs_ptr = NULL;

    template <class Net> static bool 
add_vertex_to_vs(Net &, typename Net::Node * node, typename Net::Arc*)
{
  ((Aleph::set<typename Net::Node*>*)vs_ptr)->insert(node);
  return false; // indica que debe seguir explorando
}

template <class N> class No_Res_Arc
{
public:

 bool operator () (typename N::Node *, typename N::Arc * a) const
  {
    return not a->is_residual;
  }

  bool operator () (N&, typename N::Arc * a) const
  {
    return not a->is_residual;
  }

 bool operator () (typename N::Arc * a) const
  {
    return not a->is_residual;
  }
};

  template <class Net, template <class> class Maxflow>
typename Net::Flow_Type min_cut(Net &                             net, 
                                Aleph::set<typename Net::Node*> & vs, 
                                Aleph::set<typename Net::Node*> & vt,
                                DynDlist<typename Net::Arc *> &   cuts,
                                DynDlist<typename Net::Arc *> &   cutt,
                                const bool leave_residual = false)
{
  Maxflow <Net> () (net, true); // calcula flujo máximo y red residual

  typename Net::Node * source = net.get_source();

  vs_ptr = &vs;
  depth_first_traversal <Net, Res_F<Net>> 
    (net, source, &add_vertex_to_vs);
  const size_t size_vt = net.get_num_nodes() - vs.size();
  for (Node_Iterator<Net> it(net); it.has_current() and 
         vt.size() < size_vt; it.next())
    {
      typename Net::Node * p = it.get_current();
      if (not IS_NODE_VISITED(p, Aleph::Depth_First))
        vt.insert(p);
    }     

  for (Arc_Iterator<Net, No_Res_Arc<Net> > it(net); it.has_curr(); it.next())
    {
      typename Net::Arc * arc = it.get_current();
      if (arc->flow == 0) // ¿candidato a arco de retroceso?
        if (vt.count(net.get_src_node(arc)) > 0 and // ¿de vt hacia vs?
            vs.count(net.get_tgt_node(arc)) > 0)
          {
            cutt.insert(arc);
            continue;
          }

      if (arc->flow == arc->cap) // ¿candidato a arco de cruce?
        if (vs.count(net.get_src_node(arc)) > 0 and // ¿de vs hacia vt?
            vt.count(net.get_tgt_node(arc)) > 0)
          cuts.insert(arc); 
    }

  if (not leave_residual)
    { 
      net.unmake_residual_net();
      net.unmake_super_nodes();
    }

  return source->out_flow;
}

  template <class Net, 
            template <class> class Maxflow = Heap_Preflow_Maximum_Flow>
class Min_Cut
{
public:
  typename Net::Flow_Type operator () (Net &                             net, 
                                       Aleph::set<typename Net::Node*> & vs, 
                                       Aleph::set<typename Net::Node*> & vt,
                                       DynDlist<typename Net::Arc *> &   cuts,
                                       DynDlist<typename Net::Arc *> &   cutt)
  {
    return min_cut <Net, Maxflow> (net, vs, vt, cuts, cutt);
  }
};

template <class N> class Flow_Filter
{

public:

  static typename N::Flow_Type min;

  typename N::Arc * operator () (typename N::Node *, typename N::Arc * arc) 
  { 
    return arc->cap - arc->flow >= min ? arc : NULL; 
  }

  typename N::Arc * operator () (typename N::Arc * arc) 
  { 
    return arc->cap - arc->flow >= min ? arc : NULL; 
  }
};

template <class N> typename N::Flow_Type Flow_Filter<N>::min = 0;

    template <class Net, template <class, class> class Find_Path>
bool find_aumenting_path(Net & net, Path<Net> & path, 
                         const typename Net::Flow_Type & min_slack)
{
  Flow_Filter<Net>::min     = min_slack;
  typename Net::Node * src  = net.get_source();
  typename Net::Node * sink = net.get_sink();

  if (not net.residual_net)
    net.make_residual_net();

  return Find_Path<Net, Flow_Filter<Net>> () (net, src, sink, path);
}

template <class Net, 
          template <class, class> class Find_Path = Find_Path_Breadth_First>
class Find_Aumenting_Path
{
public:
  bool operator () (Net & net, Path<Net> & path, 
                    const typename Net::Flow_Type & slack)
  {
    return find_aumenting_path <Net, Find_Path> (net, path, slack);
  }
};

    template <class Net, template <class, class> class Find_Path>
bool find_decrementing_path(Net & net, Path<Net> & path, 
                            const typename Net::Flow_Type & min_slack)
{
  Flow_Filter<Net>::min     = min_slack;
  typename Net::Node * src  = net.get_source();
  typename Net::Node * sink = net.get_sink();

  if (not net.residual_net)
    net.make_residual_net();

  return Find_Path<Net, Flow_Filter<Net> > () (net, sink, src, path);
}

template <class Net, 
          template <class, class> class Find_Path = Find_Path_Breadth_First>
class Find_Decrementing_Path
{
public:
  bool operator () (Net & net, Path<Net> & path, 
                    const typename Net::Flow_Type & slack)
  {
    return find_decrementing_path <Net, Find_Path> (net, path, slack);
  }
};

    template <class Net>
void increase_flow(Net & net, Path<Net> & path, 
                   const typename Net::Flow_Type & val)
{
  if (not net.residual_net)
    throw std::domain_error("Residual net is not created");

    // aumentar flujo de red por el camino de aumento
  for (typename Path<Net>::Iterator it(path); it.has_current_arc(); it.next())
    {
      typename Net::Arc * arc = it.get_current_arc();
      typename Net::Arc * img = (typename Net::Arc *) arc->img_arc;
      if (arc->cap - arc->flow < val)
        throw std::domain_error("minimum slack of aumenting path es"
                                " smaller than val");

      arc->flow += val; 
      img->flow -= val;
      if (not arc->is_residual)
        {
          net.increase_out_flow(net.get_src_node(arc), val);
          net.increase_in_flow(net.get_tgt_node(arc), val);
        }
      else
        {
          net.decrease_in_flow(net.get_src_node(arc), val);
          net.decrease_out_flow(net.get_tgt_node(arc), val);
        }
    }
}

    template <class Net>
void decrease_flow(Net & net, Path<Net> & path, 
                   const typename Net::Flow_Type & val)
{
  if (not net.residual_net)
    throw std::domain_error("Residual net is not created");

    // decrementar flujo de red por camino de decremento
  for (typename Path<Net>::Iterator it(path); it.has_current_arc(); it.next())
    {
      typename Net::Arc * arc = it.get_current_arc();
      typename Net::Arc * img = (typename Net::Arc *) arc->img_arc;

      if (arc->cap - arc->flow < val)
        throw std::domain_error("minimum slack of decrementing path "
                                " is smaller than val");

      arc->flow -= val; 
      img->flow += val;
      if (not arc->is_residual)
        {
          net.decrease_out_flow(get_src_node(arc), val);
          net.decrease_in_flow(net.get_tgt_node(arc), val);
        }
      else
        {
          net.increase_out_flow(get_tgt_node(arc), val);
          net.increase_in_flow(net.get_src_node(arc), val);
        }
    }
}

} // end namespace Aleph

# endif // TPL_NETGRAPH_H