MODULE GDDrawing;

IMPORT Rd, Wr, Fmt, Text, Thread; 
FROM Stdio IMPORT stderr;

IMPORT FileFmt, PSPlot, FGet, NGet, FPut, NPut;
FROM PSPlot IMPORT Color;

IMPORT GDGraph AS Graph;
IMPORT GDGeometry AS Geo;
FROM GDGraph IMPORT 
  Node, Nodes, InsertNode, DeleteNode, NodeDict, 
  Dir, Degrees, NoNodes;
FROM GDGeometry IMPORT 
  Point, Points, Interval, Weights, Axis, NAT, INT, LONG, BOOL, SIGN;

REVEAL
  T = PublicT BRANDED OBJECT
      injective: BOOL;                 (* TRUE means all node positions must be distinct. *)
      yRange: Interval;                (* Nominal range of "y" coordinates. *)
      layer: REF ARRAY OF LayerInfo;   (* Layer nodes. *)
        (*
          Invariant: all nodes of "D" have positions in the range 
          "[yRange.lo..yRange.hi]".  For each "y" in that range,
          the nodes in layer "y" are listed in "layer[y-yRange.lo]",
          sorted by abscissa. *)
          
    OVERRIDES
      init := InitMethod;
      copy := CopyMethod;
      baseNode := BaseNodeMethod;
      layerNodes := LayerNodesMethod;
      layerXRange := LayerXRangeMethod;
      nodeAt := NodeAtMethod;
      nodeRunAt := NodeRunAtMethod;
      addJoint := AddJointMethod;
      moveNode := MoveNodeMethod;
      moveNodes := MoveNodesMethod;
      moveAllNodes := MoveAllNodesMethod;
      crunch := CrunchMethod;
    END;

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

PROCEDURE InitMethod(dr: T; G: Graph.T; injective: BOOL; allocJoints: NAT := 20): T =
  BEGIN
    WITH
      nBaseV = G.nV,
      allocNodes = nBaseV + allocJoints
    DO
      dr.G := G;
      (* Create the drawing graph "D", add the base nodes to it: *)
      WITH
        D = NEW(Graph.T).init(allocNodes := allocNodes),
        posR = NEW(REF Points, allocNodes), 
        pos = SUBARRAY(posR^, 0, nBaseV),
        yBase = 0
      DO
        FOR u := 0 TO nBaseV-1 DO
          WITH v = D.addNode(G.degrees(u)), y = yBase + u DO 
            <* ASSERT u = v *>
            pos[u] := Point{0, y}
          END;
        END;
        dr.D := D;
        dr.pos := posR;
        dr.injective := injective;
      END;
      RebuildLayerTable(dr, allocLayers := allocNodes);
      RETURN dr
    END
  END InitMethod;
  
PROCEDURE CopyMethod(dr: T): T =
  BEGIN
    WITH
      nV = dr.D.nV,
      allocNodes = NUMBER(dr.pos^),
      pos = SUBARRAY(dr.pos^,0,nV),
      cp = NEW(T)
    DO
      cp.G := dr.G;
      cp.D := dr.D.copy();
      cp.pos := NEW(REF Points, allocNodes); 
      cp.injective := dr.injective;
      SUBARRAY(cp.pos^,0,nV) := pos;
      WITH allocLayers = MAX(0, dr.yRange.hi - dr.yRange.lo + 1) DO
        RebuildLayerTable(cp, allocLayers := allocLayers);
      END;
      RETURN cp
    END
  END CopyMethod;

PROCEDURE BaseNodeMethod(dr: T; u: Node; dir: Dir): Node =
  BEGIN
    <* ASSERT u < dr.D.nV *>
    WITH D = dr.D, nV = dr.G.nV DO
      WHILE u # NoNode AND u >= nV DO u := D.nbor(u, dir, 0) END;
      RETURN u
    END
  END BaseNodeMethod;

PROCEDURE LayerNodesMethod(dr: T; y: INT): REF Nodes =
  BEGIN
    WITH f = GetLayerInfo(dr, y) DO
      IF f.n = 0 THEN
        RETURN NoNodes
      ELSE
        RETURN Graph.CopyNodes(SUBARRAY(f.nodes^,0,f.n))
      END
    END
  END LayerNodesMethod;

PROCEDURE LayerXRangeMethod(dr: T; y: INT): Interval =
  BEGIN
    WITH f = GetLayerInfo(dr, y) DO
      IF f.n = 0 THEN
        RETURN Interval{lo := +1, hi := -1}
      ELSE
        WITH 
          uu = f.nodes^, 
          pos = SUBARRAY(dr.pos^,0,dr.D.nV),
          uLo = uu[0], uHi = uu[f.n-1]
        DO
          RETURN Interval{pos[uLo][0], pos[uHi][0]}
        END
      END
    END
  END LayerXRangeMethod;
  
PROCEDURE NodeAtMethod(dr: T; p: Point): Node =
  VAR lo, hi: INT;
  BEGIN
    WITH layer = GetLayerInfo(dr, p[1]), n = layer.n DO
      IF n = 0 THEN 
        RETURN NoNode
      ELSE
        WITH
          uu = SUBARRAY(layer.nodes^, 0, n),
          pos = SUBARRAY(dr.pos^, 0, dr.D.nV)
        DO
          GetNodesAtX(p[0], uu, pos, lo, hi);
          (* We'd better have at most one node there: *)
          <* ASSERT lo >= hi *>
          IF lo > hi THEN RETURN NoNode ELSE RETURN uu[lo] END
        END
      END
    END
  END NodeAtMethod;

PROCEDURE NodeRunAtMethod(dr: T; p: Point): REF Nodes =
  VAR lo, hi: INT; vvR: REF Nodes := NoNodes;
  BEGIN
    WITH layer = GetLayerInfo(dr, p[1]), n = layer.n DO
      (* Wr.PutText(stderr, "node run at " & FPt(p) & " ="); *)
      IF n > 0 THEN
        WITH
          uu = SUBARRAY(layer.nodes^, 0, n),
          pos = SUBARRAY(dr.pos^, 0, dr.D.nV)
        DO
          GetNodeRunAtX(p[0], uu, pos, lo, hi);
          IF lo <= hi THEN
            <* ASSERT pos[uu[lo]][0] <= p[0] *>
            <* ASSERT pos[uu[hi]][0] >= p[0] *>
            WITH m = hi - lo + 1 DO
              vvR := NEW(REF Nodes, m);
              vvR^ := SUBARRAY(uu, lo, m);
            END;
            (*
            FOR k := lo TO hi DO
              Wr.PutText(stderr, " " & FI(uu[k]) & FPt(pos[uu[k]]))
            END
            *)
          END
        END
      END;
      (* Wr.PutText(stderr, "\n"); *)
    END;
    RETURN vvR
  END NodeRunAtMethod;
  
PROCEDURE AddJointMethod(dr: T; p: Point): Node =
  BEGIN
    EnsureLayerExists(dr, p[1]);
    WITH
      u = dr.D.addNode(Degrees{1,0,1})
    DO
      EnsurePosExists(dr.pos, u);
      dr.pos[u] := p;
      WITH layer = dr.layer[p[1] - dr.yRange.lo] DO
        InsertNodeInLayer(u, layer, SUBARRAY(dr.pos^, 0, dr.D.nV), dr.injective)
      END;
      RETURN u
    END
  END AddJointMethod;
  
PROCEDURE MoveNodeMethod(dr: T; u: Node; p: Point) =
  BEGIN
    <* ASSERT u < dr.D.nV *>
    WITH
      pos = SUBARRAY(dr.pos^, 0, dr.D.nV),
      injective = dr.injective,
      xOld = pos[u][0] + 0, yOld = pos[u][1] + 0,
      xNew = p[0], yNew = p[1]
    DO
      (* Wr.PutText(stderr, "moving " & FI(u) & " from " & FPt(pos[u])); *)
      (* Wr.PutText(stderr, " to " & FPt(p) & "\n"); *)
      IF yOld # yNew THEN
        WITH oldLayer = dr.layer[yOld - dr.yRange.lo] DO
          DeleteNodeFromLayer(u, oldLayer, pos)
        END;
        pos[u] := p;
        EnsureLayerExists(dr, yNew);
        WITH newLayer = dr.layer[yNew - dr.yRange.lo] DO
          InsertNodeInLayer(u, newLayer, pos, injective)
        END
      ELSIF xOld # xNew THEN
        pos[u] := p;
        WITH layer = dr.layer[yNew - dr.yRange.lo] DO
          ReSortLayer(layer, pos, injective)
        END
      END;
    END
  END MoveNodeMethod;
  
PROCEDURE MoveNodesMethod(dr: T; READONLY uu: Nodes; READONLY pp: Points) =
  BEGIN
    <* ASSERT NUMBER(uu) = NUMBER(pp) *>
    WITH
      pos = SUBARRAY(dr.pos^, 0, dr.D.nV),
      injective = dr.injective
    DO
      (* First pass - delete nodes from layers. *)
      (* 
        (It would be more efficient to delete only those nodes that
        change layers, but we don't know how to identify them in the
        third pass.)
      *)
      FOR i := 0 TO LAST(uu) DO
        <* ASSERT uu[i] < dr.D.nV *>
        WITH u = uu[i], yOld = pos[u][1] + 0 DO
          (* Wr.PutText(stderr, "moving " & FI(u) & " from " & FPt(pos[u])); *)
          (* Wr.PutText(stderr, " to " & FPt(pp[i]) & "\n"); *)
          WITH oldLayer = dr.layer[yOld - dr.yRange.lo] DO
            DeleteNodeFromLayer(u, oldLayer, pos)
          END;
        END
      END;
      (* Second pass - move nodes to new positions. *)
      FOR i := 0 TO LAST(uu) DO pos[uu[i]] := pp[i] END;
      (* Third pass - reinsert nodes in the layer structure. *)
      FOR i := 0 TO LAST(uu) DO
        WITH u = uu[i] + 0, yNew = pp[i][1] + 0 DO
          EnsureLayerExists(dr, yNew);
          WITH newLayer = dr.layer[yNew - dr.yRange.lo] DO
            InsertNodeInLayer(u, newLayer, pos, injective)
          END
        END
      END;
    END
  END MoveNodesMethod;
  
PROCEDURE MoveAllNodesMethod(dr: T; READONLY pp: Points) =
  BEGIN
    WITH
      pos = SUBARRAY(dr.pos^, 0, dr.D.nV),
      injective = dr.injective
    DO
      <* ASSERT NUMBER(pp) <= NUMBER(pos) *>
      (* First pass - delete nodes from layers. *)
      (* 
        (It would be more efficient to delete only those nodes that
        change layers, but we don't know how to identify them in the
        third pass.)
      *)
      FOR u := 0 TO LAST(pp) DO
        WITH yOld = pos[u][1] + 0 DO
          (* Wr.PutText(stderr, "moving " & FI(u) & " from " & FPt(pos[u])); *)
          (* Wr.PutText(stderr, " to " & FPt(pp[u]) & "\n"); *)
          WITH oldLayer = dr.layer[yOld - dr.yRange.lo] DO
            DeleteNodeFromLayer(u, oldLayer, pos)
          END;
        END
      END;
      (* Second pass - move nodes to new positions. *)
      FOR u := 0 TO LAST(pp) DO pos[u] := pp[u] END;
      (* Third pass - reinsert nodes in the layer structure. *)
      FOR u := 0 TO LAST(pp) DO
        WITH yNew = pp[u][1] + 0 DO
          EnsureLayerExists(dr, yNew);
          WITH newLayer = dr.layer[yNew - dr.yRange.lo] DO
            InsertNodeInLayer(u, newLayer, pos, injective)
          END
        END
      END;
    END
  END MoveAllNodesMethod;
  
PROCEDURE CrunchMethod(dr: T): NodeDict =
  BEGIN
    (* Wr.PutText(stderr, "(crunching drawing..."); *)
    WITH
      nOldV = dr.D.nV + 0,
      dict = dr.D.crunch(dr.G.nV),
      nNewV = dr.D.nV + 0,
      new = SUBARRAY(dict.new^, 0, nOldV),
      old = SUBARRAY(dict.old^, 0, nNewV)
    DO
      (* Map position table: *)
      WITH 
        posOld = SUBARRAY(dr.pos^, 0, nOldV),
        allocNodes = (3*nNewV) DIV 2 + 20,
        posNewR = NEW(REF Points, allocNodes), 
        posNew = SUBARRAY(posNewR^, 0, nNewV)
      DO
        FOR uNew := 0 TO nNewV-1 DO
          posNew[uNew] := posOld[old[uNew]]
        END;
        dr.pos := posNewR
      END;
      (* Cleanup and translate layer table: *)
      TranslateLayerTable(dr, new);
      (* Wr.PutText(stderr, "done)"); *)
      RETURN dict
    END
  END CrunchMethod;
  
(* UTILITY TOOLS *)

PROCEDURE MedBaseY(dr: T): LONG =
  BEGIN
    WITH nBaseV = dr.G.nV, pos = SUBARRAY(dr.pos^, 0, nBaseV) DO
      RETURN Geo.Median(pos, 1)
    END
  END MedBaseY;
  
PROCEDURE MinBaseY(dr: T): LONG =
  BEGIN
    WITH nBaseV = dr.G.nV, pos = SUBARRAY(dr.pos^, 0, nBaseV) DO
      RETURN FLOAT(Geo.GetRange(pos, 1).lo, LONG)
    END
  END MinBaseY;
  
PROCEDURE MaxBaseY(dr: T): LONG =
  BEGIN
    WITH nBaseV = dr.G.nV, pos = SUBARRAY(dr.pos^, 0, nBaseV) DO
      RETURN FLOAT(Geo.GetRange(pos, 1).hi, LONG)
    END
  END MaxBaseY;

PROCEDURE GetNodesAtX(
    x: INT;
    READONLY uu: Nodes;
    READONLY pos: Points;
    VAR lo, hi: INT;
  ) =
  BEGIN
    WITH r = NUMBER(uu)-1 DO
      lo := 0; WHILE lo <= r AND pos[uu[lo]][0] < x DO INC(lo) END;
      hi := r; WHILE hi >= 0 AND pos[uu[hi]][0] > x DO DEC(hi) END;
    END
  END GetNodesAtX;

PROCEDURE GetNodeRunAtX(
    x: INT;
    READONLY uu: Nodes;
    READONLY pos: Points;
    VAR lo, hi: INT;
  ) =
  BEGIN
    GetNodesAtX(x, uu, pos, lo, hi);
    IF lo <= hi THEN
      WITH r = NUMBER(uu)-1 DO
        WHILE lo > 0 AND pos[uu[lo-1]][0] >= pos[uu[lo]][0] - 1 DO DEC(lo) END;
        WHILE hi < r AND pos[uu[hi+1]][0] <= pos[uu[hi]][0] + 1 DO INC(hi) END;
      END
    END
  END GetNodeRunAtX;

(* LAYER TABLE MANAGEMENT TOOLS *)

TYPE
  LayerInfo = RECORD 
      n: NAT;           (* Number of nodes in layer. *)
      nodes: REF Nodes; (* "nodes[0..n-1]" are the nodes in layer. *)
    END;

PROCEDURE GetLayerInfo(dr: T; y: INT): LayerInfo =
  BEGIN
    WITH yLo = dr.yRange.lo, yHi = dr.yRange.hi DO
      IF y < yLo OR y > yHi THEN 
        RETURN LayerInfo{n := 0, nodes := NIL}
      ELSE
        RETURN dr.layer[y-yLo]
      END
    END
  END GetLayerInfo;

PROCEDURE InsertNodeInLayer(
    u: Node; 
    VAR layer: LayerInfo; 
    READONLY pos: Points;
    injective: BOOL
  ) =
  BEGIN
    InsertNode(u, layer.nodes, layer.n);
    ReSortNodesByX(SUBARRAY(layer.nodes^, 0, layer.n), pos, injective)
  END InsertNodeInLayer;
  
PROCEDURE DeleteNodeFromLayer(u: Node; VAR layer: LayerInfo; READONLY pos: Points) =
  VAR nOld := layer.n;
  BEGIN
    DeleteNode(u, layer.nodes, layer.n);
    <* ASSERT layer.n = nOld-1 *>
    ReSortNodesByX(SUBARRAY(layer.nodes^, 0, layer.n), pos, FALSE)
  END DeleteNodeFromLayer;
  
PROCEDURE ReSortLayer(VAR layer: LayerInfo; READONLY pos: Points; injective: BOOL) =
  BEGIN
    ReSortNodesByX(SUBARRAY(layer.nodes^, 0, layer.n), pos, injective)
  END ReSortLayer;
  
PROCEDURE SortNodesByCoords(VAR uu: Nodes; READONLY pos: Points) =
  
  PROCEDURE Compare(u, v: Node): SIGN =
    BEGIN
      RETURN CompareNodesByCoords(u, v, pos)
    END Compare;
    
  BEGIN
    Graph.SortNodes(uu, Compare)
  END SortNodesByCoords;

PROCEDURE ReSortNodesByX(VAR uu: Nodes; READONLY pos: Points; injective: BOOL) =
  (*
    Reorders the nodes in "uu" by abscissa; handy after minor changes. *)
  VAR j: NAT;
  BEGIN
    (* Insertion sort. Should be good enough. *)
    FOR i := 1 TO LAST(uu) DO
      VAR u := uu[i]; xu := pos[u][0]; v: Node; xv: INT;
      BEGIN
        j := i;
        LOOP
          IF j = 0 THEN EXIT END;
          v := uu[j-1]; xv := pos[v][0];
          (* Here we enforce distinctness: *)
          IF injective THEN <* ASSERT xu # xv *> END;
          IF xu >= xv THEN EXIT END;
          uu[j] := v; DEC(j)
        END;
        IF j < i THEN uu[j] := u END;
      END
    END;
    (* Paranoid check: *)
    IF injective THEN 
      FOR i := 1 TO LAST(uu) DO <* ASSERT pos[uu[i-1]][0] < pos[uu[i]][0] *> END
    ELSE
      FOR i := 1 TO LAST(uu) DO <* ASSERT pos[uu[i-1]][0] <= pos[uu[i]][0] *> END
    END
  END ReSortNodesByX;

PROCEDURE CompareNodesByCoords(u, v: Node; READONLY pos: Points): SIGN =
  BEGIN
    WITH
      pu = pos[u], xu = pu[0], yu = pu[1],
      pv = pos[v], xv = pv[0], yv = pv[1]
    DO
      IF yu < yv THEN 
        RETURN -1
      ELSIF yu > yv THEN
        RETURN +1
      ELSIF xu < xv THEN 
        RETURN -1
      ELSIF xu > xv THEN
        RETURN +1
      ELSIF u < v THEN 
        RETURN -1
      ELSIF u > v THEN
        RETURN +1
      ELSE
        RETURN 00
      END
    END
  END CompareNodesByCoords;

PROCEDURE TranslateLayerTable(dr: T; READONLY new: Nodes) =
  BEGIN
    WITH
      yLo = dr.yRange.lo,
      yHi = dr.yRange.hi
    DO
      IF yLo <= yHi THEN
        FOR k := 0 TO yHi-yLo DO
          TranslateLayer(dr.layer[k], new)
        END
      END
    END;
  END TranslateLayerTable;

PROCEDURE TranslateLayer(VAR layer: LayerInfo; READONLY new: Nodes) =
  VAR n: NAT := 0;
  BEGIN
    IF layer.n > 0 THEN
      WITH uu = layer.nodes^ DO
        FOR kOld := 0 TO layer.n-1 DO
          WITH uOld = uu[kOld], uNew = new[uOld] DO
            IF uNew # NoNode THEN uu[n] := uNew; INC(n) END
          END
        END
      END
    END;
    layer.n := n;
    IF n = 0 THEN
      layer.nodes := NIL
    ELSIF 4*n + 20 < NUMBER(layer.nodes^) THEN
      (* Reallocate "layer.nodes^": *)
      WITH tR = NEW(REF Nodes, 2*n+10) DO
        SUBARRAY(tR^,0,n) := SUBARRAY(layer.nodes^,0,n); 
        layer.nodes := tR
      END
    END
  END TranslateLayer;
  
PROCEDURE RebuildLayerTable(dr: T; allocLayers: NAT) =
  (*
    Rebuilds the layer data structure of "dr" from the drawing graph "dr.D" 
    and the node position table "dr.pos".  Initially allow for at least
    "allocLayers" layers.  *)
  BEGIN
    WITH
      D = dr.D,
      pos = SUBARRAY(dr.pos^, 0, D.nV),
      yRange = Geo.GetRange(pos, 1),
      minLayers = MAX(0, yRange.hi - yRange.lo),
      extraLayers = MAX(0, allocLayers - minLayers),
      yLo = yRange.lo - (extraLayers DIV 2),
      yHi = MAX(yRange.hi, yRange.lo) + (extraLayers DIV 2),
      nLayers = yHi - yLo + 1
    DO
      dr.yRange := Interval{yLo, yHi};
      dr.layer := NEW(REF ARRAY OF LayerInfo, nLayers);
      <* ASSERT LAST(dr.layer^) = dr.yRange.hi - dr.yRange.lo *>
      FOR k := 0 TO nLayers-1 DO dr.layer[k] := LayerInfo{0, NIL} END;
      (* Create the drawing graph "D", add the base nodes to it: *)
      FOR u := 0 TO D.nV-1 DO
        WITH y = pos[u][1] DO
          InsertNodeInLayer(u, dr.layer[y-yLo], pos, FALSE)
        END;
      END;
    END
  END RebuildLayerTable;
  
PROCEDURE EnsureLayerExists(dr: T; y: INT) =
  (*
    Expands "dr.yRange", "dr.layerV", and "dr.layerN" if necessary
    to include layer "y". *)
  VAR yJoin: Interval; mOld: NAT;
  BEGIN
    WITH yLo = dr.yRange.lo, yHi = dr.yRange.hi DO
      IF y < yLo OR y > yHi THEN
        IF yLo > yHi THEN
          mOld := 0; yJoin := Interval{y,y}
        ELSE
          mOld := yHi - yLo + 1; yJoin := Interval{MIN(y,yLo), MAX(y,yHi)}
        END;
        WITH 
          mJoin = yJoin.hi - yJoin.lo + 1,
          mNew = MAX(2*mOld + 20, mJoin),
          yNewLo = yJoin.lo - (mNew - mJoin) DIV 2,
          yNewHi = yNewLo + mNew - 1
        DO
          ReallocLayerTable(dr, yNewLo, yNewHi)
        END
      END;
    END
  END EnsureLayerExists;

PROCEDURE ReallocLayerTable(dr: T; yNewLo, yNewHi: INT) =
  (*
    Reallocates the layer table of "dr" to span the range 
    "[yNewLo..yNewHi]" exactly. *)
  BEGIN
    WITH
      yOldLo = dr.yRange.lo + 0,
      yOldHi = dr.yRange.hi + 0,
      mNew = MAX(0, yNewHi - yNewLo + 1),
      layerNewR = NEW(REF ARRAY OF LayerInfo, mNew),
      layerNew = layerNewR^
    DO
      FOR k := 0 TO mNew-1 DO
        WITH y = yNewLo + k DO
          IF y < yOldLo OR y > yOldHi THEN
            layerNew[k] := LayerInfo{0, NIL};
          ELSE
            layerNew[k] := dr.layer[y - yOldLo];
          END
        END
      END;
      dr.layer := layerNewR;
      dr.yRange := Interval{yNewLo, yNewHi}
    END
  END ReallocLayerTable;

(* POSITION TABLE MANAGEMENT TOOLS *)

PROCEDURE EnsurePosExists(VAR posR: REF Points; u: Node) =
  (*
    Expands "posR^" if needed to include node "u". *)
  VAR mOld: NAT;
  BEGIN
    IF posR = NIL THEN mOld := 0 ELSE mOld := NUMBER(posR^) END;
    IF u >= mOld THEN
      WITH
        mNew = MAX(2*mOld + 20, u+1), 
        posNewR = NEW(REF Points, mNew)
      DO
        IF mOld > 0 THEN SUBARRAY(posNewR^,0,mOld) := posR^ END;
        posR := posNewR
      END
    END
  END EnsurePosExists;

CONST DrawingFileVersion = "2000-01-07";

PROCEDURE Read(rd: Rd.T; allocJoints: NAT := 0): T =
  VAR nBaseV: NAT;
  BEGIN
    Wr.PutText(stderr, "(drawing header...");
    FileFmt.ReadHeader(rd, "drawing", DrawingFileVersion);
    EVAL FileFmt.ReadComment(rd, '|');
    Wr.PutText(stderr, "base node count...");
    nBaseV := NGet.Int(rd, "baseNodes"); FGet.EOL(rd);
    Wr.PutText(stderr, "base graph...");
    WITH
      D = Graph.Read(rd, nBaseV + allocJoints),
      posR = NEW(REF Points, nBaseV + allocJoints), 
      pos = SUBARRAY(posR^, 0, D.nV)
    DO
      Wr.PutText(stderr, "node positions...");
      FOR v := 0 TO D.nV-1 DO 
        WITH x = FGet.Int(rd), y = FGet.Int(rd) DO
          pos[v] := Point{x,y}
        END;
        FGet.EOL(rd);
      END;
      Wr.PutText(stderr, "drawing footer...");
      FileFmt.ReadFooter(rd, "drawing");
      VerifyDegreeInvariants(D, nBaseV);
      WITH 
        G = RecoverBaseGraph(D, nBaseV),
        dr = NEW(T)
      DO
        dr.G := G;
        dr.D := D;
        dr.pos := posR;
        RebuildLayerTable(dr, allocLayers := 0);
        Wr.PutText(stderr, "done)");
        RETURN dr
      END;
    END
  END Read;
  
PROCEDURE RecoverBaseGraph(D: Graph.T; nBaseV: NAT): Graph.T =
  (*
    Extracts the base graph from a drawing graph "D", given the
    number of base nodes "nBaseV". *)
  BEGIN
    WITH G = NEW(Graph.T).init(nBaseV) DO
      (* Create the vertices of the base graph "G": *)
      <* ASSERT D.nV >= nBaseV *>
      FOR u := 0 TO nBaseV-1 DO
        WITH deg = D.degrees(u), v = G.addNode(deg) DO
          <* ASSERT u = v *>
        END
      END;
      (* Create the edges of the base graph "G": *)
      FOR u := 0 TO nBaseV-1 DO
        FOR ku := 0 TO D.degree(u,+1)-1 DO
          VAR r := u; v := D.nbor(u,+1,ku);
          BEGIN
            (* Follow polyline until a base base node: *)
            WHILE v # NoNode AND v >= nBaseV DO r := v; v := D.nbor(r,+1,0) END;
            IF v # NoNode THEN
              (* Insert corresponding edge in "G": *)
              WITH kv = D.findSlot(v,-1,r) DO G.addEdge(+1,u,ku,v,kv) END
            END
          END
        END
      END;
      RETURN G
    END
  END RecoverBaseGraph;

PROCEDURE Write(wr: Wr.T; dr: T; comment: TEXT := "") =
  BEGIN
    FileFmt.WriteHeader(wr, "drawing", DrawingFileVersion);
    IF NOT Text.Empty(comment) THEN FileFmt.WriteComment(wr, comment, '|') END;
    NPut.Int(wr, "baseNodes", dr.G.nV); FPut.EOL(wr);
    Graph.Write(wr, dr.D, comment := "");
    WITH nV = dr.D.nV, pos = SUBARRAY(dr.pos^, 0, nV) DO
      FOR v := 0 TO nV-1 DO 
        FPut.Int(wr, pos[v][0]);
        FPut.Int(wr, pos[v][1]);
        FPut.EOL(wr);
      END
    END;
    FileFmt.WriteFooter(wr, "drawing");
  END Write;
  
PROCEDURE PlotEPS(
    name: TEXT; 
    dr: T;
    label: REF ARRAY OF TEXT; 
    labelPad: NAT;            
    fontSize: REAL;         
    scale: LONG := 1.0d0;     
  ) =
  VAR eps: PSPlot.EPSFile;
  BEGIN
    WITH
      nV = dr.D.nV,
      pos = SUBARRAY(dr.pos^,0,nV),
      bb = Geo.GetBounds(pos),
      xlo = FLOAT(bb[0].lo - 1, LONG),
      xhi = FLOAT(bb[0].hi + 1, LONG),
      ylo = FLOAT(bb[1].lo - 1, LONG),
      yhi = FLOAT(bb[1].hi + 1, LONG),
      dx = (xhi - xlo),
      dy = (yhi - ylo),
      fName = name & ".eps",
      xSize = 10.0d0*scale*dx,  (* Nominal EPS figure width (mm) *)
      ySize = 10.0d0*scale*dy   (* Nominal EPS figure height (mm) *)
    DO
      <* ASSERT xSize < 1000.0d0 *>
      <* ASSERT ySize < 1000.0d0 *>
      eps := NEW(PSPlot.EPSFile).open(fName, xSize, ySize);
      Wr.PutText(stderr, "writing " & fName & " ...");
      eps.setScale(PSPlot.Axis.X, +xlo, +xhi);
      eps.setScale(PSPlot.Axis.Y, -yhi, -ylo);
      Plot(eps, dr, label, labelPad, fontSize, scale);
      eps.close()
    END;
    Wr.PutText(stderr, "\n");
  END PlotEPS;

PROCEDURE Plot(
    eps: PSPlot.File; 
    dr: T;
    label: REF ARRAY OF TEXT;
    labelPad: NAT;
    fontSize: REAL;
    scale: LONG;
  ) =

  VAR 
    (* Scale factor as "REAL" for convenience: *)
    sc: REAL := FLOAT(scale);
    
    (* Charbox dimensions (grid units). *)
    charHX := 0.0125 * fontSize;   (* Charbox half-width. *)
    charHY := 0.0250 * fontSize;   (* Charbox half-height. *)
  
  CONST 
    DY = 0.0d0;         (* Vertical offset for edge-base node attachment point. *)

    DoBaseNodeDots = FALSE;
  
  PROCEDURE PlotEdge(p: Point; pBaseNode: BOOL; q: Point; qBaseNode: BOOL) =
    BEGIN
      WITH
        dp = + DY * FLOAT(ORD(pBaseNode), LONG),
        dq = - DY * FLOAT(ORD(qBaseNode), LONG),
        xp = FLOAT(p[0], LONG), yp = - (FLOAT(p[1], LONG) + dp),
        xq = FLOAT(q[0], LONG), yq = - (FLOAT(q[1], LONG) + dq)
      DO
        IF yp < yq THEN
          eps.setLineDashed(ARRAY OF REAL{2.0*sc, 1.0*sc}, skip := 1.0*sc);
          eps.setLineColor(PSPlot.Red)
        ELSE
          eps.setLineSolid();
          eps.setLineColor(PSPlot.Black)
        END;
        eps.segment(xp, yp, xq, yq)
      END
    END PlotEdge;

  PROCEDURE PlotBaseNodeDots(p: Point; inDot, otDot: BOOL) =
    BEGIN
      WITH
        xp = FLOAT(p[0], LONG),  yp = - FLOAT(p[1], LONG)
      DO
        IF inDot THEN eps.circle(xp, yp + DY, sc*0.6) END;
        IF otDot THEN eps.circle(xp, yp - DY, sc*0.6) END;
      END;
    END PlotBaseNodeDots;

  PROCEDURE PlotBaseNode(v: BaseNode; p: Point) =
    VAR txt: TEXT;
    BEGIN
      IF label # NIL THEN txt := label[v] ELSE txt := Fmt.Int(v) END;
      WITH
        xp = FLOAT(p[0], LONG),  yp = - FLOAT(p[1], LONG),
        HX = FLOAT(charHX * FLOAT(MAX(labelPad, Text.Length(txt))), LONG),
        HY = FLOAT(charHY, LONG)
      DO
        (* eps.rectangle(xp - HX, xp + HX, yp - VY, yp + VY); *)
        eps.ellipse(xp, yp, HX, 0.0d0, 0.0d0, HY);
        eps.label(txt, xp, yp, 0.5, 0.5, size := sc*fontSize);
      END;
    END PlotBaseNode;

  BEGIN
    WITH 
      G = dr.G,
      nBaseV = G.nV,
      D = dr.D,
      nV = D.nV,
      pos = SUBARRAY(dr.pos^, 0, nV)
    DO

      (* Plot edges: *)
      eps.setLineWidth(sc*0.20);
      eps.setLineSolid();
      FOR u := 0 TO nV-1 DO
        WITH deg = dr.D.degree(u,+1) DO
          FOR k := 0 TO deg-1 DO
            WITH v = dr.D.nbor(u,+1,k) DO
              IF v # NoNode THEN
                PlotEdge(pos[u], u < nBaseV, pos[v], v < nBaseV)
              END
            END
          END
        END
      END;

      IF DoBaseNodeDots THEN
        (* Plot base node attachment dots: *)
        eps.setLineColor(PSPlot.Invisible);
        eps.setFillColor(PSPlot.Black);
        FOR v := 0 TO nBaseV-1 DO
          PlotBaseNodeDots(pos[v], G.degree(v,-1) > 0, G.degree(v,+1) > 0)
        END;
      END;

      (* Plot vertices: *)
      eps.setLineSolid();
      eps.setLineColor(PSPlot.Black);
      eps.setLineWidth(sc*0.20);
      eps.setFillColor(Color{1.000, 1.000, 0.750});
      eps.setTextColor(PSPlot.Blue);
      FOR v := 0 TO nBaseV-1 DO
        PlotBaseNode(v, pos[v])
      END;
    END    
  END Plot;

PROCEDURE NeighborMean(
    u: Node;                (* A node of "H". *)
    H: Graph.T;           (* A graph. *)
    READONLY pos: Points;   (* Coordinates for the nodes of H. *)
    default: LONG;          (* Default coordinate. *)
    wt: REF Weights;        (* Weights for the nodes of "H". *)
    axis: Axis;             (* Coordinate to average. *)
  ): LONG = 
  
  VAR wxSum: INT := 0; wSum: INT := 0; wv: INT;
  BEGIN
    FOR dir := -1 TO +1 BY 2 DO
      WITH deg = H.degree(u, dir) DO
        FOR k := 0 TO deg-1 DO
          WITH v = H.nbor(u,dir,k), xv = pos[v][axis] + 0 DO
            IF v # NoNode THEN 
              IF wt = NIL THEN wv := 1 ELSE wv := wt[v] END;
              wxSum := wxSum + xv*wv; wSum := wSum + wv
            END
          END
        END
      END
    END;
    IF wSum = 0 THEN 
      RETURN default
    ELSE
      RETURN FLOAT(wxSum, LONG)/FLOAT(wSum, LONG)
    END;
  END NeighborMean;

PROCEDURE VerifyDegreeInvariants(D: Graph.T; nBaseV: NAT) =
  BEGIN
    (* Verify that all vertices beyond the first "nBaseV" are joints: *)
    <* ASSERT D.nV >= nBaseV *>
    FOR u := nBaseV TO D.nV-1 DO
      <* ASSERT D.degrees(u) = Degrees{1,0,1} *>
    END;
  END VerifyDegreeInvariants;  
  
PROCEDURE VerifySlotCorrespondence(G, D: Graph.T) =
  VAR v, w: Node;
  BEGIN
    WITH nBaseV = G.nV DO
      FOR u := 0 TO nBaseV-1 DO
        FOR dir := -1 TO +1 BY 2 DO 
          WITH deg = G.degree(u, dir) DO
            FOR ku := 0 TO deg-1 DO
              v := D.nbor(u,dir,ku);
              w := G.nbor(u,dir,ku);
              IF v # NoNode THEN
                WHILE v # NoNode AND v >= nBaseV DO v := D.nbor(v,dir,0) END;
                <* ASSERT v = w *>
              END
            END
          END
        END
      END
    END
  END VerifySlotCorrespondence;

PROCEDURE VerifyDistinctPositions(READONLY pos: Points) =
  BEGIN
    WITH
      uu = NEW(REF Nodes, NUMBER(pos))^
    DO
      FOR i := 0 TO LAST(uu) DO uu[i] := i END;
      SortNodesByCoords(uu, pos);
      FOR i := 1 TO LAST(uu) DO
        <* ASSERT pos[uu[i]] # pos[uu[i-1]] *>
      END;
    END
  END VerifyDistinctPositions;

<*UNUSED*>
PROCEDURE FI(x: INT): TEXT = 
  BEGIN RETURN Fmt.Int(x) END FI;

<*UNUSED*>
PROCEDURE FPt(p: Point): TEXT = 
  BEGIN RETURN "(" & Fmt.Int(p[0]) & "," & Fmt.Int(p[1]) & ")" END FPt;

BEGIN
END GDDrawing.
(* Last edited on 2000-01-13 10:40:07 by stolfi *)