possibly improve the iterator for PermutationGroups

[moritz]

Currently, the iterator looks like this.

def __iter__(self):

    from sage.sets.recursively_enumerated_set import RecursivelyEnumeratedSet
    return iter(RecursivelyEnumeratedSet(seeds=[self.one()],
                successors=lambda g: (g._mul_(h) for h in self.gens())))

Perhaps it would be faster to use a strong_generating_system to iterate through the group.

[moritz]

I added an alternative suggestion __iter_alt__ in the patch. However, sometimes calculation the strong_generating_system takes much longer than actually iterating with the standard iterator above. Here are some timing experiments:

sage: p = [(i,i+1) for i in range(1,601,2)]
sage: q = [tuple(range(1+i,601,3)) for i in range(3)]
sage: A = PermutationGroup([p,q])
sage: A.cardinality()
60000
sage: time len(list(A))
CPU times: user 441 ms, sys: 58 ms, total: 499 ms
Wall time: 478 ms
60000
sage: time len(list(A.__iter_alt__()))
CPU times: user 3.39 s, sys: 690 ms, total: 4.08 s
Wall time: 4.1 s
60000
sage: P = PermutationGroup(SymmetricGroup(10).gens())
sage: time len(list(P))
CPU times: user 9.06 s, sys: 182 ms, total: 9.24 s
Wall time: 9.24 s
3628800
sage: time len(list(P.__iter_alt__()))
CPU times: user 3.64 s, sys: 144 ms, total: 3.79 s
Wall time: 3.79 s
3628800

Any suggestions? What are good other groups to do benchmarking?

---

New commits:

​74488fe

a suggestion for an alternative iterator

[moritz]

More timings:

sage: P = CoxeterGroup(["E",7], implementation="permutation")
sage: time len(list(P.iteration(algorithm="depth", tracking_words=False)))
CPU times: user 4.55 s, sys: 525 ms, total: 5.07 s
Wall time: 5.14 s
2903040
sage: time len(list(P.__iter__()))
CPU times: user 21.2 s, sys: 693 ms, total: 21.8 s
Wall time: 21.8 s
2903040
sage: time len(list(P.__iter_alt__()))
CPU times: user 4.33 s, sys: 231 ms, total: 4.56 s
Wall time: 4.58 s
2903040

[tscrim]

You can probably get some more speed by doing the initial step in the __iter__ function and then do the remaining steps with SGS being explicitly given. It is probably faster to also do

            S = SGS.pop()
            if not SGS:
                for g in S:
                    yield g
            else:
                for s in elements(self, SGS):
                    for g in S:
                        yield s._mul_(g)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3594002

directly yield in the last step

[moritz]

Replying to tscrim:

You can probably get some more speed by doing the initial step in the __iter__ function and then do the remaining steps with SGS being explicitly given.

I don't quite follow. We would still have to compute SGS in this case, won't we? This seems to be the slow part, at least in the first example above with 60000 elements.

It is probably faster to also do ...

That makes it a little faster, thanks!

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3dc4822

removing the if from the loop

[stumpc5]

+        def elements(self, SGS=None):
+            S = SGS.pop()

seems weird. I assume that you want to remove the =None ?

[moritz]

Replying to stumpc5:

+        def elements(self, SGS=None):
+            S = SGS.pop()

seems weird. I assume that you want to remove the =None ?

yes, and I just did. (Before seeing, that you did the same thing..)

---

New commits:

​7780abc

remove 'None'

[stumpc5]

I did some more, could you also merge that change in to the branch?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​07faa6b	trivial simplifications!
​7d3d68f	Merge branch 'u/stumpc5/enumerate_PermutationGroup' of git://trac.sagemath.org/sage into enumerate_PermutationGroup

[vdelecroix]

In def elements(self, SGS): the argument self is never used and you can remove it.

[stumpc5]

sage: time len(list(A.__iter_alt__()))
CPU times: user 3.39 s, sys: 690 ms, total: 4.08 s
Wall time: 4.1 s
60000

This computation takes much less time in gap directly, I thus expect that doing

SGS = self.strong_generating_system(base)

directly in gap rather than in sage would improve the iteration in this case drastically. (I just tried for a bit but failed. The code in gap directly is

#############################################################################
##
#F  ElementsStabChain(<S>)
##
InstallGlobalFunction(ElementsStabChain,function ( S )
    local   elms,               # element list, result
            stb,                # elements of the stabilizer
            pnt,                # point in the orbit of <S>
            rep;                # inverse representative for that point

    # if <S> is trivial then it is easy
    if Length(S.generators) = 0  then
        elms := [ S.identity ];

    # otherwise
    else

        # start with the empty list
        elms := [];

        # compute the elements of the stabilizer
        stb := ElementsStabChain( S.stabilizer );

        # loop over all points in the orbit
        for pnt  in S.orbit  do

           # add the corresponding coset to the set of elements
           rep := S.identity;
           while S.orbit[1] ^ rep <> pnt  do
                rep := LeftQuotient( S.transversal[pnt/rep], rep );
           od;
           UniteSet( elms, stb * rep );

        od;

   fi;

   # return the result
   return elms;
end);

[stumpc5]

Replying to vdelecroix:

In def elements(self, SGS): the argument self is never used and you can remove it.

done.

I also made some tests and think that the permutation group method strong_generating_system is much slower than necessary (see the above example where this algo is slower than the brute force). In particular, there is a lot of talking back and forth between sage and gap. Anyway, for my purposes, this is a non-issue. @Dima, do you know how to easily make this method fast?

---

New commits:

​285524e

removed self from elements arguments

[tscrim]

Dima, see comment:18.

[stumpc5]

posted question and code at ​https://groups.google.com/d/topic/sage-devel/wivuFBq7n7k/discussion.

[stumpc5]

def SGS_old(G): # taken from strong_generating_system
    sgs = []
    stab = G
    base = G.base()
    for j in base:
        sgs.append(stab.transversals(j))
        stab = stab.stabilizer(j)
    return sgs

def SGS_new(G, gap_output=True):
    S = G._gap_().StabChain()
    gap.execute(gap_cosets)
    cosets = S.CosetsStabChain()
    if gap_output:
        return cosets                           #  <-------

    elt_class = G._element_class()
    gap2sage = lambda elt: elt_class(elt, G, check=False)
    return [ [ gap2sage(elt) for elt in coset]  #  <-------
             for coset in cosets ]              #  <-------

gap_cosets = """CosetsStabChain := function ( S )
    local   cosets,             # element list, result
            new_coset,          #
            pnt,                # point in the orbit of <S>
            rep;                # inverse representative for that point

    # if <S> is trivial then it is easy
    if Length(S.generators) = 0  then
        cosets := [ [S.identity] ];

    # otherwise
    else

        # compute the elements of the stabilizer
        cosets := CosetsStabChain( S.stabilizer );

        # loop over all points in the orbit
        new_coset := [];
        for pnt  in S.orbit  do

            # add the corresponding coset to the set of elements
            rep := S.identity;
            while S.orbit[1] ^ rep <> pnt  do
                 rep := LeftQuotient( S.transversal[pnt/rep], rep );
            od;
            Add( new_coset, rep );
        od;
        Add( cosets, new_coset );
   fi;

   # return the result
   return cosets;
end;
"""

[stumpc5]

The new method is much faster than the old, but then slowed down by casting the output to sage:

sage: sage: p = [(i,i+1) for i in range(1,601,2)]
....: sage: q = [tuple(range(1+i,601,3)) for i in range(3)]
....: sage: A = PermutationGroup([p,q])

sage: %time XX = SGS_old(A)
CPU times: user 3.67 s, sys: 664 ms, total: 4.34 s
Wall time: 4.45 s

sage: %time XX = SGS_new(A, True)
CPU times: user 4.94 ms, sys: 712 µs, total: 5.65 ms
Wall time: 14 ms

sage: %time XX = SGS_new(A, False)
CPU times: user 2.48 s, sys: 116 ms, total: 2.59 s
Wall time: 2.61 s

[vdelecroix]

You should be using libgap and not gap. The communication would be much faster.

[vdelecroix]

And why not using StrongGeneratorsStabChain that is builtin in GAP?

[stumpc5]

You should be using libgap and not gap

Could you provide a code snippet in the above example of how to do that?

And why not using StrongGeneratorsStabChain? that is builtin in GAP?

Because this provides a list of strong generators for the stabilizer chain, but not a list of coset reprs:

gap> P := Group([(1,2),(1,2,3,4)]);
Group([ (1,2), (1,2,3,4) ])
gap> S := StabChain(P);
<stabilizer chain record, Base [ 1, 2, 3 ], Orbit length 4, Size: 24>
gap> CosetsStabChain(S);
[ [ () ], [ (), (3,4) ], [ (), (2,4,3), (2,3,4) ], [ (), (1,4)(2,3), (1,2)(3,4), (1,3)(2,4) ] ]
gap> StrongGeneratorsStabChain(S);
[ (3,4), (2,3,4), (1,2)(3,4), (1,4)(2,3) ]

[stumpc5]

And why not using...

Let me add that the gap code above is simply taken from the ElementsStabChain code from gap (see PageSource(ElementsStabChain)) and returning the cosets rather than the elements. If there was a builtin function that does this, I would be surprised if it is not used there (though it might be possible for speed reasons).

[dimpase]

Historically, little attention was paid to the problem at hand, and indeed, iteration over all the elements of a group G, rather than doing something more intelligent, is almost surely an indication of something done without enough thought.

E.g.,

• would knowing the coset representatives of a subgroup H<G and generators for H suffice?
• would knowing the words in generators rather than permutations suffice?
• can an adavntage be taken of a short ​base for G?

[dimpase]

Replying to vdelecroix:

You should be using libgap and not gap. The communication would be much faster.

Yeah, if the iterator at hand uses (pexpect) GAP, merely converting it to libGAP would be quite a speedup, I suppose.

[stumpc5]

Well, let me wrap up the discussion so for:

• Moritz implemented an iteration using a short base and a stabilizer chain which are both implemented in gap (using the Schreier-Sims algorithm), see the attached branch. This already speeded the iteration in most situations drastically.

• Permutation group elements are faster multiplied in sage than in gap, so we want to do it there.

• The remaining bottleneck is the computation of the strong generating system for which I provided code in gap. This is now also really fast, except for the translation of the gap permutations to sage. I still need help to get this last issue disappear.

[dimpase]

GAP has a function Elements() which enumerates the elements of a group, did you try it directly? Probably you're reinventing the wheel...

I cannot believe that permutation multiplication in (lib)GAP is slower than in Sage. You are probably doing lots of slow conversions between Sage and pexpect GAP.

[stumpc5]

Before comparing any of these functions, we should have a way to fast translate gap permutations to sage permutations. Otherwise, this would even much more so spoil the use of the gap function Elements which then results in all permutations to be translated from gap to sage.

Here are some preliminary timings for the Coxeter group of type E7:

In gap, the iteration needs 9.2 secs:

gap> P := Group([(1,64)(3,12)(8,19)(15,23)(16,20)(21,26)(25,30)(27,33)(29,35)(31,34)(32,40)(36,41)(37,44)(38,43)(45,47)(49,52)(62,63)(66,75)(71,82)(78,86)(79,83)(84,89)(88,93)(90,96)(92,98)(94,97)(95,103)(99,104)(100,107)(101,106)(108,110)(112,115)(125,126), (2,65)(4,9)(8,15)(13,18)(14,24)(16,21)(19,23)(20,26)(22,28)(25,31)(27,36)(30,34)(33,41)(50,55)(54,56)(57,59)(58,60)(67,72)(71,78)(76,81)(77,87)(79,84)(82,86)(83,89)(85,91)(88,94)(90,99)(93,97)(96,104)(113,118)(117,119)(120,122)(121,123), (1,12)(3,66)(4,8)(9,15)(13,16)(14,25)(18,21)(22,27)(24,31)(28,36)(35,39)(40,42)(43,46)(44,48)(47,51)(52,53)(61,62)(64,75)(67,71)(72,78)(76,79)(77,88)(81,84)(85,90)(87,94)(91,99)(98,102)(103,105)(106,109)(107,111)(110,114)(115,116)(124,125), (2,9)(3,8)(4,67)(5,13)(11,14)(12,19)(17,22)(21,29)(26,35)(31,32)(34,40)(36,38)(41,43)(48,50)(51,54)(53,57)(60,61)(65,72)(66,71)(68,76)(74,77)(75,82)(80,85)(84,92)(89,98)(94,95)(97,103)(99,101)(104,106)(111,113)(114,117)(116,120)(123,124), (4,13)(5,68)(6,11)(8,16)(9,18)(10,17)(15,21)(19,20)(23,26)(32,37)(38,45)(40,44)(42,48)(43,47)(46,51)(57,58)(59,60)(67,76)(69,74)(71,79)(72,81)(73,80)(78,84)(82,83)(86,89)(95,100)(101,108)(103,107)(105,111)(106,110)(109,114)(120,121)(122,123), (5,11)(6,69)(7,10)(13,14)(16,25)(18,24)(20,30)(21,31)(26,34)(29,32)(35,40)(39,42)(45,49)(47,52)(51,53)(54,57)(56,59)(68,74)(70,73)(76,77)(79,88)(81,87)(83,93)(84,94)(89,97)(92,95)(98,103)(102,105)(108,112)(110,115)(114,116)(117,120)(119,122), (6,10)(7,70)(11,17)(14,22)(24,28)(25,27)(30,33)(31,36)(32,38)(34,41)(37,45)(40,43)(42,46)(44,47)(48,51)(50,54)(55,56)(69,73)(74,80)(77,85)(87,91)(88,90)(93,96)(94,99)(95,101)(97,104)(100,108)(103,106)(105,109)(107,110)(111,114)(113,117)(118,119)]);
<permutation group with 7 generators>
gap> Length(Elements(P)); time;
2903040
9255

while in sage, it only takes 4.3 secs of which 0.7 secs are spend constructing the strong generating system:

sage: P = CoxeterGroup(["E",7], implementation="permutation")
sage: %time X = P.strong_generating_system(P.base())
CPU times: user 632 ms, sys: 71.6 ms, total: 703 ms
Wall time: 718 ms

sage: P = CoxeterGroup(["E",7], implementation="permutation")
sage: %time len(list(P.__iter_alt__()))
CPU times: user 4.11 s, sys: 208 ms, total: 4.32 s
Wall time: 4.34 s
2903040

So I expect this to go down below 4 sec if we speed the strong generating system.

[dimpase]

Do you know that GAP actually has iterators? With P as in comment 31:

gap> si:=IteratorStabChain(StabChain(P));
<iterator>
gap> l:=[];;
gap> for i in si do Add(l,i); od;
gap> Length(l);

is about 30% faster than Elements(P) - the latter actually will also sort the elements, perhaps that's why.

Moreover, in for i loop about 60% of time is taken by Add(). (I measured this by replacing Add() by a counter increment).

[stumpc5]

Oh, I forgot that Elements sorts the output, thanks for the reminder and for improving the timing comparisons:

gap> si:=IteratorStabChain(StabChain(P));; time;
4
gap> i:=0;;
gap> for j in si do i:=i+1; od; time;
2608
gap> i;
2903040

vs.

sage: i = 0
sage: %time for _ in P.__iter_alt__(): i+=1
CPU times: user 1.97 s, sys: 80.2 ms, total: 2.05 s
Wall time: 2.08 s
sage: i
2903040

(Is this now the timing you think is optimal in gap?) If this is the case, sage is still faster, and, again, 0.7secs are overhead that can be remomved from the sage timing, pushing it to roughly 1.5secs.

[dimpase]

I wonder whether the SGSs in these tests use the same base.

I must say that I also don't understand the recursive code in the branch on the ticket. I suspect it won't work correctly if SGS's elements are not involutions. (E.g. if the group is cyclic, then there should be a loop over the powers of the only generator in SGS, no?)

[stumpc5]

The algorithm in sage now is the very same as the one in gap, there is no structural difference. For a general description, see for example ​https://math.stackexchange.com/a/1706006/234820.

In particular, both use the same base and the same stabilizer chain.

I suspect it won't work correctly if SGS's elements are not involutions.

This is not correct if I understand it correctly. In particular, the elements in the SGS are usually not only involutions, for example

sage: S = PermutationGroup([(1,2),(1,2,3)])
sage: S.strong_generating_system()
[[(), (1,2,3), (1,3,2)], [(), (2,3)], [()]]

and the SGS is not a generating set for the group itself. Another example of a cyclic group is

sage: P = PermutationGroup([(1,2,3,4)])
sage: for elt in P.__iter_alt__(): print elt
()
(1,4,3,2)
(1,3)(2,4)
(1,2,3,4)

and the stabilizer chain is trivial (i.e., the chain of cosets is a single coset of the trivial subgroup).

[dimpase]

ok, so these are for g in S loops that iterate over subgroups S, right?

What iterator is there? Do you want to use your iterator there, explicitly?

[stumpc5]

Sorry, I am slightly lost now in what you are asking -- let me try to answer anyway.

Let S be a permutation group in sage. The current iterator for g in S iterates by a brute force algorithm through the elements keeping the complete group in memory to check if an element was seen before.

The new algorithm reimplements the algorithm used in gap, and indeed takes the important ingredients from gap. gap computes a chain e = S1 in S2 in ... in Sk = S of subgroups (obtained using the Schreier-Sims algorithm) for which the cosets of Si inside S{i+1} are particularly easy to compute. This yields a sequence of cosets [e]=C1, C2, ..., Ck such that every element in S has a unique representation g = g1*...*gk for gi in Ci.

In sage, the new algorithm is to

• directly compute in gap the list C1, C2, ..., Ck
• translate the elements in all these lists to sage (this takes too long and should be improved)
• produce an iterator for all elements in S by iterating over all products of elements in these lists.

Let me know if you meant something else...

[dimpase]

Replying to stumpc5:

Let me know if you meant something else...

I meant

​lines 844 and ​848

in your patch, both say for g in S: -- how do these iterations work? (I imagine S is a subgroup in the SGS chain).

[stumpc5]

in your patch, both say for g in S: -- how do these iterations work? (I imagine S is a subgroup in the SGS chain)

No, S = SGS.pop() so S = Ci is the next layer of coset representatives. Let us go through the code with the notation as above:

        def elements(SGS):
            S = SGS.pop()
            if not SGS:
                for g in S:
                    yield g
            else:
                for s in elements(SGS):
                    for g in S:
                        yield s._mul_(g)

Given such a sequence SGS = C1,C2,...,Ci, elements(SGS) iterates over all elements in Si. This is done by recursively doing the computation for C1,C2,..,C{i-1} and then yielding all elements in S{i-1} multiplied with all elements in S=Ci.

[dimpase]

Well, by right, SGS is a generating set. So your naming scheme (no other documentation :-)) is misleading...

[stumpc5]

SGS is a strong generating system taken from the sage naming scheme. The later might be misleading, though, since it is not a generating set but a sequence of coset reprs...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​8f70da5

speeded strong_generating_system

[stumpc5]

with the new commit, the timings from the beginning become

sage: p = [(i,i+1) for i in range(1,601,2)]
....: q = [tuple(range(1+i,601,3)) for i in range(3)]
....: A = PermutationGroup([p,q])
....: A.cardinality()
....: 
60000
sage: time len(list(A))
CPU times: user 458 ms, sys: 56.2 ms, total: 514 ms
Wall time: 512 ms
60000
sage: time len(list(A.__iter_alt__()))
CPU times: user 2.97 s, sys: 167 ms, total: 3.13 s
Wall time: 3.23 s
60000
sage: time x = A.strong_generating_system(A.base())
CPU times: user 4.51 s, sys: 396 ms, total: 4.9 s
Wall time: 5.08 s
sage: time x = A.strong_generating_system(implementation="gap")
CPU times: user 2.67 s, sys: 80.3 ms, total: 2.75 s
Wall time: 2.77 s

so the iterations are equally fast, but the translation from gap permutations to sage permutations takes 2.5 secs.

[dimpase]

Do you have a problem using libGAP rather than GAP?

[stumpc5]

You should be using libgap and not gap

Could you provide a code snippet in the above example of how to do that?

No, I just do now know how to do it, see my above comment from comment:25 .

[dimpase]

Replying to stumpc5:

You should be using libgap and not gap

Could you provide a code snippet in the above example of how to do that?

No, I just do now know how to do it, see my above comment from comment:25 .

You can have a look at ​CossidentePenttilaGraph() for the use of libgap.function_factory (so that you can basically call your own GAP function via libgap).

So you can change your code as follows:

gap_cosets = libgap.function_factory("""function ( S0 )
    local   CosetsStabChain;
    CosetsStabChain := function(S)  # for the recursive call
    local   cosets,             # element list, result
            new_coset,          #
            pnt,                # point in the orbit of <S>
            rep;                # inverse representative for that point

    # if <S> is trivial then it is easy
    if Length(S.generators) = 0  then
        cosets := [ [S.identity] ];

    # otherwise
    else

        # compute the elements of the stabilizer
        cosets := CosetsStabChain( S.stabilizer );

        # loop over all points in the orbit
        new_coset := [];
        for pnt  in S.orbit  do

            # add the corresponding coset to the set of elements
            rep := S.identity;
            while S.orbit[1] ^ rep <> pnt  do
                 rep := LeftQuotient( S.transversal[pnt/rep], rep );
            od;
            Add( new_coset, rep );
        od;
        Add( cosets, new_coset );
   fi;

   # return the result
   return cosets;
   end;
   return CosetsStabChain(S0);
end;""")

Now one can do

sage: G=libgap.AlternatingGroup(4)
sage: S=G.StabChain()
sage: gap_cosets(S)
[ [ () ], [ (), (2,4,3), (2,3,4) ], [ (), (1,4,3), (1,3,4), (1,2,4) ] ]

another useful thing is syntax like G=libgap(PermutationGroup([(1,2)])) to convert from Sage's group (implemented with pexpect GAP, I gather) to libgap.

[stumpc5]

Thank you, Dima! I did these changes together with tons of small changes that had a huge impact on the running time.

Beside this, I arrived at the following, and wonder if this is okay doing: Given I have

• a PermutationGroupElement "sample" and
• a list "new_perm" representing a permutation in one-line notation in the *same permutation group* as "sample"

It seems that, oddly, we do not have a fast way of generating new permutations from these lists, despite the fact that we can multiply permutations very fast. I thus provide a new PermutationGroupElement method _generate_new that basically takes a copy from the given PermutationGroupElement and only overwrites its underlying c list. The code is basically copied from _mul_.

An example would be

sage: S = PermutationGroup(SymmetricGroup(4).gens())
sage: sample = S([2,3,4,1])
sage: lst1 = []
sage: lst2 = [2,1]
sage: lst3 = [4,3,2,1]
sage: pi1 = sample._generate_new(lst1); pi1
()
sage: pi2 = sample._generate_new(lst2); pi2
(1,2)
sage: pi3 = sample._generate_new(lst3); pi3
(1,4)(2,3)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​f5a48dd

new permutation group element method _generate_new, improved speed for strong generating system

[stumpc5]

And here is the new timings:

sage: sage: p = [(i,i+1) for i in range(1,601,2)]
....: ....: q = [tuple(range(1+i,601,3)) for i in range(3)]
....: ....: A = PermutationGroup([p,q])
....: ....: A.cardinality()
....: 
60000
sage: sage: time len(list(A))
CPU times: user 514 ms, sys: 52.2 ms, total: 566 ms
Wall time: 535 ms
60000
sage: sage: time len(list(A.__iter_alt__()))
CPU times: user 588 ms, sys: 91 ms, total: 679 ms
Wall time: 669 ms
60000
sage: sage: time len(list(A))
CPU times: user 469 ms, sys: 72.5 ms, total: 542 ms
Wall time: 529 ms
60000
sage: sage: time len(list(A.__iter_alt__()))
CPU times: user 375 ms, sys: 52.3 ms, total: 427 ms
Wall time: 419 ms
60000
sage: sage: time x = A.strong_generating_system(A.base())
CPU times: user 3.82 s, sys: 784 ms, total: 4.6 s
Wall time: 4.69 s
sage: sage: time x = A.strong_generating_system(implementation="gap")
CPU times: user 281 ms, sys: 127 µs, total: 281 ms
Wall time: 280 ms

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​762d318

added doc for _generate_new

[stumpc5]

I set this to "needs review" for a technical review of the code where I still hope for input for improvements (which will not only be relevant for my particular purpose!).

Once this is settled, I propose to use the now PermutationGroup iterator per default (and I would also be interested in examples where the old iterator is still faster).

[stumpc5]

Comparing the gap and the sage iterators from comment:33, we now have:

gap> si:=IteratorStabChain(StabChain(P));
<iterator>
gap> i:=0;;
gap> for j in si do i:=i+1; od; time;
2523
gap> i;
2903040

vs.

sage: sage: i = 0
....: sage: %time for _ in P.__iter_alt__(): i+=1
....: 
CPU times: user 1.49 s, sys: 28.2 ms, total: 1.52 s
Wall time: 1.59 s
sage: i
2903040

[dimpase]

The right thing to do (not on this ticket, I guess it's a nontrivial amount of work) is to convert the whole sage/groups/perm_gps to libGAP.

This ticket shows that the performance improvements are big.

[stumpc5]

This ticket shows that the performance improvements are big

Let me remark that the biggest performance boost came from generating new permutation group elements from old using the new _generate_new method. The use of libGAP also played a role, though.

[dimpase]

1) Regarding the timing (GAP iterator vs your iterator). Am I right that they both operate with (lib)GAP data? (i.e. that mul call on permutations in __iter_alt__, is it actually in (lib)GAP?)

2) I don't understand the hack one._generate_new(libgap.ListPerm?(elt).sage()). I thought that libgap has a fast method to convert GAP permutations to Sage. E.g.

sage: G=libgap.AlternatingGroup(4)
sage: gg=G.GeneratorsOfGroup()
sage: gg.sage()
[(1,2,3), (2,3,4)]
sage: type(g.sage()[0])
<type 'sage.groups.perm_gps.permgroup_element.PermutationGroupElement'>
sage: type(g[0])
<type 'sage.libs.gap.element.GapElement_Permutation'>

[stumpc5]

i.e. that mul call on permutations in iter_alt, is it actually in (lib)GAP?

No, this is pure cython and uses PermutationGroupElement._mul_(left, _right), no (lib)GAP involved here.

Since both algorithms are the same, and the strong generating set in the Sage version is taken through libGAP, I suspect that the time advantage of Sage really comes from this superfast cython multiplication of permutation group elements.

I thought that libgap has a fast method to convert GAP permutations to Sage

This is not the case: libgap uses GapElement_Permutation.sage which returns PermutationGroupElement(libgap.ListPerm(self).sage()), so it goes through the ordinary permutation group element constructor that is very slow because of tons of checks. This is why I wrote _generate_new that just uses the data from a given element and overwrites its c array containing the permutation:

sage: P = PermutationGroup([(1,2),(1,2,3,4)])
sage: one = P.one()
sage: l = [4,3,2,1]
sage: %timeit P(l)
1000 loops, best of 3: 579 µs per loop
sage: %timeit one._generate_new(l)
The slowest run took 39.93 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 203 ns per loop

[dimpase]

Replying to stumpc5:

I thought that libgap has a fast method to convert GAP permutations to Sage

This is not the case: libgap uses GapElement_Permutation.sage which returns PermutationGroupElement(libgap.ListPerm(self).sage()), so it goes through the ordinary permutation group element constructor that is very slow because of tons of checks. This is why I wrote _generate_new that just uses the data from a given element and overwrites its c array containing the permutation:

Shouldn't then this be used to improve the code in src/sage/libs/gap/? Cause that's what it does, you provide a better implementation for <libgap perm>.sage(), right?

[stumpc5]

Not really, because my method (A) makes use of the fact that one knows the permutation group the element lives in a priori, and (B) does not check the input whatsoever for correctness (and might result in segfaults if not used properly).

The first is an issue because <libgap perm> does not know which parent a gap permutation lives in (and neither does gap iirc), while the second might be neglected in case we expect that we always get something reasonable fromo gap.

[dimpase]

OK, I see. Still, libgap should be using the low-level C array with GAP permutation data, rather than converting GAP list -> Sage list -> PermutationGroupElement? (creating, if needed, a parent permutation group along the way, no?)

I don't understand whether this involves calling pexpect GAP or not.

[stumpc5]

I don't understand whether this involves calling pexpect GAP or not.

What involved pexpect GAP? This creation of a PermutationGroupElement? Then the answer is no, there is no gap object attached by default to a PermutationGroupElement.

[dimpase]

Replying to stumpc5:

I don't understand whether this involves calling pexpect GAP or not.

What involved pexpect GAP? This creation of a PermutationGroupElement? Then the answer is no, there is no gap object attached by default to a PermutationGroupElement.

OK, thanks - it's not immediately clear with all that methods using _gap_ in that class.

Anyhow, I don't understand why libGAP creates PermutationGroupElement rather than a Permutation. It would be natural to me to have this as the default, mostly as it does not involve constructing a parent (GAP permutations don't really have parents in GAP, they are all in every SymmetricGroup() of degree at least LargestMovedPoint() of the permutation. Then your hack can be used for fast libgap->Sage conversion of permutations.

[stumpc5]

It would be natural to me to have this as the default

I agree, this would get down the speed quite a bit

sage: P = PermutationGroup([(1,2),(1,2,3,4,5,6,7,8)])
sage: l = [8,7,6,5,4,3,2,1]
sage: %timeit P(l, check=False)
The slowest run took 7.58 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 44.3 µs per loop
sage: %timeit Permutation(l, check_input=False)
The slowest run took 7.85 times longer than the fastest. This could mean that an intermediate result is being cached.
100000 loops, best of 3: 6.74 µs per loop

but my hack would not work in this case because it uses the c structure underlying PermutationGroupElement which doesn't exist (caution: haven't rechecked!) for Permutation, and this hack is possibly still an order of magnitude faster

sage: one = P.one()
sage: %timeit one._generate_new(l)
The slowest run took 128.70 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 241 ns per loop

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​0986e67	cosmetic change to SGS
​952b292	yippie -- I got the strong generating system computation really fast now by directly converting the gap list of c ints to an <int*>!

[stumpc5]

Okay, now everything is as fast as I was expecting it to be possible! I now do not create a sage list from the libgap list and then fill a <int*> with this data, but create the latter directly from the libgap list entries. This removes all the overhead!

Here are the new timings:

sage: p = [(i,i+1) for i in range(1,601,2)]
....: q = [tuple(range(1+i,601,3)) for i in range(3)]
....: A = PermutationGroup([p,q])
....: A.cardinality()
....: 
60000
sage: i=0
sage: %time for _ in A.__iter__(): i+=1
CPU times: user 512 ms, sys: 40.3 ms, total: 552 ms
Wall time: 543 ms
sage: i
60000
sage: i=0
sage: %time for _ in A.__iter_alt__(): i+=1
CPU times: user 132 ms, sys: 23.7 ms, total: 156 ms
Wall time: 137 ms
sage: i
60000
sage: %time X = A.strong_generating_system(A.base())
CPU times: user 3.5 s, sys: 675 ms, total: 4.17 s
Wall time: 4.22 s
sage: [ len(x) for x in X ]
[600, 100]
sage: %time Y = A.strong_generating_system(implementation="gap")
CPU times: user 38.3 ms, sys: 129 µs, total: 38.5 ms
Wall time: 38.2 ms
sage: [ len(y) for y in Y ]
[1, 100, 600]

[stumpc5]

I consider it wrong to define (as I did) the method _generate_perm(self, old) for libGAP_List (handing over a PermutationGroupElement old. I think it should rather exist a method _generate_new_libgap(self, gaplst for PermutationGroupElement that takes as input a libGAP_List gaplst.

I was not able to do that due to import errors and my ignorance how to properly cast the objects into the correct types, but this seems only a matter of proper coding in cython. @Travis or @Dima, could you help with this?

I think once this is settled, the code is ready for review!

@Moritz: rather incredible speedup compared to your original call of strong_generating_system, hm?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​d49798c

removed todo

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​2be32af

removed todo

[tscrim]

I can probably try it late tomorrow or the day after (I promised Darij I would review #25326 first). So Dima, let me know if you are going to work on this tomorrow.

[moritz]

Replying to stumpc5:

@Moritz: rather incredible speedup compared to your original call of strong_generating_system, hm?

Yes, absolutely great! As I see it, the new __iter_alt_ method is now faster then __iter__ in all cases, right? Therefore I propose to replace __iter__ with the __iter_alt__ before puuting this ticket on review.

[stumpc5]

No, not quite. The new __iter_alt__ is still slower for small groups (size ~ <500 elts). I still propose to use ___iter_alt__ per default, but leave the option to use the brute force iteration if desired (in case you want to loop through many small groups.

The issue is that if the group is really small, waiting even for a tiny information from libgap is slower than exploring the group directly.

[dimpase]

Replying to tscrim:

I can probably try it late tomorrow or the day after (I promised Darij I would review #25326 first). So Dima, let me know if you are going to work on this tomorrow.

I suggest that the documentation says that libGAP (rather than just GAP) is used. That's all for me for the time being (swamped with other things).

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​4482025

Merge branch 'develop' into t/25477/enumerate_PermutationGroup

[dimpase]

What is the current state here?

[stumpc5]

Replying to stumpc5:

I consider it wrong to define (as I did) the method _generate_perm(self, old) for libGAP_List (handing over a PermutationGroupElement old. I think it should rather exist a method _generate_new_libgap(self, gaplst for PermutationGroupElement that takes as input a libGAP_List gaplst.

I was not able to do that due to import errors and my ignorance how to properly cast the objects into the correct types, but this seems only a matter of proper coding in cython. @Travis or @Dima, could you help with this?

I am still waiting for someone to have a look at this...

In principle, I think this implementation is a strong improvement to permutation groups in Sage, and since I am not very good in optimizing in cython, there might a few simple improvements left to be found. I don't think it is ready for a positive review, but it certainly is for a review of the technical details.

[dimpase]

Replying to stumpc5:

Replying to stumpc5:

I consider it wrong to define (as I did) the method _generate_perm(self, old) for libGAP_List (handing over a PermutationGroupElement old. I think it should rather exist a method _generate_new_libgap(self, gaplst for PermutationGroupElement that takes as input a libGAP_List gaplst.

I was not able to do that due to import errors and my ignorance how to properly cast the objects into the correct types, but this seems only a matter of proper coding in cython. @Travis or @Dima, could you help with this?

I am still waiting for someone to have a look at this...

Could you post a branch with this import problem? It'd be much quicker to see what's going on this way? (it need not become the new branch on this ticket, certainly)

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​85a01c5

not working attempt to move the permutation generation from libgap_list to PermutationGroupElement

[stumpc5]

Could you post a branch with this import problem? It'd be much quicker to see what's going on this way?

I had deleted that attempt, but it was basically what I just pushed. My problem was that I didn't see how to handle libGAP_INT_INTOBJ, libGAP_ELM_LIST inside permgroup_element.pyx.

---

New commits:

​85a01c5

not working attempt to move the permutation generation from libgap_list to PermutationGroupElement

[dimpase]

for starters, I think you need

+++ b/src/sage/groups/perm_gps/permgroup_element.pyx
@@ -77,7 +77,7 @@ from sage.sets.finite_enumerated_set import FiniteEnumeratedSet
 import sage.structure.coerce as coerce
 from sage.structure.richcmp cimport richcmp_not_equal, rich_to_bool
 
-from sage.libs.gap.element cimport libGAP_Obj, libGAP_INT_INTOBJ, libGAP_ELM_LIST
+from sage.libs.gap.gap_includes cimport libGAP_Obj, libGAP_INT_INTOBJ, libGAP_ELM_LIST
 
 import operator
 

after this one gets

[sagelib-8.3.beta5] Error compiling Cython file:
[sagelib-8.3.beta5] ------------------------------------------------------------
[sagelib-8.3.beta5] ...
[sagelib-8.3.beta5]         for i from vn <= i < old.n:
[sagelib-8.3.beta5]             new.perm[i] = i
[sagelib-8.3.beta5]         return new
[sagelib-8.3.beta5] 
[sagelib-8.3.beta5]     cpdef _generate_new_GAP(old, self):
[sagelib-8.3.beta5]         cdef libGAP_Obj obj = self.value
[sagelib-8.3.beta5]                                  ^
[sagelib-8.3.beta5] ------------------------------------------------------------
[sagelib-8.3.beta5] 
[sagelib-8.3.beta5] sage/groups/perm_gps/permgroup_element.pyx:894:34: Cannot convert Python object to 'libGAP_Obj'
[sagelib-8.3.beta5] Traceback (most recent call last):
[sagelib-8.3.beta5]   File "/home/dima/sagetrac-mirror/local/lib/python2.7/site-packages/Cython/Build/Dependencies.py", line 1164, in cythonize_one_helper
[sagelib-8.3.beta5]     return cythonize_one(*m)
[sagelib-8.3.beta5]   File "/home/dima/sagetrac-mirror/local/lib/python2.7/site-packages/Cython/Build/Dependencies.py", line 1146, in cythonize_one
[sagelib-8.3.beta5]     raise CompileError(None, pyx_file)
[sagelib-8.3.beta5] CompileError: sage/groups/perm_gps/permgroup_element.pyx

which hopefully makes sense to you.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​49b68c4

here is the next import crash

[stumpc5]

okay, here is the next one that doesn't. It now compiles but crashes on startup with the following import error:

ImportError: /home/stumpc5/Programs/sage/local/lib/python2.7/site-packages/sage/groups/perm_gps/permgroup_element.so: undefined symbol: libGAP_ElmListFuncs

[dimpase]

You have this line:

cdef GapElement_List self = <GapElement_List> self_tmp

which does not make sense to me. Is it mean to be a C cast? Something else?

[stumpc5]

Yes, it is meant to be a C cast: self_tmp (sorry for the stupid names, I would make these nicer once it works) arrives as a python object because it is a cpdef, so to access its "C part", I needed to do this (actually without the <GapElement_List> cast, but only the cdef GapElement_List self = self_tmp iirc.

[dimpase]

Replying to stumpc5:

Yes, it is meant to be a C cast: self_tmp (sorry for the stupid names, I would make these nicer once it works) arrives as a python object because it is a cpdef,

I would be very surprised seeing a C cast magically converting a Python object into a GAP one.

so to access its "C part", I needed to do this (actually without the <GapElement_List> cast, but only the cdef GapElement_List self = self_tmp iirc.

[stumpc5]

I would be very surprised seeing a C cast magically converting a Python object into a GAP one.

Well, please see that the very same code in sage.libs.gap.element.GapElement_List._generate_perm works, and it doesn't without the cast cdef PermutationGroupElement tmp = <PermutationGroupElement> old. I thought this is because it arrives as a python object PermutationGroupElement (which doesn't have access to the cython data structures) and I cast it into a cython object.

[dimpase]

Does one really need to make it cpdef? Perhaps start with cdef, make a cpdef wrapper then?

[dimpase]

now I am confused. The current branch works for me.

[stumpc5]

now I am confused. The current branch works for me.

compiling, or also running sage? I got the import error as in comment:80 only when running sage.

Does one really need to make it cpdef? Perhaps start with cdef, make a cpdef wrapper then?

I could try this, but it seems so simple: I have a cpdef function on class A that takes an element of class B as second argument, and I want to move this to a cpdef function on class B that takes an element of class A as second argument...

[dimpase]

OK, more precisely, all is fine on OSX, but on Linux Sage crashes at startup, with the import error as you posted!

from the crash report:

...
/home/dima/sagetrac-mirror/local/lib/python2.7/site-packages/sage/groups/perm_gps/permgroup.py in <module>()
    131 #  Distributed under the terms of the GNU General Public License (GPL)
...
    145 from sage.interfaces.gap import gap, GapElement
--> 146 from sage.groups.perm_gps.permgroup_element import PermutationGroupElement, standardize_generator
        global sage.groups.perm_gps.permgroup_element = undefined
...

ImportError: /home/dima/sagetrac-mirror/local/lib/python2.7/site-packages/sage/groups/perm_gps/permgroup_element.so: undefined symbol: libGAP_ElmListFuncs

[dimpase]

And this is not a surprise:

$ ldd local/lib/python2.7/site-packages/sage/groups/perm_gps/permgroup_element.so
        linux-vdso.so.1 =>  (0x00007fffe53f0000)
        libpython2.7.so.1.0 => /home/dima/sagetrac-mirror/local/lib/libpython2.7.so.1.0 (0x00007c9a2c87c000)
        libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007c9a2c65f000)
        libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007c9a2c295000)
        libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007c9a2c091000)
        libutil.so.1 => /lib/x86_64-linux-gnu/libutil.so.1 (0x00007c9a2be8e000)
        libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007c9a2bb85000)
        /lib64/ld-linux-x86-64.so.2 (0x00007c9a2cedf000)

it is not linked with libgap, that's why.

[stumpc5]

I have no idea what that means! (Or rather how to fix this.)

[dimpase]

Compare this with the libgap.so extension:

$ ldd local/lib/python2.7/site-packages/sage/libs/gap/libgap.so 
        linux-vdso.so.1 =>  (0x00007ffe0d3fd000)
        libgap.so.4 => /home/dima/sagetrac-mirror/local/lib/libgap.so.4 (0x00007f4959fa5000)
        libpython2.7.so.1.0 => /home/dima/sagetrac-mirror/local/lib/libpython2.7.so.1.0 (0x00007f4959b85000)
        libpthread.so.0 => /lib/x86_64-linux-gnu/libpthread.so.0 (0x00007f4959968000)
        libc.so.6 => /lib/x86_64-linux-gnu/libc.so.6 (0x00007f495959e000)
        libm.so.6 => /lib/x86_64-linux-gnu/libm.so.6 (0x00007f4959295000)
        libdl.so.2 => /lib/x86_64-linux-gnu/libdl.so.2 (0x00007f4959091000)
        libutil.so.1 => /lib/x86_64-linux-gnu/libutil.so.1 (0x00007f4958e8e000)
        /lib64/ld-linux-x86-64.so.2 (0x00007f495ac26000)

you see there libgap.so.4, so it can find the symbols there...

Surely it's some declaration missing somewhere, and that's all. Or, perhaps, the corresponding GAP function should somehow be blessed in libgap extension.

[stumpc5]

Does this mean I should rather not move the method and keep the current version?

[dimpase]

No, this just means we have to ask people who know Cython better. It must be totally trivial to fix.

I have an idea how it is done via setup.py, but here it's done in some other way, and cannot see any trace anywhere of how and why libgap.so extension gets linked to the correct dynamic library.

[dimpase]

I've asked here: ​https://groups.google.com/d/msg/sage-devel/_d28OHHsnE4/WnYr1L6NCAAJ

[dimpase]

here is a fix that makes it work:

$ git diff
diff --git a/src/sage/groups/perm_gps/permgroup_element.pxd b/src/sage/groups/perm_gps/permgroup_element.pxd
index 08a928b..e261780 100644
--- a/src/sage/groups/perm_gps/permgroup_element.pxd
+++ b/src/sage/groups/perm_gps/permgroup_element.pxd
@@ -1,3 +1,4 @@
+# distutils: libraries = gap
 from sage.structure.element cimport MultiplicativeGroupElement, MonoidElement, Element
 from sage.structure.list_clone cimport ClonableIntArray

(at least with this I can start sage all right)

PS. Tests pass, too.

[jdemeyer]

Replying to dimpase:

here is a fix that makes it work:

$ git diff
diff --git a/src/sage/groups/perm_gps/permgroup_element.pxd b/src/sage/groups/perm_gps/permgroup_element.pxd
index 08a928b..e261780 100644
--- a/src/sage/groups/perm_gps/permgroup_element.pxd
+++ b/src/sage/groups/perm_gps/permgroup_element.pxd
@@ -1,3 +1,4 @@
+# distutils: libraries = gap
 from sage.structure.element cimport MultiplicativeGroupElement, MonoidElement, Element
 from sage.structure.list_clone cimport ClonableIntArray

That fix shouldn't be needed (I created #25608 for that).

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​6654f75

moved things around, ready for review

[stumpc5]

Okay, the method move is now finished, setting it back to review. Still open to see (maybe @tscrim?) if there is a "simple" way to further speed _generate_new and _generate_new_GAP in permgroup_element.

[jdemeyer]

Adding # distutils: gap to src/sage/groups/perm_gps/permgroup_element.pxd is not correct since permgroup_element.pxd has nothing to do with libGAP. Only the implementation requires libGAP, so the distutils declaration should be in the .pyx file.

But really, you should just use #25608 instead and not add # distutils: gap.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​70fbfff

Merge branch 'develop' into t/25477/enumerate_PermutationGroup

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​44b9ee0	Further clean up of libGAP declarations
​99fa0d0	Merge branch 't/25608/further_clean_up_of_libgap_declarations' into t/25477/enumerate_PermutationGroup
​f99aee1	removed gap import in permgroup_element

[jdemeyer]

1. You should use range() instead of this out-dated syntax:

   for i from 0 <= i < vn:
   

2. Are you certain that cdef int vn = lst.__len__() won't overflow? As far as I know, both Python and libGAP support lists with over 2³² elements. Better use Py_ssize_t, which is anyway the standard CPython type for lengths.

3. Why did you write lst.__len__() instead of len(lst)?

4. Why are the new methods cpdef instead of def? They are not called from Cython, so that's only a pessimization.

5. Why cdef PermutationGroupElement tmp = <PermutationGroupElement> self?

6. This is dangerous: <GapElement_List> lst_in since you don't know that lst_in is a GapElement_List. Either do an explicit isinstance() check or use <GapElement_List?>lst_in.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​c9b495e

fixed issues raised by jdemeyer

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1be187e

removed old artefact in code

[stumpc5]

Replying to jdemeyer:

1. You should use range() instead of this out-dated syntax:

done

2. Are you certain that cdef int vn = lst.__len__() won't overflow? As far as I know, libGAP supports lists with over 2^32^ elements. Better use Py_ssize_t, which is anyway the standard CPython type for lengths.

done

3. Why did you write lst.__len__() instead of len(lst)?

doesn't the second call the first? I thought that I'd rather use the first then.

4. Why are the new methods cpdef instead of def? They are not called from Cython, so that's only a pessimization.

my impression is that these are now the fastest way to create new permutations in sage, so I want to propagate their use, especially in critical cythonized places.

5. Why cdef PermutationGroupElement tmp = <PermutationGroupElement> self?

removed

6. This is dangerous: <GapElement_List> lst_in since you don't know that lst_in is a GapElement_List. Either do an explicit isinstance() check or use <GapElement_List?>lst_in.

I would like to do cpdef _generate_new_GAP(self, GapElement_List lst_in) but do not know how

---

New commits:

​c9b495e	fixed issues raised by jdemeyer
​1be187e	removed old artefact in code

[dimpase]

I don't understand why you're moving things back and forth - don't you want to provide an optimised replacement for GapElement_Permutation in sage.libs.gap.element ?

[stumpc5]

No, I want to obtain an optimized way to create a PermutationGroupElement from a GapElement_Permutation given that I have given a permutation group element of the same parent. This thus became a method for PermutationGroupElement with argument GapElement_List. In particular, I have this method now for python lists and for gap list.

As a byproduct, one may call this from the one of the parent if only the parent is given but not the element.

See that this later is indeed slower:

sage: P = PermutationGroup([(1,2),(1,2,3,4,5,6,7,8)])
sage: one = P.one()
sage: perm = libgap.eval('[]')
sage: %timeit one._generate_new_GAP(perm)
The slowest run took 138.45 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 174 ns per loop
sage: %timeit P.one()._generate_new_GAP(perm)
The slowest run took 77.00 times longer than the fastest. This could mean that an intermediate result is being cached.
1000000 loops, best of 3: 220 ns per loop

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​82244f3

Merge branch 'develop' into t/25477/enumerate_PermutationGroup

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1d1e70c

updated and extended doctests

[stumpc5]

Again ready for review!

6. This is dangerous: <GapElement_List> lst_in since you don't know that lst_in is a GapElement_List. Either do an explicit isinstance() check or use <GapElement_List?>lst_in.

I would like to do cpdef _generate_new_GAP(self, GapElement_List lst_in) but do not know how

this is still open.

[tscrim]

So I don't understand the comment about iterations. So I removed it. However, if SGS is not consistent across different computers, we have a problem with doctests. The RES is guaranteed to have a fixed iteration order, so I can only assume the SGS does not have such a guarantee. I changed the doctests to what I needed for them to pass, but I am guessing they will fail for you. If they pass for you and SGS is consistent, then the comment should have been removed anyways.

I made a few other tweaks.

---

New commits:

​83d485a

Reviewer changes.

[stumpc5]

Thanks Travis!

-            G = libgap.Group(self.gens())
+            G = libgap(self) #libgap.Group(self.gens())

I chose the first because it was slightly faster. Also I wanted to keep

-            # the following also works but is much slower:
-            #gap2sage = lambda elt: one._generate_new(libgap.ListPerm(elt).sage())
-            #gap2sage = lambda elt: libgap.ListPerm(elt)._generate_perm(one)

to see the other but slower ways to do this.

-            in the two algorithms, and is not even guaranteed to always
-            be the same for one of the two.

This now depends on gap behaviour, also I prefer not to guarantee that I will not change the output when changing the implementation. Anyways, if you prefer to remove it, that's fine too.

[tscrim]

Replying to stumpc5:

Thanks Travis!

-            G = libgap.Group(self.gens())
+            G = libgap(self) #libgap.Group(self.gens())

I chose the first because it was slightly faster.

I get the latter is faster:

sage: A = PermutationGroup([(1,2),(1,2,3,4,5,6,7,8,9)])
sage: %timeit libgap.Group(A.gens())
The slowest run took 13.32 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 53.5 µs per loop
sage: %timeit libgap(A)
The slowest run took 4.21 times longer than the fastest. This could mean that an intermediate result is being cached.
10000 loops, best of 3: 35.6 µs per loop

Also I wanted to keep

-            # the following also works but is much slower:
-            #gap2sage = lambda elt: one._generate_new(libgap.ListPerm(elt).sage())
-            #gap2sage = lambda elt: libgap.ListPerm(elt)._generate_perm(one)

to see the other but slower ways to do this.

Commented out code gets removed quite quickly by other people and looks bad to put slower code. The data is in the history if someone wants to look it up.

-            in the two algorithms, and is not even guaranteed to always
-            be the same for one of the two.

This now depends on gap behaviour, also I prefer not to guarantee that I will not change the output when changing the implementation. Anyways, if you prefer to remove it, that's fine too.

If a particular version of GAP is consistent, than no comment is necessary (a number of doctests typically change when GAP upgrades). We have never guaranteed an output order for group iteration, so I don't see any point in mentioning it.

[stumpc5]

Okay, let's see which doctests have to be undated (given the new ordering of the permgroup iteration).

[stumpc5]

I updated the permutation group iteration order in judson-abstract-algebra and cc'ed Rob Beezer her as requested there.

---

New commits:

​63d057d

fixed doctests in judson-abstract-algebra

[stumpc5]

Currently, the slow iteration using RecursivelyEnumeratedSet uses depth first search with generators being the group generators. Which of the following would you prefer:

1. Depth or breadth first search ?
2. Generators or generators plus their inverses?

I would think that bfs would be a more natural ordering, especially for small groups that you actually look at and see how the Cayley graph is propagated. Using also inverses is another question and I don't actually have a preference. Sometimes generating sets for Cayley graphs are expected to also contain inverses, sometimes not.

[dimpase]

while you are at it, could you fix the doc for PermutationGroupElement? It appears that what you get from sage: PermutationGroupElement? has not made it into the corresponding rst file, e.g. it's not possible to see an example with an explicit parent in the reference manual (in section ​http://doc.sagemath.org/html/en/reference/groups/sage/groups/perm_gps/permgroup_element.html)

Also, an example with an explicit parent out to be given in Examples, not in Tests only.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a2dd202	updated doctests for generic graph
​596dc64	fixed some more doctests

[stumpc5]

I fixed a few more doctests, not done yet. Let me mention that I find it questionable (and risky) that so many doctests fail due to this change of iteration...

[dimpase]

to me, this rather shows that make doctests depend upon an order which need to be preserved is not a good idea in the first place...

[stumpc5]

yes, and even more is true: some doctest pick out some unspecified element

         sage: SGA.group()
         Symmetric group of order 4! as a permutation group
         sage: SGA.an_element()
-        () + 2*(1,2) + 4*(1,2,3,4)
+        () + (1,2,3,4) + 2*(1,3)(2,4) + 3*(1,4)(2,3)

do stuff with it

 
             sage: G = SymmetricGroup(4).algebra(QQ)
             sage: G.retract_plain(G.an_element(), 3)
-            () + 2*(1,2)
+            ()

which might be boring for other elements (such as here), or even pick special elements with intended properties from it

     sage: D = DihedralGroup(5)
     sage: elements = D.list(); elements
-    [(), (1,5)(2,4), (1,2,3,4,5), (1,4)(2,3), (1,3,5,2,4), (2,5)(3,4),
-     (1,3)(4,5), (1,5,4,3,2), (1,4,2,5,3), (1,2)(3,5)]
-
-~~~~~~~~~~~~~~~~~~~~~~ ::
-
-    sage: rotate = elements[2]
-    sage: flip = elements[3]

[dimpase]

I would have rewritten the test

sage: SGA.an_element()
...

as

sage: SGA.an_element() in SGA
True

etc...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​a9fb42b

fixed some more doctests

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​83d485a	Reviewer changes.
​a6d7793	Merge branch 'u/tscrim/iteration_permutation_groups-25477' into t/25477/enumerate_PermutationGroup

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​55d5c64

implemented breadth and depth first search iteration

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​9504bc2

fixed more doctests

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​97ce629

more doctests

[tscrim]

Replying to dimpase:

I would have rewritten the test

sage: SGA.an_element()
...

as

sage: SGA.an_element() in SGA
True

etc...

Well, the generic element is usually generic enough to be useful. This test only works when testing _an_element_. It is used elsewhere to show how the element changes.

[tscrim]

Replying to stumpc5:

I fixed a few more doctests, not done yet. Let me mention that I find it questionable (and risky) that so many doctests fail due to this change of iteration...

I agree there is some risk, but a good review should be checking that it is just an order difference. So generically, I would say your concern is misplaced. However, I am very worried there is an order dependency on the machine, which means we have write more annoying doctests (i.e., calling sorted a lot).

[tscrim]

Replying to stumpc5:

Currently, the slow iteration using RecursivelyEnumeratedSet uses depth first search with generators being the group generators. Which of the following would you prefer:

1. Depth or breadth first search ?
2. Generators or generators plus their inverses?

I would think that bfs would be a more natural ordering

I would just let it run its default, but if I had to choose, I would say BFS as well.

[tscrim]

Also one little thing, but could you break the long output line in sage/categories/finite_dimensional_lie_algebras_with_basis.py (like how it was before).

[stumpc5]

Replying to tscrim:

I would just let it run its default, but if I had to choose, I would say BFS as well.

Well, we now have the three options SGS, BFS, DFS. I think that should be okay now...

[tscrim]

Replying to stumpc5:

Replying to tscrim:

I would just let it run its default, but if I had to choose, I would say BFS as well.

Well, we now have the three options SGS, BFS, DFS. I think that should be okay now...

A better way would be to pass options to RES.

[rbeezer]

Dear Christian,

Thanks very much for the cc on the changes to Judson's AATA. I'll be able to do a diff and even more quickly discern what has changed. I'm waiting on 8.3 to do a new update for the book, so this may miss that cut-off, but that should not be a crisis by any means. Thanks everyone for your work keeping this code working well.

Rob

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b586d0c

added examples of perm gp elts into documentation

[stumpc5]

Replying to dimpase:

while you are at it, could you fix the doc for PermutationGroupElement?

done.

[tscrim]

I know, I need to finish a review of this (including doing a rebase). This might also conflict with #26045. This is now been able to become high on my todo list.

[tscrim]

I have rebased it to 8.4.beta1 and over #26045. However, let's let this run on some patchbots to see if there is any machine dependency (in addition, if some of the other people here can do so). If there are no doctest failures, then we are safe to assume the SGS iterator is not machine dependent. With that, I would be okay setting a positive review.

---

New commits:

​19f82d1	various fixes for py3 in combinat and root_systems
​7456780	pyflakes cleanup in one file
​3c25ec3	fix sorting key for affine root systems
​4a2c244	Merge branch 'u/chapoton/26045' in 8.4.b1
​79c25c9	Merge branch 'u/stumpc5/iteration_permutation_groups-25477' of git://trac.sagemath.org/sage into public/group_theory/improve_iterator_perm_gps-25477
​ca2d1b3	Fixing doctests due to bad rebase and some doc tweaks.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1b00bf2

merged

[stumpc5]

@tscrim: what to do with these doctest failures?

[git]

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

​5af3439	Fix the erroneous fix of some corner cases
​138f85b	Turn some functions of polynomial.hilbert to methods of ETuple
​e900222	Hilbert series: Remove unnecessary bound checks, simplify code branches
​f0e5c85	Add documentation to new ETuple methods
​2d29b9b	Using slices, custom median, and formatting/P#P8 tweaks to Hilbert code.
​07ff2e1	Reverting use of get_exp for __getitem__.
​7b2ed50	Merge branch 'u/tscrim/hilbert_functions-26243' of git://trac.sagemath.org/sage into u/tscrim/hilbert_functions-26243
​5637e07	Fixing INPUT:: to INPUT: in hilbert.pyx.
​10f3921	Merge branch 'public/group_theory/improve_iterator_perm_gps-25477' of git://trac.sagemath.org/sage into public/group_theory/improve_iterator_perm_gps-25477
​f226bc3	Fixing trivial doctest failures.

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​dddbcf3

Fixing trivial doctest failures.

[tscrim]

Since #26225 is now in Sage, those doctest should essentially be sorted (or given by a fixed ordering). So this should fix those. Also you were missing a blank line in one of the book tests.

(Sorry about the noise, I was not based on develop. Fixed in the last forced push.)

[stumpc5]

Only

File "src/sage/tests/books/judson-abstract-algebra/permute-sage.py", line 268, in sage.tests.books.judson-abstract-algebra.permute-sage
Failed example:
    H.list()
Expected:
    [(), (1,3), (4,5), (1,3)(4,5)]
Got:
    [(), (4,5), (1,3), (1,3)(4,5)]

seems to be left. What kind of sorting do you want there?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1a3a69a

Fix the ordering of one more test.

[tscrim]

Replying to stumpc5:

What kind of sorting do you want there?

I suspect that is just a leftover from the test not being run. I doubt it is machine-dependent, so the fix I just pushed should take care of it.

[stumpc5]

okay, this seems to be ready to go!

[vbraun]

Can you have a look at the patchbot output: Python3 incompatible code inserted on 1 non-empty lines. Also, tests fail

sage -t --long src/sage/geometry/polyhedron/ppl_lattice_polytope.py  # 1 doctest failed

[tscrim]

Replying to vbraun:

Can you have a look at the patchbot output: Python3 incompatible code inserted on 1 non-empty lines.

That is libGAP code on that line (comes from the libgap.function_factory("""), so it is not Python3 incompatible.

Also, tests fail

sage -t --long src/sage/geometry/polyhedron/ppl_lattice_polytope.py  # 1 doctest failed

I think I forgot to run the --long when I checked this locally (and so I concluded that it was a spurious patchbot fault). However, I did confirm this locally. It seems like the group is too large to pass to libgap:

sage: libgap(aut)  # the group from the test in question
python2: libgap.c:184: libgap_get_input: Assertion `strlen(libGAP_stdin_buffer) < length' failed.

It has 1152 generators! However gap handles this fine.

So I just reverted the first change on comment:115, and that seems to fix the problem. I can see how that is more safe as I am surprised libgap is passing things by strings. Which might account for the fact that @stumpc5 and me got different timing results.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​1a85742

Using libgap.Group(self.gens()) to avoid input string limits via libgap(self).

[vbraun]

Its the Weyl group of F4...

libgap is passing by strings but there is an assumption that it fits into a certain buffer. This could be fixed...

[dimpase]

Why on Earth one generates all the 1152 elements of W(F_4)? This looks like something really weird in the implementation...

[tscrim]

Replying to vbraun:

Its the Weyl group of F4...

Right, but the computation uses all of the elements as generators, which is what makes it blowout the buffer.

libgap is passing by strings but there is an assumption that it fits into a certain buffer. This could be fixed...

I agree, but this is the less invasive (and easier) change. Ideally a PermutationGroup and its elements would directly construct the object in libgap using the C data structures.

[tscrim]

That change annoyingly causes the output order to change... So, should we fix the output order (again) or revert 1a85742 and do a more proper fix of the libgap interface?

[dimpase]

What exactly is happening in the latest commit? Do you mean to say that before it libgap(bla), for bla a permutation group, caused all the elements of bla to be listed?! If so, this was certainly a bug.

[tscrim]

No, I am saying that libgap(bla) and libgap.Group(bla.gens()) results in a different iteration order (at least in Sage) of the group elements.

[dimpase]

Replying to tscrim:

So I don't understand the comment about iterations. So I removed it. However, if SGS is not consistent across different computers, we have a problem with doctests. The RES is guaranteed to have a fixed iteration order, so I can only assume the SGS does not have such a guarantee. I changed the doctests to what I needed for them to pass, but I am guessing they will fail for you. If they pass for you and SGS is consistent, then the comment should have been removed anyways.

SGS is often computed with a randomisation involved, and so one cannot guarantee the same chain of subgroups across different runs on the same machine. If you fix the base to be used, you'd have less variance, but still there could be randomisation involved in choosing the ordering of the base elements.

We should not tweak algorithms in order to make testing easier, we should tweak tests, anyway.

What is RES?

[tscrim]

RES = RecursivelyEnumeratedSet.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​b52770e	Merge branch 'public/group_theory/improve_iterator_perm_gps-25477' of git://trac.sagemath.org/sage into public/group_theory/improve_iterator_perm_gps-25477
​916cea6	Fixing doctests due to change in output order.

[tscrim]

I've just updated the doctests to reflect the change in ordering (or remove the dependence on the ordering). So if the patchbots come back green, I think we can set this back to a positive review (that should be reasonable enough confirmation that the doctest outputs are not machine-dependent).

[tscrim]

One green patchbot and one almost green. I don't see how the error in computing the volume of a polytope could be caused by this ticket, so I suspect it is an unrelated failure.

[dimpase]

Let us get this in. The polyhedral computation stuff done, as in the doctest in question, over RDF, is dodgy anyway (more so as it is done on PPC64, of all platforms), and I don't see how it can have anything to do with this ticket.

[vbraun]

Merge conflict

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​ca255bb

Merge branch 'develop' into public/group_theory/improve_iterator_perm_gps-25477

[tscrim]

Trivial merge conflict.

[vbraun]

File "src/sage/rings/number_field/galois_group.py", line 137, in sage.rings.number_field.galois_group.GaloisGroup_v1.group
Failed example:
    list(P)                       # optional - database_gap
Expected:
    [(), (1,2), (1,2,3), (2,3), (1,3,2), (1,3)]
Got:
    [(), (1,2,3), (1,3,2), (2,3), (1,2), (1,3)]
**********************************************************************
1 item had failures:
   1 of   5 in sage.rings.number_field.galois_group.GaloisGroup_v1.group
    [142 tests, 1 failure, 3.05 s]
----------------------------------------------------------------------
sage -t --long src/sage/rings/number_field/galois_group.py  # 1 doctest failed

[dimpase]

Just sort these doctest results with something stable and platform-independent...

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​648e61c

Fixing database_gap doctest.

[tscrim]

It actually is not machine-dependent, it was just not tested by the patchbots because it is marked with an optional gap package to make sure there were no others (well, there was one also marked with magma, but I cannot test that one).

[dimpase]

I don't think this is a right fix. Any new version of GAP might mess this up... How about

sage: set(list(P))
{(), (2,3), (1,2), (1,2,3), (1,3,2), (1,3)}

or

sage: sorted(list(P))
[(), (2,3), (1,2), (1,2,3), (1,3,2), (1,3)]

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​6c63409

Using sorted on doctest.

[tscrim]

I'm not a fan of decorating doctests that don't really need it, but I've added it.

[vbraun]

On 32-bit:

**********************************************************************
File "src/sage/geometry/polyhedron/library.py", line 360, in sage.geometry.polyhedron.library.Polytopes.icosahedron
Failed example:
    ico.volume()
Expected:
    2.1816949907715726
Got:
    2.181694990771572
**********************************************************************
1 item had failures:
   1 of  11 in sage.geometry.polyhedron.library.Polytopes.icosahedron
    [211 tests, 1 failure, 118.70 s]
----------------------------------------------------------------------
sage -t --long src/sage/geometry/polyhedron/library.py  # 1 doctest failed

[dimpase]

oh, and that's why: the icosahedron generator contains

       verts = [p(v) for p in AlternatingGroup(3) for v in pts]
       return Polyhedron(vertices=verts, base_ring=base_ring, backend=backend)

how about we just add ... to that doctest, like 2.18169499... ?

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​693ada7

Adding rel tol to polyhedron/library.py test.

[tscrim]

I added a rel tol in case there ends up being other slight numerical differences in the future.

[dimpase]

Replying to tscrim:

I added a rel tol in case there ends up being other slight numerical differences in the future.

The digits after 2.18169499 are not right, for the exact evaluation gives 2.18169499062....

sage: (5/12*sqrt(5) + 5/4).n()
2.18169499062491

That's why I'd prefer truncation as I suggested.

[tscrim]

I think that is part of the point of the doctest, that it is different. If you insist, I will change it to truncation, but I also believe that suggests a much larger numerical instability than is actually present.

[dimpase]

IMHO it is not an instability, it is an artefact of naive computing with finite precision, i.e. accumulation of a rounding error.

[git]

Branch pushed to git repo; I updated commit sha1. New commits:

​3ec5057

Changing to truncated version.

[tscrim]

Okay, that is not quite the right terminology. It still suggests that there is significantly more variation in the doctest than what it really has. However, I changed it because I want this to get in and don't want it cut from the next stable release.

[dimpase]

Replying to tscrim:

Okay, that is not quite the right terminology. It still suggests that there is significantly more variation in the doctest than what it really has. However, I changed it because I want this to get in and don't want it cut from the next stable release.

It rather illustrates that RDF is not a field, and that it isn't even commutative... There could be more variation coming from a yet another CPU used. E.g. if you have a better FPU then you might get closer to the correct value, which would again break the doctest if it had too much numerical noise left.