Add AbsorptionExpectedSteps

Author: mtorpey

doc/attr.xml

@@ -1138,7 +1138,8 @@ gap> ArticulationPoints(D);
 
     Strongly connected components are indexed in the order given by
     <Ref Attr="DigraphStronglyConnectedComponents"/>. See also
-    <Ref Oper="DigraphRandomWalk"/>.
+    <Ref Oper="DigraphRandomWalk"/> and
+    <Ref Attr="DigraphAbsorptionExpectedSteps"/>.
     <P/>
 <Example><![CDATA[
 gap> gr := Digraph([[2, 3, 4], [3], [2], []]);
@@ -1152,6 +1153,38 @@ gap> DigraphAbsorptionProbabilities(gr);
 </ManSection>
 <#/GAPDoc>
 
+<#GAPDoc Label="DigraphAbsorptionExpectedSteps">
+<ManSection>
+  <Attr Name="DigraphAbsorptionExpectedSteps" Arg="digraph"/>
+  <Returns>A list of rational numbers.</Returns>
+  <Description>
+    A random walk of unbounded length through a finite digraph may pass through
+    several different strongly connected components (SCCs).  However, with
+    probability 1 it must eventually reach an SCC which it can never leave,
+    because the SCC has no out-edges leading to other SCCs.  When this happens,
+    we say the walk has been <E>absorbed</E>.
+    <P/>
+    
+    <C>DigraphAbsorptionExpectedSteps</C> returns a list <C>L</C> with length
+    equal to the number of vertices of <A>digraph</A>, where <C>L[v]</C>
+    is the expected number of steps a random walk starting at vertex <C>v</C>
+    should take before it is absorbed – that is, the expected number of steps to
+    reach an SCC that it can never leave.
+    <P/>
+
+    See also <Ref Oper="DigraphRandomWalk"/> and
+    <Ref Attr="DigraphAbsorptionProbabilities"/>.
+    <P/>
+<Example><![CDATA[
+gap> gr := Digraph([[2], [3, 4], [], [2]]);
+<immutable digraph with 4 vertices, 4 edges>
+gap> DigraphAbsorptionExpectedSteps(gr);
+[ 4, 3, 0, 4 ]
+]]></Example>
+  </Description>
+</ManSection>
+<#/GAPDoc>
+
 <#GAPDoc Label="StrongOrientation">
 <ManSection>
   <Oper Name="StrongOrientation" Arg="D"/>

doc/z-chap4.xml

@@ -70,6 +70,7 @@
     <#Include Label="DigraphShortestPath">
     <#Include Label="DigraphRandomWalk">
     <#Include Label="DigraphAbsorptionProbabilities">
+    <#Include Label="DigraphAbsorptionExpectedSteps">
     <#Include Label="Dominators">
     <#Include Label="DominatorTree">
     <#Include Label="IteratorOfPaths">

gap/attr.gd

@@ -51,7 +51,9 @@ DeclareAttribute("DigraphOddGirth", IsDigraph);
 DeclareAttribute("DigraphUndirectedGirth", IsDigraph);
 DeclareAttribute("ArticulationPoints", IsDigraph);
 DeclareSynonymAttr("CutVertices", ArticulationPoints);
+DeclareAttribute("DIGRAPHS_AbsorbingMarkovChain", IsDigraph);
 DeclareAttribute("DigraphAbsorptionProbabilities", IsDigraph);
+DeclareAttribute("DigraphAbsorptionExpectedSteps", IsDigraph);
 
 DeclareAttribute("DigraphAllSimpleCircuits", IsDigraph);
 DeclareAttribute("DigraphLongestSimpleCircuit", IsDigraph);

gap/attr.gi

@@ -192,26 +192,16 @@ function(I, subtracted_set, size_bound)
 end
 );
 
-InstallMethod(DigraphAbsorptionProbabilities,
+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, N, chances_of_absorption, output, comp_no,
-        c;
-  # No vertices: trivial matrix
-  if DigraphNrVertices(D) = 0 then
-    return [];
-  fi;
-
+        w_comp, sink_comp_index, j, fundamental_mat;
   scc := DigraphStronglyConnectedComponents(D);
 
-  # One SCC only: all vertices stay in it with probability 1.
-  if Length(scc.comps) = 1 then
-    return ListWithIdenticalEntries(DigraphNrVertices(D), [1]);
-  fi;
-
   # 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);
@@ -229,6 +219,14 @@ function(D)
   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);
@@ -257,21 +255,57 @@ function(D)
     od;
   od;
 
-  # Compute the chances as t tends to infinity using formula (I - Q)^-1 * R
-  N := Inverse(IdentityMat(nr_transients) - transient_mat);
-  chances_of_absorption := N * absorption_mat;
+  # "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);
 
-  # Rows are vertices, columns are SCCs
+  # 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 above
+  # 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
@@ -283,6 +317,28 @@ function(D)
   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,

tst/standard/attr.tst

@@ -2883,23 +2883,6 @@ gap> DigraphStronglyConnectedComponents(gr).comps;  # this ordering is used
 [ [ 2, 3 ], [ 4 ], [ 1 ] ]
 gap> DigraphAbsorptionProbabilities(gr);
 [ [ 2/3, 1/3, 0 ], [ 1, 0, 0 ], [ 1, 0, 0 ], [ 0, 1, 0 ] ]
-gap> soccer := Digraph([[7, 2, 3, 5, 1, 4],  # Motivating example:
->                       [7, 2, 3, 5, 5, 4],  # game of 'Soccer Dice'
->                       [7, 5, 5, 5, 5, 1],
->                       [7, 7, 2, 5, 5, 5],
->                       [6, 6, 7, 7, 7, 4],
->                       [], []]);;
-gap> DigraphStronglyConnectedComponents(soccer).comps;  # this ordering is used
-[ [ 7 ], [ 6 ], [ 1, 2, 3, 5, 4 ] ]
-gap> DigraphAbsorptionProbabilities(soccer) =
->     [[3473 / 4493, 1020 / 4493, 0],
->      [3365 / 4493, 1128 / 4493, 0],
->      [3211 / 4493, 1282 / 4493, 0],
->      [3471 / 4493, 1022 / 4493, 0],
->      [2825 / 4493, 1668 / 4493, 0],
->      [0, 1, 0],
->      [1, 0, 0]];
-true
 gap> DigraphAbsorptionProbabilities(EmptyDigraph(0));
 [  ]
 gap> DigraphAbsorptionProbabilities(EmptyDigraph(1));
@@ -2908,6 +2891,8 @@ gap> DigraphAbsorptionProbabilities(CompleteDigraph(5));
 [ [ 1 ], [ 1 ], [ 1 ], [ 1 ], [ 1 ] ]
 gap> DigraphAbsorptionProbabilities(ChainDigraph(4));
 [ [ 0, 0, 0, 1 ], [ 0, 0, 0, 1 ], [ 0, 0, 0, 1 ], [ 0, 0, 0, 1 ] ]
+gap> DigraphAbsorptionProbabilities(CycleDigraph(5));
+[ [ 1 ], [ 1 ], [ 1 ], [ 1 ], [ 1 ] ]
 gap> gr := ChainDigraph(250);;
 gap> probs := DigraphAbsorptionProbabilities(gr);;
 gap> scc := DigraphStronglyConnectedComponents(gr);;
@@ -2917,6 +2902,53 @@ gap> ForAll(probs,  # all zeros except for the sink
 >                       comp -> v[comp] = 0
 >                       or (v[comp] = 1 and scc.id[comp] = sink)));
 true
+gap> gr := Digraph([[1], []]);;
+gap> DigraphAbsorptionProbabilities(gr);
+[ [ 1, 0 ], [ 0, 1 ] ]
+
+# DigraphAbsorptionExpectedSteps
+gap> DigraphAbsorptionExpectedSteps(Digraph([[2], [3], [2]]));
+[ 1, 0, 0 ]
+gap> DigraphAbsorptionExpectedSteps(Digraph([]));
+[  ]
+gap> DigraphAbsorptionExpectedSteps(Digraph([[1]]));
+[ 0 ]
+gap> DigraphAbsorptionExpectedSteps(Digraph([[1, 2], [2]]));
+[ 2, 0 ]
+gap> DigraphAbsorptionExpectedSteps(ChainDigraph(5));
+[ 4, 3, 2, 1, 0 ]
+gap> DigraphAbsorptionExpectedSteps(CompleteDigraph(5));
+[ 0, 0, 0, 0, 0 ]
+
+# Absorption motivating example: game of 'Soccer Dice'
+gap> soccer := Digraph([[7, 2, 3, 5, 1, 4],
+>                       [7, 2, 3, 5, 5, 4],
+>                       [7, 5, 5, 5, 5, 1],
+>                       [7, 7, 2, 5, 5, 5],
+>                       [6, 6, 7, 7, 7, 4],
+>                       [], []]);;
+gap> DigraphStronglyConnectedComponents(soccer).comps;  # this ordering is used
+[ [ 7 ], [ 6 ], [ 1, 2, 3, 5, 4 ] ]
+gap> DigraphAbsorptionProbabilities(soccer) =
+>     [[3473 / 4493, 1020 / 4493, 0],
+>      [3365 / 4493, 1128 / 4493, 0],
+>      [3211 / 4493, 1282 / 4493, 0],
+>      [3471 / 4493, 1022 / 4493, 0],
+>      [2825 / 4493, 1668 / 4493, 0],
+>      [0, 1, 0],
+>      [1, 0, 0]];
+true
+gap> DigraphAbsorptionExpectedSteps(soccer);
+[ 13026/4493, 11868/4493, 10716/4493, 9510/4493, 6078/4493, 0, 0 ]
+
+# Absorption motivating example: a flea on the vertices of a dodecahedron.
+gap> gr := Digraph("dodecahedral");;
+gap> gr := DigraphRemoveEdges(gr, List(OutNeighbours(gr)[1], v -> [1, v]));
+<immutable digraph with 20 vertices, 57 edges>
+gap> DigraphAbsorptionExpectedSteps(gr)[1];
+0
+gap> DigraphAbsorptionExpectedSteps(gr)[2];
+19
 
 # Unbind local variables, auto-generated by etc/tst-unbind-local-vars.py
 gap> Unbind(A);