Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have a function that generates lists:

f[n_] := Join @@ Permutations /@ PadLeft@IntegerPartitions[n];

For example

f[3]
{{0, 0, 3}, {0, 3, 0}, {3, 0, 0}, {0, 2, 1}, {0, 1, 2}, 
 {2, 0, 1}, {2, 1, 0}, {1, 0, 2}, {1, 2, 0}, {1, 1, 1}}

f[3] returns 10 lists of three elements, which I want to classify this according to certain laws;

{0, 0, 3}, {0, 3, 0} and {3, 0, 0} belong one class with 3 members because RotateLeft can transform any one of these into any of the others.

Similarly, {0,1,2}, {1,2,0} and {2,0,1} belong one class with 3 members.

{0, 2, 1}, {2, 1, 0} and {1, 0, 2} belong one class with 3 members.

{1, 1, 1} belongs to a class of which it is the only member.

So the result I want is

enter image description here

I have code that will generate the desired result

m = AbsoluteTime[];
f[n_Integer] := Join @@ Permutations /@ PadLeft@IntegerPartitions[n];
f2[list_] := Module[{n = 1, m = 0},
  While[n <= Length[data1], 
   If[Or @@ (SameQ[data1[[n, 1]], #] & /@ 
       NestList[RotateLeft, list, Length[list] - 1]), data1[[n, 2]]++;
     m++; Break[]];
   n++];
  If[m == 0, AppendTo[data1, {list, 1}]]]
data1 = {{First@f[7], 1}};
f2 /@ Drop[f[7], 1];
Length[f[7]]
data1 // Length
AbsoluteTime[] - m

but it's too slow. When n = 7, it takes about 11 seconds, but I want to work with n is the range from 14 to 20. how can I improve efficiency the efficiency of my code?

share|improve this question
    
Notice, that there is unscoped data1 in f2 definition. –  Kuba Oct 25 '13 at 8:59
    
Hey Chenminqi, I just wanted to ask if my (rather long) answer was of any use to you. Nice to see you are still active on the site :) –  Jacob Akkerboom Nov 5 at 12:21

2 Answers 2

up vote 5 down vote accepted

I believe I'm missing some basic function that can easily reproduce the test here. But that's a detail, the main thing is to use Tally with the second argument specification:

test[x_List] := Sort@Array[RotateLeft[x, #] &, Length[x]]

Tally[f[3], test[#1] == test[#2] &] // MatrixForm

$\left( \begin{array}{cc} \{0,0,3\} & 3 \\ \{0,2,1\} & 3 \\ \{0,1,2\} & 3 \\ \{1,1,1\} & 1 \\ \end{array} \right)$

It's not the optimal way, for n = 10 it takes about 5 sec on my old pc. I will think about improvements.

share|improve this answer
    
@user10193 Please hold on with an accept a little, it is a good habit in order to not discourage others. I'm 100% sure that better answers will appear here. :) –  Kuba Oct 25 '13 at 8:57
    
Tkx very much.I considered about tally but I just know use Tally[list] but Tally[list, test].so thank you. –  Chenminqi Oct 25 '13 at 9:01
    
@Kuba do you mean First@Sort in test? Now you end up comparing lists of lists, whereas you only need to compare one representative. I agree there should probably be some trick to do this step :) –  Jacob Akkerboom Oct 25 '13 at 10:58
    
@JacobAkkerboom No I didn't, but good idea with First it may help with long lists. I was thinking about something like built in testing if two lists are related with an odd/even permutation. But I must admit my math skils are not existant in this area. I don't even know if what I'm thinking about is called even permutation in english :). I was trying to avoid to many Sorting with MemberQ[test[#],#2]& as the second argument but it seems it is longer. –  Kuba Oct 25 '13 at 11:10
    
@Kuba it is an interesting problem. Odd/even permutations are a big thing, so its probably the right word ;). If you can tackle it like that I'd be curious. For me the main problem becomes finding that the "lowest subsequence" of length 3 of {0,1,2,0,0,3} is {0,0,3} and not {0,1,2}. Hrm :P –  Jacob Akkerboom Oct 25 '13 at 11:38

Introduction

The key to concept to tackling this problem is rotational symmetry. Using this, it turns out that we can find out how many elements your classes have relatively easily (in terms of memory and CPU usage, it does take some math). Actually finding and storing the elements that represent the classes is more of a programming challenge than a mathematical one, and this is what will cost most computational resources. The algorithm I present finds exactly one representative of each class, so in this way the problem is nicely cut into two parts.

Theory

Now let's do some group theory. Throughout the text of this answer, let $n$ be the n of your question, which is the integer for which we want to describe such classes. In my code I will call this variable in for input. It turns out the classes you describe are orbits of the group $\mathbb{Z}_n$ acting on the set of weak compositions of $n$. Let's give the set of weak compositions a name

$$X_n = \lbrace \vec x \in (\mathbb{Z}_n)^n | \sum_{k=1}^n x_k = n \rbrace$$

We wish to know how many elements there are in any particular class a.k.a. orbit. The section on wikipedia article that explains about orbits gives a result that we can use here, in its subsection Orbit-stabilizer theorem and Burnside's lemma. The result states

$$|G.x| = |G|\; /\; |G_x|$$

On the LHS is the quantity we are after, i.e. the number of elements of the orbit of $x$. Formally, $G_x$ is called the stabiliser of x, a notion which is also explained in the article. However, we can forget about stabilisers and look directly at $|G_x|$. $|G_x|$ is the number of elements of $G$ for which $g.x =x$. For our group $G = \mathbb{Z}_n$ this is something quite intuitive. For example for $n = 6$ and

$$x=(1,1,1,1,1,1)$$

we see that all elements of $G$ map x to x. So $|G_x| = n = 6$. So for this particular $x$, $|G.x| = |G|\; /\; |G_x| = n/n = 1$. Of course it was quite clear that the class to which $x$ belongs has only one element. Now consider

$$x=(0,0,3,0,0,3)$$

the class to which this $x$ belongs has 3 elements, because there are 2 ways we can left shift this onto itself. We can shift by 0 or 3 spots. So in this case $|G.x| = 6/2 = 3$. In this last example let us say that $x$ has period 3, because the smallest repeating subsequence that comprises $x$ is $(0,0,3)$ which has length 3. It should be clear that if we know the period of a sequence $x$, it is easy to find the length of its orbit.

There is another advantage to considering the periods of sequences. We can partition the set of orbits by their periods. Then finding all the classes of which the elements have a certain period can be transformed into the original problem. If for $n = 6$, we wish to find all sequences with period 3, we simply have to solve the same problem for $n = 3$. That is, for $n=3$ we have the sequences

$$ (1,1,1)\quad (0,1,2) \quad (0,2,1)\quad (0,0,3)$$

and so for $n=6$ we have the corresponding sequences

$$ (1,1,1,1,1,1)\quad (0,1,2,0,1,2) \quad (0,2,1,0,1,2)\quad (0,0,3,0,0,3)$$

But both sequences that consist only of 1's have period 1 and not 3. So in the case for $n=3$, we have to consider the trivial subproblem for $n=1$. To get a proper one-to-one correspondence between such sequences, we should only consider the aperiodic elements in the subproblems. Here, by aperiodic I mean that an element of $x \in X_n$ is aperiodic if the period of $x$ is $n$.

One last nice bit of theory is the following. We can also derive how many orbits there are. It is nice to know how many orbits there will be, as this will tell us how long our arrays will have to be. We will first only consider the case where $n$ is prime, which will lead to a useful result when we work with aperiodic sequences. For $n$ prime, we can easily relate the number of orbits to $|X_n|$. Let $Y_n$ be the set of orbits. As every orbit contains $n$ elements, as $n$ is prime, we have

$$|Y_n| = |X_n|\; /\; n $$

Now recall that $X_n$ is the set of weak compositions of $n$. In the section Number of compositions of this wiki article, we find that the number of weak compositions is $\binom{2 n -1}{n -1}$ (use the formula with $k=n$). So we have

$$ |Y_n| = \binom{2 n -1}{n -1}\; /\; n $$

For $n$ which are not prime, partition $X_n$ by period and count the elements in each of the parts. Maybe now the following function explains itself. The idea is that in later code there will be a big array main, consisting of sequences of length $n$ (in in the code). In this array the elements are ordered by period. Maybe the first f[1] will have period $n$ and the elements f[1]+1;;f[2] will have period $n/2$. We have that outSize is the length of array main and that subListLength gives the length of the parts (f[i_]:=subListLength[[i]]).

ClearAll[mInitialization, divs, divDivisors, divLen, terms, bins, 
outSize]
mInitialization[in_] :=
  (
   divs = Divisors[in];

   divDivisors =
    Outer[If[# == #2, False, Divisible[#, #2]] &, divs, divs];

   divLen = Length@divs;

   terms = ConstantArray[0, divLen];

   bins = Binomial[2 # - 1, #] & /@ divs;

   Do[
    terms[[iter]] = 
     bins[[iter]] - 
      Sum[If[divDivisors[[iter, kkkk]], terms[[kkkk]], 0], {kkkk, 
        divLen}]

    ,
    {iter, divLen}

    ];

   subListLengths = Reverse[terms/divs];

   outSize = Total[subListLengths];

   );

Compiled function

I will not go into the details of finding the representations of the classes. There is a little bit more interesting math here and maybe I will post that later.

First let's make a tidy overview of our variables. The variables in hConsts are basically values that will inlined in the compiled function.

hConsts = Hold[walkedBack, walkedForward, maxDepth];
hLocals = 
  Hold[test, depth, direction, total, iters, (*activeGap,*) 
   distanceToGap, initVal, candidate, storage, storeCount];
ClearAll @@ hLocals;
ClearAll @@ hConsts;
walkedBack = 0;
walkedForward = 1;
maxDepth := in - 1

Here is the initialization of the compiled function we will use. It is not actually used as a function, but it is inlined, as it will be part of a compiled function.

(*OwnValues of initialization are used*)
initialization :=
 (
  iters = Table[0, {in - 1}];
  iters[[in - 1]] = in;
  total = 0;

  depth = maxDepth;

  test = True;

  (*activeGap=0;*)
  direction = walkedForward;
  distanceToGap = 0;

  storage = ConstantArray[0, {outSize, in}];

  storeCount = 0;
  )

Here is the code that is used to store a sequence in our array. This code is also inlined in the compiled function.

store[] := 
  (
   Increment[storeCount];
   Do[
    storage[[storeCount, kkkk + 1]] = iters[[kkkk]]
    ,
    {kkkk, in - 1}
    ];
   );

This is the main code of the compiled function.

codezzz :=
 (
  While[
   test
   ,
   If[
    depth < maxDepth
    ,
    If[
     direction == walkedBack
     ,


     If[
      total + iters[[depth]] <= 
       in - 2 (*we know a sequence cannot end in a 0, 
      so we need space for at least a 1. This argument should be \
expanded.*)
      ,

      Increment[iters[[depth]]];
      total += iters[[depth]];
      Increment[depth];
      direction = walkedForward;
      (*activeGap=depth;*)
      distanceToGap = 0;

      (*after we have walked back, we always start with a 0, 
      unless depth \[Equal] maxDepth*)

      ,
      Decrement[depth];
      direction = walkedBack;
      (*as we added 1 before prematurely*)

      If[
       depth == 0
       ,
       test = False;
       ,
       total = total - iters[[depth]];
       ]
      ]

     ,
     (*if cont: walkedForward*)

     If[
      distanceToGap > 0,
      initVal = iters[[distanceToGap]],
      initVal = 0;
      ];


     If[
      total + initVal <= in - 1
      ,
      (*initialization*)
      iters[[depth]] = initVal;

      Increment[depth];
      direction = walkedForward;
      Increment[distanceToGap];
      total = total + initVal
      ,
      Decrement[depth];
      direction = walkedBack;
      total -= iters[[depth]]

      ]


     ]

    ,

    If[
     distanceToGap == 0,

     iters[[depth]] = in - total;
     store[]
     ,

     candidate = in - total;

     If[
      candidate > iters[[distanceToGap]]
      ,
      iters[[depth]] = candidate;
      store[]
      ]

     ];
    Decrement[depth];
    direction = walkedBack;
    (*activeGap=0;*)
    total -= iters[[depth]]

    ]

   ];
  storage
  )

This code creates the compiled function

cfu =
  ReleaseHold[
     Unevaluated[
      Compile[{{in, _Integer}, {outSize, _Integer}},
       Block[
        #1
        ,
        #2
        ]
       , CompilationTarget :> "C"
       ]
      ]
     ] & @@
   {
    List @@@ Hold[Evaluate[hLocals]]
    ,
    (
      (
        Hold[
          initialization;
          codezzz
          ] /. Join[OwnValues[codezzz], OwnValues[initialization]]
        ) /. Apply[Join, List @@ OwnValues /@ hConsts]
      ) /. DownValues[store]
    };

The compiled function finds only aperiodic sequences. We need more functions! This function handles cacheing, which is especially useful as we need to solve subproblems. The W in cfuW is for wrapper.

<<Developer`
cached[x_] := False;
ClearAll[cfuW];
cfuW[inSub_ /; inSub > 2, outSize_, isAperiodicSize_] :=

  If[isAperiodicSize,
   If[
    TrueQ@cached[inSub],
    cache[inSub],
    {
      cache[inSub] = cfu[inSub, outSize],
      cached[inSub] = True
      } // First
    ],
   cfu[inSub, outSize]
   ];
cfuW[1, 1, True] := {{1}} // ToPackedArray;
cfuW[2, 1, True] := {{0, 2}} // ToPackedArray;

Another wrapper to localize some variables in the aperiodic case

findAperiodicComps[inSub_] :=
 Block[
  {in, outSize, divs, terms, subListLengths},
  in = inSub;
  mInitialization[in]; 
  cfuW[in, Last[terms]/Last[divs], True]
  ]

Finally, the function that deals with the general case.

ClearAll[findComps]
findComps[inSub_] :=
 Block[
  {in, main, rest, spanLists, divs, terms, subListLengths}
  ,
  in = inSub;
  mInitialization[in]; (*sets outSize etc*)

  main = ConstantArray[0, {outSize, in}];
  main[[;; Last[terms]/Last[divs]]] = 
   cfuW[in, Last[terms]/Last[divs], True];
  rest =
   cfuW[#, #2, True] &
    @@@
    Reverse[
     Delete[
      Transpose[
       {
        divs,
        terms/divs
        }
       ]
      ,
      {{-1}}
      ]
     ];
  spanLists =
   Transpose[
    {1 + Delete[#, -1], Delete[#, 1]} &@
     Accumulate@Reverse[terms/divs]
    ];
  Function[
    main[[Span @@ #3]] = ArrayPad[
      #
      ,
      {0, {0, in - #2}}
      ,
      "Periodic"
      ]
    ] @@@
   Transpose[
    {rest, Reverse[Delete[divs, -1]], spanLists}

    ];
  {main, divs, subListLengths}
  ]

Functions to conform to OPs format

The format required by the OP forces the unpacking of arrays. But it is useful to show how get the data in this format.

take[result_, terms_, subListLengths_, a_ ;; b_] :=
 Transpose[
  {
   result[[a ;; b]]
   ,
   Flatten[
     ConstantArray[#, #2] & @@@ 
      Transpose[{Reverse@terms, subListLengths}], 1][[a ;; b]]
   }
  ]

take[result_, terms_, subListLengths_, All] :=
 Transpose[
  {
   result
   ,
   Flatten[
    ConstantArray[#, #2] & @@@ 
     Transpose[{Reverse@terms, subListLengths}], 1]
   }
  ]

takeSimple[in_] := take[##, All] & @@ findComps[in]

Result

We now have

takeSimple[3] // MatrixForm

enter image description here

Which has a different ordering than both the OPs code and Kuba's code, but sure enough is correct.

Comparison set up

f[n_]:=Join@@Permutations/@PadLeft@IntegerPartitions[n];
test[x_List]:=Sort@Array[RotateLeft[x,#]&,Length[x]]
representation[x_List] := 
 First@test[x];
kuba[n_]:=Tally[f[n],test[#1]==test[#2]&]

Actual comparison

n = 10;
(resKubaShort = kuba[n][[All, 1]]) // Timing // First
(resJacobShort = findComps[n][[1]]) // Timing // First
 2.231080  
 0.004269
SameQ @@ Sort /@ {representation /@ resKubaShort, resJacobShort}
True

This also indicates that all the results in resJacobShort are canonical. For what its worth.

A timing for higher n would be

n = 15;
(resJacobShort = findComps[n][[1]]) // Timing // First
3.690403

I think I can also still run it for n=16. But for higher n, we run out of memory. Memory usage in this code does not display much CS skill. We store integers that cannot be larger than n as 64 bit integers. No effort is made to for example store part of the list on the hard drive to load it on demand. Maybe someday..

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.