Add graphs from group actions graph_from_group_action([::Type{T},] G, L, act, adj) where {T <: Union{Directed, Undirected}} and graph_from_group_action([::Type{T},] Omega::GSet, adj) where {T <: Union{Directed, Undirected}}

[long-zm123]

graph_from_group_action([::Type{T},] G, L, act, adj) where {T <: Union{Directed, Undirected}}
graph_from_group_action([::Type{T},] Omega::GSet, adj) where {T <: Union{Directed, Undirected}}

Return the graph on the vertices 1:length(L) with a directed edge (i, j)
if and only if adj(L[i], L[j]) returns true.

The group G must act on the list L via the action function act, i.e.
act(x, g) is the image of x in L under g in G, and L must be invariant
under this action. The elements of L must be pairwise distinct. The
adjacency predicate adj must be invariant under the action of G, i.e.
adj(x, y) == adj(act(x, g), act(y, g)) holds for all x, y in L and all
g in G.

If a G-set Omega is passed instead of G, L, act, then the group, the
vertex list and the action function are taken from Omega; the call is
equivalent to graph_from_group_action(T, acting_group(Omega), collect(Omega), action_function(Omega), adj). The vertices of the resulting graph are the
elements of collect(Omega) in their given order.

If the first argument is omitted, or is Directed, then the returned graph is
directed; if it is Undirected, then the returned graph is undirected and its
edges are obtained by symmetrizing the directed ones.

Examples

julia> G = symmetric_group(5);

julia> L = [[i, j] for i in 1:5 for j in i+1:5];

julia> g = graph_from_group_action(G, L, on_sets,
           (x, y) -> isempty(intersect(x, y)));

julia> n_vertices(g), n_edges(g)
(10, 30)

julia> gu = graph_from_group_action(Undirected, G, L, on_sets,
           (x, y) -> isempty(intersect(x, y)));

julia> n_vertices(gu), n_edges(gu)
(10, 15)

julia> Omega = gset(symmetric_group(5), on_sets, [[1, 2]]);

julia> g = graph_from_group_action(Omega, (x, y) -> isempty(intersect(x, y)));

julia> n_vertices(g), n_edges(g)
(10, 30)

[fieker]

@mjrodgers and @ThomasBreuer will have a (quick) look at this.

[mjrodgers]

Hi thanks for contributing this! I notice that this is sort of done without using a lot of the functionality we have already. For example, you seem to convert the group action to permutations on [1,2,...n] and then compute orbits on pairs by constructing the stabilizers. But we can get these orbits naturally using the existing GSet functionality like so:

function orbits_on_pairs(Omega::GSet)
  acthom = action_homomorphism(Omega)
  H = image(acthom)
  d = degree(H)
  OmegaPairs = gset(H, [[i,j] for i in 1:d, j in 1:d])
  return orbits(OmegaPairs)
end

this avoids any complications from exotic action functions. And then it is easy to decide which orbits get used as edges in the graph. Can you use this to make your code simpler? Is the performance still acceptable using this (and would it still do everything you want)?

(I guess you would want to eliminate the pairs [i,i] from the orbit calculation, but this is simple to do.)

[long-zm123]

Hi thanks for contributing this! I notice that this is sort of done without using a lot of the functionality we have already. For example, you seem to convert the group action to permutations on [1,2,...n] and then compute orbits on pairs by constructing the stabilizers. But we can get these orbits naturally using the existing GSet functionality like so:

function orbits_on_pairs(Omega::GSet)
  acthom = action_homomorphism(Omega)
  H = image(acthom)
  d = degree(H)
  OmegaPairs = gset(H, [[i,j] for i in 1:d, j in 1:d])
  return orbits(OmegaPairs)
end

this avoids any complications from exotic action functions. And then it is easy to decide which orbits get used as edges in the graph. Can you use this to make your code simpler? Is the performance still acceptable using this (and would it still do everything you want)?

(I guess you would want to eliminate the pairs [i,i] from the orbit calculation, but this is simple to do.)

Thanks for the review. We carefully evaluated the orbital (pair-orbit) approach and benchmarked it against the current implementation. The two are mathematically equivalent and produce identical edge sets, but we propose keeping the current implementation, because the design goal is to scale to larger graphs, which conflicts with the O(d²) bound of the orbital approach.

• Complexity: the orbital approach must materialize all d² ordered pairs and then compute their orbits — O(d²·m) time and O(d²) memory. The current implementation computes the same orbits via point-stabilizer orbits of representatives (orbit–stabilizer duality), without enumerating pairs: O(d·m + E) time and O(d + E) memory. Benchmark on PSL(2,41) with 574 vertices: current implementation 0.09 s / a few MiB; orbital approach ≈1.35 s / 268 MiB. Since the cost grows quadratically, the orbital approach is infeasible around d ≈ 10⁴, while the current one still handles large sparse graphs. This is also the reason GAP's Digraphs/GRAPE use representatives + stabilizer orbits + Schreier-vector translation; our implementation is the Julia port of that algorithm.
• adj calls are identical: both approaches call adj once per orbit (16 orbits in the benchmark); the orbital approach does not reduce the number of calls.
• Generality: the orbital approach relies on action_homomorphism, which requires the points and the action function to be GAP-convertible. The current implementation evaluates act directly in Julia; only the stabilizer step goes through GAP, and it is guarded (every generator must fix the representative, otherwise we fall back to a pure-Julia Schreier computation). In testing, GAP's Stabilizer with a Julia action function returned an incorrect subgroup in the left-coset case; our guard detects and recovers, while the orbital path has no such check.
• Details: image(acthom) returns the tuple (H, emb), so it should be image(acthom)[1]; and the diagonal pairs [i,i] should not be removed from the seeds, since that changes loop semantics — they should be left for adj to decide.

We therefore suggest keeping the current implementation. If pair-orbits are useful as a standalone feature, we can add orbits_on_pairs as a separate utility without changing the main path.

[HechtiDerLachs]

Briefly touched upon this in triage. @mjrodgers has some more comments which will probably follow, soon. @ThomasBreuer and him are still looking into this.

[ThomasBreuer]

@long-zm123 Concerning your comment

for exotic point types and Julia action functions the GAP stabilizer computation can be unreliable

Could you please show an example where this problem occurs? I have no idea how to construct such a situation.

[long-zm123]

Hi @ThomasBreuer, thanks for the question — here is a concrete, reproducible example.

Minimal example. Let G = S₄, H = Stab(1) ≅ S₃, and let G act on the left cosets of H by left multiplication, written as a Julia function:

using Oscar

G = symmetric_group(4)
H, _ = stabilizer(G, 1)                  # H = Stab(1), |H| = 6
c = left_coset(H, G([2, 1, 3, 4]))       # the left coset (1,2)·H
act = (x, g) -> g * x                    # left multiplication, a Julia function

st, _ = stabilizer(G, c, act)            # goes to GAP's Stabilizer
println("order of returned stabilizer: ", order(st))
for h in gens(st)
    println(h, "  fixes c? ", act(c, h) == c)
end

Output (Oscar 1.9.0-DEV, GAP 4.16):

order of returned stabilizer: 6
(2,4,3)  fixes c? false
(3,4)    fixes c? true

The true stabilizer of c = (1,2)·H is (1,2)·H·(1,2)⁻¹, which has order 6. GAP's Stabilizer also returns a subgroup of order 6, but one of its generators, (2,4,3), does not fix c at all — it sends c to the coset (1,2,4,3)·H. So the returned subgroup is not the stabilizer of c.

Why this happens. GAP's group-action machinery assumes a right action, μ(μ(x,g),h) = μ(x, g·h), and does not verify it — from the GAP Reference Manual (Chapter 41): "GAP does not test whether the acting function μ satisfies the conditions for a group operation but silently assumes that it does. (If it does not, results are unpredictable.)" Left multiplication act(x,g) = g·x is a left action (it satisfies μ(μ(x,g),h) = μ(x, h·g)), so GAP's orbit/stabilizer algorithms, which compose the action with the group product, silently compute the wrong object. As a control, right multiplication act(x,g) = x·g (a genuine right action) on the right cosets gives correct stabilizers in every case (0 bad generators).

This is not a rare corner case. For G = S₆, H = Stab({1,2}) (|H| = 48) acting by left multiplication on its 15 left cosets, GAP's Stabilizer returned generators that fail to fix the point for 12 of the 15 cosets (20 bad generators in total).

Why the guard in the PR is there. The implementation does not blindly trust the GAP result: for every returned generator hg it checks act(verts[r], hg) == verts[r], and if any check fails it discards the GAP generators and falls back to pure-Julia Schreier generators of the stabilizer (computed from the orbit BFS tree), which always fix the point. On the S₆ example above, the guarded implementation produces exactly the edge set of a brute-force reference (135 edges), so the final graph is correct even though the raw GAP stabilizer is wrong.

(For completeness: the alternative action_homomorphism-based pair-orbit approach suggested earlier is affected in the same way, since it routes the Julia action function through the same GAP machinery — on the left-coset G-set of the minimal example, action_homomorphism disagrees with act on 6 generator evaluations, while the right-action control agrees on all of them.)

[ThomasBreuer]

@long-zm123 Left multiplication (x, g) -> g * x does in general not define a group action.
For a group action f, the condition f(x, g1 * g2) = f(f(x, g1), g2) must be satisfied.