Saf zero subspaces iterator

.gitignore

@@ -9,5 +9,5 @@ doc/build/
 dist/
 surface_dynamics/flat_surfaces/origamis/origami_dense.c
 surface_dynamics/interval_exchanges/lyapunov_exponents.c
-surface_dynamics/interval_exchanges/iet_family.cpp
+surface_dynamics/interval_exchanges/iet_family_pyx.cpp
 surface_dynamics/interval_exchanges/integer_iet.c

doc/news/saf-zero.rst

@@ -0,0 +1,4 @@
+**Added:**
+
+* iterator over SAF zero subspaces
+* interface with libintervalxt for Boshernitzan criterion

doc/source/interval_exchanges/iet_family_pyx.rst

@@ -0,0 +1,4 @@
+Interval exchange families (Cython backend)
+===========================================
+
+.. automodule:: surface_dynamics.interval_exchanges.iet_family_pyx

doc/source/interval_exchanges/index.rst

@@ -23,4 +23,5 @@ Interval exchange transformations
    iet
    integer_iet
    iet_family
+   iet_family_pyx
    flip_sequence

doc/source/misc/additive_multivariate_generating_series.rst

@@ -0,0 +1,4 @@
+Additive multivariate generating series
+=======================================
+
+.. automodule:: surface_dynamics.misc.additive_multivariate_generating_series

doc/source/misc/even_permutation.rst

@@ -0,0 +1,4 @@
+Even permutation
+================
+
+.. automodule:: surface_dynamics.misc.even_permutation

doc/source/misc/factored_denominator.rst

@@ -0,0 +1,4 @@
+Factored denominator
+====================
+
+.. automodule:: surface_dynamics.misc.factored_denominator

doc/source/misc/generalized_multiple_zeta_values.rst

@@ -0,0 +1,4 @@
+Generalized multiple zeta values
+================================
+
+.. automodule:: surface_dynamics.misc.generalized_multiple_zeta_values

doc/source/misc/index.rst

@@ -4,6 +4,18 @@ Miscellaneous
 .. toctree::
    :maxdepth: 3
 
+   additive_multivariate_generating_series
    constellation
+   even_permutation
+   factored_denominator
+   generalized_multiple_zeta_values
    group_representation
+   iscc
+   latex
+   linalg
+   multiplicative_multivariate_generating_series
+   permutation
    plane_tree
+   ppl_utils
+   sql_db
+   statistics

doc/source/misc/iscc.rst

@@ -0,0 +1,4 @@
+iscc interface
+==============
+
+.. automodule:: surface_dynamics.misc.iscc

doc/source/misc/latex.rst

@@ -0,0 +1,4 @@
+Latex
+=====
+
+.. automodule:: surface_dynamics.misc.latex

doc/source/misc/linalg.rst

@@ -0,0 +1,4 @@
+Constellations
+==============
+
+.. automodule:: surface_dynamics.misc.linalg

doc/source/misc/multiplicative_multivariate_generating_series.rst

@@ -0,0 +1,4 @@
+Multiplicative multivariate generating series
+=============================================
+
+.. automodule:: surface_dynamics.misc.multiplicative_multivariate_generating_series

doc/source/misc/permutation.rst

@@ -0,0 +1,4 @@
+Permutation
+===========
+
+.. automodule:: surface_dynamics.misc.permutation

doc/source/misc/ppl_utils.rst

@@ -0,0 +1,4 @@
+ppl utilities
+=============
+
+.. automodule:: surface_dynamics.misc.ppl_utils

doc/source/misc/sql_db.rst

@@ -0,0 +1,4 @@
+sqlite database interface
+=========================
+
+.. automodule:: surface_dynamics.misc.sql_db

doc/source/misc/statistics.rst

@@ -0,0 +1,4 @@
+Statistics
+==========
+
+.. automodule:: surface_dynamics.misc.statistics

setup.py

@@ -70,10 +70,10 @@
         },
 
     # build iet_family only if pplpy is available
-    'iet_family': {
-        'name': 'surface_dynamics.interval_exchanges.iet_family',
+    'iet_family_pyx': {
+        'name': 'surface_dynamics.interval_exchanges.iet_family_pyx',
         'dir': 'interval_exchanges',
-        'sources': ['iet_family.pyx'],
+        'sources': ['iet_family_pyx.pyx'],
         'headers': [],
         'condition': WITH_PPL
         }

surface_dynamics/flat_surfaces/abelian_strata.py

@@ -1247,11 +1247,10 @@ def permutation_representative(self, left_degree=None, reduced=True, alphabet=No
 
         if reduced:
             from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
-            p = ReducedPermutationIET([l0, l1])
-
+            p = ReducedPermutationIET([l0, l1], reduced=True)
         else:
             from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
-            p = LabelledPermutationIET([l0, l1])
+            p = LabelledPermutationIET([l0, l1], reduced=False)
 
         p.alphabet(alphabet)
         return p
@@ -1283,9 +1282,9 @@ def lyapunov_exponents_approx(self, **kargs):
     # TODO
     # def carea(self)
 
-    def random_standard_permutation(self, nsteps=64):
+    def random_permutation(self, reduced=True, nsteps=None):
         r"""
-        Perform a random walk on rauzy diagram stopped on a standard permutation.
+        Return a random permutation obtained via a random walk in the extended rauzy diagram.
 
         INPUT:
 
@@ -1301,6 +1300,64 @@ def random_standard_permutation(self, nsteps=64):
 
         - bot rauzy move (proba 9/20)
 
+        EXAMPLES:
+
+            sage: from surface_dynamics import *
+
+            sage: C = AbelianStratum(10).hyperelliptic_component()
+            sage: p = C.random_permutation(); p   # random
+            0 1 2 3 4 5 6 7 8 9 10 11
+            11 10 9 8 7 6 5 4 3 2 1 0
+            sage: p.stratum_component()
+            H_6(10)^hyp
+
+
+        TESTS::
+
+            sage: from surface_dynamics import *
+            sage: from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
+            sage: from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
+            sage: for A in [AbelianStratum(6), AbelianStratum(1,1,1,1), AbelianStratum(3,3)]:
+            ....:     p = C.random_permutation(reduced=True)
+            ....:     assert isinstance(p, ReducedPermutationIET), C
+            ....:     assert p._labels is None, C
+            ....:     p = C.random_permutation(reduced=False)
+            ....:     assert isinstance(p, LabelledPermutationIET), C
+            ....:     assert p._labels is not None, C
+        """
+        import sage.misc.prandom as prandom
+
+        p = self.permutation_representative(reduced=True)
+        if nsteps is None:
+            nsteps = 64 * self.stratum().dimension()
+
+        for _ in range(nsteps):
+            rd = prandom.random()
+            if rd < 0.1:   # (inplace) symmetric with proba 1/10
+                p._inversed_twin()
+                p._reversed_twin()
+            elif rd < .55: # (inplace) rauzy move top with proba 9/20
+                p._move(1, 0, 1, p._twin[0][0])
+            else:          # (inplace) rauzy move bot with proba 9/20
+                p._move(0, 0, 0, p._twin[1][0])
+
+        p._check()
+        if reduced:
+            return p
+        if not reduced:
+            from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
+            return LabelledPermutationIET([p[0], p[1]], reduced=False)
+
+    def random_standard_permutation(self, reduced=True, nsteps=None):
+        r"""
+        Return a random standard permutation representative obtained by performing a random
+        walk in the Rauzy diagram.
+
+        INPUT:
+
+        - ``nsteps`` - integer or None - perform nsteps and then stops as soon
+          as a Strebel differential is found.
+
         EXAMPLES:
 
             sage: from surface_dynamics import *
@@ -1336,34 +1393,7 @@ def random_standard_permutation(self, nsteps=64):
             sage: p.stratum_component()
             H_17(32)^odd
         """
-        import sage.misc.prandom as prandom
-
-        p = self.permutation_representative()
-        if nsteps is None:
-            nsteps = 64 * self.stratum().dimension()
-        d = len(p)-1
-
-        for _ in range(nsteps):
-            rd = prandom.random()
-            if rd < 0.1:   # (inplace) symmetric with proba 1/10
-                p._inversed_twin()
-                p._reversed_twin()
-            elif rd < .55: # (inplace) rauzy move top with proba 9/20
-                p._move(1, 0, 1, p._twin[0][0])
-            else:          # (inplace) rauzy move bot with proba 9/20
-                p._move(0, 0, 0, p._twin[1][0])
-
-        while not p.is_standard():
-            rd = prandom.random()
-            if rd < 0.1:   # (inplace) symmetric with proba 1/10
-                p._inversed_twin()
-                p._reversed_twin()
-            elif rd < .55: # (inplace) rauzy move top with proba 9/20
-                p._move(1, 0, 1, p._twin[0][0])
-            else:          # (inplace) rauzy move bot with proba 9/20
-                p._move(0, 0, 0, p._twin[1][0])
-
-        return p
+        return self.random_permutation(reduced, nsteps).to_standard()
 
     def rauzy_diagram(self, *args, **kwds):
         r"""
@@ -2192,6 +2222,20 @@ def permutation_representative(self, left_degree=None, reduced=True, alphabet=No
             1
             sage: p.rauzy_diagram()
             Rauzy diagram with 20 permutations
+
+        TESTS::
+
+            sage: from surface_dynamics import AbelianStratum
+            sage: from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
+            sage: from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
+            sage: for A in [AbelianStratum(6), AbelianStratum(1,1,1,1), AbelianStratum(3,3)]:
+            ....:     for C in A.components():
+            ....:         p = C.permutation_representative(reduced=True)
+            ....:         assert isinstance(p, ReducedPermutationIET)
+            ....:         assert p._labels is None, C
+            ....:         p = C.permutation_representative(reduced=False)
+            ....:         assert isinstance(p, LabelledPermutationIET)
+            ....:         assert p._labels is not None, C
         """
         g = self._stratum.genus()
         n = self._stratum.nb_fake_zeros()
@@ -2238,15 +2282,16 @@ def permutation_representative(self, left_degree=None, reduced=True, alphabet=No
 
         if reduced:
             from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
-            p = ReducedPermutationIET([l0, l1])
-
+            p = ReducedPermutationIET([l0, l1], reduced=True)
         else:
             from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
-            p = LabelledPermutationIET([l0, l1])
+            p = LabelledPermutationIET([l0, l1], reduced=False)
         if alphabet is not None:
             p.alphabet(alphabet)
         elif relabel:
             p.alphabet(range(len(p)))
+
+        p._check()
         return p
 
     def rauzy_class_cardinality(self, left_degree=None, reduced=True):
@@ -2410,10 +2455,10 @@ def standard_permutations(self, reduced=True):
 
             if reduced:
                 from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
-                p = ReducedPermutationIET([l0, l1])
+                p = ReducedPermutationIET([l0, l1], reduced=True)
             else:
                 from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
-                p = LabelledPermutationIET([l0, l1])
+                p = LabelledPermutationIET([l0, l1], reduced=False)
 
             return [p]
 
@@ -2861,11 +2906,10 @@ def permutation_representative(self, left_degree=None, reduced=True, alphabet=No
 
         if reduced:
             from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
-            p = ReducedPermutationIET([l0, l1])
-
+            p = ReducedPermutationIET([l0, l1], reduced=True)
         else:
             from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
-            p = LabelledPermutationIET([l0, l1])
+            p = LabelledPermutationIET([l0, l1], reduced=False)
 
         if alphabet is not None:
             p.alphabet(alphabet)
@@ -3223,11 +3267,10 @@ def permutation_representative(self, left_degree=None, reduced=True, alphabet=No
 
         if reduced:
             from surface_dynamics.interval_exchanges.reduced import ReducedPermutationIET
-            p = ReducedPermutationIET([l0, l1])
-
+            p = ReducedPermutationIET([l0, l1], reduced=True)
         else:
             from surface_dynamics.interval_exchanges.labelled import LabelledPermutationIET
-            p = LabelledPermutationIET([l0, l1])
+            p = LabelledPermutationIET([l0, l1], reduced=False)
 
         if alphabet is not None:
             p.alphabet(alphabet)