Re: Counting
- To: mathgroup at smc.vnet.net
- Subject: [mg114977] Re: Counting
- From: Leonid Shifrin <lshifr at gmail.com>
- Date: Sat, 25 Dec 2010 02:34:29 -0500 (EST)
- References: <201012200540.AAA22760@smc.vnet.net>
Yaroslav,
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
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