Counting
Yaroslav Bulatov
1, 1, 0, 0, 0, 0, 0, 0, 0} that are not equivalent under the symmetry
of DihedralGroup[16]. In other words, count the ways of assigning
those integers to vertices of a 4 dimensional cube.
This takes about a minute in another systme using "OrbitsDomain" command. My
Mathematica approach is below, however it doesn't finish within 10
minutes, any advice how to make it tractable?
nonequivalentPermutations[lst_, group_] := (
removeEquivalent[{}] := {};
removeEquivalent[list_] := (
Sow[First[list]];
equivalents = Permute[First[list], #] & /@ GroupElements[group];
DeleteCases[list, Alternatives @@ equivalents]
);
reaped = Reap@FixedPoint[removeEquivalent, Permutations@lst];
reaped[[2, 1]] // Length
);
nonequivalentPermutations[{2, 2, 2, 2, 2, 2, 2, 1, 1, 0, 0, 0, 0, 0,
0, 0}, DihedralGroup[16]]
Leonid Shifrin
this is a follow - up to my previous post. I think that the code I just
posted grossly overestimates
the number of non-equivalent permutations, because I forgot about the
non-abelian nature of
permutations and my method explores only a part of the group space. I did
not figure out the
rigorous way to solve this problem (perhaps one should look more carefully
at the group structure),
but we can at least modify the code to give an upper bound on the number of
non-equivalent
permutations:
1. Modify the code for getUniquePermutatations, to apply the longest cycles
first
Clear[getUniquePermutatationsAlt];
getUniquePermutatationsAlt[permutations_, actions : {__Cycles}] /;
Length[permutations] > 0 :=
With[{rlen = Range[Length[permutations[[1]]]]},
Fold[reducePermutationsC, permutations,
Map[Permute[rlen, #] &, SortBy[actions, Length@#[[1]] &]]]]
2. Use this function inside FixedPoint:
In[80]:= FixedPoint[ getUniquePermutatationsAlt[#,actions]&,perms]//Timin=
g
Out[80]=
{5.016,{{0,0,0,0,0,0,0,2,1,2,2,2,2,1,2,2},{0,0,0,0,0,0,0,2,1,2,2,2,1,2,2,2}=
,{0,0,0,0,0,0,0,2,1,2,2,1,2,2,2,2},{0,0,0,0,0,0,0,2,1,2,1,2,2,2,2,2},{0,0,0=
,0,0,0,0,2,1,1,2,2,2,2,2,2},{0,0,0,0,0,0,0,1,2,2,2,2,2,2,2,1},{0,0,0,0,0,0,=
0,1,2,2,2,2,2,2,1,2},{0,0,0,0,0,0,0,1,2,2,2,2,2,1,2,2},{0,0,0,0,0,0,0,1,2,2=
,2,2,1,2,2,2},{0,0,0,0,0,0,0,1,2,2,2,1,2,2,2,2},{0,0,0,0,0,0,0,1,2,2,1,2,2,=
2,2,2},{0,0,0,0,0,0,0,1,2,1,2,2,2,2,2,2},{0,0,0,0,0,0,0,1,1,2,2,2,2,2,2,2}}=
}
I don't have a proof with this method that none of the above is really
equivalent to others, alas. I have
some (empirical) evidence that this may be the case, but obviously better
arguments are needed here.
That is, if I understood the problem at all correctly, which is something =
I
am also starting to doubt.
Regards,
Leonid
On Mon, Dec 20, 2010 at 8:40 AM, Yaroslav Bulatov <yaros...@gmail.com>wro=
te:
Leonid Shifrin
This wasn't easy, but this was fun. You better have M8 and some C compiler
installed
(if you lack the latter, remove the CompilationTarget->"C" option where it
is present, but
the code will be 2-4 times slower).
(* These first 2 functions courtesy of Norbert Pozar *)
Clear[extractPositionFromSparseArray, positionExtr];
extractPositionFromSparseArray[
HoldPattern[SparseArray[u___]]] := {u}[[4, 2, 2]]
positionExtr[x_List, n_] :=
Flatten@extractPositionFromSparseArray[
SparseArray[Unitize[x - n], Automatic, 1]]
Clear[less];
less =
Compile[{{fst, _Integer, 1}, {sec, _Integer, 1}},
Module[{i = 1},
While[fst[[i]] == sec[[i]], i++];
fst[[i]] < sec[[i]]]];
Clear[positionsOfSameInSorted ];
positionsOfSameInSorted =
Compile[{{fsorted, _Integer, 2}, {ssorted, _Integer, 2}},
Module[{i, j, pctr,
positions =
Table[{0, 0}, {Max[Length[fsorted], Length[ssorted]]}]},
For[i = 1; j = 1; pctr = 0,
i <= Length[fsorted] && j <= Length[ssorted],
If[fsorted[[i]] == ssorted[[j]],
positions[[++pctr]] = {i, j};
i++;
j++,
(* else *)
If[less[fsorted[[i]], ssorted[[j]]],
i++,
(* else *)
j++
]]];
Take[positions, pctr]
], CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}];
Clear[positionsOfSame];
positionsOfSame =
Compile[{{fst, _Integer, 2}, {sec, _Integer, 2}},
Module[{ord1, ord2, ps, psrt},
ord1 = Ordering[fst];
ord2 = Ordering[sec];
psrt = positionsOfSameInSorted[fst[[ord1]], sec[[ord2]]];
If[Length[psrt] != 0,
ps = Transpose[ psrt];
Transpose[{ord1[[ps[[1]]]], ord2[[ps[[2]]]]}],
(* else *)
{{0}}
]
], CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}];
Clear[permuteNTimes ];
permuteNTimes =
Compile[{{list, _Integer, 2}, {permutation, _Integer,
1}, {n, _Integer}},
Nest[#[[All, permutation]] &, list, n]];
Clear[reducePermutationsC];
reducePermutationsC =
Compile[{{list, _Integer, 2}, {permutation, _Integer, 1}},
Module[{perms = list, n = 1, pos = Table[-1, {Length[list]}],
permuted = list, psame = Table[{-1, -1}, {Length[list]}],
done = False},
While[! done,
permuted = permuteNTimes[perms , permutation, n++];
psame = positionsOfSame[perms , permuted];
If[psame != {{0}},
pos =
positionExtr[UnitStep[(#[[2]] - #[[1]]) &@Transpose[psame]], 0];
If[pos != {},
perms = Delete[perms, List /@ pos],
(* else *)
done = True
],
(* else *)
done = True
]];
perms
], {{positionExtr[_, _], _Integer, 1}},
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}, CompilationTarget -> "C"];
Clear[getUniquePermutatations];
getUniquePermutatations[permutations_, actions : {__Cycles}] /;
Length[permutations] > 0 :=
With[{rlen = Range[Length[permutations[[1]]]]},
Fold[reducePermutationsC, permutations,
Map[Permute[rlen, #] &, actions]]]
Now, let's see. The warm-up:
In[17]:= getUniquePermutatations[Permutations[Range[8]],
DeleteCases[GroupElements@DihedralGroup[8],Cycles[{}]]]//Short//Timing
Out[17]=
{0.25,{{1,2,3,5,6,8,7,4},{1,2,3,5,7,4,6,8},{1,2,3,5,7,4,8,6},<<3255>>,{8,6,5,4,3,1,2,7},{8,6,5,4,3,2,1,7}}}
The real deal:
In[18]:= (perms = Permutations[{2,2,2,2,2,2,2,1,1,0,0,0,0,0,0,0}])//Length
Out[18]= 411840
In[26]:= (un = getUniquePermutatations[perms,
DeleteCases[GroupElements@DihedralGroup[16],Cycles[{}]]])//Short//Timing
Out[26]=
{5.469,{{0,0,0,2,0,2,2,0,0,2,2,2,2,1,1,0},<<15183>>,{0,0,0,0,0,0,0,1,1,2,2,2,2,2,2,2}}}
So, about 6 seconds on my machine. If need faster, could probably
parallelize, I did not bother. Hope that
the code is correct. Some simple checks like say this one:
In[27]:= pm =
Permute[#,Cycles[{{2,16},{3,15},{4,14},{5,13},{6,12},{7,11},{8,10}}]]&/@un;
In[28]:= positionsOfSame[un,pm]
Out[28]=
{{11886,11886},{9790,9790},{8547,8547},{7889,7889},{5970,5970},{4821,4821},{4142,4142},{3798,3798},{1333,1333},{105,105}}
suggest that it is: the only elements that are the same in the result and
result permuted according to
some group element are those that are invariant under the transform (have
same positions).
I'd say, it is a good showcase for a new Compile.
Regards and Happy Holidays,
Leonid
Leonid Shifrin
Yaroslav,
Looks like I have produced a chain of incorrect solutions. Here is a
modified one, which
I hope is correct. At least, it is more clear what is going on there. The
idea is that I permute
all the list of permutations with a single group element (transformation) at
once, this can
be done very efficiently with the list[[All, permutation]] idiom. Then, I
use an auxiliary function
positionsOfSame, which finds pairs of positions of same permutations in the
original list
and the permuted one, so for example {{1,3},{2,4}} would mean that first
permutation when
transformed became third one in the original list, and second became fourth.
Then, I use
these pairs as a graph specification (edges), to find connected components -
they will give
the positions in the original list of all permutations which are equivalent
to each other under
a given transformation. I pick the first element (position) in each
component, and delete
from the permutation list all elements (permutations) which are at remaining
positions in
these connected components. This is done for each group element
(transformation), one
after another, with Fold. The result now seems to be independent on the
order in which the
transformations are applied. Starting with the longest cycles is an
optimization, since the
connected components are larger then (of the size of more or less the length
of the cycle),
which means that long permutation list is quickly reduced in size and other
transformations
work with smaller lists and are then faster. Here is the code:
(* Find connected components and a list of positions in the permutation
list, at which positions to
delete elements *)
Clear[getTree, listSplit, gatherBy, getComponents,
getPositionsToDelete];
getTree = Compile[{{pairs, _Integer, 2}},
Module[{t = 0, i = 0, j = 0, xl = 0, yl = 0, k = 0,
len = Max[pairs], dad = Table[0, {Max[pairs]}]},
For[k = 1, k <= Length[pairs], k++, xl = i = pairs[[k, 1]];
yl = j = pairs[[k, 2]];
If[xl == yl, Continue[]];
While[dad[[i]] > 0, i = dad[[i]]];
While[dad[[j]] > 0, j = dad[[j]]];
While[dad[[xl]] > 0, t = xl; xl = dad[[xl]]; dad[[t]] = i];
While[dad[[yl]] > 0, t = yl; yl = dad[[yl]]; dad[[t]] = j];
If[i != j, If[dad[[j]] <= dad[[i]], dad[[j]] += dad[[i]] - 1;
dad[[i]] = j;,(*else*)(dad[[i]] += dad[[j]] - 1);
dad[[j]] = i];];];
For[k = 1, k <= len, k++, If[dad[[k]] <= 0, dad[[k]] = k]];
dad]];
listSplit[x_List, lengths_List] :=
MapThread[Take[x, {##}] &, {Most[#], Rest[#] - 1}] &@
Accumulate[Prepend[lengths, 1]];
gatherBy[lst_List, flst_List] :=
listSplit[lst[[Ordering[flst]]], (Sort@Tally[flst])[[All, 2]]];
getComponents[pairs_] :=
With[{dad = getTree[pairs]},
gatherBy[Range[Length[dad]], FixedPoint[dad[[#]] &, dad]]];
(* Delete all elements except first, in each component *)
getPositionsToDelete[pairs_] :=
Flatten[#[[All, 2 ;; -1]]] &@getComponents[pairs];
(* Functions to find correspondence between positions of same permutations
in the original
and transformed list of permutations *)
(* Main functions *)
Clear[reducePermutationsC];
reducePermutationsC =
Compile[{{list, _Integer, 2}, {permutation, _Integer, 1}},
Module[{perms = list, pos = Table[-1, {Length[list]}],
permuted = list, psame = Table[{-1, -1}, {Length[list]}]},
permuted = perms[[All, permutation]];
psame = positionsOfSame[perms , permuted];
If[psame != {{0}},
pos = getPositionsToDelete[psame];
If[pos != {}, perms = Delete[perms, List /@ pos]]];
perms],
{{getPositionsToDelete[_], _Integer, 1}},
CompilationOptions -> {"InlineCompiledFunctions" -> True,
"InlineExternalDefinitions" -> True}, CompilationTarget -> "C"];
Clear[nonequivalentPermutations];
nonequivalentPermutations[list_, group_] :=
With[{perms = Permutations[list],
actions = DeleteCases[GroupElements@group, Cycles[{}]]},
With[{rlen = Range[Length[perms[[1]]]]},
Fold[reducePermutationsC, perms,
Map[Permute[rlen, #] &, SortBy[actions, Length@#[[1]] &]]]]];
The syntax is now the same as in your solution. Here are examples of use:
In[179]:= nonequivalentPermutations[{2,1,1,0},DihedralGroup[4]]
Out[179]= {{2,1,1,0},{2,1,0,1}}
In[180]:=
nonequivalentPermutations[{2,2,2,2,2,2,2,1,1,1,1,1,1,0,0,0},DihedralGroup[16]]//Short//Timing
Out[180]=
{5.437,{{2,2,2,2,2,2,2,1,1,1,1,1,1,0,0,0},<<30098>>,{2,1,2,0,2,1,2,0,2,1,1,2,1,1,2,0}}}
The resulting number of permutations is somewhat close to an estimate
(30030) given to you by someone at
StackOverflow. The connected component functionality is my port of
implementation of weight-balancing path-compression union-find algorithm
taken straight from the Sedgewick's book on algorithms in C, with some minor
changes as needed for adoption to Mathematica. I also heavily optimized it
for speed. The reason I did not compile it to C is that the bottleneck there
is actually in the gatherBy function (rather than getTree), which can not
alas be Compiled since it produces non-tensor structures (list of sublists
of non-constant length). If it could, we could perhaps get another 30-40 %
speed increase for the entire solution. But even this "home-brewed" gatherBy
is about 4 times faster than the built-in GatherBy for this particular
problem.
Anyway, I really hope that now this is a correct solution, and this is then
my last post on this topic :)
Hope this helps.
Regards,
Maxim
Conceptually, the code that you want is Length@GroupOrbits[group,
Permutations@lst, Permute]. This isn't going to be very efficient
though (I believe it will be improved in future versions of
Mathematica). Here's a brute force approach:
In[1]:= Module[{gens, verts, edges},
gens = PermutationList /@ GroupGenerators@DihedralGroup[16];
verts =
Permutations@{2, 2, 2, 2, 2, 2, 2, 1, 1, 0, 0, 0, 0, 0, 0, 0};
edges = Join @@ (Transpose@{verts, verts[[All, #]]} &) /@ gens;
Length@ConnectedComponents@Graph[Rule @@@ Union@edges]] // Timing
Out[1]= {10.578, 12940}
Maxim Rytin
m...@inbox.ru