Anforderungen  |   Konzepte  |   Entwurf  |   Entwicklung  |   Qualitätssicherung  |   Lebenszyklus  |   Steuerung
 
 
 
 


SSL attr.gi   Sprache: unbekannt

 
#############################################################################
##
##  attr.gi
##  Copyright (C) 2014-21                                James D. Mitchell
##
##  Licensing information can be found in the README file of this package.
##
#############################################################################
##

InstallMethod(DigraphNrVertices, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep], DIGRAPH_NR_VERTICES);

InstallGlobalFunction(OutNeighbors, OutNeighbours);

# The next method is (yet another) DFS which simultaneously computes:
# 1. *articulation points* as described in
#   https://www.eecs.wsu.edu/~holder/courses/CptS223/spr08/slides/graphapps.pdf
# 2. *bridges* as described in https://stackoverflow.com/q/28917290/
#   (this is a minor adaption of the algorithm described in point 1).
# 3. a *strong orientation* as alluded to somewhere on the internet that I can
#    no longer find. It's essentially just "orient every edge in the DFS tree
#    away from the root, and every other edge (back edges) from the node with
#    higher `pre` value to the one with lower `pre` value (i.e. they point
#    backwards from later nodes in the DFS to earlier ones). If the graph is
#    bridgeless, then it is guaranteed that the orientation of the last
#    sentence is strongly connected."

BindGlobal("DIGRAPHS_ArticulationPointsBridgesStrongOrientation",
function(D)
  local N, copy, articulation_points, bridges, orientation, nbs, counter, pre,
        low, nr_children, stack, u, v, i, w, connected;

  N := DigraphNrVertices(D);

  if HasIsConnectedDigraph(D) and not IsConnectedDigraph(D) then
    # not connected, no articulation points, no bridges, no strong orientation
    return [false, [], [], fail];
  elif N < 2 then
    # connected, no articulation points (removing 0 or 1 nodes does not make
    # the graph disconnected), no bridges, strong orientation (since
    # the digraph with 0 nodes is strongly connected).
    return [true, [], [], D];
  elif not IsSymmetricDigraph(D) then
    copy := DigraphSymmetricClosure(DigraphMutableCopyIfMutable(D));
    MakeImmutable(copy);
  else
    copy := D;
  fi;

  # outputs
  articulation_points := [];
  bridges             := [];
  orientation         := List([1 .. N], x -> BlistList([1 .. N], []));

  # Get out-neighbours once, to avoid repeated copying for mutable digraphs.
  nbs := OutNeighbours(copy);

  # number of nodes encountered in the search so far
  counter := 0;

  # the order in which the nodes are visited, -1 indicates "not yet visited".
  pre := ListWithIdenticalEntries(N, -1);

  # low[i] is the lowest value in pre currently reachable from node i.
  low := [];

  # nr_children of node 1, for articulation points the root node (1) is an
  # articulation point if and only if it has at least 2 children.
  nr_children := 0;

  stack := Stack();
  u := 1;
  v := 1;
  i := 0;

  repeat
    if pre[v] <> -1 then
      # backtracking
      i := Pop(stack);
      v := Pop(stack);
      u := Pop(stack);
      w := nbs[v][i];

      if v <> 1 and low[w] >= pre[v] then
        AddSet(articulation_points, v);
      fi;
      if low[w] = pre[w] then
        Add(bridges, [v, w]);
      fi;
      if low[w] < low[v] then
        low[v] := low[w];
      fi;
    else
      # diving - part 1
      counter := counter + 1;
      pre[v]  := counter;
      low[v]  := counter;
    fi;
    i := PositionProperty(nbs[v], w -> w <> v, i);
    while i <> fail do
      w := nbs[v][i];
      if pre[w] <> -1 then
        # v -> w is a back edge
        if w <> u and pre[w] < low[v] then
          low[v] := pre[w];
        fi;
        orientation[v][w] := not orientation[w][v];
        i := PositionProperty(nbs[v], w -> w <> v, i);
      else
        # diving - part 2
        if v = 1 then
          nr_children := nr_children + 1;
        fi;
        orientation[v][w] := true;
        Push(stack, u);
        Push(stack, v);
        Push(stack, i);
        u := v;
        v := w;
        i := 0;
        break;
      fi;
    od;
  until Size(stack) = 0;

  if counter = DigraphNrVertices(D) then
    connected := true;
    if nr_children > 1 then
      # The `AddSet` is not really needed, and could just be `Add`, since 1
      # should not already be in articulation_points. But let's use AddSet
      # to keep the output sorted.
      AddSet(articulation_points, 1);
    fi;
    if not IsEmpty(bridges) then
      orientation := fail;
    else
      orientation := DigraphByAdjacencyMatrix(DigraphMutabilityFilter(D),
                                              orientation);
    fi;
  else
    connected           := false;
    articulation_points := [];
    bridges             := [];
    orientation         := fail;
  fi;
  if IsImmutableDigraph(D) then
    SetIsConnectedDigraph(D, connected);
    SetArticulationPoints(D, articulation_points);
    SetBridges(D, bridges);
    if IsSymmetricDigraph(D) then
      SetStrongOrientationAttr(D, orientation);
    fi;
  fi;
  return [connected, articulation_points, bridges, orientation];
end);

InstallMethod(ArticulationPoints, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
D -> DIGRAPHS_ArticulationPointsBridgesStrongOrientation(D)[2]);

InstallMethod(Bridges, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
D -> DIGRAPHS_ArticulationPointsBridgesStrongOrientation(D)[3]);

InstallMethodThatReturnsDigraph(StrongOrientation,
"for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  if not IsSymmetricDigraph(D) then
    ErrorNoReturn("not yet implemented");
  fi;
  return DIGRAPHS_ArticulationPointsBridgesStrongOrientation(D)[4];
end);

# Utility function to calculate the maximal independent sets (as BLists) of a
# subgraph induced by removing a set of vertices.
BindGlobal("DIGRAPHS_MaximalIndependentSetsSubtractedSet",
function(I, subtracted_set, size_bound)
  local induced_mis, temp, i;
  # First remove all vertices in the subtracted set from each MIS
  temp := List(I, i -> DifferenceBlist(i, subtracted_set));
  induced_mis := [];
  # Then remove any sets which are no longer maximal
  # Sort in decreasing size
  Sort(temp, {x, y} -> SizeBlist(x) > SizeBlist(y));
  # Then check elements from back to front for if they are a subset
  for i in temp do
    if SizeBlist(i) <= size_bound and
        ForAll(induced_mis, x -> not IsSubsetBlist(x, i)) then
      Add(induced_mis, i);
    fi;
  od;
  return induced_mis;
end
);

InstallMethod(DIGRAPHS_AbsorbingMarkovChain,
"for a digraph",
[IsDigraph],
function(D)
  # Helper for DigraphAbsorptionProbabilities and DigraphAbsorptionExpectedSteps
  local scc, is_sink_comp, transient_vertices, i, comp, v, sink_comps,
        nr_transients, transient_mat, absorption_mat, neighbours, chance, w,
        w_comp, sink_comp_index, j, fundamental_mat;
  scc := DigraphStronglyConnectedComponents(D);

  # Find the "sink components" (components from which there is no escape)
  # We could use QuotientDigraph and DigraphSinks here, but this avoids copying.
  is_sink_comp := ListWithIdenticalEntries(Length(scc.comps), true);
  transient_vertices := [];
  for i in [1 .. Length(scc.comps)] do
    comp := scc.comps[i];
    for v in comp do
      if ForAny(OutNeighboursOfVertex(D, v), w -> not w in comp) then
        is_sink_comp[i] := false;
        Append(transient_vertices, comp);
        break;
      fi;
    od;
  od;
  sink_comps := Positions(is_sink_comp, true);
  nr_transients := Length(transient_vertices);

  # If no transient vertices, then we can't return anything interesting.
  if nr_transients = 0 then
    return rec(transient_vertices := transient_vertices,
               fundamental_mat := [],
               sink_comps := sink_comps,
               absorption_mat := []);
  fi;

  # transient_mat[i][j] is the chance of going
  # from transient vertex i to transient vertex j
  transient_mat := NullMat(nr_transients, nr_transients);

  # absorption_mat[i][j] is the chance of going
  # from transient vertex i to sink component j
  absorption_mat := NullMat(nr_transients, Length(sink_comps));

  for i in [1 .. Length(transient_vertices)] do
    v := transient_vertices[i];
    neighbours := OutNeighboursOfVertex(D, v);
    chance := 1 / Length(neighbours);  # chance of following each edge
    for w in neighbours do
      w_comp := scc.id[w];
      if is_sink_comp[w_comp] then
        # w is in a sink component
        sink_comp_index := Position(sink_comps, w_comp);
        Assert(1, w in scc.comps[sink_comps[sink_comp_index]]);
        absorption_mat[i][sink_comp_index] :=
            absorption_mat[i][sink_comp_index] + chance;
      else
        # w is a transient vertex
        j := Position(transient_vertices, w);
        transient_mat[i][j] := transient_mat[i][j] + chance;
      fi;
    od;
  od;

  # "Fundamental matrix": mat[i][j] is the expected number of visits to each
  # transient vertex j given a random walk starting at transient vertex i.
  # Compute this using formula (I - Q)^-1
  fundamental_mat := Inverse(IdentityMat(nr_transients) - transient_mat);

  # Return all info needed for further computation in DigraphAbsorption* methods
  return rec(transient_vertices := transient_vertices,
             fundamental_mat := fundamental_mat,
             sink_comps := sink_comps,
             absorption_mat := absorption_mat);
end);

InstallMethod(DigraphAbsorptionProbabilities,
"for a digraph",
[IsDigraph],
function(D)
  local scc, markov, chances_of_absorption, output, sink_comps, comp_no, v,
        transient_vertices, i, j, c;
  # Strongly connected components are an important part of definition
  scc := DigraphStronglyConnectedComponents(D);

  # One SCC only (or none): all vertices stay in it with probability 1.
  if Length(scc.comps) <= 1 then
    return ListWithIdenticalEntries(DigraphNrVertices(D), [1]);
  fi;

  # Get data about absorbing Markov chain represented by this digraph
  markov := DIGRAPHS_AbsorbingMarkovChain(D);

  # Calculate chances of absorption
  # Rows are transient vertices, columns are absorbing SCCs
  if Length(markov.transient_vertices) = 0 then
    chances_of_absorption := [];
  else
    chances_of_absorption := markov.fundamental_mat * markov.absorption_mat;
  fi;

  # Convert to output format
  # Rows are all vertices, columns are all SCCs
  output := NullMat(DigraphNrVertices(D), Length(scc.comps));

  # Non-transient vertices stay in their own SCC with probability 1
  sink_comps := markov.sink_comps;
  for comp_no in sink_comps do
    for v in scc.comps[comp_no] do
      output[v][comp_no] := 1;
    od;
  od;

  # Transient vertices have chances as calculated in Markov code
  transient_vertices := markov.transient_vertices;
  for i in [1 .. Length(transient_vertices)] do
    v := transient_vertices[i];
    for j in [1 .. Length(sink_comps)] do
      c := sink_comps[j];
      output[v][c] := chances_of_absorption[i][j];
    od;
  od;

  return output;
end);

InstallMethod(DigraphAbsorptionExpectedSteps,
"for a digraph",
[IsDigraph],
function(D)
  local markov, N, transient_expected_steps, out, i;
  if DigraphNrVertices(D) = 0 then
    return [];
  fi;
  markov := DIGRAPHS_AbsorbingMarkovChain(D);
  N := markov.fundamental_mat;
  if Length(N) = 0 then  # no transient vertices
    transient_expected_steps := [];
  else
    transient_expected_steps := N * ListWithIdenticalEntries(Length(N), 1);
  fi;
  out := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
  for i in [1 .. Length(transient_expected_steps)] do
    out[markov.transient_vertices[i]] := transient_expected_steps[i];
  od;
  return out;
end);

BindGlobal("DIGRAPHS_ChromaticNumberLawler",
function(D)
  local n, vertices, subset_colours, s, S, i, I, subset_iter, x,
  mis, subset_mis;
  n := DigraphNrVertices(D);
  vertices := List(DigraphVertices(D));
  # Store all the Maximal Independent Sets, which can later be used for
  # calculating the maximal independent sets of induced subgraphs.
  mis := DigraphMaximalIndependentSets(D);
  # Convert each MIS to a Blist
  mis := List(mis, x -> BlistList(vertices, x));
  # Store current best colouring for each subset
  subset_colours := ListWithIdenticalEntries(2 ^ n, infinity);
  # Empty set can be colouring with only one colour.
  subset_colours[1] := 0;
  # Iterator for blist subsets.
  subset_iter := ListWithIdenticalEntries(n, [false, true]);
  subset_iter := IteratorOfCartesianProduct2(subset_iter);
  # Skip the first one, which should be the empty set.
  S := NextIterator(subset_iter);
  # Iterate over all vertex subsets.
  for S in subset_iter do
    # Cartesian iterator is ascending lexicographically, but we want reverse
    # lexicographic ordering. We treat this as iterating over the complement.
    # Index the current subset that is being iterated over.
    s := 1;
    for x in [1 .. n] do
      # Need to negate, as we are iterating over the complement.
      if not S[x] then
        s := s + 2 ^ (x - 1);
      fi;
    od;
    # Iterate over the maximal independent sets of V[S]
    subset_mis :=
      DIGRAPHS_MaximalIndependentSetsSubtractedSet(mis, S, infinity);
    # Flip the list, as we now need the actual set.
    FlipBlist(S);
    for I in subset_mis do
        # Calculate S \ I. This is destructive, but is undone.
        SubtractBlist(S, I);
        # Index S \ I
        i := 1;
        for x in [1 .. n] do
          if S[x] then
            i := i + 2 ^ (x - 1);
          fi;
        od;
        # The chromatic number of this subset is the minimum value of all
        # the maximal independent subsets of D[S].
        subset_colours[s] := Minimum(subset_colours[s], subset_colours[i] + 1);
        # Undo the changes to the subset.
        UniteBlist(S, I);
    od;
  od;
  return subset_colours[2 ^ n];
end
);

BindGlobal("DIGRAPHS_UnderThreeColourable",
function(D)
  local nr;
  nr := DigraphNrVertices(D);
  if DigraphHasLoops(D) then
    ErrorNoReturn("the argument <D> must be a digraph with no loops,");
  elif nr = 0 then
    return 0;  # chromatic number = 0 iff <D> has 0 verts
  elif IsNullDigraph(D) then
    return 1;  # chromatic number = 1 iff <D> has >= 1 verts & no edges
  elif IsBipartiteDigraph(D) then
    return 2;  # chromatic number = 2 iff <D> has >= 2 verts & is bipartite
               # <D> has at least 2 vertices at this stage
  elif DigraphColouring(D, 3) <> fail then
               # Check if there is a 3 colouring
    return 3;
  fi;
  return infinity;
end
);

BindGlobal("DIGRAPHS_ChromaticNumberByskov",
function(D)
  local n, a, vertices, subset_colours, S, i, j, I, subset_iter,
  index_subsets, vertex_blist, k, MIS;

  n := DigraphNrVertices(D);
  vertices := DigraphVertices(D);
  vertex_blist := BlistList(vertices, vertices);
  # Store all the Maximal Independent Sets, which can later be used for
  # calculating the maximal independent sets of induced subgraphs.
  MIS := DigraphMaximalIndependentSets(D);
  # Convert each MIS to a Blist
  MIS := List(MIS, x -> BlistList(vertices, x));
  # Store current best colouring for each subset
  subset_colours := ListWithIdenticalEntries(2 ^ n, infinity);
  # Empty set is 0 colourable
  subset_colours[1] := 0;
  # Function to index the subsets of the vertices of D
  index_subsets := function(subset)
    local x, index;
    index := 1;
    for x in [1 .. n] do
      if subset[x] then
        index := index + 2 ^ (x - 1);
      fi;
    od;
    return index;
  end;
  # Iterate over vertex subsets
  subset_iter := IteratorOfCombinations(vertices);
  # Skip the first one, which should be the empty set
  S := NextIterator(subset_iter);
  Assert(1, IsEmpty(S), "First set from iterator should be the empty set");
  # First find the 3 colourable subgraphs of D
  for S in subset_iter do
    a := DIGRAPHS_UnderThreeColourable(InducedSubdigraph(D, S));
    S := BlistList(vertices, S);
    i := index_subsets(S);
    # Mark this as three or less colourable if it is.
    subset_colours[i] := Minimum(a, subset_colours[i]);
  od;
  # Process 4 colourable subgraphs
  for I in MIS do
    SubtractBlist(vertex_blist, I);
    # Iterate over all subsets of V(D) \ I as blists
    # This is done by taking the cartesian product of n copies of [true, false]
    # or [true] if the vertex is in I. The [true] is used as each element will
    # be flipped to get reverse lexicographic ordering.
    subset_iter := EmptyPlist(n);
    for i in [1 .. n] do
      if I[i] then
        subset_iter[i] := [true];
      else
        subset_iter[i] := [true, false];
      fi;
    od;
    subset_iter := IteratorOfCartesianProduct2(subset_iter);
    # Skip the empty set.
    NextIterator(subset_iter);
    for S in subset_iter do
      FlipBlist(S);
      i := index_subsets(S);
      if subset_colours[i] = 3 then
        # Index union of S and I
        j := index_subsets(UnionBlist(S, I));
        subset_colours[j] := Minimum(subset_colours[j], 4);
      fi;
    od;
    # Undo the changes made.
    UniteBlist(vertex_blist, I);
  od;
  # Iterate over vertex subset blists.
  subset_iter := ListWithIdenticalEntries(n, [true, false]);
  subset_iter := IteratorOfCartesianProduct2(subset_iter);
  # Skip the first one, which should be the empty set
  S := NextIterator(subset_iter);
  for S in subset_iter do
    # Cartesian iteratator goes in lexicographic order, but we want reverse.
    FlipBlist(S);
    # Index the current subset that is being iterated over
    i := index_subsets(S);
    if 4 <= subset_colours[i] and subset_colours[i] < infinity then
      k := SizeBlist(S) / subset_colours[i];
      # Iterate over the maximal independent sets of D[V \ S]
      for I in DIGRAPHS_MaximalIndependentSetsSubtractedSet(MIS, S, k) do
        # Index S union I
        j := index_subsets(UnionBlist(S, I));
        subset_colours[j] := Minimum(subset_colours[j], subset_colours[i] + 1);
      od;
    fi;
  od;
  return subset_colours[2 ^ n];
end
);

BindGlobal("DIGRAPHS_ChromaticNumberZykov",
function(D)
  local nr, ZykovReduce, chrom;
  nr := DigraphNrVertices(D);
  # Recursive function call
  ZykovReduce := function(D)
    local nr, D_contract, vertices, v, x, y, i, j, adjacent;
    nr := DigraphNrVertices(D);
    # Update upper bound if possible.
    chrom := Minimum(nr, chrom);
    # Leaf nodes are either complete graphs or cliques that have size equal to
    # the current upper bound. The chromatic number is then the smallest clique
    # found.
    # Cliques finder arguments:
    # digraph = D - The graph
    # hook = fail - hook is not required
    # user_param = [] - user_param is a list as hook is fail
    # limit = 1 - We only need one clique
    # include = exclude = [] - We check all vertices
    # max = false - This clique need not be maximal
    # size = chrom - We want a clique the size of our upper bound
    # reps = true - As we only care about the existence of the clique,
    # we can instead search for representatives which is more efficient.
    if not IsCompleteDigraph(D) and IsEmpty(CliquesFinder(D, fail, [], 1, [],
                                                          [], false, chrom,
                                                          true)) then
      # Sort vertices by degree, so that higher degree vertices are picked first
      # Picking higher degree vertices will make it more likely a clique will
      # form in one of the modified graphs, which will terminate the recursion.
      vertices := DigraphWelshPowellOrder(D);
      # Get the adjacency function
      adjacent := DigraphAdjacencyFunction(D);
      # Choose two non-adjacent vertices x, y
      for i in [1 .. nr] do
        x := vertices[i];
        # Search for the first vertex not adjacent to all others
        # This is guaranteed to exist as D is not the complete graph.
        if OutDegreeOfVertex(D, x) < nr - 1 then
          # Now search for a non-adjacent vertex, prioritising higher degree
          # ones
          for j in [i + 1 .. nr] do
            y := vertices[j];
            if not adjacent(x, y) then
              break;
            fi;
          od;
          break;
        fi;
      od;
      Assert(1, x <> y, "x and y must be different");
      # Colour the vertex contraction.
      # A contraction of a graph effectively merges two non adjacent vertices
      # into a single new vertex with the edges merged.
      # We merge y into x, keeping x.

      # We could potentially use quotient digraph here, but the increased
      # generality might cause this to get slower.
      D_contract := DigraphMutableCopy(D);
      for v in vertices do
         # Iterate over all vertices that are not x or y
         if v = x or v = y then
           continue;
         fi;
         # Add any edge that involves y, but not already x to avoid duplication.
         if adjacent(v, y) and not adjacent(v, x) then
            DigraphAddEdge(D_contract, x, v);
            DigraphAddEdge(D_contract, v, x);
         fi;
      od;
      DigraphRemoveVertex(D_contract, y);
      ZykovReduce(D_contract);
      # Colour the edge addition
      # This just adds symmetric edges between x and y;
      DigraphAddEdge(D, [x, y]);
      DigraphAddEdge(D, [y, x]);
      ZykovReduce(D);
      # Undo changes to the graph
      DigraphRemoveEdge(D, [x, y]);
      DigraphRemoveEdge(D, [y, x]);
    fi;
  end;
  # Algorithm requires an undirected graph without multiple edges.
  D := DigraphMutableCopy(D);
  D := DigraphRemoveAllMultipleEdges(D);
  D := DigraphSymmetricClosure(D);
  # Use greedy colouring as an upper bound
  chrom := RankOfTransformation(DigraphGreedyColouring(D), nr);
  ZykovReduce(D);
  return chrom;
end
);

BindGlobal("DIGRAPHS_ChromaticNumberChristofides",
function(D)
  local nr, I, n, T, b, unprocessed, i, v_without_t, j, u, min_occurrences,
  cur_occurrences, chrom, colouring, stack, vertices;

  nr := DigraphNrVertices(D);
  vertices := List(DigraphVertices(D));
  # Initialise the required variables.
  # Calculate all maximal independent sets of D.
  I := DigraphMaximalIndependentSets(D);
  # Convert each MIS into a BList
  I := List(I, i -> BlistList(vertices, i));
  # Upper bound for chromatic number.
  chrom := nr;
  # Set of vertices of D not in the current subgraph at level n.
  T := ListWithIdenticalEntries(nr, false);
  # Current search level of the subgraph tree.
  n := 0;
  # The maximal independent sets of V \ T at level n.
  b := [ListWithIdenticalEntries(nr, false)];
  # Number of unprocessed MIS's of V \ T from level 1 to n
  unprocessed := ListWithIdenticalEntries(nr, 0);
  # Would be jth colour class of the chromatic colouring of G.
  colouring := List([1 .. nr], i -> BlistList(vertices, [i]));
  # Stores current unprocessed MIS's of V \ T at level 1 to level n
  stack := [];
  # Now perform the search.
  repeat
    # Step 2
    if n < chrom then
      # Step 3
      # If V = T then we've reached a null subgraph
      if SizeBlist(T) = nr then
        chrom := n;
        SubtractBlist(T, b[n + 1]);
        for i in [1 .. chrom] do
          colouring[i] := b[i];
          # TODO set colouring attribute
        od;
      else
        # Step 4
        # Compute the maximal independent sets of V \ T
        v_without_t := DIGRAPHS_MaximalIndependentSetsSubtractedSet(I, T,
                                                                    infinity);
        # Step 5
        # Pick u in V \ T such that u is in the fewest maximal independent sets.
        u := -1;
        min_occurrences := infinity;
        # Flip T to get V \ T
        FlipBlist(T);
        # Convert to list to iterate over the vertices.
        for i in ListBlist(vertices, T) do
          # Count how many times this vertex appears in a MIS
          cur_occurrences := Number(v_without_t, j -> j[i]);
          if cur_occurrences < min_occurrences then
            min_occurrences := cur_occurrences;
            u := i;
          fi;
        od;
        # Revert changes to T
        FlipBlist(T);
        Assert(1, u <> -1, "Vertex must be picked");
        # Remove maximal independent sets not containing u.
        v_without_t := Filtered(v_without_t, x -> x[u]);
        # Add these MISs to the stack
        Append(stack, v_without_t);
        # Search has moved one level deeper
        n := n + 1;
        unprocessed[n] := Length(v_without_t);
      fi;
    else
      # if n >= g then T = T \ b[n]
      # This exceeds the current best bound, so stop search.
      SubtractBlist(T, b[n + 1]);
    fi;
    # Step 6
    while n <> 0 do
      # step 7
      if unprocessed[n] = 0 then
        n := n - 1;
        SubtractBlist(T, b[n + 1]);
      else
        # Step 8
        # take an element from the top of the stack
        i := Remove(stack);
        unprocessed[n] := unprocessed[n] - 1;
        b[n + 1] := i;
        UniteBlist(T, i);
        break;
      fi;
    od;
  until n = 0;
  return chrom;
end
);

InstallMethod(ChromaticNumber, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local nr, comps, upper, chrom, tmp_comps, tmp_upper, n, comp, bound, clique,
  c, i, greedy_bound, brooks_bound;
  nr := DigraphNrVertices(D);

  if DigraphHasLoops(D) then
    ErrorNoReturn("the argument <D> must be a digraph with no loops,");
  elif nr = 0 then
    return 0;  # chromatic number = 0 iff <D> has 0 verts
  elif IsNullDigraph(D) then
    return 1;  # chromatic number = 1 iff <D> has >= 1 verts & no edges
  elif IsBipartiteDigraph(D) then
    return 2;  # chromatic number = 2 iff <D> has >= 2 verts & is bipartite
               # <D> has at least 2 vertices at this stage
  fi;

  # Value option for alternative dispatch
  if ValueOption("lawler") <> fail then
    return DIGRAPHS_ChromaticNumberLawler(D);
  elif ValueOption("byskov") <> fail then
    return DIGRAPHS_ChromaticNumberByskov(D);
  elif ValueOption("zykov") <> fail then
    return DIGRAPHS_ChromaticNumberZykov(D);
  elif ValueOption("christofides") <> fail then
    return DIGRAPHS_ChromaticNumberChristofides(D);
  fi;

  # The chromatic number of <D> is at least 3 and at most nr
  D := DigraphMutableCopy(D);
  D := DigraphRemoveAllMultipleEdges(D);
  D := DigraphSymmetricClosure(D);
  MakeImmutable(D);

  if IsCompleteDigraph(D) then
    # chromatic number = nr iff <D> has >= 2 verts & this cond.
    return nr;
  elif nr = 4 then
    # if nr = 4, then 3 is only remaining possible chromatic number
    return 3;
  elif 2 * nr = DigraphNrEdges(D)
      and IsRegularDigraph(D) and Length(OutNeighboursOfVertex(D, 1)) = 2 then
    # <D> is an odd-length cycle graph
    return 3;
  fi;

  # The chromatic number of <D> is at least 3 and at most nr - 1

  # The variable <chrom> is the current best known lower bound for the
  # chromatic number of <D>.
  chrom := 3;

  # Prepare a list of connected components of D whose chromatic number we
  # do not yet know.
  if IsConnectedDigraph(D) then
    comps := [D];
    greedy_bound := RankOfTransformation(DigraphGreedyColouring(D), nr);
    brooks_bound := Maximum(OutDegrees(D));  # Brooks' theorem
    upper := [Minimum(greedy_bound, brooks_bound)];
    chrom := Maximum(CliqueNumber(D), chrom);
  else
    tmp_comps := [];
    tmp_upper := [];
    for comp in DigraphConnectedComponents(D).comps do
      n := Length(comp);
      if chrom < n then
        # If chrom >= n, then we can colour the vertices of comp using any n of
        # the required (at least) chrom colours, and we do not have to consider
        # comp.

        # Note that n > chrom >= 3 and so comp is not null, so no need to check
        # for that.
        comp := InducedSubdigraph(DigraphMutableCopy(D), comp);
        if IsCompleteDigraph(comp) then
          # Since n > chrom, this is an improved lower bound for the overall
          # chromatic number.
          chrom := n;
        elif not IsBipartiteDigraph(comp) then
          # If comp is bipartite, then its chromatic number is 2, and, since
          # the chromatic number of D is >= 3, this component can be
          # ignored.
          greedy_bound := RankOfTransformation(DigraphGreedyColouring(comp),
                                               DigraphNrVertices(comp));
          # Don't need to take odd cycles into account for Brooks' theorem,
          # since they are 3-colourable and the chromatic number of D is >= 3.
          brooks_bound := Maximum(OutDegrees(comp));
          bound := Minimum(greedy_bound, brooks_bound);
          if bound > chrom then
            # If bound <= chrom, then comp can be coloured by at most chrom
            # colours, and so we can ignore comp.
            clique := CliqueNumber(comp);
            if clique = bound then
              # The chromatic number of this component is known, and it can be
              # ignored, and clique = bound > chrom, and so clique is an
              # improved lower bound for the chromatic number of D.
              chrom := clique;
            else
              Add(tmp_comps, comp);
              Add(tmp_upper, bound);
              if clique > chrom then
                chrom := clique;
              fi;
            fi;
          fi;
        fi;
      fi;
    od;

    # Remove the irrelevant components since we have a possibly improved value
    # of chrom.
    comps := [];
    upper := [];

    for i in [1 .. Length(tmp_comps)] do
      if chrom < DigraphNrVertices(tmp_comps[i]) and chrom < tmp_upper[i] then
        Add(comps, tmp_comps[i]);
        Add(upper, tmp_upper[i]);
      fi;
    od;

    # Sort by size, since smaller components are easier to colour
    SortParallel(comps, upper, {x, y} -> Size(x) < Size(y));
  fi;
  for i in [1 .. Length(comps)] do
    # <c> is the current best upper bound for the chromatic number of comps[i]
    c := upper[i];
    while c > chrom and DigraphColouring(comps[i], c - 1) <> fail do
      c := c - 1;
    od;
    if c > chrom then
      chrom := c;
    fi;
  od;
  return chrom;
end);

#
# The following method is currently useless, as the OutNeighbours are computed
# and set whenever a digraph is created.  It could be reinstated later if we
# decide to allow digraphs to exist without known OutNeighbours.
#

# InstallMethod(OutNeighbours,
# "for a digraph with representative out neighbours and group",
# [IsDigraph and HasRepresentativeOutNeighbours and HasDigraphGroup],
# function(D)
#   local gens, sch, reps, out, trace, word, i, w;
#
#   gens := GeneratorsOfGroup(DigraphGroup(D));
#   sch  := DigraphSchreierVector(D);
#   reps := RepresentativeOutNeighbours(D);
#
#   out  := EmptyPlist(DigraphNrVertices(D));
#
#   for i in [1 .. Length(sch)] do
#     if sch[i] < 0 then
#       out[i] := reps[-sch[i]];
#     fi;
#
#     trace  := DIGRAPHS_TraceSchreierVector(gens, sch, i);
#     out[i] := out[trace.representative];
#     word   := trace.word;
#     for w in word do
#        out[i] := OnTuples(out[i], gens[w]);
#     od;
#   od;

#   return out;
# end);

InstallMethod(DigraphAdjacencyFunction, "for a digraph by out-neighbours",
[IsDigraph], D -> {u, v} -> IsDigraphEdge(D, u, v));

InstallMethod(AsTransformation, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  if not IsFunctionalDigraph(D) then
    return fail;
  fi;
  return Transformation(Concatenation(OutNeighbours(D)));
end);

InstallMethod(DigraphNrEdges, "for a digraph", [IsDigraphByOutNeighboursRep],
function(D)
  local m;
  m := DIGRAPH_NREDGES(D);
  if IsImmutableDigraph(D) then
    SetIsEmptyDigraph(D, m = 0);
  fi;
  return m;
end);

InstallMethod(DigraphNrAdjacencies, "for a digraph",
[IsDigraphByOutNeighboursRep], DIGRAPH_NRADJACENCIES);

InstallMethod(DigraphNrAdjacenciesWithoutLoops, "for a digraph",
[IsDigraphByOutNeighboursRep], DIGRAPH_NRADJACENCIESWITHOUTLOOPS);

InstallMethod(DigraphNrLoops,
"for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local i, j, sum;
    sum := 0;
  if HasDigraphHasLoops(D) and not DigraphHasLoops(D) then
    return 0;
  else
    for i in DigraphVertices(D) do
      for j in OutNeighbours(D)[i] do
        if i = j then
          sum := sum + 1;
        fi;
      od;
    od;
  fi;
  if IsImmutableDigraph(D) then
    SetDigraphHasLoops(D, sum <> 0);
  fi;
  return sum;
end);

InstallMethod(DigraphNrLoops,
"for a digraph that knows its adjacency matrix",
[IsDigraphByOutNeighboursRep and HasAdjacencyMatrix],
function(D)
  local A, i;
  A := AdjacencyMatrix(D);
  return Sum(DigraphVertices(D), i -> A[i][i]);
end);

InstallMethod(DigraphEdges, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local out, adj, nr, i, j;
  out := EmptyPlist(DigraphNrEdges(D));
  adj := OutNeighbours(D);
  nr := 0;

  for i in DigraphVertices(D) do
    for j in adj[i] do
      nr := nr + 1;
      out[nr] := [i, j];
    od;
  od;
  return out;
end);

# Hash related attributes

InstallMethod(DigraphHash, "for a digraph", [IsDigraph], DIGRAPH_HASH);

# To make built in Orbit function use DigraphHash
InstallMethod(SparseIntKey, "for an object and digraph",
[IsObject, IsDigraph],
{coll, D} -> DigraphHash
);

# To make orb package use DigraphHash
InstallMethod(ChooseHashFunction, "for a digraph and positive integer",
[IsDigraph, IsInt],
# TODO: (reiniscirpons) remove the rank when new semigroups version gets
# released.
2,
{D, hashlen} -> rec(func := {x, data} -> 1 + (DigraphHash(x) mod data[1]),
                    data := [hashlen])
);

# attributes for digraphs . . .

InstallMethod(AsGraph, "for a digraph", [IsDigraph], Graph);

InstallMethod(DigraphVertices, "for a digraph", [IsDigraph],
D -> [1 .. DigraphNrVertices(D)]);

InstallMethod(DigraphRange,
"for an immutable digraph by out-neighbours",
[IsDigraphByOutNeighboursRep and IsImmutableDigraph],
function(D)
  if not IsBound(D!.DigraphRange) then
    DIGRAPH_SOURCE_RANGE(D);
    SetDigraphSource(D, D!.DigraphSource);
  fi;
  return D!.DigraphRange;
end);

InstallMethod(DigraphRange,
"for a mutable digraph by out-neighbours",
[IsDigraphByOutNeighboursRep and IsMutableDigraph],
D -> DIGRAPH_SOURCE_RANGE(D).DigraphRange);

InstallMethod(DigraphSource,
"for an immutable digraph by out-neighbours",
[IsDigraphByOutNeighboursRep and IsImmutableDigraph],
function(D)
  if not IsBound(D!.DigraphSource) then
    DIGRAPH_SOURCE_RANGE(D);
    SetDigraphRange(D, D!.DigraphRange);
  fi;
  return D!.DigraphSource;
end);

InstallMethod(DigraphSource,
"for a mutable digraph by out-neighbours",
[IsDigraphByOutNeighboursRep and IsMutableDigraph],
D -> DIGRAPH_SOURCE_RANGE(D).DigraphSource);

InstallMethod(InNeighbours, "for a digraph", [IsDigraphByOutNeighboursRep],
D -> DIGRAPH_IN_OUT_NBS(OutNeighbours(D)));

InstallMethod(AdjacencyMatrix, "for a digraph", [IsDigraphByOutNeighboursRep],
ADJACENCY_MATRIX);

InstallMethod(BooleanAdjacencyMatrix, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local n, nbs, mat, i, j;
  n := DigraphNrVertices(D);
  nbs := OutNeighbours(D);
  mat := List(DigraphVertices(D), x -> BlistList([1 .. n], []));
  for i in DigraphVertices(D) do
    for j in nbs[i] do
      mat[i][j] := true;
    od;
  od;
  return mat;
end);

InstallMethod(DigraphShortestDistances, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local vertices, data, sum, distances, v, u;
  if HasDIGRAPHS_ConnectivityData(D) then
    vertices := DigraphVertices(D);
    data := DIGRAPHS_ConnectivityData(D);
    sum := 0;
    for v in vertices do
      if IsBound(data[v]) then
        sum := sum + 1;
      fi;
    od;
    if sum > Int(0.9 * DigraphNrVertices(D))
        or (HasDigraphGroup(D) and not IsTrivial(DigraphGroup(D)))  then
      # adjust the constant 0.9 and possibly make a decision based on
      # how big the group is
      distances := [];
      for u in vertices do
        distances[u] := [];
        for v in vertices do
          distances[u][v] := DigraphShortestDistance(D, u, v);
        od;
      od;
      return distances;
    fi;
  fi;
  return DIGRAPH_SHORTEST_DIST(D);
end);

# returns the vertices (i.e. numbers) of <D> ordered so that there are no
# edges from <out[j]> to <out[i]> for all <i> greater than <j>.

InstallMethod(DigraphTopologicalSort, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
D -> DIGRAPH_TOPO_SORT(OutNeighbours(D)));

InstallMethod(DigraphStronglyConnectedComponents,
"for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local verts;

  if HasIsAcyclicDigraph(D) and IsAcyclicDigraph(D) then
    verts := DigraphVertices(D);
    return rec(comps := List(verts, x -> [x]), id := verts * 1);

  elif HasIsStronglyConnectedDigraph(D)
      and IsStronglyConnectedDigraph(D) then
    verts := DigraphVertices(D);
    return rec(comps := [verts * 1], id := verts * 0 + 1);
  fi;

  return GABOW_SCC(OutNeighbours(D));
end);

InstallMethod(DigraphNrStronglyConnectedComponents, "for a digraph",
[IsDigraph],
D -> Length(DigraphStronglyConnectedComponents(D).comps));

InstallMethod(DigraphConnectedComponents, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
DIGRAPH_CONNECTED_COMPONENTS);

InstallMethod(DigraphNrConnectedComponents, "for a digraph",
[IsDigraph],
D -> Length(DigraphConnectedComponents(D).comps));

InstallImmediateMethod(DigraphNrConnectedComponents,
"for a digraph with known connected components",
IsDigraph and HasDigraphConnectedComponents, 0,
D -> Length(DigraphConnectedComponents(D).comps));

InstallMethod(OutDegrees, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local adj, degs, i;
  adj := OutNeighbours(D);
  degs := EmptyPlist(DigraphNrVertices(D));
  for i in DigraphVertices(D) do
    degs[i] := Length(adj[i]);
  od;
  return degs;
end);

InstallMethod(InDegrees, "for a digraph with in neighbours",
[IsDigraph and HasInNeighbours],
2,  # to beat the method for IsDigraphByOutNeighboursRep
function(D)
  local inn, degs, i;
  inn := InNeighbours(D);
  degs := EmptyPlist(DigraphNrVertices(D));
  for i in DigraphVertices(D) do
    degs[i] := Length(inn[i]);
  od;
  return degs;
end);

InstallMethod(InDegrees, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local adj, degs, x, i;
  adj := OutNeighbours(D);
  degs := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
  for x in adj do
    for i in x do
      degs[i] := degs[i] + 1;
    od;
  od;
  return degs;
end);

InstallMethod(OutDegreeSequence, "for a digraph", [IsDigraph],
function(D)
  D := ShallowCopy(OutDegrees(D));
  Sort(D, {a, b} -> b < a);
  return D;
  # return SortedList(OutDegrees(D), {a, b} -> b < a);
end);

InstallMethod(OutDegreeSequence,
"for a digraph by out-neighbours with known digraph group",
[IsDigraphByOutNeighboursRep and HasDigraphGroup],
function(D)
  local out, adj, orbs, orb;
  out := [];
  adj := OutNeighbours(D);
  orbs := DigraphOrbits(D);
  for orb in orbs do
    Append(out, [1 .. Length(orb)] * 0 + Length(adj[orb[1]]));
  od;
  Sort(out, {a, b} -> b < a);
  return out;
end);

InstallMethod(OutDegreeSet, "for a digraph", [IsDigraph],
D -> Set(ShallowCopy(OutDegrees(D))));

InstallMethod(InDegreeSequence, "for a digraph", [IsDigraph],
function(D)
  D := ShallowCopy(InDegrees(D));
  Sort(D, {a, b} -> b < a);
  return D;
  # return SortedList(OutDegrees(D), {a, b} -> b < a);
end);

InstallMethod(InDegreeSequence,
"for a digraph with known digraph group and in-neighbours",
[IsDigraph and HasDigraphGroup and HasInNeighbours],
function(D)
  local out, adj, orbs, orb;
  out := [];
  adj := InNeighbours(D);
  orbs := DigraphOrbits(D);
  for orb in orbs do
    Append(out, [1 .. Length(orb)] * 0 + Length(adj[orb[1]]));
  od;
  Sort(out, {a, b} -> b < a);
  return out;
end);

InstallMethod(InDegreeSet, "for a digraph", [IsDigraph],
D -> Set(ShallowCopy(InDegrees(D))));

InstallMethod(DigraphSources, "for a digraph with in-degrees",
[IsDigraph and HasInDegrees], 3,
function(D)
  local degs;
  degs := InDegrees(D);
  return Filtered(DigraphVertices(D), x -> degs[x] = 0);
end);

InstallMethod(DigraphSources, "for a digraph with in-neighbours",
[IsDigraph and HasInNeighbours],
2,  # to beat the method for IsDigraphByOutNeighboursRep
function(D)
  local inn, sources, count, i;

  inn := InNeighbours(D);
  sources := EmptyPlist(DigraphNrVertices(D));
  count := 0;
  for i in DigraphVertices(D) do
    if IsEmpty(inn[i]) then
      count := count + 1;
      sources[count] := i;
    fi;
  od;
  ShrinkAllocationPlist(sources);
  return sources;
end);

InstallMethod(DigraphSources, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local out, seen, tmp, next, v;
  out  := OutNeighbours(D);
  seen := BlistList(DigraphVertices(D), []);
  for next in out do
    for v in next do
      seen[v] := true;
    od;
  od;
  # FIXME use FlipBlist (when available)
  tmp  := BlistList(DigraphVertices(D), DigraphVertices(D));
  SubtractBlist(tmp, seen);
  return ListBlist(DigraphVertices(D), tmp);
end);

InstallMethod(DigraphSinks, "for a digraph with out-degrees",
[IsDigraph and HasOutDegrees],
2,  # to beat the method for IsDigraphByOutNeighboursRep
function(D)
  local degs;
  degs := OutDegrees(D);
  return Filtered(DigraphVertices(D), x -> degs[x] = 0);
end);

InstallMethod(DigraphSinks, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local out, sinks, count, i;

  out   := OutNeighbours(D);
  sinks := [];
  count := 0;
  for i in DigraphVertices(D) do
    if IsEmpty(out[i]) then
      count := count + 1;
      sinks[count] := i;
    fi;
  od;
  return sinks;
end);

InstallMethod(DigraphPeriod, "for a digraph", [IsDigraphByOutNeighboursRep],
function(D)
  local comps, out, deg, nrvisited, period, stack, len, depth, current,
        olddepth, i;

  if HasIsAcyclicDigraph(D) and IsAcyclicDigraph(D) then
    return 0;
  fi;

  comps := DigraphStronglyConnectedComponents(D).comps;
  out := OutNeighbours(D);
  deg := OutDegrees(D);

  nrvisited := [1 .. Length(DigraphVertices(D))] * 0;
  period := 0;

  for i in [1 .. Length(comps)] do
    stack := [comps[i][1]];
    len := 1;
    depth := EmptyPlist(Length(DigraphVertices(D)));
    depth[comps[i][1]] := 1;
    while len <> 0 do
      current := stack[len];
      if nrvisited[current] = deg[current] then
        len := len - 1;
      else
        nrvisited[current] := nrvisited[current] + 1;
        len := len + 1;
        stack[len] := out[current][nrvisited[current]];
        olddepth := depth[current];
        if IsBound(depth[stack[len]]) then
          period := GcdInt(period, depth[stack[len]] - olddepth - 1);
          if period = 1 then
            return period;
          fi;
        else
          depth[stack[len]] := olddepth + 1;
        fi;
      fi;
    od;
  od;

  if period = 0 and IsImmutableDigraph(D) then
    SetIsAcyclicDigraph(D, true);
  fi;

  return period;
end);

InstallMethod(DIGRAPHS_ConnectivityData, "for a digraph", [IsDigraph],
D -> EmptyPlist(0));

BindGlobal("DIGRAPH_ConnectivityDataForVertex",
function(D, v)
  local data, out_nbs, record, orbnum, reps, i, next, laynum, localGirth,
        layers, sum, localParameters, nprev, nhere, nnext, lnum, localDiameter,
        layerNumbers, x, y, tree, expand, stab, edges, edge;

  data := DIGRAPHS_ConnectivityData(D);
  if IsBound(data[v]) then
    return data[v];
  fi;

  expand := false;
  if HasDigraphGroup(D) then
    stab            := DigraphStabilizer(D, v);
    record          := DIGRAPHS_Orbits(stab, DigraphVertices(D));
    orbnum          := record.lookup;
    reps            := List(record.orbits, Representative);
    if Length(record.orbits) < DigraphNrVertices(D) then
        expand := true;
    fi;
  else
    orbnum          := [1 .. DigraphNrVertices(D)];
    reps            := [1 .. DigraphNrVertices(D)];
  fi;
  out_nbs         := OutNeighbours(D);
  i               := 1;
  next            := [orbnum[v]];
  laynum          := ListWithIdenticalEntries(Length(reps), 0);
  laynum[next[1]] := 1;
  tree            := ListWithIdenticalEntries(DigraphNrVertices(D), 0);
  tree[v]         := -1;
  localGirth      := -1;
  layers          := [next];
  sum             := 1;
  localParameters := [];

  # localDiameter is the length of the longest shortest path starting at v
  #
  # localParameters is a list of 3-tuples [a_{i - 1}, b_{i - 1}, c_{i - 1}] for
  # each i between 1 and localDiameter where c_i (respectively a_i and b_i) is
  # the number of vertices at distance i − 1 (respectively i and i + 1) from v
  # that are adjacent to a vertex w at distance i from v.

  # <tree> gives a shortest path spanning tree rooted at <v> and is used by
  # the ShortestPathSpanningTree method.

  while Length(next) > 0 do
    next := [];
    for x in layers[i] do
      nprev := 0;
      nhere := 0;
      nnext := 0;
      for y in out_nbs[reps[x]] do
        lnum := laynum[orbnum[y]];
        if i > 1 and lnum = i - 1 then
          nprev := nprev + 1;
        elif lnum = i then
          nhere := nhere + 1;
        elif lnum = i + 1 then
          nnext := nnext + 1;
        elif lnum = 0 then
          AddSet(next, orbnum[y]);
          nnext := nnext + 1;
          laynum[orbnum[y]] := i + 1;
          if not expand then
            tree[y] := reps[x];
          else
            edges := Orbit(stab, [reps[x], y], OnTuples);
            for edge in edges do
              if tree[edge[2]] = 0 or tree[edge[2]] > edge[1] then
                tree[edge[2]] := edge[1];
              fi;
            od;
          fi;
        fi;
      od;
      if (localGirth = -1 or localGirth = 2 * i - 1) and nprev > 1 then
        localGirth := 2 * (i - 1);
      fi;
      if localGirth = -1 and nhere > 0 then
        localGirth := 2 * i - 1;
      fi;
      if not IsBound(localParameters[i]) then
         localParameters[i] := [nprev, nhere, nnext];
      else
         if nprev <> localParameters[i][1] then
            localParameters[i][1] := -1;
         fi;
         if nhere <> localParameters[i][2] then
            localParameters[i][2] := -1;
         fi;
         if nnext <> localParameters[i][3] then
            localParameters[i][3] := -1;
         fi;
      fi;
    od;
    if Length(next) > 0 then
      i := i + 1;
      layers[i] := next;
      sum := sum + Length(next);
    fi;
  od;
  if sum = Length(reps) then
     localDiameter := Length(layers) - 1;
  else
     localDiameter := -1;
  fi;

  layerNumbers := [];
  for i in DigraphVertices(D) do
     layerNumbers[i] := laynum[orbnum[i]];
  od;
  data[v] := rec(layerNumbers    := layerNumbers,
                 localDiameter   := localDiameter,
                 localGirth      := localGirth,
                 localParameters := localParameters,
                 layers          := layers,
                 spstree         := tree);
  return data[v];
end);

BindGlobal("DIGRAPHS_DiameterAndUndirectedGirth",
function(D)
  local outer_reps, diameter, girth, v, record, localGirth,
        localDiameter, i;

  # This function attempts to find the diameter and undirected girth of a given
  # D, using its DigraphGroup.  For some digraphs, the main algorithm will
  # not produce a sensible answer, so there are checks at the start and end to
  # alter the answer for the diameter/girth if necessary.  This function is
  # called, if appropriate, by DigraphDiameter and DigraphUndirectedGirth.

  if DigraphHasNoVertices(D) and IsImmutableDigraph(D) then
    SetDigraphDiameter(D, fail);
    SetDigraphUndirectedGirth(D, infinity);
    return rec(diameter := fail, girth := infinity);
  fi;

  # TODO improve this, really check if the complexity is better with the group
  # or without, or if the group is not known, but the number of vertices makes
  # the usual algorithm impossible.

  outer_reps := DigraphOrbitReps(D);
  diameter   := 0;
  girth      := infinity;

  for i in [1 .. Length(outer_reps)] do
    v := outer_reps[i];
    record     := DIGRAPH_ConnectivityDataForVertex(D, v);
    localGirth := record.localGirth;
    localDiameter := record.localDiameter;

    if localDiameter > diameter then
      diameter := localDiameter;
    fi;

    if localGirth <> -1 and localGirth < girth then
      girth := localGirth;
    fi;
  od;

  # Checks to ensure both components are valid
  if not IsStronglyConnectedDigraph(D) then
    diameter := fail;
  fi;
  if DigraphHasLoops(D) then
    girth := 1;
  elif IsMultiDigraph(D) then
    girth := 2;
  fi;

  if IsImmutableDigraph(D) then
    SetDigraphDiameter(D, diameter);
    SetDigraphUndirectedGirth(D, girth);
  fi;
  return rec(diameter := diameter, girth := girth);
end);

InstallMethod(DigraphDiameter, "for a digraph", [IsDigraph],
function(D)
  if not IsStronglyConnectedDigraph(D) then
    # Diameter undefined
    return fail;
  elif HasDigraphGroup(D) and Size(DigraphGroup(D)) > 1 then
    # Use the group to calculate the diameter
    return DIGRAPHS_DiameterAndUndirectedGirth(D).diameter;
  fi;
  # Use the C function
  return DIGRAPH_DIAMETER(D);
end);

InstallMethod(DigraphUndirectedGirth, "for a digraph", [IsDigraph],
function(D)
  # This is only defined on undirected graphs (i.e. symmetric digraphs)
  if not IsSymmetricDigraph(D) then
    ErrorNoReturn("the argument <D> must be a symmetric digraph,");
  elif DigraphHasLoops(D) then
    # A loop is a cycle of length 1
    return 1;
  elif IsMultiDigraph(D) then
    # A pair of multiple edges is a cycle of length 2
    return 2;
  fi;
  # Otherwise D is simple
  return DIGRAPHS_DiameterAndUndirectedGirth(D).girth;
end);

InstallMethod(DigraphGirth, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local verts, girth, out, dist, i, j;
  if DigraphHasLoops(D) then
    return 1;
  fi;
  # Only consider one vertex from each orbit
  if HasDigraphGroup(D) and not IsTrivial(DigraphGroup(D)) then
    verts := DigraphOrbitReps(D);
  else
    verts := DigraphVertices(D);
  fi;
  girth := infinity;
  out := OutNeighbours(D);
  for i in verts do
    for j in out[i] do
      dist := DigraphShortestDistance(D, j, i);
      # distance [j,i] + 1 equals the cycle length
      if dist <> fail and dist + 1 < girth then
        girth := dist + 1;
        if girth = 2 then
          return girth;
        fi;
      fi;
    od;
  od;
  return girth;
end);

InstallMethod(DigraphOddGirth, "for a digraph",
[IsDigraph],
function(digraph)
  local comps, girth, oddgirth, A, B, gr, k, comp;

  if IsAcyclicDigraph(digraph) then
    return infinity;
  elif IsOddInt(DigraphGirth(digraph)) then
    # No need to check girth isn't infinity, as we have
    # that digraph is not acyclic.
    return DigraphGirth(digraph);
  fi;
  comps := DigraphStronglyConnectedComponents(digraph).comps;
  if Length(comps) > 1 and IsMutableDigraph(digraph) then
    # Necessary because InducedSubdigraph alters mutable args
    digraph := DigraphImmutableCopy(digraph);
  fi;
  oddgirth := infinity;
  for comp in comps do
    if Length(comps) > 1 then  # i.e. if not IsStronglyConnectedDigraph(digraph)
      gr := InducedSubdigraph(digraph, comp);
    else
      gr := digraph;
      # If digraph is strongly connected, then we needn't
      # induce the subdigraph of its strongly connected comp.
    fi;
    if not IsAcyclicDigraph(gr) then
      girth := DigraphGirth(gr);
      if IsOddInt(girth) and girth < oddgirth then
        oddgirth := girth;
      else
        k := girth - 1;
        A := AdjacencyMatrix(gr) ^ k;
        B := AdjacencyMatrix(gr) ^ 2;
        while k <= DigraphNrEdges(gr) + 2 and k < DigraphNrVertices(gr)
            and k < oddgirth - 2 do
          A := A * B;
          k := k + 2;
          if Trace(A) <> 0 then
            # It suffices to find the trace as the entries of A are positive.
            oddgirth := k;
          fi;
        od;
      fi;
    fi;
  od;

  return oddgirth;
end);

InstallMethod(DigraphLongestSimpleCircuit, "for a digraph",
[IsDigraph],
function(D)
  local circs, lens, max;
  if IsAcyclicDigraph(D) then
    return fail;
  fi;
  circs := DigraphAllSimpleCircuits(D);
  lens := List(circs, Length);
  max := Maximum(lens);
  return circs[Position(lens, max)];
end);

InstallMethod(DigraphAllSimpleCircuits, "for a digraph by out-neighbours",
[IsDigraphByOutNeighboursRep],
function(D)
  local UNBLOCK, CIRCUIT, out, stack, endofstack, C, scc, n, blocked, B,
  c_comp, comp, s, loops, i, d_labels, c_labels;

  if IsEmptyDigraph(D) then
    return [];
  fi;

  UNBLOCK := function(u)
    local w;
    blocked[u] := false;
    while not IsEmpty(B[u]) do
      w := B[u][1];
      Remove(B[u], 1);
      if blocked[w] then
        UNBLOCK(w);
      fi;
    od;
  end;

  CIRCUIT := function(v, component)
    local f, buffer, dummy, w;

    f := false;
    endofstack := endofstack + 1;
    stack[endofstack] := v;
    blocked[v] := true;

    for w in OutNeighboursOfVertex(component, v) do
      if w = 1 then
        buffer := stack{[1 .. endofstack]};
        Add(out, DigraphVertexLabels(component){buffer});
        f := true;
      elif blocked[w] = false then
        dummy := CIRCUIT(w, component);
        if dummy then
          f := true;
        fi;
      fi;
    od;

    if f then
      UNBLOCK(v);
    else
      for w in OutNeighboursOfVertex(component, v) do
        if not w in B[w] then
          Add(B[w], v);
        fi;
      od;
    fi;

    endofstack := endofstack - 1;
    return f;
  end;

  out := [];
  stack := [];
  endofstack := 0;

  # Reduce the D, remove loops, and store the correct vertex labels
  C := DigraphRemoveLoops(ReducedDigraph(DigraphMutableCopyIfMutable(D)));
  MakeImmutable(C);
  if HaveVertexLabelsBeenAssigned(D)
      and DigraphVertexLabels(D) <> DigraphVertices(D) then
    # We require the labels of the digraph <C> to be the original nodes in <D>
    # (this is used in CIRCUIT above). If <D> has other vertex labels, then
    # these are copied into <C> by the functions above, and aren't then
    # necessarily the original nodes in <D>.
    d_labels := DigraphVertexLabels(D);
    c_labels := List(DigraphVertexLabels(C), x -> Position(d_labels, x));
    SetDigraphVertexLabels(C, c_labels);
  fi;

  # Strongly connected components of the reduced graph
  scc := DigraphStronglyConnectedComponents(C);

  # B and blocked only need to be as long as the longest connected component
  n := Maximum(List(scc.comps, Length));
  blocked := BlistList([1 .. n], []);
  B := List([1 .. n], x -> []);

  # Perform algorithm once per connected component of the whole digraph
  for c_comp in scc.comps do
    n := Length(c_comp);
    if n = 1 then
      continue;
    fi;
    c_comp := InducedSubdigraph(C, c_comp);  # C is definitely immutable
    comp := c_comp;
    s := 1;
    while s < n do
      if s <> 1 then
        comp := InducedSubdigraph(c_comp, [s .. n]);
        comp := InducedSubdigraph(comp,
                                  DigraphStronglyConnectedComponent(comp, 1));
      fi;

      if not IsEmptyDigraph(comp) then
        # TODO would it be faster/better to create blocked as a new BlistList?
        # Are these things already going to be initialised anyway?
        for i in DigraphVertices(comp) do
          blocked[i] := false;
          B[i] := [];
        od;
        CIRCUIT(1, comp);
      fi;
      s := s + 1;
    od;
  od;
  loops := List(DigraphLoops(D), x -> [x]);
  return Concatenation(loops, out);
end);

# Compute all undirected simple circuits by filtering the output
# of DigraphAllSimpleCircuits
InstallMethod(DigraphAllUndirectedSimpleCircuits, "for a digraph", [IsDigraph],
function(D)
  local digraph, cycles, keep, cycle, cycleRev;
  digraph := DigraphSymmetricClosure(
      DigraphRemoveAllMultipleEdges(DigraphMutableCopyIfMutable(D)));
  cycles := Filtered(DigraphAllSimpleCircuits(digraph),
    c -> Length(c) > 2 or Length(c) = 1);
  keep := HashSet();
  for cycle in cycles do
    cycleRev := [cycle[1]];
    Append(cycleRev, Reversed(cycle{[2 .. Length(cycle)]}));
    if (not cycle in keep and not cycleRev in keep) or Length(cycle) = 1 then
      AddSet(keep, cycle);
    fi;
  od;
  return Set(keep);
end);

# Compute all chordless cycles for a given symmetric digraph
# Algorithm based on https://arxiv.org/pdf/1404.7610
InstallMethod(DigraphAllChordlessCyclesOfMaximalLength,
"for a digraph and an integer", [IsDigraph, IsInt],
function(D, maxLength)
  local BlockNeighbours, UnblockNeighbours,
  Triplets, CCExtension, digraph, temp, T, C, blocked, triple;

    if IsEmptyDigraph(D) then
        return [];
    fi;

    BlockNeighbours := function(digraph, v, blocked)
        local u;
        for u in OutNeighboursOfVertex(digraph, v) do
            blocked[u] := blocked[u] + 1;
        od;
        return blocked;
    end;

    UnblockNeighbours := function(digraph, v, blocked)
        local u;
        for u in OutNeighboursOfVertex(digraph, v) do
            if blocked[u] > 0 then
                blocked[u] := blocked[u] - 1;
            fi;
        od;
        return blocked;
    end;

    # Computes all possible triplets
    Triplets := function(digraph)
        local T, C, u, pair, x, y, labels;
        T := [];
        C := [];
        for u in DigraphVertices(digraph) do
            for pair in Combinations(OutNeighboursOfVertex(digraph, u), 2) do
                x := pair[1];
                y := pair[2];
                labels := DigraphVertexLabels(digraph);
                if labels[u] < labels[x] and labels[x] < labels[y] then
                    if not IsDigraphEdge(digraph, x, y) then
                        Add(T, [x, u, y]);
                    else
                        Add(C, [x, u, y]);
                    fi;
                elif labels[u] < labels[y] and labels[y] < labels[x] then
                    if not IsDigraphEdge(digraph, x, y) then
                        Add(T, [y, u, x]);
                    else
                        Add(C, [y, u, x]);
                    fi;
                fi;
            od;
        od;
        return [T, C];
    end;

    # Extends a given chordless path if possible
    CCExtension := function(digraph, path, C, key, blocked)
        local v, extendedPath, data;
        blocked := BlockNeighbours(digraph, path[Length(path)], blocked);
        for v in OutNeighboursOfVertex(digraph, path[Length(path)]) do
            if DigraphVertexLabel(digraph, v) > key and blocked[v] = 1
              and Length(path) < maxLength then
                extendedPath := Concatenation(path, [v]);
                if IsDigraphEdge(digraph, v, path[1]) then
                    Add(C, extendedPath);
                else
                    data := CCExtension(digraph, extendedPath, C, key, blocked);
                    C := data[1];
                    blocked := data[2];
                fi;
            fi;
        od;
        blocked := UnblockNeighbours(digraph, path[Length(path)], blocked);
        return [C, blocked];
    end;

    digraph := DigraphMutableCopy(D);
    DigraphSymmetricClosure(DigraphRemoveLoops(
                 DigraphRemoveAllMultipleEdges(digraph)));
    MakeImmutable(digraph);
    SetDigraphVertexLabels(digraph,
                 Reversed(DigraphDegeneracyOrdering(digraph)));

    temp := Triplets(digraph);
    T := temp[1];
    C := temp[2];
    blocked := List(DigraphVertices(digraph), i -> 0);
    while T <> [] do
        triple := Remove(T);
        blocked := BlockNeighbours(digraph, triple[2], blocked);
        temp := CCExtension(digraph, triple, C,
                              DigraphVertexLabel(digraph, triple[2]), blocked);
        C := temp[1];
        blocked := temp[2];
        blocked := UnblockNeighbours(digraph, triple[2], blocked);
    od;
    return C;
end);

InstallMethod(DigraphAllChordlessCycles, "for a digraph",
[IsDigraph],
D -> DigraphAllChordlessCyclesOfMaximalLength(D, INTOBJ_MAX));

# Compute for a given rotation system the facial walks
InstallMethod(FacialWalks, "for a digraph and a list",
[IsDigraph, IsDenseList],
function(D, rotationSystem)

    local FacialWalk, facialWalks, remEdges, cycle;

    if not IsEulerianDigraph(D) then
        Error("the 1st argument (digraph <D>) must be Eulerian, but it is not");
    fi;

    if Length(rotationSystem) <> DigraphNrVertices(D) then
        Error("the 2nd argument (list <rotationSystem>) is not a rotation ",
              "system for the 1st argument (digraph <D>), expected a ",
              "dense list of length ", DigraphNrVertices(D),
              "but found dense list of length ", Length(rotationSystem));
    fi;

    if Difference(Union(rotationSystem), DigraphVertices(D))
       <> [] then
        Error("the 2nd argument (dense list <rotationSystem>) is not ",
              "a rotation system for the 1st argument (digraph <D>), ",
              "expected the union to be ", DigraphVertices(D), " but found ",
              Union(rotationSystem));
    fi;

    # computes a facial cycles starting with the edge 'startEdge'
    FacialWalk := function(rotationSystem, startEdge)
        local startVertex, preVertex, actVertex, cycle, nextVertex, pos;

        startVertex := startEdge[1];
        actVertex := startEdge[2];
        preVertex := startVertex;

        cycle := [startVertex, actVertex];

        nextVertex := 0;  # just an initialization
        while true do
            pos := Position(rotationSystem[actVertex], preVertex);

            if pos < Length(rotationSystem[actVertex]) then
                nextVertex := rotationSystem[actVertex][pos + 1];
            else
                nextVertex := rotationSystem[actVertex][1];
            fi;
            if nextVertex <> startEdge[2] or actVertex <> startVertex then
                Add(cycle, nextVertex);
                Remove(remEdges, Position(remEdges, [preVertex, actVertex]));
                preVertex := actVertex;
                actVertex := nextVertex;
            else
                break;
            fi;
        od;
        Remove(remEdges, Position(remEdges, [preVertex, startVertex]));
        # Remove the last vertex, otherwise otherwise
        # the start vertex is contained twice
        Remove(cycle);
        return cycle;
    end;

    D := DigraphRemoveLoops(DigraphRemoveAllMultipleEdges(
        DigraphMutableCopyIfMutable(D)));

    facialWalks := [];
    remEdges := ShallowCopy(DigraphEdges(D));

    while remEdges <> [] do
        cycle := FacialWalk(rotationSystem, remEdges[1]);
        Add(facialWalks, cycle);
    od;
    return facialWalks;
end);

# Computes the minimal cyclic edge cut of connected cubic graphs with at
# least 8 vertices based on the paper "An Algorithm for Cyclic Edge Connectivity
# of Cubic Graphs"
InstallMethod(MinimalCyclicEdgeCut, "for a digraph", [IsDigraph],
function(digraph)
  local girth, cut, cutsize, vertex, treedepth, paths, treev, treew,
  minimal_cycle, v, w, FullTree, FindCut, FindPath, e;

  # Compute a full tree of the given depth centered at the given vertex
  FullTree := function(digraph, vertex, depth)
    local result, i, lastTree, node;
    result := Set([vertex]);
    for i in [1 .. depth] do
      lastTree := ShallowCopy(result);
      for node in lastTree do
          UniteSet(result, OutNeighboursOfVertex(digraph, node));
      od;
    od;
    return result;
  end;

  # Compute a maximal amount of paths between the two distinct sets of
  # nodes treev and treew by using a variation of the Ford-Fulkerson algorithm
  FindPath := function(digraph, treev, treew)
    local digraphCopy, pathlist, newPathFound, node1, node2, path, i;

    digraphCopy := DigraphMutableCopy(digraph);
    pathlist := [];
    newPathFound := true;
    while newPathFound do
      newPathFound := false;
      for node1 in treev do
        for node2 in treew do
          path := DigraphShortestPath(digraphCopy, node1, node2);

          if path <> fail then
            for i in [1 .. (Length(path[1]) - 1)] do
              # remove edges corresponding to the current path
              digraphCopy := DigraphRemoveEdge(digraphCopy,
                [path[1][i], path[1][i + 1]]);
              # add backward edges, if they are missing
              if not [path[1][i + 1], path[1][i]] in
                    DigraphEdges(digraphCopy) then
                  digraphCopy := DigraphAddEdge(digraphCopy,
                    [path[1][i + 1], path[1][i]]);
              fi;
            od;
            Append(pathlist, [path[1]]);
            newPathFound := true;
            break;
          fi;
        od;
      od;
    od;
    return pathlist;
  end;

  # Finds a cut that disconnects the node sets treev and treew in
  # the given digraph by using the paths found in FindPath
  FindCut := function(digraph, treev, treew, allPaths)
    local pathByEdges, path, edgeList, i, cutEdges, edge, component1,
    component2, edgeSet, digraphCopy, permutations, nodeSet,
    v, pathInducedSubgraph;

    # Convert the paths into the list of corresponding edges
    pathByEdges := [];
    for path in allPaths do
      edgeList := [];
      for i in [1 .. Length(path) - 1] do
          Append(edgeList, [[path[i], path[i + 1]]]);
      od;
      Append(pathByEdges, [edgeList]);
    od;

    nodeSet := Set([]);
    for path in allPaths do
      for v in path do
        AddSet(nodeSet, v);
      od;
    od;

    edgeSet := Set([]);
    for v in DigraphVertices(digraph) do
      for w in OutNeighboursOfVertex(digraph, v) do
        UniteSet(edgeSet, [[v, w]]);
        UniteSet(edgeSet, [[w, v]]);
      od;
    od;

    # We can find a cut by removing one specific edge from every path,
--> --------------------

--> maximum size reached

--> --------------------

[ zur Elbe Produktseite wechseln0.54Quellennavigators  Analyse erneut starten  ]

                                                                                                                                                                                                                                                                                                                                                                                                     


Neuigkeiten

     Aktuelles
     Motto des Tages

Software

     Produkte
     Quellcodebibliothek

Aktivitäten

     Artikel über Sicherheit
     Anleitung zur Aktivierung von SSL

Muße

     Gedichte
     Musik
     Bilder

Jenseits des Üblichen ....

Besucherstatistik

Besucherstatistik

Monitoring

Montastic status badge