Re: Counting
- To: mathgroup at smc.vnet.net
- Subject: [mg114979] Re: Counting
- From: Leonid Shifrin <lshifr at gmail.com>
- Date: Sun, 26 Dec 2010 04:01:06 -0500 (EST)
- References: <201012200540.AAA22760@smc.vnet.net>
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 *)
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}];
(* 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,
Leonid
On Mon, Dec 20, 2010 at 8:40 AM, Yaroslav Bulatov <yaroslavvb at gmail.com>wrote:
> I'd like to count the number of permutations of {2, 2, 2, 2, 2, 2, 2,
> 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]]
>
>
- References:
- Counting
- From: Yaroslav Bulatov <yaroslavvb@gmail.com>
- Counting