utils: add stabilize_period and matric_formula_period

graph_tools/utils.py

@@ -3,7 +3,9 @@ def _init_():
     Singleton, OneOfKind, \
     UnionFind, DynamicLRU, \
     max_, select, prod, split_on_none, powerset, \
-    matrix_to_formula
+    matrix_to_formula, \
+    stabilize_period, matrix_formula_period, \
+    jit
 
   from itertools import chain, combinations
   from functools import reduce
@@ -12,6 +14,14 @@ def _init_():
   import weakref
   import inspect
 
+
+  try:
+    from numba import njit as jit
+  except:
+    print("Numba not found!")
+    def jit(fn): return fn
+
+
   def push_tests(cls):
     mod = inspect.getmodule(cls)
     if not hasattr(mod, "__test__"): setattr(mod, "__test__", {})
@@ -152,6 +162,8 @@ def _init_():
   def powerset(iterable):
       """Enumerate all subsets of `iterable`.
 
+      Vraci tuply!!! TODO
+
       EXAMPLE:
       >>> list(powerset([1,2,3]))
       [(), (1,), (2,), (3,), (1, 2), (1, 3), (2, 3), (1, 2, 3)]
@@ -161,6 +173,39 @@ def _init_():
       return chain.from_iterable(combinations(s, r) for r in range(len(s)+1))
 
 
+  def stabilize_period(X, vectors=False, shift=0):
+    return matrix_formula_period(
+      X['finalize'], X['step_matrix'], X['initial_vector'], vectors, shift)
+
+
+  def matrix_formula_period(l, A, r, vectors=False, shift=0):
+
+    from sage.all import matrix
+
+    assert all( x >= 0 for x in l.dict().values() )
+    assert all( x >= 0 for x in A.dict().values() )
+    assert all( x >= 0 for x in r.dict().values() )
+
+    def to_zero_one(M):
+      for (r, c), v in M.dict().items():
+        assert v > 0
+        M[r, c] = 1
+      return M
+
+    l = to_zero_one(matrix(l))
+    A = to_zero_one(matrix(A))
+    r = to_zero_one(matrix(r))
+
+    seen = []
+    while True:
+      for i in range(len(seen)):
+        if seen[i] == r:
+          return (i + shift,
+            [None]*shift + (seen if vectors else [min(1, (l*v)[0, 0]) for v in seen ]))
+
+      seen.append(r)
+      r = to_zero_one(A*r)
+
   def matrix_to_formula(l, A, r, k=None, base_ring=None):
     """Tranform matrix formula l A^k r to simple formula.