Is this Sedgewick code correct?

Question

I am solving an optimization problem in which, among other things, I must maximize flow networks. I implemented a c++ code based flow-maximization algorithm based in the following java code that appears in the book of Sedgewick "Algorithms in Java, Third Edition, Part 5: Graph Algorithms", which maximizes a network flow using a vertex-based PREFLOW-push algorithm:

class NetworkMaxFlow
{ private Network G; private int s, t;
  private int[] h, wt;
  private void initheights()
  NetworkMaxFlow(Network G, int s, int t)
  { this.G = G; this.s = s; this.t = t;
    wt = new int[G.V()]; h = new int[G.V()];
    initheights();
    intGQ gQ = new intGQ(G.V());
    gQ.put(s); wt[t] = -(wt[s] = Edge.M*G.V());
    while (!gQ.empty())
    { int v = gQ.get();
      AdjList A = G.getAdjList(v);
      for (Edge e = A.beg(); !A.end(); e = A.nxt())
        { int w = e.other(v), cap = e.capRto(w);
          int P = cap < wt[v] ? cap : wt[v];
          if (P > 0 && v == s || h[v] == h[w]+1) // first observation (see below)
            { e.addflowRto(w, P);
              wt[v] -= P; wt[w] += P;
              if ((w != s) && (w != t)) gQ.put(w); // enqueue w if it is not source or sink
            }
        }
      if (v != s && v != t && wt[v] > 0) // why check v != t if t never enter the queue?
        { h[v]++; gQ.put(v); }
    }
  }
}

My implementation, based on that code, fails to maximize the following network After execution, the the resulting flow is as follows With this flow, the flow value is 8, but the maximum is 9, as indicated by the flow of the figure below According to my understanding, the algorithm is consistent with the explanation of the book. However, I see two strange things

1. There is no explicit preflow phase from source. It is included in the while and executed first and only once when the predicate P > 0 && v == s is true. Maybe this was done to shorten the code
2. According to my understanding, and discourse of the book, the sink never enters the queue. However, when the height is increased, the code checks that v != t. Any reason for this?

This is an excerpt from my implementation of this algorithm in C++

template <class Net, class Q_Type> typename Net::Flow_Type
generic_preflow_vertex_push_maximum_flow(Net & net)
{
  init_height_in_nodes(net); // breadth first traverse from sink to
                 // source. Nodes are labeled with their
                 // minimal distance (in nodes) to sink
  auto source = net.get_source();
  auto sink   = net.get_sink();

  using Itor = __Net_Iterator<Net>;
  Q_Type q; // generic queue (can be fifo, heap or random) of active nodes

  // preflow: floods all nodes connected to the source 
  for (Itor it(source); it.has_curr(); it.next()) 
    {
      auto arc  = it.get_curr();  
      arc->flow = arc->cap; // saturate arc to its maximum 
      auto tgt = net.get_tgt_node(arc);
      put_in_active_queue(q, tgt);
      assert(node_height<Net>(source) == node_height<Net>(tgt) + 1);
      assert(not is_residual<Net>(source, arc));
    }

  while (not q.is_empty()) // while there are active nodes
    {
      auto src = get_from_active_queue(q);
      auto excess = net.get_in_flow(src) - net.get_out_flow(src);

      for (Itor it(src); it.has_curr(); it.next()) 
        {
          auto arc = it.get_curr();
          auto tgt = net.get_connected_node(arc, src);

          if (node_height<Net>(src) != node_height<Net>(tgt) + 1)
            continue; // this arc is not eligible

          typename Net::Flow_Type flow_to_push;
          if (is_residual<Net>(src, arc))
            {
              flow_to_push = std::min(arc->flow, excess);
              arc->flow -= flow_to_push;
            }
          else
            {
              flow_to_push = std::min(arc->cap - arc->flow, excess);
              arc->flow += flow_to_push;
            }

          excess -= flow_to_push;
          if (tgt != sink and tgt != source)
            put_in_active_queue(q, tgt);
        }

    if (excess > 0) // src still active?
      { 
        node_height<Net>(src)++;
        put_in_active_queue(q, src);
      }
  }

  return net.flow_value(); // sum of all outing flow from source
}

¿Someone find any logical inconsistency between my code and the code of Sedgewick? I have the impression that my code (and perhaps also the Sedgewick) is not properly handling the increases in heights. But I do not manage to understand why

I show a detailed execution trace with the network that fails to maximize (the trace start from the first q.get() of while. The values in parentheses are the values of the heights. IN is the incoming flow to the node. OUT the outcoming one.

As example, the line

    4104 (2) --> 0 (1) pushing 1 from 4104 toward 0

refers the eligible arc 4104-->0. The node 4104 has height 2 and the node 0 has height 1. The expression "pushing 1" means that 1 unity of flow is pushed toward the target node (0). The line ================ separates each queue extraction. The queue is FIFO and its state is printed at the end of each processing.

Note that many times zero flow units are pushed or reduceed, but the destination node becomes active.

This is the execution trace

Initial Queue = 4104 4105 4106 4107 4108

Active node 4104 Height = 2 IN = 1 OUT = 0
    4104 (2) --> source (3) not eligible
    4104 (2) --> 0 (1) pushing 1 from 4104 toward 0
    4104 (2) --> 1 (1) pushing 0 from 4104 toward 1
    4104 (2) --> 2 (1) pushing 0 from 4104 toward 2
    4104 (2) --> 4 (1) pushing 0 from 4104 toward 4
    Excess = 0
    Queue = 4105 4106 4107 4108 0 1 2 4
================
Active node 4105 Height = 2 IN = 3 OUT = 0
    4105 (2) --> source (3) not eligible
    4105 (2) --> 1 (1) pushing 1 from 4105 toward 1
    4105 (2) --> 4 (1) pushing 1 from 4105 toward 4
    4105 (2) --> 6 (1) pushing 1 from 4105 toward 6
    Excess = 0
    Queue = 4106 4107 4108 0 1 2 4 6
================
Active node 4106 Height = 2 IN = 1 OUT = 0
    4106 (2) --> source (3) not eligible
    4106 (2) --> 1 (1) pushing 1 from 4106 toward 1
    4106 (2) --> 5 (1) pushing 0 from 4106 toward 5
    Excess = 0
    Queue = 4107 4108 0 1 2 4 6 5
================
Active node 4107 Height = 2 IN = 1 OUT = 0
    4107 (2) --> source (3) not eligible
    4107 (2) --> 1 (1) pushing 1 from 4107 toward 1
    4107 (2) --> 2 (1) pushing 0 from 4107 toward 2
    4107 (2) --> 3 (1) pushing 0 from 4107 toward 3
    4107 (2) --> 4 (1) pushing 0 from 4107 toward 4
    4107 (2) --> 6 (1) pushing 0 from 4107 toward 6
    Excess = 0
    Queue = 4108 0 1 2 4 6 5 3
================
Active node 4108 Height = 2 IN = 3 OUT = 0
    4108 (2) --> source (3) not eligible
    4108 (2) --> 1 (1) pushing 1 from 4108 toward 1
    4108 (2) --> 2 (1) pushing 1 from 4108 toward 2
    4108 (2) --> 4 (1) pushing 1 from 4108 toward 4
    4108 (2) --> 5 (1) pushing 0 from 4108 toward 5
    4108 (2) --> 6 (1) pushing 0 from 4108 toward 6
    Excess = 0
    Queue = 0 1 2 4 6 5 3
================
Active node 0 Height = 1 IN = 1 OUT = 0
    0 (1) --> sink (0) pushing 1 from 0 toward sink
    0 (1) --> 4104 (2) not eligible
    Excess = 0
    Queue = 1 2 4 6 5 3
================
Active node 1 Height = 1 IN = 4 OUT = 0
    1 (1) --> sink (0) pushing 2 from 1 toward sink
    1 (1) --> 4105 (2) not eligible
    1 (1) --> 4106 (2) not eligible
    1 (1) --> 4107 (2) not eligible
    1 (1) --> 4108 (2) not eligible
    Excess = 2    1 goes back onto queue with label 2
    Queue = 2 4 6 5 3 1
================
Active node 2 Height = 1 IN = 1 OUT = 0
    2 (1) --> sink (0) pushing 1 from 2 toward sink
    2 (1) --> 4108 (2) not eligible
    Excess = 0
    Queue = 4 6 5 3 1
================
Active node 4 Height = 1 IN = 2 OUT = 0
    4 (1) --> sink (0) pushing 2 from 4 toward sink
    4 (1) --> 4105 (2) not eligible
    4 (1) --> 4108 (2) not eligible
    Excess = 0
    Queue = 6 5 3 1
================
Active node 6 Height = 1 IN = 1 OUT = 0
    6 (1) --> sink (0) pushing 1 from 6 toward sink
    6 (1) --> 4105 (2) not eligible
    Excess = 0
    Queue = 5 3 1
================
Active node 5 Height = 1 IN = 0 OUT = 0
    5 (1) --> sink (0) pushing 0 from 5 toward sink
    Excess = 0
    Queue = 3 1
================
Active node 3 Height = 1 IN = 0 OUT = 0
    3 (1) --> sink (0) pushing 0 from 3 toward sink
    Excess = 0
    Queue = 1
================
Active node 1 Height = 2 IN = 4 OUT = 2
    1 (2) --> 4105 (2) not eligible
    1 (2) --> 4106 (2) not eligible
    1 (2) --> 4107 (2) not eligible
    1 (2) --> 4108 (2) not eligible
    Excess = 2    1 goes back onto queue with label 3
    Queue = 1
================
Active node 1 Height = 3 IN = 4 OUT = 2
    1 (3) --> 4105 (2) Reducing 1 from 1 toward 4105
    1 (3) --> 4106 (2) Reducing 1 from 1 toward 4106
    1 (3) --> 4107 (2) Reducing 0 from 1 toward 4107
    1 (3) --> 4108 (2) Reducing 0 from 1 toward 4108
    Excess = 0
    Queue = 4105 4106 4107 4108
================
Active node 4105 Height = 2 IN = 3 OUT = 2
    4105 (2) --> source (3) not eligible
    4105 (2) --> 1 (3) not eligible
    Excess = 1    4105 goes back onto queue with label 3
    Queue = 4106 4107 4108 4105
================
Active node 4106 Height = 2 IN = 1 OUT = 0
    4106 (2) --> source (3) not eligible
    4106 (2) --> 1 (3) not eligible
    4106 (2) --> 5 (1) pushing 1 from 4106 toward 5
    Excess = 0
    Queue = 4107 4108 4105 5
================
Active node 4107 Height = 2 IN = 1 OUT = 1
    4107 (2) --> source (3) not eligible
    4107 (2) --> 2 (1) pushing 0 from 4107 toward 2
    4107 (2) --> 3 (1) pushing 0 from 4107 toward 3
    4107 (2) --> 4 (1) pushing 0 from 4107 toward 4
    4107 (2) --> 6 (1) pushing 0 from 4107 toward 6
    Excess = 0
    Queue = 4108 4105 5 2 3 4 6
================
Active node 4108 Height = 2 IN = 3 OUT = 3
    4108 (2) --> source (3) not eligible
    4108 (2) --> 5 (1) pushing 0 from 4108 toward 5
    4108 (2) --> 6 (1) pushing 0 from 4108 toward 6
    Excess = 0
    Queue = 4105 5 2 3 4 6
================
Active node 4105 Height = 3 IN = 3 OUT = 2
    4105 (3) --> source (3) not eligible
    4105 (3) --> 1 (3) not eligible
    Excess = 1    4105 goes back onto queue with label 4
    Queue = 5 2 3 4 6 4105
================
Active node 5 Height = 1 IN = 1 OUT = 0
    5 (1) --> sink (0) pushing 1 from 5 toward sink
    5 (1) --> 4106 (2) not eligible
    Excess = 0
    Queue = 2 3 4 6 4105
================
Active node 2 Height = 1 IN = 1 OUT = 1
    2 (1) --> sink (0) pushing 0 from 2 toward sink
    2 (1) --> 4108 (2) not eligible
    Excess = 0
    Queue = 3 4 6 4105
================
Active node 3 Height = 1 IN = 0 OUT = 0
    3 (1) --> sink (0) pushing 0 from 3 toward sink
    Excess = 0
    Queue = 4 6 4105
================
Active node 4 Height = 1 IN = 2 OUT = 2
    4 (1) --> 4105 (4) not eligible
    4 (1) --> 4108 (2) not eligible
    Excess = 0
    Queue = 6 4105
================
Active node 6 Height = 1 IN = 1 OUT = 1
    6 (1) --> sink (0) pushing 0 from 6 toward sink
    6 (1) --> 4105 (4) not eligible
    Excess = 0
    Queue = 4105
================
Active node 4105 Height = 4 IN = 3 OUT = 2
    4105 (4) --> source (3) Reducing 1 from 4105 toward source
    4105 (4) --> 1 (3) pushing 0 from 4105 toward 1
    Excess = 0
    Queue = 1
================
Active node 1 Height = 3 IN = 2 OUT = 2
    1 (3) --> 4107 (2) Reducing 0 from 1 toward 4107
    1 (3) --> 4108 (2) Reducing 0 from 1 toward 4108
    Excess = 0
    Queue = 4107 4108
================
Active node 4107 Height = 2 IN = 1 OUT = 1
    4107 (2) --> source (3) not eligible
    4107 (2) --> 2 (1) pushing 0 from 4107 toward 2
    4107 (2) --> 3 (1) pushing 0 from 4107 toward 3
    4107 (2) --> 4 (1) pushing 0 from 4107 toward 4
    4107 (2) --> 6 (1) pushing 0 from 4107 toward 6
    Excess = 0
    Queue = 4108 2 3 4 6
================
Active node 4108 Height = 2 IN = 3 OUT = 3
    4108 (2) --> source (3) not eligible
    4108 (2) --> 5 (1) pushing 0 from 4108 toward 5
    4108 (2) --> 6 (1) pushing 0 from 4108 toward 6
    Excess = 0
    Queue = 2 3 4 6 5
================
Active node 2 Height = 1 IN = 1 OUT = 1
    2 (1) --> sink (0) pushing 0 from 2 toward sink
    2 (1) --> 4108 (2) not eligible
    Excess = 0
    Queue = 3 4 6 5
================
Active node 3 Height = 1 IN = 0 OUT = 0
    3 (1) --> sink (0) pushing 0 from 3 toward sink
    Excess = 0
    Queue = 4 6 5
================
Active node 4 Height = 1 IN = 2 OUT = 2
    4 (1) --> 4105 (4) not eligible
    4 (1) --> 4108 (2) not eligible
    Excess = 0
    Queue = 6 5
================
Active node 6 Height = 1 IN = 1 OUT = 1
    6 (1) --> sink (0) pushing 0 from 6 toward sink
    6 (1) --> 4105 (4) not eligible
    Excess = 0
    Queue = 5
================
Active node 5 Height = 1 IN = 1 OUT = 1
    5 (1) --> sink (0) pushing 0 from 5 toward sink
    5 (1) --> 4106 (2) not eligible
    Excess = 0
    Queue =

Albumin

You may want to change the title. I don't think it reflects the quality of the post - typically "is my code correct?" titles have very poor questions – Plaque

Thanks for your suggestion. I'll think about a new title. Anyway if you help me with a suggestion I would thank you even more – Albumin

This is a better stack exchange website for posts like this.codereview.stackexchange.com – Gayelord

i think you should find a smaller counter example. Just write a generator and a simpler max flow implementation. As your problem is with the value of the flow and not with capacities any max-flow algorithm will work. – Galenism

Answer 1

Your questions

Preflow

There is no explicit preflow phase from source. It is included in the while and executed first and only once when the predicate P > 0 && v == s is true. Maybe this was done to shorten the code

Yes, I think that was done to compact the code. Due to the assignment wt[s] = Edge.M*G.V() at the beginning the source has enough "virtual" excess in order to fuel the preflow in the first iteration of the algorithm. If you wanted to do the same trick in your implementation, you would have to inflate the excess variable (that you calculate on the fly instead of storing it in an array like Sedgewick) during the first iteration, when you encounter the source node. But you don't have to do it that way, your explicit preflooding seems just fine and probably even makes things more readable.

I also suspect that in the condition P > 0 && v == s || h[v] == h[w]+1 the implicit operator precedence (P > 0 && v == s) || h[v] == h[w]+1 was not intended. The way it is written, the check P > 0 is only performed for the source case. Not checking it in the other cases won't hurt because executing the body with P == 0 will just push a 0 excess to other nodes (which doesn't do anything) and unnecessarily add those other nodes to the queue, just for them being removed again immediately. I've seen that you implemented it the same way. That's OK, even though a slight waste of computational time. Probably P > 0 && (v == s || h[v] == h[w]+1) was really intended.

Sink in the queue

According to my understanding, and discourse of the book, the sink never enters the queue. However, when the height is increased, the code checks that v != t. Any reason for this?

I agree, this is weird. The only way the sink can enter the queue is by being the source at the same time (in a graph with only one node and no edges). However, in that case the condition v != s already avoids the infinite loop, no extra condition necessary.

Wrong results

After execution, the resulting flow is as follows [...] With this flow, the flow value is 8, but the maximum is 9

Initial heights

You are getting wrong results because your initial heights are wrong. I cannot compare Sedgewick's initialization code with yours because your post doesn't specify either. But I learn from your log file that you are starting with a height of 3 for source. That is clearly against the conditions for height functions: The source has to start out with a height equal to the number of nodes in your graph and will keep it during the whole algorithm.

In your case the height of source is too close to the height of its downstream neighbors. This causes the downstream neighbors to push some flow back to the source quite soon, before trying hard enough to let it flow further downstream. They are supposed to do that only if there is no way to let it flow down to the sink.

Broken graph?

What further concerns me is that the edge [4104 -> 1] seems to disappear. While it is mentioned during the processing of the node 4104, it never shows up during the processing of the node 1. All other incoming edges are mentioned, be it with regard to "not eligible", "pushing" or "reducing". I'd expect it to appear as one of the three types, just like the others. But then again, I don't know where exactly you put those logging statements. Still, it leaves a bit of a concern and I thought I'd mention it in case you still have trouble after fixing your heights.

Hassle

thanks a lot for having taken the time to analyze answer in detail my question. You're right. The error is in the way of initializing the heights. These are initialized with the minimal distance in arcs from the sink. This way is suggested by Sedgewick in the function initheights() as well as others authors. However, the exception is for the source, which never must be active. So, your correction, to init source in |V| is fine. Now the implementation computes correctly the maximum flow. – Albumin

I do not think Sedgewick was wrong, but I think his text could mislead. He puts as exercise to implement initheights() using breadth-first search from the sink. His code does not explicitly initialize the source and in his discourse he deals with the source sometimes, I would say, not clearly. For example, in his corollary about the height values he says "The height of the source never changes, and it is initially no greater than V". That is correct but I think he must have been more explicit. Arjuna et al, for example, in their algorithm very explicitly initialize the source height in |V| – Albumin

Hopefully, this question remains as a lesson. Again, thank you for your answer – Albumin

I just got my hands on his text and I do think he is wrong: He proves the correctness of the algorithm and uses just the other properties of h (without using h(s)=|V|). And that's obviously not enough, otherwise your implementation would work. So his proof has to have a flaw somewhere. I see the flaw in his proof of Property 22.11, where he doesn't distinguish properly between "edges in the residual network" and "eligible edges in the res. network". He says that h(u) was increased because there were no eligible edges, but then he uses that fact without the eligible limitation. – Hassle