MODULE MaxFlowProbs EXPORTS Main;

(* Generates test cases for unit-cost directed max-flow problems *)
(* Created April 1997 by J. Stolfi  *)
(* See the copyright notice at the end of this file. *)

IMPORT Wr, Rd, Thread, Fmt, ASCII, Math, XRandom;
IMPORT ParseParams, Process, Text;
IMPORT DirGraph;
FROM DirGraph IMPORT
  VertexNum, VertexNums, EdgeNum, EdgeNums,
  Edge, Edges, Cost, Costs, Mark, Marks, Count, Counts;
FROM XRandom IMPORT Seed;
FROM Stdio IMPORT stdout, stderr;

<* FATAL Thread.Alerted, Wr.Failure *>

CONST
  MaxVertices = 32*256*256;
  MaxEdges = 16*MaxVertices;
  Version = "1.1";

TYPE
  LONG = LONGREAL;
  BOOL = BOOLEAN;
  
  Kind = {Rete};
  
  Graph = DirGraph.T;
  
  SubstSpec = RECORD
      fraction: LONG;     (* Fraction of marked edges to replace. *)
      marks: SET OF Mark; (* Replace all edges with these marks. *)
      kind: Kind;         (* Block type. *)
      NV: Count;          (* Vertex count. *)
      NE: Count;          (* Edge count. *)
      W: Count;           (* Max flow. *)
      L: Count;           (* Mean flow length. *)
      scaleUppers: BOOL;  (* TRUE if unreplaced upercase edges must be scaled by "W". *)
      comment: TEXT;      (* Comments specific to this graph. *)
    END;
    
  SubstSpecs = ARRAY OF SubstSpec;
  
  Options = RECORD
      spec: REF SubstSpecs;   (* Sequence of graphs to substitute. *)
      omit: SET OF Mark;      (* Edges to omit from the output. *)
      seed: XRandom.Seed;     (* Random number generator seed. *)
      dontScramble: BOOLEAN;  (* TRUE prevents the random renumbering of vertices. *)
      reverse: BOOLEAN;       (* TRUE reverses the default vertex numbering. *)
      comment: TEXT;          (* General comment text. *)
    END;
    
PROCEDURE DoIt() =
  VAR G: Graph;
  BEGIN
    WITH
      o = GetOptions(),
      rnd = NEW(XRandom.T).init(seed := o.seed)
    DO
      G := MakeUnitFlowGraph(o, rnd);
      IF o.dontScramble THEN
        Wr.PutText(stderr, "Preserving raw vertex numbering...\n")
      ELSIF o.reverse THEN
        Wr.PutText(stderr, "Reversing raw vertex numbering...\n");
        DirGraph.ComplementVertexNums(G)
      ELSE
        Wr.PutText(stderr, "Scrambling vertex numbering...\n");
        DirGraph.PermuteVertexNums(G, rnd)
      END;
      Wr.PutText(stderr, "Writing graph...\n");
      DirGraph.Write(stdout, G, vBase := 1);
      Wr.PutText(stderr, "Done.\n");
    END;
  END DoIt;
  
PROCEDURE GetOptions (): Options =
  <* FATAL Thread.Alerted, Wr.Failure *>
  VAR o: Options;
      rs: REF SubstSpecs := NEW(REF SubstSpecs, MaxSubsts);
      nSubsts: CARDINAL := 0;
  CONST
    MaxSubsts = 20;
  BEGIN
     WITH pp = NEW(ParseParams.T).init(stderr) DO
      TRY
        o.comment := "Generated by stolfi/MaxFlowProbs-" & Version;
        
        pp.getKeyword("-seed");
        o.seed := pp.getNextInt(FIRST(Seed), LAST(Seed));
        o.comment := o.comment & "\n" & "  -seed " & Fmt.Int(o.seed);

        IF  pp.keywordPresent("-dontScramble") THEN
          o.dontScramble := TRUE; o.reverse := FALSE;
          o.comment := o.comment & "\n" & "  -dontScramble"
        ELSIF pp.keywordPresent("-reverse") THEN
          o.dontScramble := FALSE; o.reverse := TRUE;
          o.comment := o.comment & "\n" & "  -reverse"
        ELSE
          o.reverse := FALSE; o.dontScramble := FALSE
        END;
        
        IF pp.keywordPresent("-writeOnly") THEN
          o.omit := ASCII.Letters;
          WITH letters = pp.getNext() DO 
            FOR i := 0 TO Text.Length(letters) - 1 DO
              WITH c = Text.GetChar(letters, i) DO
                IF NOT c IN ASCII.Letters THEN 
                  pp.error("invalid mark \"" & Text.FromChar(c) & "\"")
                END;
                IF NOT (c IN o.omit) THEN
                  pp.error("repeated mark \"" & Text.FromChar(c) & "\"")
                END;
                o.omit := o.omit - SET OF Mark{c}
              END
            END;
            IF o.omit = ASCII.Letters THEN
              pp.error("\"-writeOnly\" with empty set of marks")
            END;
            o.comment := o.comment & " -writeOnly " & letters
          END;
        ELSE
          o.omit := DirGraph.NoMark
        END;
        
        TRY
          FOR i := 0 TO MaxSubsts-1 DO
            WITH spec = ParseSubstSpecs(pp) DO
              rs[nSubsts] := spec; INC(nSubsts)
            END
          END
        EXCEPT Rd.EndOfFile => (*OK*)
        END;
        o.spec := NEW(REF SubstSpecs, nSubsts);
        o.spec^ := SUBARRAY(rs^, 0, nSubsts);

        pp.finish();                                       
      EXCEPT                                                            
      | ParseParams.Error =>                                              
          Wr.PutText(stderr, "Usage: MaxFlowProbs \\\n");
          Wr.PutText(stderr, "  -seed NUM \\\n");
          Wr.PutText(stderr, "  [ -dontScramble | -reverse ] \\\n");
          Wr.PutText(stderr, "  [ -writeOnly MARKS ] \\\n");
          Wr.PutText(stderr, "  SUBSTSPECS... \n");
          Wr.PutText(stderr, "\n");
          Wr.PutText(stderr, "where SUBSTSPECS is\n");  
          Wr.PutText(stderr, "  -subst MARKS \\n");
          Wr.PutText(stderr, "    -by { Rete | ...} \\\n");
          Wr.PutText(stderr, "    -NV NUM  -NE NUM  -W WIDTH  -L LENGTH \n");
          Process.Exit (1);
      END;
    END;
    RETURN o
  END GetOptions;
  
PROCEDURE ParseSubstSpecs(pp: ParseParams.T): SubstSpec 
  RAISES {Rd.EndOfFile, ParseParams.Error} =
  (*
    Parses the next "-spec SPECS" group of options from the command line.
    Raises "Rd.EndOfFile" if none. *)

  PROCEDURE MatchKeyword(key: TEXT) RAISES {ParseParams.Error} =
    BEGIN
      IF NOT pp.testNext(key) THEN 
        pp.error("expecting \"" & key & "\" keyword")
      END
    END MatchKeyword;

  VAR spec: SubstSpec;
  BEGIN
    IF NOT pp.keywordPresent("-subst") THEN RAISE Rd.EndOfFile END;
    spec.comment := "  -subst ";
    
    spec.fraction := pp.getNextLongReal(0.0d0, 1.0d0);
    spec.comment := spec.comment & " " &
      Fmt.LongReal(spec.fraction, style := Fmt.Style.Fix, prec := 6);
    
    WITH letters = pp.getNext() DO 
      spec.marks := DirGraph.NoMark;
      spec.scaleUppers := FALSE;
      FOR i := 0 TO Text.Length(letters) - 1 DO
        WITH c = Text.GetChar(letters, i) DO
          IF NOT c IN ASCII.Letters THEN 
            pp.error("invalid mark \"" & Text.FromChar(c) & "\"")
          END;
          IF c IN spec.marks THEN
            pp.error("repeated mark \"" & Text.FromChar(c) & "\"")
          END;
          spec.marks := spec.marks + SET OF Mark{c};
          IF c IN ASCII.Uppers THEN 
            spec.scaleUppers := TRUE
          END
        END
      END;
      IF spec.marks = DirGraph.NoMark THEN
        pp.error("\"-subst\" applied to empty set of marks")
      END;
      spec.comment := spec.comment & " " & letters
    END;
    
    MatchKeyword("-by");
    spec.comment := spec.comment & " -by ";
    WITH kt = pp.getNext() DO
      IF Text.Equal(kt, "Rete") THEN
        spec.kind := Kind.Rete
      ELSE
        pp.error("bad graph family name \"" & kt & "\"");
      END;
      spec.comment :=  spec.comment & kt
    END;

    MatchKeyword("-NV");
    spec.NV := pp.getNextInt(1, MaxVertices);
    spec.comment := spec.comment & "\n" & "  -NV " & Fmt.Int(spec.NV);

    MatchKeyword("-NE");
    WITH
      maxE = MAX(1, ROUND(MIN(FLOAT(MaxEdges, LONG), Choose2(spec.NV) - 1.0d0)))
    DO
      spec.NE := pp.getNextInt(1, maxE)
    END;
    spec.comment := spec.comment & " -NE " & Fmt.Int(spec.NE);

    MatchKeyword("-W");
    spec.W := pp.getNextInt(1, MAX(1, spec.NE DIV 2));
    spec.comment := spec.comment & "\n" & "  -W " & Fmt.Int(spec.W);

    MatchKeyword("-L");
    spec.L := pp.getNextInt(1, spec.NE DIV spec.W);
    spec.comment := spec.comment & " -L " & Fmt.Int(spec.L);
    
    IF spec.scaleUppers AND spec.fraction < 1.0d0 AND spec.W # 1 THEN
      pp.error("partial replacement of uppercase edges with non-unit \"-W\"");
    END;

    RETURN spec
  END ParseSubstSpecs;

PROCEDURE MakeUnitFlowGraph(READONLY o: Options; rnd: XRandom.T): Graph =
  VAR G: Graph;
  BEGIN
    Wr.PutText(stderr, "Generating trivial graph...\n");
    G := MakeTrivialGraph(o.comment);
    FOR ig := 0 TO LAST(o.spec^) DO
      Wr.PutText(stderr, "Expanding graph...\n");
      G := ExpandGraph(G, o.spec[ig], rnd)
    END;
    RETURN G
  END MakeUnitFlowGraph;
  
PROCEDURE MakeTrivialGraph(cmt: TEXT): Graph =
  BEGIN
    WITH
      eRef = NEW(REF Edges, 1), e = eRef^,
      mRef = NEW(REF Marks, 1), m = mRef^,
      cRef = NEW(REF Costs, 1), c = cRef^
    DO
      e[0] := Edge{0, 1};
      m[0] := 'T';
      c[0] := 1;
      RETURN Graph{
        NV := 2,
        NE := 1,
        s := 0,
        t := 1,
        e := eRef,
        m := mRef,
        c := cRef,
        class := "unit-flow",
        cmt := cmt & "\n--- trivial unit flow graph ---"
      }
    END
  END MakeTrivialGraph;
  
PROCEDURE ExpandGraph(
    READONLY G: Graph; 
    READONLY spec: SubstSpec;
    rnd: XRandom.T;
  ): Graph =
  (*
    Returns a graph "X" that is the result of replacing every selected
    edge of "G" (the `skeleton' graph) by an instance of the graph
    defined by "spec" (a `block').  The block's source and sink are
    identified with the edge's origin and destination; the other block
    vertices are new.  Unselected edges are preserved.
    
    The expansion preserves the relative ordering of the skeleton
    vertices.  Moreover, all the new vertices of each substituted block
    will be numbered consecutively, with numbers intermediate between
    its source and sink.
    
    Every edge of the original graph will become a set of
    consecutive edges in the expanded graph.  The ordering of these
    edge sets in the edge array is the same as the original edges.
    
    The block-building procedures generally require that the source be
    smaller than the sink.  Therefore, if a marked edge of "G" goes
    from a high vertex to a low vertex, we reverse its direction,
    replace it by a block, and then complement all vertices in the 
    latter.
  *)

  BEGIN (* ExpandGraph *)
    WITH
      (* Parameters of "G": *)
      NEG = NUMBER(G.e^),
      (* Correspondence between "G" and expanded graph: *)
      tb = ComputeExpTables(G, spec, rnd),
      (* Parameters of expanded graph "X": *)
      NVX = tb.NV,
      NEX = tb.NE,
      (* Elements of expanded graph "X": *)
      eRef = NEW(REF Edges, NEX), e = eRef^,
      mRef = NEW(REF Marks, NEX), m = mRef^,
      cRef = NEW(REF Costs, NEX), c = cRef^
    DO
      FOR ie := 0 TO NEG-1 DO
        WITH
          (* Endpoints, mark, and cost of edge "ie" in "G": *)
          uvOld = G.e[ie], 
          uOld = uvOld[0],
          vOld = uvOld[1],
          (* Endpoints in "X" of block substituted for edge "ie": *)
          u = tb.vMap[uOld],
          v = tb.vMap[vOld],
          (* Index range of "X" edges and new "X" vertices *)
          (* corresponding to edge "ie" of "G": *)
          eLo = tb.eLo[ie], eHi = tb.eHi[ie],
          vLo = tb.vLo[ie], vHi = tb.vHi[ie],
          (* Subset of edges, marks, and costs of substituted block: *)
          es = SUBARRAY(e, eLo, eHi - eLo),
          ms = SUBARRAY(m, eLo, eHi - eLo),
          cs = SUBARRAY(c, eLo, eHi - eLo)
        DO
          IF tb.exp[ie] THEN
            <* ASSERT eHi = eLo + spec.NE *>
            <* ASSERT vHi = vLo + spec.NV - 2 *>
            ExpandEdge(G.m[ie], G.c[ie], spec, rnd, u, v, vLo, vHi, es, ms, cs);
          ELSE 
            <* ASSERT eHi = eLo + 1 *>
            <* ASSERT vHi <= vLo *>
            CopyEdge(G.m[ie], G.c[ie], spec, u, v, es[0], ms[0], cs[0])
          END;
        END;
      END;
      RETURN Graph{
        NV := NVX,
        NE := NEX,
        e := eRef,
        s := tb.vMap[G.s],
        t := tb.vMap[G.t],
        c := cRef,
        m := mRef,
        cmt := G.cmt & "\n--- substituted with ---\n" & spec.comment,
        class := "unit-flow"
      }
    END
  END ExpandGraph;
  
PROCEDURE ExpandEdge(
    mOld: Mark;               (* Mark of original edge. *)
    cOld: Cost;               (* Cost of original edge. *)
    READONLY spec: SubstSpec; (* Substitution parameters. *)
    rnd: XRandom.T;           (* Randomness source. *)
    u, v: VertexNum;          (* New source and sink of block. *)
    vLo, vHi: VertexNum;      (* Range of new internal vertices. *)
    (*OUT*)
    VAR e: Edges;             (* New edges. *)
    VAR m: Marks;             (* New edge marks. *)
    VAR c: Costs;             (* New edge costs. *)
  ) =
  BEGIN 
    AddSubGraph(
      spec.kind,
      W := spec.W, L := spec.L, 
      s := MIN(u,v), t := MAX(u,v),
      vLo := vLo, vHi := vHi,
      rnd := rnd,
      (*OUT*)
      e := e, m := m, c := c
    );
    (* Reverse direction of edges if the orignal edge was decreasing: *)
    IF u > v THEN 
      DirGraph.ReverseEdges(e)
    END;
    (* Scale edge costs of block by original edge cost: *)
    FOR je := 0 TO LAST(e) DO 
      c[je] := c[je] * cOld
    END;
    (* If original edge was not part of the maximum flow, this entire block isn't: *)
    IF mOld IN ASCII.Lowers THEN
      FOR je := 0 TO LAST(e) DO 
        m[je] := ASCII.Lower[m[je]]
      END;
    END
  END ExpandEdge;
  
PROCEDURE CopyEdge(
    mOld: Mark;                (* Mark of original edge. *)
    cOld: Cost;                (* Cost of original edge. *)
    READONLY spec: SubstSpec;  (* Substitution parameters. *)
    u, v: VertexNum;           (* Endpoints of new edge. *)
    (*OUT*)
    VAR e: Edge;               (* New edge. *)
    VAR m: Mark;               (* New edge mark. *)
    VAR c: Cost;               (* New edge cost. *)
  ) =
  BEGIN
    (* Just copy the edge "ie" *)
    e := Edge{u,v}; 
    m := mOld; 
    IF mOld IN ASCII.Uppers AND spec.scaleUppers THEN
      (* Edge was part of maximum flow, and the latter got scaled by "spec.W": *)
      c := cOld * spec.W
    ELSE
      c := cOld
    END
  END CopyEdge; 

TYPE 
  ExpTables = RECORD (* Mapping between "G" and "X" *)
      NV: Count;                (* Number of vertices in expanded graph "X" *)
      NE: Count;                (* Number of edges in expanded graph "X" *)
      vMap: REF VertexNums;     (* Maps "G" vertex nums to "X" vertex nums. *)
      exp: REF ARRAY OF BOOL;   (* "exp[i] = TRUE" means expand edge "e[i]". *)
      eLo, eHi: REF EdgeNums;   (* Map "G" edge nums to "X" range of edge nums. *)
      vLo, vHi: REF VertexNums; (* Map "G" edge nums to "X" range of vertex nums. *)
    END;

PROCEDURE ComputeExpTables(
    READONLY G: Graph; 
    READONLY spec: SubstSpec;
    rnd: XRandom.T;
  ): ExpTables =
  (*
    Returns a set of tables "tb" that describes the correspondence
    between elements of the skeleton graph  "G" and those of 
    the expanded graph "X".

    Specifically, 

      * Vertex number "iv" of "G" corresponds to vertex number
      "tb.vMap[iv]" of the expanded graph "X".

      * The expansion of edge "G.e[ie]" results in edges "X.e[lo..hi-1]"
      of "X", where "lo = tb.eLo[ie]" and "hi = tb.eHi[ie]".

      * The expansion of edge "G.e[ie]" results may result in one or
      more `new' vertices of "X" (which do not correspond to any
      vertex of "G").  These vertices are consecutively numbered
      "[lo..hi-1]", where "lo = tb.vLo[ie]" and "hi = tb.eHi[ie]".

    If the block has no edges, or no new nodes, then the corresponding
    table entries have "hi = lo".

    The procedure assumes that every selected edge of "G" becomes a
    block with "spec.NV" vertices (including its source and sink) and
    "spec.NE" edges, while unselected edges are just copied.
  *)
  VAR ek: EdgeNum; vk: VertexNum;
  BEGIN
    WITH
      NEX = NEW(REF Count)^,
      NVX = NEW(REF Count)^,
      nInsert = NEW(REF Counts, G.NV)^,
      vMapRef = NEW(REF VertexNums, G.NV), vMap = vMapRef^,
      expRef = NEW(REF ARRAY OF BOOL, G.NE), exp = expRef^,
      eLoRef = NEW(REF EdgeNums, G.NE), eLo = eLoRef^,
      eHiRef = NEW(REF EdgeNums, G.NE), eHi = eHiRef^,
      vLoRef = NEW(REF VertexNums, G.NE), vLo = vLoRef^,
      vHiRef = NEW(REF VertexNums, G.NE), vHi = vHiRef^
    DO
      (* 
        First, choose the edges to be expanded,
        marking them in "tb.exp[ie]". 
        Also store in "nInsert[iv]" the number of new "X" vertices
        that will be inserted just after vertex number "iv" of "G".
        Also, for all "ie" in "[0..G.NE-1]", define "eLo[ie]"
        and "eHi[ie]" as defined above. Also sets "tb.vLo[ie]" 
        and "tb.vHi[ie]", but relative to the number of the "X" 
        vertex corresponding to the lower endpoint of that edge.
      *)
      FOR iv := 0 TO G.NV-1 DO nInsert[iv] := 0 END;
      ek := 0;
      FOR ie := 0 TO G.NE-1 DO
        (*
          At this point, "ek" is the number of "X" edges allocated
          to edges "G.e[0..ie-1]".  For each vertex number "k" in "G",
          "nInsert[k]" is the number of new "X" vertices produced by 
          those edges that will be inserted just after vertex "k".
        *)
        WITH 
          e = G.e[ie], 
          m = G.m[ie], 
          iv = MIN(e[0], e[1])
        DO
          eLo[ie] := ek;
          vLo[ie] := nInsert[iv];
          exp[ie] := m IN spec.marks 
            AND (spec.fraction = 1.0d0 OR spec.fraction >= rnd.longreal(0.0d0, 1.0d0));
          IF exp[ie] THEN
            ek := ek + spec.NE; 
            nInsert[iv] := nInsert[iv] + (spec.NV-2)
          ELSE
            ek := ek + 1;
          END;
          vHi[ie] := nInsert[iv];
          eHi[ie] := ek;
        END
      END;
      NEX := ek;

      (* 
        Now, accumulate all the entries in "nInsert" to obtain
        the new vertex nums "vMap":
      *)
      vk := 0;
      FOR iv := 0 TO G.NV-1 DO
        vMap[iv] := vk;
        vk := vk + 1 + nInsert[iv]
      END;
      NVX := vk;

      (*
        Finally, convert all entries of "vLo" and "(vHi-1)" from
        relative to absolute vertex numbers: 
      *)
      FOR ie := 0 TO G.NE-1 DO
        WITH 
          e = G.e[ie], 
          iv = MIN(e[0], e[1])
        DO
          vLo[ie] := vLo[ie] + vMap[iv] + 1;
          vHi[ie] := vHi[ie] + vMap[iv] + 1;
        END
      END;

      RETURN ExpTables{
        NV := NVX, NE := NEX,
        vMap := vMapRef, 
        exp := expRef,
        eLo := eLoRef, eHi := eHiRef,
        vLo := vLoRef, vHi := vHiRef
      }
    END;
  END ComputeExpTables;

(******************************************************************************)
(* SUBGRAPH CONSTRUCTION ROUTINES *********************************************)
(******************************************************************************)

(* 
  Each procedure in this section creates a subgraph of a specific type. *)

PROCEDURE AddSubGraph(
    kind: Kind;              (* Graph type *)
    s, t: CARDINAL;          (* The source and sink vertices. *)
    W: CARDINAL;             (* Maximum flow. *)
    L: CARDINAL;             (* Mean flow length. *)
    vLo, vHi: VertexNum;     (* New vertices (excl. "s" and "t") are "[vlo..vHi-1]". *)
    rnd: XRandom.T;          (* Randomness source *)
    (*OUT*) VAR e: Edges;    (* Edges. *)
    (*OUT*) VAR m: Marks;    (* Edge marks. *)
    (*OUT*) VAR c: Costs;    (* Edge capacities. *)
  ) =
  (*
    Creates a subgraph of the given "kind".
    
    The resulting subgraph has a total of "NE = NUMBER(e)" edges and 
    "NV = 2 + vHi - vLo" vertices.  The latter include a specified 
    `source' "s < vLo", a specified `sink' "t >= vHi", and zero or
    more `internal' vertices numbered "[vLo..vHi-1]".

    The subgraph is determined by its `kind' (a general family of
    graphs), the number of vertices "NV", the number of edges "NE", and
    two `shape' parameters "W >= 0" (the `width', usually the subgraph's
    maximum flow) and "L > 0" (its `length', usually the mean flow length).
    Within these parameters, the graph is randomly selected using the
    bit source "rnd".

    The edges of the subgraph will be stored in "e", and their
    capacities in "c".  The marks "m" depend on the graph type.
  *)
  BEGIN
    CASE kind OF 
    | Kind.Rete => AddReteSubGraph(
        W := W, L := L, 
        s := s, t := t,
        vLo := vLo, vHi := vHi,
        rnd := rnd,
        (*OUT*)
        e := e, m := m, c := c
      );
    END
  END AddSubGraph;

PROCEDURE AddReteSubGraph(
    s, t: CARDINAL;          (* The source and sink vertices. *)
    W: CARDINAL;             (* Maximum flow. *)
    L: CARDINAL;             (* Mean flow length. *)
    vLo, vHi: VertexNum;     (* New vertices (excl. "s" and "t") are "[vlo..vHi-1]". *)
    rnd: XRandom.T;          (* Randomness source *)
    (*OUT*) VAR e: Edges;    (* Edges. *)
    (*OUT*) VAR m: Marks;    (* Edge marks. *)
    (*OUT*) VAR c: Costs;    (* Edge capacities (all 1 for now). *)
  ) =
  (*
    Creates a "Rete" subgraph with source "s", sink "t", with maximum
    flow "W > 0" and mean flow length "L > 0".
    
    The subgraph "G" built by this procedure consists of two kinds of
    edges, `red' and `green'.
    
    The red edges make up "W" paths of length "L" from "s" to "t".
    The paths are edge- and vertex-disjoint, except for "s" and "t".
    The vertex numbers increase along each path; Moreover, if "i < j"
    then all vertices in path "i" are less than those in path "j".
    
    The green edges connect only internal vertices, in "[vLo..vHi]".
    Every green edge has destination less than source, so the green
    subgraph is acyclic and directed `against the grain' of the red
    subgraph.
    
    Green edges are classified  as `RR', `RN', `NR', or `NN', depending
    on whether their origin and destination belong to a red path (`R') or 
    not (`N'). A green edge is called `short' if it has both ends on
    the same red path, and is `long' otherwise.  
    
    The optimum flow value is obviously "W": the red paths demonstrate
    the lower bound, and a trivial cut at "s" or "t" proves the upper
    bound.
    
    Moreover, the red paths are the *unique* optimum flow.  This is
    equivalent to saying that the graph obtained by reversing the
    direction of all red edges is acyclic, which is obviously true.
    
    If "L > 1", the procedure requires that "NV >= W*(L-1) + 2", 
    "NE >= W*L", and "NE <= 2*W + (NV-2:2)", where "(n:k)" denotes
    the binomial coefficient: "n!/k!/(n-k)" (or 0 if "k>n").
    All edges will have unit capacity.
    
    If "L = 1", then "NE" must be at least "1" and at most "1 +
    (NV-2:2)". The red subgraph will be a single edge
    with capacity "W" connecting "s" and "t". All other edges
    will have unit capacity.
    
    On output, every edge "e[i]" will have cost "c[i] = 1".
    Red edges will be marked 'T', and green edges will be
    marked 'f'. 
  *)
  
  CONST
    RedMark = 'T';
    GreenMark = 'f';
  
  TYPE PathRange = ARRAY [0..1] OF VertexNum;
  
  PROCEDURE SelectRedPathNodes(): REF ARRAY OF PathRange =
    (*
      Returns a vector "pr" such that the internal vertices
      of red path number "k" are "[lo..hi-1]" where
      "lo = pr[k][0]" and "hi = pr[k][1]".
      
      The remaining (green) vertices are distributed more or 
      less uniformly between the paths.
    *)
    VAR remGV: Count; (* Green vertices yet to be assigned *)
    BEGIN
      <* ASSERT W > 0 *>
      <* ASSERT L > 0 *>
      WITH
        NV = vHi - vLo + 2,
        rpr = NEW(REF ARRAY OF PathRange, W), pr = rpr^,
        NGV = NV - 2 - W * (L - 1)
      DO
        (* Enough vertices for all red paths? *)
        <* ASSERT NV - 2 >= W * (L - 1) *>
        
        IF W = 1 THEN
          WITH vMd = (vLo + vHi) DIV 2 DO
            pr[0][0] := vMd - (L-1) DIV 2;
            pr[0][1] := pr[0][0] + L - 1;
          END
        ELSE
          <* ASSERT L > 1 *>
          pr[0][0] := vLo;
          pr[0][1] := pr[0][0] + L - 1;
          remGV := NGV;
          FOR k := 1 TO W-1 DO
            WITH g = remGV DIV (W-k) DO
              pr[k][0] := pr[k-1][1] + g;
              pr[k][1] := pr[k][0] + L - 1;
              <* ASSERT pr[k][1] <= vHi *>
              remGV := remGV - g
            END
          END;
          <* ASSERT remGV = 0 *>
        END;
        RETURN rpr
      END
    END SelectRedPathNodes;
    
  TYPE PathIndex = INTEGER;
  
  PROCEDURE AssignNodesToRedPaths(
      READONLY pr: ARRAY OF PathRange
    ): REF ARRAY OF PathIndex =
    (* 
      Returns a vector "px" such that "px[i]" is "k" if
      vertex "vLo + i" belongs to red path "k"; or "-1" if 
      it doesn't belong to any red path.
    *)
    BEGIN
      WITH
        NV = vHi - vLo + 2,
        rpx = NEW(REF ARRAY OF PathIndex, NV - 2), px = rpx^
      DO
        FOR v := 0 TO NV-3 DO px[v] := -1 END;
        FOR k := 0 TO W-1 DO
          FOR v := pr[k][0] TO pr[k][1]-1 DO
            px[v - vLo] := k
          END
        END;
        RETURN rpx
      END 
    END AssignNodesToRedPaths;
    
  VAR 
    NRE: CARDINAL;        (* Number of distinct red edges required. *)
    NGE: CARDINAL;        (* Number of distinct green edges required. *)
    ek: CARDINAL := 0;    (* Edge counter. *)
  
  PROCEDURE AddRedEdges(
      READONLY pr: ARRAY OF PathRange;
      <*UNUSED*> READONLY px: ARRAY OF PathIndex;
    ) =
    (* 
      Adds the "NRE" red edges.
    *)
    VAR ctXRE, ctIRE: Count := 0;
    BEGIN
      <* ASSERT L > 0 *>
      <* ASSERT W > 0 *>
      <* ASSERT NUMBER(e) - ek >= NRE *>
      IF L = 1 THEN
        <* ASSERT NRE = 1 *>
        e[ek] := Edge{s,t}; c[ek] := W; m[ek] := RedMark;
        INC(ek); INC(ctXRE)
      ELSE
        <* ASSERT vLo < vHi *>
        <* ASSERT NRE = W*L *>
        WITH
          NV = vHi - vLo + 2,
          (* Split "NE" among the red, short green, and long green edges: *)
          NIRE = MAX(0, W*(L-2)), (* Internal red edges. *) 
          NXRE = NRE - NIRE       (* External red edges, attached to "s" and/or "t" *)
        DO
          (* Enough vertices for the red edges? *)
          <* ASSERT NV - 2 >= W * (L - 1) *>

          FOR k := 0 TO W-1 DO
            WITH u = pr[k][0], v = pr[k][1]-1 DO
              e[ek] := Edge{s, u}; c[ek] := 1; m[ek] := RedMark;
              INC(ek); INC(ctXRE);
              FOR t := u TO v - 1 DO
                e[ek] := Edge{t, t+1}; c[ek] := 1; m[ek] := RedMark;
                INC(ek); INC(ctIRE)
              END;
              e[ek] := Edge{v, t}; c[ek] := 1; m[ek] := RedMark;
              INC(ek); INC(ctXRE);
            END
          END;
          <* ASSERT ctXRE = NXRE *>
          <* ASSERT ctIRE = NIRE *>
        END;
      END;
    END AddRedEdges;

  PROCEDURE AddGreenEdges(
      READONLY pr: ARRAY OF PathRange;
      <*UNUSED*> READONLY px: ARRAY OF PathIndex;
    ) =
    (* 
      Adds the "NGE" green edges. 
      Biased towards `perverse shortcuts'. *)
    
    PROCEDURE AddLongRREdges(
        VAR ctGE: Count;   (* Green edges already added *)
        VAR remGE: LONG;   (* Potential green edges remaining *)
      ) =
      (*
        Tries to add long RR edges, biased towards `perverse shortcuts'. *)
      BEGIN
        FOR d := 0 TO 2*(L-2) DO
          (* Enumerate long RR edges that create an "s-t" path of length "d+3". *)
          (* Any such path must consist of 1 + r edges from one red path, *)
          (* the green edge, and d - r + 1 edges from a lower-numbered red path, *)
          (* where both "r" and "d-r" are between 0 and "L-2". *)
          FOR t := 1 TO W-1 DO
            FOR r := MAX(0, d - (L-2)) TO MIN(L-2, d) DO
               FOR ku := 0 TO W-1-t DO
                WITH kv = ku + t DO
                  WITH 
                    v = pr[kv][0] + r,
                    u = pr[ku][1] - 1 - (d - r)
                  DO
                    <* ASSERT u < v *>
                    <* ASSERT v < pr[kv][1] *>
                    <* ASSERT u >= pr[ku][0] *>
                    WITH
                      prob = 0.75d0 + 0.25d0 * FLOAT(NGE-ctGE, LONG)/remGE
                    DO
                      IF rnd.longreal(0.0d0, 1.0d0) <= prob THEN 
                        e[ek] := Edge{v, u}; c[ek] := 1; 
                        m[ek] := GreenMark;
                        INC(ek); INC(ctGE);
                        IF ctGE >= NGE THEN RETURN END;
                      END;
                      remGE := remGE - 1.0d0
                    END
                  END
                END
              END
            END
          END;
        END;
      END AddLongRREdges;
    
    PROCEDURE AddRNAndNREdges(
        VAR ctGE: Count;   (* Green edges already added *)
        VAR remGE: LONG;   (* Potential green edges remaining *)
      ) =
      (*
        Tries to add RN and NR edges, biased towards shortcuts. *)
      VAR uLo, uHi, vR: VertexNum;
      BEGIN
        WITH
          NRV = W * (L-1),
          NGV = vHi - vLo - NRV,
          u = NEW(REF ARRAY OF VertexNum, NGV)^,
          k = NEW(REF ARRAY OF CARDINAL, NGV)^,
          pw = 1.0d0/FLOAT(W, LONG)
        DO
          (* Set "u[i]" to the "i"th green vertex, *)
          (* and "k[i]" to the number of red paths below it: *)
          VAR nG: CARDINAL := 0;
          BEGIN
            FOR kG := 0 TO W DO
              IF kG = 0 THEN uLo := vLo ELSE uLo := pr[kG-1][1] END;
              IF kG = W THEN uHi := vHi ELSE uHi := pr[kG][0] END;
              FOR uG := uLo TO uHi-1 DO
                u[nG] := uG; k[nG] := kG; INC(nG)
              END
            END;
            <* ASSERT nG = NGV *>
          END;
          (* Scramble the green vertices: *)
          FOR i := 0 TO LAST(u)-1 DO
            WITH j = rnd.integer(i+1, LAST(u)) DO
              VAR t := u[i]; BEGIN u[i] := u[j]; u[j] := t END;
              VAR t := k[i]; BEGIN k[i] := k[j]; k[j] := t END;
            END
          END;
          (* Scan green vertices: *)
          FOR i := 0 TO LAST(u) DO
            WITH uG = u[i], kG = k[i] DO
              (* Scan red vertices: *)
              FOR r := 0 TO L-2 DO
                FOR kR := 0 TO W-1 DO
                  IF kR < kG THEN 
                    vR := pr[kR][1] - 1 - r;
                    <* ASSERT vR < uG *>
                  ELSE
                    vR := pr[kR][0] + r;
                    <* ASSERT vR > uG *>
                  END;
                  WITH 
                    prob = pw + (1.0d0 - pw) * FLOAT(NGE-ctGE, LONG)/remGE
                  DO
                    IF rnd.longreal(0.0d0, 1.0d0) <= prob THEN 
                      IF uG < vR THEN 
                        e[ek] := Edge{vR, uG}
                      ELSE 
                        e[ek] := Edge{uG, vR}
                      END;
                      c[ek] := 1; m[ek] := GreenMark;
                      INC(ek); INC(ctGE);
                      IF ctGE >= NGE THEN RETURN END;
                    END;
                    remGE := remGE - 1.0d0
                  END;
                END
              END
            END
          END
        END
      END AddRNAndNREdges;
      
    PROCEDURE AddNNEdges(
        VAR ctGE: Count;   (* Green edges already added *)
        VAR remGE: LONG;   (* Potential green edges remaining *)
      ) =
      (*
        Tries to add NN edges, at random. *)
      VAR uLo, uHi: VertexNum;
          prob: LONG;
      BEGIN
        WITH
          NRV = W*(L-1),
          NGV = vHi - vLo - NRV,
          u = NEW(REF ARRAY OF VertexNum, NGV)^,
          pw = 1.0d0/FLOAT(W, LONG)
        DO
          (* Set "u[i]" to the "i"th green vertex, *)
          VAR nG: CARDINAL := 0;
          BEGIN
            FOR kG := 0 TO W DO
              IF kG = 0 THEN uLo := vLo ELSE uLo := pr[kG-1][1] END;
              IF kG = W THEN uHi := vHi ELSE uHi := pr[kG][0] END;
              FOR uG := uLo TO uHi-1 DO
                u[nG] := uG; INC(nG)
              END
            END;
            <* ASSERT nG = NGV *>
          END;
          (* Scramble the green vertices: *)
          FOR i := 0 TO LAST(u)-1 DO
            WITH j = rnd.integer(i+1, LAST(u)) DO
              VAR t := u[i]; BEGIN u[i] := u[j]; u[j] := t END;
            END
          END;
          (* Scan pairs of green vertices: *)
          FOR d := NUMBER(u)-1 TO 1 BY -1 DO
            FOR i := 0 TO LAST(u) - d DO
              WITH uG = u[i], vG = u[i+d] DO
                prob := pw + (1.0d0 - pw) * FLOAT(NGE-ctGE, LONG)/remGE;
                IF rnd.longreal(0.0d0, 1.0d0) <= prob THEN 
                  IF uG < vG THEN 
                    e[ek] := Edge{vG, uG}
                  ELSE 
                    e[ek] := Edge{uG, vG}
                  END;
                  c[ek] := 1; m[ek] := GreenMark;
                  INC(ek); INC(ctGE);
                  IF ctGE >= NGE THEN RETURN END;
                END;
                remGE := remGE - 1.0d0
              END
            END
          END
        END
      END AddNNEdges;
      
    PROCEDURE AddShortRREdges(
        VAR ctGE: Count;   (* Green edges already added *)
        VAR remGE: LONG;   (* Potential green edges remaining *)
      ) =
      (*
        Tries to add the short RR edges. *)
      BEGIN
        WITH
          pw = 1.0d0/FLOAT(W, LONG)
        DO
          FOR d := 1 TO (L-2) DO
            (* Enumerate short RR edges that junp back by "d" edges along. *)
            (* some red path. *)
            FOR r := d TO L-2 DO
              FOR k := 0 TO W-1 DO
                WITH 
                  v = pr[k][0] + r,
                  u = v - d
                DO
                  <* ASSERT u < v *>
                  <* ASSERT v < pr[k][1] *>
                  <* ASSERT u >= pr[k][0] *>
                  WITH
                    prob = pw + (1.0d0 - pw) * FLOAT(NGE-ctGE, LONG)/remGE
                  DO
                    IF rnd.longreal(0.0d0, 1.0d0) <= prob THEN 
                      e[ek] := Edge{v, u}; c[ek] := 1; 
                      m[ek] := GreenMark;
                      INC(ek); INC(ctGE);
                      IF ctGE >= NGE THEN RETURN END;
                    END;
                    remGE := remGE - 1.0d0
                  END
                END
              END
            END
          END
        END;
      END AddShortRREdges;
    
    VAR ctGE: Count;  (* Number of green edges already added. *)
        remGE: LONG;  (* Number of remaining chances to pick a green edge. *)
    BEGIN
      <* ASSERT L > 0 *>
      <* ASSERT W > 0 *>
      <* ASSERT NUMBER(e) - ek >= NGE *>
      (* Enumerate the decreasing pairs in order of decreasing desirability.*)
      WITH
        NV = vHi - vLo + 2
      DO
        ctGE := 0;
        remGE := Choose2(NV-2); (* Max possible number of green edges. *)

        (* Enough vertices for all green edges? *)
        <* ASSERT FLOAT(NGE, LONG) <= remGE *>

        IF ctGE >= NGE THEN RETURN END;
        AddLongRREdges(ctGE, remGE);
        IF ctGE >= NGE THEN RETURN END;
        AddRNAndNREdges(ctGE, remGE);
        IF ctGE >= NGE THEN RETURN END;
        AddNNEdges(ctGE, remGE);
        IF ctGE >= NGE THEN RETURN END;
        AddShortRREdges(ctGE, remGE);
        <* ASSERT ctGE = NGE *>
        <* ASSERT remGE >= 0.0d0 *>
      END     
    END AddGreenEdges;

  BEGIN (* AddReteSubGraph *)
    <* ASSERT W > 0 *>
    <* ASSERT L > 0 *>
    WITH
      NV = vHi - vLo + 2,
      NE = NUMBER(e)
    DO
      (* Compute number of red and green edges required: *)
      IF L = 1 THEN NRE := 1 ELSE NRE := W*L END;
      NGE := MAX(0, NE - NRE);
      WITH
        minV = FLOAT(W * (L - 1) + 2, LONG),
        maxV = FLOAT(MaxVertices, LONG),
        minE = FLOAT(NRE, LONG),
        maxE = FLOAT(NRE, LONG) + Choose2(NV-2)
      DO
        Require("NV", NV, minV, maxV);
        Require("NE", NE, minE, maxE);
        WITH
          pr = SelectRedPathNodes()^,
          px = AssignNodesToRedPaths(pr)^
        DO
          ek := 0;
          AddRedEdges(pr, px);
          AddGreenEdges(pr, px);
          <* ASSERT ek = NE *>
        END
      END
    END;
  END AddReteSubGraph;
  
(******************************************************************************)
(* UTILITY ROUTINES ***********************************************************)
(******************************************************************************)

PROCEDURE Choose2(n: CARDINAL): LONG = 
  BEGIN
    WITH fn = FLOAT(n, LONG) DO
      RETURN fn*(fn - 1.0d0)/2.0d0
    END;
  END Choose2;
  
<*UNUSED*>
PROCEDURE RB(rnd: XRandom.T): BOOLEAN =
  BEGIN
    RETURN rnd.boolean()
  END RB;

<*UNUSED*>
PROCEDURE RP(
    rnd: XRandom.T; 
    n: CARDINAL; 
    lo, hi: CARDINAL; 
    VAR v: ARRAY OF CARDINAL;
  ) =
  (*
    Puts in "v[0..n-1]" a random sequence of "n" distinct
    elements from [lo..hi].
  *)
  BEGIN
    RS(rnd, n, lo, hi, v);
    Scramble(rnd, v);
  END RP;
  
PROCEDURE Scramble(rnd: XRandom.T; VAR v: ARRAY OF CARDINAL) =
  (* Permutes "v" in random order *)
  VAR t: CARDINAL;
  BEGIN
    FOR k := LAST(v) TO 1 BY -1 DO 
      WITH i = rnd.integer(0, k) DO 
        IF i # k THEN
          t := v[i]; v[i] := v[k]; v[k] := t
        END
      END
    END
  END Scramble;
  
PROCEDURE RS(
    rnd: XRandom.T; 
    n: CARDINAL; 
    lo, hi: CARDINAL; 
    VAR v: ARRAY OF CARDINAL;
  ) =
  (*
    Puts in "v[0..n-1]" a random subset of "n" distinct
    elements from [lo..hi], in increasing order.
  *)
  VAR t, e: INTEGER;
  BEGIN
    <* ASSERT hi >= lo + n - 1 *>
    FOR W := 0 TO n-1 DO
      e := rnd.integer(lo + W, hi);
      t := W;
      WHILE t > 0 AND v[t-1] >= e DO v[t] := v[t-1]; t := t-1; e := e - 1 END;
      v[t] := e
    END
  END RS;
  
<*UNUSED*>
PROCEDURE SQRT(x: REAL): REAL =
  BEGIN
    RETURN FLOAT(Math.sqrt(FLOAT(x, LONG)))
  END SQRT;

<*UNUSED*>
PROCEDURE RR(rnd: XRandom.T): REAL = 
  BEGIN
    RETURN rnd.real()
  END RR;

PROCEDURE Require(what: TEXT; x: INTEGER; lo, hi: LONG) =
  BEGIN
    IF FLOAT(x, LONG) < lo OR FLOAT(x, LONG) > hi THEN
      Wr.PutText(stderr, "Invalid parameter: ");
      Wr.PutText(stderr, what);
      Wr.PutText(stderr, " is ");
      Wr.PutText(stderr, Fmt.Int(x));
      Wr.PutText(stderr, ", should be in [");
      Wr.PutText(stderr, Fmt.LongReal(lo, style := Fmt.Style.Fix, prec := 0));
      Wr.PutText(stderr, "..");
      Wr.PutText(stderr, Fmt.LongReal(hi, style := Fmt.Style.Fix, prec := 0));
      Wr.PutText(stderr, "]\n");
      Wr.Flush(stderr);
      Process.Exit(1)
    END
  END Require;

BEGIN
  DoIt()
END MaxFlowProbs.