Find group isomorphisms in Mathematica

findGroupIsomorphisms[
  group1_,
  group2_,
  max_: \[Infinity]
  ] :=
 Module[{ng1, ng2, reverse, g1, g2, order, gs1, go1, gm1, gm2,
   os1, os2, inng1, abelianQ, targets, isomorphisms, w1,
   mt1, \[CurlyPhi]m1, \[CurlyPhi], homomorphismQ, t, ts},
  (*
   * Choose the domain group with fewest generators
   *)
  If[GroupOrder[group1] != GroupOrder[group2],
   Return[{}];
   ];
  {ng1, ng2} = Length[GroupGenerators[#]] & /@ {group1, group2};
  reverse = ng2 < ng1;
  If[reverse,
   {g1, g2} = {group2, group1};
   {ng1, ng2} = {ng2, ng1},
   {g1, g2} = {group1, group2}
   ];
  (*
   * Do some quick checks for isomorphism
   *)
  order = GroupOrder[g1];
  {gm1, gm2} = GroupElements /@ {g1, g2};
  {os1, os2} = Map[PermutationOrder, {gm1, gm2}, {2}];
  If[Sort[Tally[os1]] != Sort[Tally[os2]],
   Return[{}]
   ];
  (*
   * Pick possible targets in g2
   *)
  gs1 = GroupGenerators[g1];
  go1 = PermutationOrder /@ gs1;
  targets = Table[
    Pick[gm2, Thread[os2 == n]],
    {n, go1}
    ];
  targets = Tuples[targets];
  (*
   * List the inner automorphisms
   *)
  inng1 = Outer[
    GroupElementPosition[
      g1, #1\[PermutationProduct]#2\[PermutationProduct]\
InversePermutation[#1]] &,
    gm1, gs1
    ];
  inng1 = Union[inng1];
  abelianQ = Length[inng1] == 1;
  (*
   *
  Here I should do something much more efficient for abelian groups...
\
   *)
  (* Stub *)
  (*
   *
  Test every possible mapping of the generators
   *)
  w1 = GroupElementToWord[g1, #] & /@ gm1;
  mt1 = GroupMultiplicationTable[g1];
  isomorphisms = {};
  While[targets =!= {},
   t = targets[[-1]]; targets = Drop[targets, -1];
   \[CurlyPhi] = With[{t = t},
     Function[i, GroupElementFromWord[PermutationGroup[t], w1[[i]]]]
     ];
   ts = Map[\[CurlyPhi], inng1, {2}];
   targets = Complement[targets, ts];
   If[GroupOrder[PermutationGroup[t]] != order,
    Continue[]
    ];
   \[CurlyPhi]m1 =
    GroupElementFromWord[PermutationGroup[t], #] & /@ w1;
   homomorphismQ =
    Map[\[CurlyPhi], mt1, {2}] ==
     Outer[PermutationProduct, \[CurlyPhi]m1, \[CurlyPhi]m1];
   If[homomorphismQ,
    (*
     * We found some isomorphisms!
     *)
    isomorphisms = Join[isomorphisms, ts]
    ];
   If[Length[isomorphisms] >= max,
    Break[]
    ]
   ];
  (*
   * Return the isomorphisms as rule lists
   *)
  isomorphisms =
   Take[isomorphisms, Min[max, Length[isomorphisms]]];
  isomorphisms = {gs1, #} & /@ isomorphisms;
  If[reverse, isomorphisms = Reverse /@ isomorphisms];
  Apply[Rule, Transpose /@ isomorphisms, {2}]
  ]


isomorphicGroupsQ[group1_, group2_] :=
 MatchQ[findGroupIsomorphisms[group1, group2, 1], {_}]