#line 311 "/home/travis/build/felix-lang/felix/src/packages/trees.fdoc"
  // Directed Cyclic graph
  
  include "std/datatype/dlist";
  include "std/datatype/partition";
  
  class DiGraph[V,E with Str[V], Str[E]] // V,E labels for graph parts
  {
    // vertices are stored in an array, so they're identified
    // by their slot number 0 origin
    typedef digraph_t = (vertices: darray[vertex_t], nedges: int);
    ctor digraph_t () => (vertices= #darray[vertex_t], nedges=0);
  
    // x index implicit, the edge source
    // y index is the edge destination
    typedef edge_t = (elabel:E, x:int,y:int, weight:double);
    typedef vertex_t = (vlabel:V, outedges: list[edge_t]);
  
    fun len (d:digraph_t) => d.vertices.len;
  
    virtual fun default_vlabel: 1 -> V;
    virtual fun default_elabel: 1 -> E;
    fun default_vertex () => (vlabel = #default_vlabel, outedges = Empty[edge_t]);
  
    // Add an isolated vertex
    // If the vertex is already in the graph,
    // this routine just replaces the label
    // this allows adding out of order vertices
    // and adding vertices implicitly by adding edges
    proc add_vertex (d:&digraph_t, v:V, x:int)
    {
      while x >= d*.vertices.len.int call push_back (d*.vertices, #default_vertex);
      var pv: &V = (d*.vertices,x.size).unsafe_get_ref.vlabel;
      pv <- v;
    }
  
    proc add_weighted_edge (d:&digraph_t, x:int, y:int, elab:E, weight:double)
    {
      while x >= d*.vertices.len.int call add_vertex (d,#default_vlabel,d*.vertices.len.int);
      while y >= d*.vertices.len.int call add_vertex (d,#default_vlabel,d*.vertices.len.int);
      var pedges : &list[edge_t] = (d*.vertices,x.size).unsafe_get_ref.outedges;
      pedges <- (elabel=elab,x=x,y=y,weight=weight) ! *pedges;
      d.nedges.pre_incr;
    }
  
    proc add_edge (d:&digraph_t, x:int, y:int, elab:E) =>
      add_weighted_edge (d,x,y,elab,1.0)
    ;
  
    // add and edge and its reverse edge, distinct labels
    proc add_weighted_edge_pair (d:&digraph_t, x:int, y:int, felab:E, relab:E, weight:double)
    {
      add_weighted_edge(d,x,y,felab, weight);
      add_weighted_edge(d,y,x,relab, weight);
    }
  
    proc add_edge_pair (d:&digraph_t, x:int, y:int, felab:E, relab:E) =>
      add_weighted_edge_pair (d,x,y,felab,relab,1.0)
    ;
  
    // add and edge and its reverse edge, same label
    // use for undirected graph
    proc add_edge_pair (d:&digraph_t, x:int, y:int, elab:E)
    {
      add_edge(d,x,y,elab);
      add_edge(d,y,x,elab);
    }
  
  
    fun dump_digraph (d:digraph_t) : string =
    {
      var out = "";
      reserve (&out,10000);
      var x = 0;
      for vertex in d.vertices do
        out += x.str + " " + vertex.vlabel.str + "\n";
        for edge in vertex.outedges do
          out += "  " + edge.x.str + "->" + edge.y.str + " " +
            edge.elabel.str +
            if edge.weight != 1.0 then " "+edge.weight.str else "" endif +
            "\n"
          ;
        done
      ++x;
      done
      return out;
    }
  
    union Vstate = Undiscovered | Discovered | Processed;
  
    typedef digraph_visitor_processing_t =
    (
      process_vertex_early: digraph_t -> int -> 0,
      process_vertex_late: digraph_t -> int -> 0,
      process_edge: digraph_t -> int * int -> 0
    );
  
    proc dflt_pve (g:digraph_t) (x:int) {};
    proc dflt_pvl (g:digraph_t) (x:int) {};
    proc dflt_pe (g:digraph_t) (x:int, y:int) {};
  
    // default visitor does nothing
    ctor digraph_visitor_processing_t () => (
      process_vertex_early= dflt_pve,
      process_vertex_late= dflt_pvl,
      process_edge= dflt_pe
    );
  
    interface mutable_collection_t[T] {
       add: T -> 0;
       remove: 1 -> opt[T];
    }
  
    gen iterator[T] (x:mutable_collection_t[T]) () : opt[T] => x.remove ();
  
    object gstack_t[T] () implements mutable_collection_t[T] = {
      open DList[T];
      var d = dlist_t();
      method proc add (x:T) => push_back (&d,x);
      method gen remove () => pop_back (&d);
    }
  
    object gqueue_t[T] () implements mutable_collection_t[T] = {
      open DList[T];
      var d = dlist_t();
      method proc add (x:T) => push_back (&d,x);
      method gen remove () => pop_front (&d);
    }
  
    proc iter
      (var pending:mutable_collection_t[int])
      (d:digraph_t) (startv:int)
      (p:digraph_visitor_processing_t)
    {
      var state = varray[Vstate] (bound=d.len,default=Undiscovered);
      pending.add startv;
      set (state,startv,Discovered);
      //var parent = -1;
      for v in pending do // all vertex indices in queue
        p.process_vertex_early d v;
        set (state,v,Processed);
        for edge in d.vertices.v.outedges do
          var y = edge.y;
          p.process_edge d (v, y);
          match state.y do
          | #Undiscovered =>
            pending.add y;
            set (state,y,Discovered);
            //parent = v;
          | _ => ;
          done
        done
        p.process_vertex_late d v;
      done // vertices
    }
  
    proc breadth_first_iter (d:digraph_t) (startv:int) (p:digraph_visitor_processing_t) =>
      iter #gqueue_t[int] d startv p
    ;
  
    proc depth_first_iter (d:digraph_t) (startv:int) (p:digraph_visitor_processing_t) =>
      iter #gstack_t[int] d startv p
    ;
  
    // This routine returns a list of vertices from startv to fin, inclusive ..
    // not a list of edges.
    gen find_shortest_unweighted_path (d:digraph_t) (startv:int, fin:int) : opt[list[int]] =
    {
      if startv == fin return Some (list(startv));
  
      open DList[int];
      var state = varray[Vstate] (bound=d.len,default=Undiscovered);
      var parents = varray[int] (bound=d.len,default= -1);
      var q = queue_t();
      enqueue &q startv;
      set (state,startv,Discovered);
      set(parents,startv,-1);
      for v in &q // all vertex indices in queue
        for edge in d.vertices.v.outedges do
          var y = edge.y;
          if y == fin do
            var path = Empty[int];
            set(parents,y,v);
            while y != startv do
              path = Cons (y,path);
              y = parents.y;
            done
            path = Cons (y,path);
            return Some path;
          else
            match state.y do
            | #Undiscovered =>
              enqueue &q y;
              set (state,y,Discovered);
              set(parents,y,v);
            | _ => ;
            done
          done
        done
      return None[list[int]];
    }
  
    // find minimum spanning tree
    // Prim's algorithm, enhanced as in Skiena
    // only returns list of vertices from starting point
    gen prim (d:digraph_t) (startv:int) : list[int * int] =
    {
      var INF=DINFINITY;
      var intree = varray[bool] (bound=d.len, default=false);
      var distance = varray[double] (bound=d.len, default=INF);
      var fromv = varray[int] (bound=d.len, default= -1);
      var span = Empty[int * int];
      var src = -1;
      var v = startv;
      while not intree.v do
        set(intree,v,true);
        for edge in d.vertices.v.outedges do
          var w = edge.y;
          var weight = edge.weight;
          if distance.w > weight and not intree.w do
            set(distance,w,weight);
            set(fromv,w,v);
          done
        done
  
        // find closest out of tree vertex
        var dist = INF;
        src = -1;
        for var i in 0 upto intree.len.int - 1 do
          if not intree.i and dist > distance.i do
            dist = distance.i;
            v = i;
            src = fromv.i;
          done // not in tree
        done // each vertex i
        // v is set to closest out of tree vertex and
        // src to the vertex it comes from
        // if there is one, otherwise v is unchanged and so remains in tree
        // and src stays at -1
        if src != -1 do span = Cons ( (src,v), span); done
      done // each v not in tree
      return rev span;
    }
  
  }
  
  instance DiGraph[string, string]
  {
    fun default_vlabel () => "Unlabelled Vertex";
    fun default_elabel () => "Unlabelled Edge";
  }