## Gap 4 script based on paper GRS1997b: Recognition of actions on pairs
## II: with applications to groups acting on graphs.
##
## pairsII1.g: solution of unordered pairs problem using work on orbitals

OrbitalActionsUnorderedPairs := function(G, n)
# G acts on [1..n]

    local output,     # list containing output information, result
          trans,      # factored inverse transversal for Gw in G
          Gw,         # Stabilizer of 1 in G
          GwOrbs,     # orbits of Gw on [1..n]; non-singleton Gw-orbits
          Stabs,      # stabilizers Gw,w' for an elt w' of each elt of GwOrbs
          index2blks, # list of all Gw-blocks found whose stabilizer has
                      # index 2 in Gw
          Delta,      # elt of GwOrbs
          ww,         # favourite point in Delta
          Gww,        # Stabilizer of ww in Gw
          trans1,     # factored inverse transversal for Gww in Gw
          Hblk,       # generating set for block in (Gw, Delta) which we will
                      # factor out by
          Hblks,      # Translates of block generated by Hblk
          rep,        # maps Delta to [1..Length(Hblks)]
          elt,        # maps [1..Length(Hblks)] to elts of members of Hblk 
          GwImgGens,  # images of gens of Gw under map given by rep
          GwImg,      # isomorphic to group induced by Gw on Hblks
          H,          # stabilizer of Hblk, but as subgroup of GwImg
          found,      # boolean used to test whether H is normal in Gw
          ind2blks,   # used to construct the blocks of GwImg/H 
          basis,      # gradually extended to a basis of GwImg/H
          trans2,     # factored inverse transversal of 1 in GwImg/H
          Iblks,      # system of translates of Iblk in GwImg/H
          Iblk,       # block in GwImg/H generated by 1 and basis
          newIblk,    # used to construct new Iblk from old
          newind2blks, # used to construct new ind2blks from old
          D,          # used to construct action on cosets of a subgroup 
                      # of index 2 in Gw (the stab of an elt of index2blks)
          Gcos,       # action of G on these cosets
          out,        # result of calling an orbital actions routine
          transc,     # factored inverse transversal for Gcos_1 in Gcos
    # Miscellaneous localised variables
          orbit,pnt, gen, gpnt, img, g,i,h, z, img2, gz,b,pr, x, gx;

    # Find inverse factored transversal
    orbit  := [1];
    trans := [];
    trans[1] := ();
    for pnt in orbit do
      for gen in  GeneratorsOfGroup(G) do
        if not IsBound( trans[ pnt/gen] ) then
          Add( orbit, pnt/gen );
          trans[ pnt/gen ] := gen;
        fi;
      od;
    od;

    output := [];

    # Index 1 case
    out := OrbitalActionsNonSelfPaired(G, n);
    Info(InfoOperation, 1, "UOP2core: (index 1) ",Length(out),
              " actions found");

    for z in out do
      for x in z[2] do
        gx := [];
        img := x;
        while img <> 1 do
          Add(gx, trans[img]);
          img := img ^ trans[img];
        od;
        img := z[1];
        for g in Reversed(gx) do
          img := List(img, x -> x / g);
        od;
        img := Set(img);
        Add(output, Set([z[1], img]));
      od;
    od;



    # Index 2 case - find subgroups of index 2
    Gw := Stabilizer(G, 1);
    GwOrbs := Orbits(Gw, [1..n]);

    GwOrbs := Filtered(GwOrbs, x -> Length(x) > 1 );
    Stabs := [];

    index2blks := [];

    for Delta in GwOrbs do
      ww := Delta[1];
      Gww := Stabilizer(Gw, ww);

      # Find inverse factored transversal
      orbit  := [ww];
      trans1 := [];
      trans1[1] := ();
      for pnt in orbit do
	for gen in  GeneratorsOfGroup(Gw) do
	  if not IsBound( trans1[ pnt/gen] ) then
	    Add( orbit, pnt/gen );
	    trans1[ pnt/gen ] := gen;
	  fi;
	od;
      od;
      
      Hblk := [ww];
      for pnt in Delta do
	gpnt := [];
	img := pnt;
	while img <> ww do
	  Add(gpnt, trans1[img]);
	  img := img ^ trans1[img];
	od;
	img := pnt;
	for g in Reversed(gpnt) do
	  img := img / g;
	od;
	Add(Hblk, img);
      od;
      
      # We want the orbit of ww under the normal closure of the
      # stabilizer of the smallest Gw-block containing Hblk.
      repeat
      
	# Find action on translates of Hblk; find point stabilizer H, 
	# decide whether H is normal
	Hblks := Blocks(Gw, Delta, Set(Hblk));
	rep := [];
	elt := [];
	for b in Hblks do
	  Add(elt, b[1]);
	  for pnt in b do
	    rep[pnt] := Length(elt);
	  od;
	od;
	GwImgGens := List(GeneratorsOfGroup(Gw), gen ->
		  PermList(List(elt, x -> rep[x^gen])));
	GwImg := Group(GwImgGens, ());
      
	H := Stabilizer(GwImg, rep[ww]);
      
	found := false;
	i := 1;
	while not found and i <= Length(GeneratorsOfGroup(H)) do
	  h := GeneratorsOfGroup(H)[i];
	  z := First([1..Length(Hblks)], x -> x^h <> x);
	  if z <> fail then 
	    found := true;
	  else 
	    i := i+1; 
	  fi;
	od;
      
	if found then
	  # add to Hblk as H is not normal yet
	  # note that z^h <> z
	  # add to Hblk the elt ww^{ghg^{-1}} where z=ww^g
	  img := elt[z^h];
	  img2 := elt[z];
	  while img2 <> ww do
	    img := img ^ trans1[img2];
	    img2 := img2 ^ trans1[img2];
	  od;
	  Add(Hblk, img);
	fi;

      until not found;
  
      # We should now have H = Gww(Gw)^2, and GwImg is an el.ab. 2-gp.
  
      if Length(Hblks) = 1 then
        ind2blks := [];
      elif Length(Hblks) = 2 then
        ind2blks := [[1]];
      else
	orbit := [1];
	trans2 := [()];
	for pnt in orbit do
	  for gen in GeneratorsOfGroup(GwImg) do
	    if not IsBound( trans2[ pnt/gen] ) then
	      Add( orbit, pnt/gen );
	      trans2[ pnt/gen ] := gen;
	    fi;
	  od;
	od;

	basis := [2,3];

	Iblks := Blocks(GwImg, [1..Length(Hblks)], Concatenation([1],basis));

	# each block should have size 4
	Iblk := First(Iblks, x -> 1 in x);

	ind2blks := [[1,2],[1,3],[1,First(Iblk, x -> not x in [1,2,3])]];

	while Length(Iblks) > 1 do

	  # extend basis and Iblk (block generated by basis + 1)
	  z :=  First(Iblks, x -> not 1 in x)[1];
	  gz := [];
	  img := z;
	  while img <> 1 do
	    Add(gz, trans2[img]);
	    img := img ^ trans2[img];
	  od;
	  gz := Reversed(gz);

	  Add(basis, z);
	  Iblks := Blocks(GwImg, [1..Length(Hblks)], Concatenation([1],basis));
	  newIblk := First(Iblks, x -> 1 in x);
	  newind2blks := [Iblk];
	  for b in ind2blks do
	    # b has index 4 in newIblk, index 2 in Iblk. Need the translates
	    # of b that lie in newIblk
	    img := b;
	    for g in gz do
	      img := List(img, x -> x / g);
	    od;
	    img := Set(img);
	    Add(newind2blks, Union(b, img));
	    Add(newind2blks, Union(b, Difference(newIblk, Union(Iblk, img))));
	  od;
	  Iblk := newIblk;
	  ind2blks := newind2blks;
	od;
      fi;

      # Now project the blocks of index 2 back into Delta.
      ind2blks := List(ind2blks, b -> Set(Flat(Hblks{b})));

      # Add to main list, checking for duplicates by making sure each
      # one cannot have been found before.
      for b in ind2blks do
        if ForAll(Stabs, H ->
               ForAny(GeneratorsOfGroup(H), h -> not (OnSets(b, h) = b))) 
            then
          Add(index2blks, b);
        fi;
      od;

      Add(Stabs, Gww);      

    od;
    Info(InfoOperation, 1, "UOP2core: (index 2) ",Length(index2blks),
             " subgroups of index 2 found");

    # Now run an orbital test on each action found
    D := Cartesian ([1..n], [true, false]);

    for b in index2blks do

      Gcos := Operation(G, D, function( d, g )
                   local res, pnt, z, tr, gen;
                   res := [];
                   res[1] := d[1] ^ g;
                   pnt := d[1];
                   z := b[1];
                   tr := [];
                   while pnt <> 1 do
                     Add(tr, trans[pnt]);
                     pnt := pnt ^ trans[ pnt ];
                   od;
                   tr := Reversed(tr);
                   for gen in tr do
                     z := z / gen;
                   od;
                   z := z ^ g;
                   pnt := res[1];
                   while pnt <> 1 do
                     z := z ^ trans[pnt];
                     pnt := pnt ^ trans[pnt];
                   od;
                   if z in b then
                     res[2] := d[2];
                   else
                     res[2] := not d[2];
                   fi;
                   return res;
                   end                  );

      # NB x = 2 suffices (by construction) as the base point is taken as 1
      out := OrbitalActionsOneSelfPaired(Gcos, 2*n, 2);
      Info(InfoOperation, 1, "UOP2core: (index 2) ",Length(out),
              " actions found");
      # Find inverse factored transversal
      orbit  := [1];
      transc := [];
      transc[1] := ();
      for pnt in orbit do
	for gen in  GeneratorsOfGroup(Gcos) do
	  if not IsBound( transc[ pnt/gen] ) then
	    Add( orbit, pnt/gen );
	    transc[ pnt/gen ] := gen;
	  fi;
	od;
      od;

      for z in out do
        for x in z[2] do
          gx := [];
          img := x;
          while img <> 1 do
            Add(gx, transc[img]);
            img := img ^ transc[img];
          od;
          img := z[1];
          for g in Reversed(gx) do
            img := List(img, x -> x / g);
          od;
          img := Set(img);
          Add(output, Set(List([z[1], img], b -> Set(List(b, i ->
                          QuoInt(i+1,2))))));
        od;
      od;
    od;

    return output;

end;