Permutation group: faster use of stabilizer chain from gap

[Christian Stump]

In https://trac.sagemath.org/ticket/25477, we (this is mainly Moritz Firsching) implement a new iterator for permutation groups using stabilizer chains. Indeed, the multiplication of permutations seems to be faster than in gap so we hope this method to be really fast in the end...

The bottleneck now is that we call

    self.strong_generating_system(base)

for a permutation group self and this is done doing tons of back and forth translations between sage and gap using the pexpect interface (if I understand it correctly).

I now have reimplemented the computation of the strong generating system directly in gap, and now the remaining bottleneck is to turn the gap permutations back to sage permutations. So my question is

    What is the currently fastest way to turn gap permutations into sage permutations?

Here is the code:

#############################

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;
"""
#############################

and here is the test case to look at:

#############################

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

#############################

so the actual strong generating system is computed very fast in gap in SGS_new(A, True) while casting the output into sage in SGS_new(A, False) then takes forever.

Thanks for your help -- I would expect this fast iteration will help not help my code to be faster!

    Christian

[Simon Brandhorst]

Use libgap:

sage: def SGS_libgap(G, gap_output=True):
....:     S = libgap(G).StabChain()
....:     libgap.eval(gap_cosets)
....:     cosets = S.CosetsStabChain()
....:     if gap_output:
....:         return cosets
....:    

sage: %time XX = SGS_libgap(A,True)
CPU times: user 31.7 ms, sys: 59 µs, total: 31.8 ms
Wall time: 31.9 ms
sage: %time XX = SGS_libgap(A,False)
CPU times: user 31.7 ms, sys: 59 µs, total: 31.8 ms
Wall time: 31.8 ms


sage: %time XX = SGS_old(A)
CPU times: user 3.43 s, sys: 535 ms, total: 3.96 s
Wall time: 4.03 s
sage: %time XX = SGS_new(A,True)
CPU times: user 3.39 ms, sys: 2.03 ms, total: 5.41 ms
Wall time: 7.98 ms
sage: %time XX = SGS_new(A,False)
CPU times: user 2.44 s, sys: 101 ms, total: 2.54 s
Wall time: 2.59 s

[Christian Stump]

Hi,

Thanks for the reply! The ticket has evolved drastically since my post here (and indeed uses libgap now), the bottleneck now is that I have a libgap list of libgap integers (which are indeed C integers) and need to use these in cython, which I currently fail to do.

Anyway, I do not see how your answer would improve the situation even in the beginning because I see that the previous

SGS_new(A,True)

is faster than

SGS_libgap(A,True)


and you do not provide a way to then turn the libgap lists into sage lists, and this takes all the time in the computation. Or am I overseeing something and I can get this output directly into sage?

Here is the new code

def SGS_new(self, gap_output=True):
    gap_cosets = libgap.function_factory("""function (G)
        local   S0,
                CosetsStabChain;

        S0 := StabChain(G);
        CosetsStabChain := function(S)  # for the recursive call

        local   cosets,             # element list, result

                new_coset,          # new coset computed along the way

                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;""")
    G = libgap.Group(self.gens())
    cosets = gap_cosets(G)
    # the following case from the gap permutation elt to a
    # sage permutation is much too slow -- stumpc5, 2018-06-05
    one = self.one()
    gap2sage = lambda elt: one._generate_new(libgap.ListPerm(elt).sage())
    if gap_output == 1:
        return cosets
    elif gap_output == 2:
        return [[elt for elt in coset] for coset in cosets]
    elif gap_output == 3:
        return [[[ i.sage() for i in libgap.ListPerm(elt)] for elt in coset] for coset in cosets]
    elif gap_output == 0:
        gap2sage = lambda elt: libgap.ListPerm(elt).sage()
        return [ [ gap2sage(elt) for elt in coset] for coset in cosets ]

As you can see

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()

sage: %time XX =  SGS_new(A,1)
CPU times: user 58.4 ms, sys: 60 µs, total: 58.4 ms
Wall time: 58.5 ms
sage: %time XX =  SGS_new(A,2)
CPU times: user 48.1 ms, sys: 75 µs, total: 48.2 ms
Wall time: 48.1 ms
sage: %time XX =  SGS_new(A,3)
CPU times: user 296 ms, sys: 3.95 ms, total: 300 ms
Wall time: 299 ms
sage: %time XX =  SGS_new(A,0)
CPU times: user 264 ms, sys: 79 µs, total: 264 ms
Wall time: 263 ms

As you can see, all the time is spent in translating thousands of gap ints into sage ints (and lists of these, but that seems fast). Help is still appreciated!

[Christian Stump]

Okay, I think I got all the bottlenecks resolved now -- I'd appreciate further comments and reviews on the ticket (https://trac.sagemath.org/ticket/25477), especially from someone familiar with the libgap interface.