Speed up orbits of permutation groups on integers

[mjrodgers]

This aims to address issue #3287 by letting GAP know when we are in certain special cases which have optimized GAP methods for computing orbits.

Computing the orbit of a single point seems to still be slower than calling GAP directly by about a factor of 2, but this is a big improvement to the previous case (where it was about a 6x slowdown).

[mjrodgers]

@ThomasBreuer If you can take a look. At the moment I am only mapping ^ to OnPoints, I have some of the other conversions commented out since I don't think they are things that are supported yet.

[fingolfin]

What about other places passing an action function to GAP?

[mjrodgers]

What about other places passing an action function to GAP?

I don't think there are any, when I checked. We don't currently experience any slowdown like this when calling stabilizer, for example. (I did a project-wide search for GapObj(action_function and this was the only place it was used.)

[ThomasBreuer]

I think we have to distinguish two situations:

• If the elements of Omega can be translated to pure GAP objects and if we have a Julia action function such that the GapObjs of the group generators act on the GAP objects in question via this function then we have a chance to hit an efficient GAP method for the computations, and we can first translate the input to GAP and in the end carry back the result to Julia.
  This is the case for permutation groups acting on positive integers or vectors or sets of them.
  (For the action on sets and vectors of integers, we have to convert also omega to a pure GAP object.)
• In general we have some group acting on some objects via some function, for example a permutation group of degree n acting on polynomials in n variables by permuting the variables.
  (For examples, see test/Groups/gsets.jl).
  Here the old implementation uses that GAP allows us to act with Oscar group elements on Oscar objects, via the Julia function wrapped into a GAP function.

[ThomasBreuer]

Do we have a better idea than keeping the generic orbit method and adding methods for those special cases that "can be translated to GAP and back", similar to the approach that was taken for stabilizer in #4281?

[fingolfin]

Tests fail.

Also what this really is missing is some systematic benchmark tests, so that we can verify each change for how it affects performance.

Also, could you please share your example that where this PR makes it go from a 6x to a 2x slow down?

[mjrodgers]

I've added a new test to specifically catch the error (which was previously thrown outside of a @test statement). This case seems to throw an error when we use acts = GapObj(gens(G); recursive=true), but not using acts = GapObj(gens(G)); on the other hand, we get an error computing the orbits of symmetric_group(5) unless we use recursive=true... @ThomasBreuer any easy fix for this?

[mjrodgers]

Working to put together some systematic benchmarks, but for this specific case we have

julia> G = symmetric_group(5)
Sym(5)

julia> @benchmark Vector{Int}(GAP.Globals.Orbit(G.X,1))
BenchmarkTools.Trial: 10000 samples with 10 evaluations.
 Range (min … max):  1.700 μs …  4.329 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     1.829 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   1.844 μs ± 86.960 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

            ▁▄▆█▇▅▆
  ▂▂▂▂▂▃▃▄▇▇███████▇▆▅▅▅▄▄▅▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▂▂▂▂▂▂▂ ▃
  1.7 μs         Histogram: frequency by time        2.21 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 28.

With this pull request:

julia> @benchmark collect(orbit(G,1))
BenchmarkTools.Trial: 10000 samples with 8 evaluations.
 Range (min … max):  3.688 μs …   8.714 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     4.146 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   4.167 μs ± 195.265 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

                   ▄▇▄█▇▆▆▆▇▃▂
  ▁▁▁▂▂▂▃▃▂▂▂▂▂▂▃▅█████████████▅▇▅▄▃▃▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  3.69 μs         Histogram: frequency by time        4.89 μs <

 Memory estimate: 3.92 KiB, allocs estimate: 91.

This 2x slowdown is consistent when using larger symmetric_groups as well.

[fingolfin]

Suggested change

	acts = GapObj(gens(G); recursive=true)
	acts = [GapObj(g) for in gens(G)]

but even this will fail if g cannot be converted to a GAP object...

[ThomasBreuer]

Sorry if my comment above ("I think we have to distinguish two situations.") was not clear enough.

What I wanted to say is that in the current setup, we need first of all one generic orbit function. Up to now, we have one that delegates to GAP but keeps the Julia objects, that is, there is a GAP list of Julia group elements acting with a Julia function on Julia objects.
This method can deal with all cases but is expected to be slow.
(For this generic situation, an orbit method written in pure Julia is expected to be faster.)

In addition to that method, we can install more orbit methods. They can delegate to GAP or be written in Julia.
In the former case, they can decide to translate their arguments entirely to GAP, and then to translate the result back to Julia; this applies in particular to the situations for which GAP has special methods, such as the action of permutations on integers or (nested) lists or sets of them. (In pull request #4281, this has been done for stabilizer, it would look similar for orbit.) It is still not clear whether also in these situations, an orbit written in Julia will be faster than delegating to GAP.

[fingolfin]

This still fails all its tests. @mjrodgers what's the status there, are you working on it?

Regarding the benchmarks: the first example actually suffers from inefficient code for converting a GAP list of integers into a Vector{Int}, so the discrepancy is actually even larger if one considers that this conversion currently is very suboptimal (see oscar-system/GAP.jl#777 for stuck WIP on improving this -- perhaps @ThomasBreuer will adopt this at some point? he already took care of the other conversion direction)

Also type instability hurts this micro benchmark. To wit:

julia> f(obj::GapObj) = Int[obj[i] for i in 1:length(obj)]  # better conversion code
f (generic function with 1 method)

julia> G = symmetric_group(5)
Sym(5)

julia> @benchmark Vector{Int}(GAP.Globals.Orbit(G.X,1))
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.361 μs …  15.579 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.440 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.463 μs ± 178.038 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

     ▂▆▇▇█▆▅▄▃▁
  ▁▂▆███████████▇▅▅▄▃▃▂▂▂▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  2.36 μs         Histogram: frequency by time         2.9 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 28.

julia> @benchmark Vector{Int}(GAP.Globals.Orbit(G.X,1)::GapObj)
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.199 μs …  12.412 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.296 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.315 μs ± 164.926 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▃▆█▇▅▂
  ▁▂▄▇▇▇███████▇▅▅▄▃▃▃▃▂▂▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  2.2 μs          Histogram: frequency by time        2.75 μs <

 Memory estimate: 1.09 KiB, allocs estimate: 28.

julia> @benchmark f(GAP.Globals.Orbit(G.X,1))
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.130 μs …  12.282 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.204 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.220 μs ± 149.329 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

         ▂▆▁█ ▄ ▁
  ▂▃▃▇▅▇▅██████▆█▄▅▄▃▃▂▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂▂▁▂▂▂▂▂▂ ▃
  2.13 μs         Histogram: frequency by time        2.56 μs <

 Memory estimate: 720 bytes, allocs estimate: 21.

julia> @benchmark f(GAP.Globals.Orbit(G.X,1)::GapObj)
BenchmarkTools.Trial: 10000 samples with 9 evaluations.
 Range (min … max):  2.116 μs …  14.546 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     2.199 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   2.215 μs ± 181.758 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

              █▁▁▁▆
  ▁▁▂▂▅▃▄▄█▇▇██████▆▅▅▄▆▂▂▂▃▂▂▂▃▂▂▂▂▂▁▁▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂
  2.12 μs         Histogram: frequency by time        2.45 μs <

 Memory estimate: 720 bytes, allocs estimate: 21.

Of course this is minor and it equally affects the Oscar wrapper for GAP's Orbit.

For that, as @ThomasBreuer already pointed out, and as we discussed a few days ago, you'll probably need a different strategy:

E.g. it only makes sense to use the "the underlying GAP objects of the group generators" if the given group actually is a GAP group (which the method signature already ensures), and the objects being acted on are "native" GAP objects -- which right now probably means only for Int.

To unwrap the generators it is probably more efficient to use GAPWrap.GeneratorsOfGroup(GapObj(G)). Indeed (note the use of $G instead of G to get accurate results):

julia> G = symmetric_group(5);

julia> @benchmark GapObj(gens($G); recursive=true)
BenchmarkTools.Trial: 10000 samples with 209 evaluations.
 Range (min … max):  363.239 ns …  27.510 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     400.718 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   409.789 ns ± 284.250 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

            ▄▅█▆▅▅▂▂▂
  ▂▃▅▅▆▄▄▄▅███████████▇█▇▇▇▆▆▅▄▃▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▃
  363 ns           Histogram: frequency by time          513 ns <

 Memory estimate: 192 bytes, allocs estimate: 6.

julia> @benchmark Oscar.GAPWrap.GeneratorsOfGroup(GapObj($G))
BenchmarkTools.Trial: 10000 samples with 998 evaluations.
 Range (min … max):  14.946 ns … 35.363 ns  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     15.030 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   15.167 ns ±  0.742 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▃▇█▅▅ ▄▂ ▂ ▃ ▁▁▂▁                                           ▂
  █████▁████▁█████████▇▁▇███▁▆▇▆▆▆▁▇▆▅▆▁▆▅▅▆▁▅▄▄▅▄▁▁▃▅▅▁▄▄▅▅▄ █
  14.9 ns      Histogram: log(frequency) by time        17 ns <

 Memory estimate: 0 bytes, allocs estimate: 0.

[fingolfin]

Next I wonder how expensive storing the orbit as an attribute is -- this potentially allocates a fresh Dict which can be costly. You could find out by temporarily disabling the set_attribute! and measuring again

[mjrodgers]

This still fails all its tests. @mjrodgers what's the status there, are you working on it?

Working with Thomas on it, this pull request might end up being superseded by the solution that Thomas describes in his previous comment; but I wanted to at least include this additional test here, since it can help track an extra "weird" case that might need some extra care.

[ThomasBreuer]

Shall I make a new pull request that adds some special orbit methods, as I had sketched above, or shall I commit to the current pull request?

[mjrodgers]

I think it is good to add to the current pull request, if that is easy.

[ThomasBreuer]

I just tried to push to this pull request, but somehow this was not successful; I got the message that I can open a new pull request at https://github.com/oscar-system/Oscar.jl/pull/new/GSet_speedup.
(Last time when I tried to push to someone else's pull request, this worked.)

[joschmitt]

I just tried to push to this pull request, but somehow this was not successful; I got the message that I can open a new pull request at https://github.com/oscar-system/Oscar.jl/pull/new/GSet_speedup. (Last time when I tried to push to someone else's pull request, this worked.)

The branch is in @mjrodgers' fork and not in oscar-system/Oscar.jl. Maybe something got confused there?

[mjrodgers]

Strange, I do have it marked so that "Maintainers are allowed to edit this pull request", which says that "users with write access to oscar-system/Oscar.jl can add new commits to your GSet_speedup branch"... Let me check if there is something else I need to do.

…

On 27. Nov 2024, at 2:00 PM, Thomas Breuer ***@***.***> wrote: I just tried to push to this pull request, but somehow this was not successful; I got the message that I can open a new pull request at https://github.com/oscar-system/Oscar.jl/pull/new/GSet_speedup. (Last time when I tried to push to someone else's pull request, this worked.) — Reply to this email directly, view it on GitHub <#4307 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AANWIWN6FYGUE2G3FRR6KET2CW67DAVCNFSM6AAAAABRWQDT5GVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKMBTHAYTSNRSGQ>. You are receiving this because you were mentioned.

[mjrodgers]

I just tried to push to this pull request, but somehow this was not successful; I got the message that I can open a new pull request at https://github.com/oscar-system/Oscar.jl/pull/new/GSet_speedup. (Last time when I tried to push to someone else's pull request, this worked.)

I think maybe you need to push directly to my branch instead of to the pull request, Lars Kastner verified that he can push there.

[ThomasBreuer]

Now it worked. Thanks for the hints.

[codecov]

Codecov Report

Attention: Patch coverage is 96.96970% with 1 line in your changes missing coverage. Please review.

Project coverage is 84.40%. Comparing base (651e595) to head (70189c0).
Report is 20 commits behind head on master.

Files with missing lines	Patch %	Lines
src/Groups/gsets.jl	95.83%	1 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##           master    #4307      +/-   ##
==========================================
+ Coverage   84.36%   84.40%   +0.04%     
==========================================
  Files         652      656       +4     
  Lines       86955    87246     +291     
==========================================
+ Hits        73360    73641     +281     
- Misses      13595    13605      +10     

Files with missing lines	Coverage Δ
src/Groups/action.jl	97.45% <100.00%> (+0.16%)	⬆️
src/Groups/gsets.jl	92.59% <95.83%> (+0.23%)	⬆️

... and 44 files with indirect coverage changes

[fingolfin]

@mjrodgers ping didn't you want to push some updates yesterday?

[mjrodgers]

@fingolfin I have some general timing tests and comparisons between our different orbit methods that I am working on, but nothing to push at the moment (and I don't know that it needs to be tied to this specific pull request),

[lgoettgens]

This last commit somehow broke history. Iirc you didn't touch 133 files in this PR

[mjrodgers]

@lgoettgens No, I pulled the changes Thomas made to my local machine, and then everything was strange. I tried to revert my commit, but still everything is not working. I'm trying to figure out what went wrong (actually I'm okay if I don't figure out what went wrong as long as I can fix it).

That was very strange, when I pulled those it basically asked me if I wanted to merge my own branch into itself, and then all those files were changed (no idea where those changes came from). Then when I tried to revert it, they were still all changed. I think it's sorted now.

[mjrodgers]

@ThomasBreuer I am noticing some inconsistent behavior in my timing tests, for example if I run

julia> G = symmetric_group(1000)
Sym(1000)

julia> Omega = gset(G)
G-set of
  Sym(1000)
  with seeds 1:1000

julia> @benchmark order(stabilizer(Omega, Set([1, 2]))[1])
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):   77.958 μs … 270.000 ms  ┊ GC (min … max):  0.00% … 94.68%
 Time  (median):      89.000 μs               ┊ GC (median):     0.00%
 Time  (mean ± σ):   122.405 μs ±   2.738 ms  ┊ GC (mean ± σ):  25.49% ±  1.60%

                        ▄▆█▆▆▃
  ▂▁▂▂▁▂▂▂▃▃▄▄▄▃▃▂▂▂▃▄▆████████▆▅▄▃▃▃▃▃▃▃▃▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▁▂ ▃
  78 μs            Histogram: frequency by time          105 μs <

 Memory estimate: 186.04 KiB, allocs estimate: 770.

but if I try

julia> G = symmetric_group(1000)
Sym(1000)

julia> Omega = gset(G, [[1, 2]])  # action on ordered pairs
G-set of
  Sym(1000)
  with seeds [[1, 2]]

julia> @benchmark order(stabilizer(Omega, [1, 2])[1])

then it takes quite a long time (I'm not sure it ever actually finishes). Maybe we are somehow not triggering the specialized code in this case? (We can talk about this tomorrow in person, but I figured I would put it here in the pull request also.)

[ThomasBreuer]

Yes, the special stabilizer methods for a group are not called by the current delegation from G-sets to groups.
I think the right solution is a generic delegation through which special stabilizer methods for groups get called.
I will make a proposal for that.

[mjrodgers]

I think this is ready to merge, we will add more timing tests in a further pull request but we have some basics, we have greatly sped up both orbit and stabilizer computations in most of the basic special cases,