Services & ResourcesWolfram ForumsMathGroup Archive

[Date Index] [Thread Index] [Author Index]

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>