25
$\begingroup$

This input:

Permutations[Range[12]]

Results in this (error) output:

Permutations::outsize: 
  The result of evaluating Permutations[{1,2,3,4,5,6,7,8,9,10,11,12}] 
  would be a packed array with 5748019200 elements, but the number of elements 
  in a packed array must be a machine integer. >>

That number (5748019200) is interesting, because it's exactly 12 times 12! (that's factorial, not exclamation point)

Presumably, Mathematica is trying to store all 12! length 12 lists in a single monolithic array. I can imagine this failing.

Usually, Mathematica shields me from these types of problems. For example, I had no trouble calculating 12*12!.

My intention is to Select some elements from this list, so I don't need to have every permutation in memory at once.

Question: Is there a different way to generate the permutations that avoids this problem?

$\endgroup$
3
  • 1
    $\begingroup$ Run the same code on a 64-bit machine ;-) Seriously though, it is possible to iterate over permutations but I'm not sure if the algorithm is implemented directly in Mathematica. $\endgroup$
    – David Z
    Commented Feb 4, 2012 at 3:28
  • 3
    $\begingroup$ @David It doesn't work on a 64-bit machine either. Harold: Apparently the maximum size of a packed array is 2^31-1. This is a huge array: a 2^31 - 1-element packed array of machine integers (i.e. assuming the most efficient storage possible) would take up 8 GB of contiguous memory. You need to have a lot of memory in your computer to be able to store such an array (definitely much more than 8 GB because of the contiguous memory block requirement) $\endgroup$
    – Szabolcs
    Commented Feb 4, 2012 at 20:13
  • 3
    $\begingroup$ @Szabolcs 8 GB is quite a bit of memory, but modern operating systems, on top of which the Mathematica kernel runs, allocate memory in page-sized chunks (most commonly 4 kilobytes each), and (virtual) memory region which appears contiguous to the Mathematica kernel doesn't need any special attention on the operating system level. It just requires lots of RAM, or swap space. Contiguous memory problem existed roughly up to mid-nineties, but not really afterwards. $\endgroup$
    – kirma
    Commented Oct 2, 2014 at 16:18

7 Answers 7

19
$\begingroup$

Combinatorica` has the function NextPermutation which allows you to iterate over the permutations. There may be ways of generating a smaller subset if you have more information about what you are looking for.

$\endgroup$
7
  • 7
    $\begingroup$ On that note: the algorithms in Combinatorica are based on old FORTRAN routines discussed in this book; OP might want to take a look at the book and see what other strategies might be appropriate for his circumstances. $\endgroup$ Commented Feb 4, 2012 at 3:48
  • $\begingroup$ @J.M., that book looks to be a tremendous resource. I'm glad I asked this question now. $\endgroup$
    – Harold
    Commented Feb 4, 2012 at 4:52
  • 2
    $\begingroup$ @Harold note, the functionality in Combinatorica` is being slowly folded into the kernel functions. So, loading it will give a warning message to that effect, but NextPermutation doesn't (yet?) seem to have been moved over. $\endgroup$
    – rcollyer
    Commented Feb 4, 2012 at 4:55
  • $\begingroup$ @J.M. Thanks a lot for the link to that book. Looks excellent $\endgroup$
    – user1066
    Commented Feb 4, 2012 at 9:34
  • 1
    $\begingroup$ @rcollyer NextPermutation has in a way already been moved over: you can use GroupElements to emulate it. $\endgroup$ Commented Oct 1, 2014 at 13:20
27
$\begingroup$

Since Mathematica 8 it is possible generate the elements of any group one by one with GroupElements. Here's for example a randomly chosen element of the permutation group on 20 elements:

GroupElements[SymmetricGroup[20], {10^6 + 1}]
{Cycles[{{11, 13, 19}, {12, 18, 17, 16}, {15, 20}}]}

The result is immediately; there's no need to build up the full group of $20! \approx 2.4 \cdot 10^{18}$ elements. We can use this to scan over the permutations of a given list without generating them all at once:

ScanPermutations[function_, list_] /; Length[list] <= 8 :=
  Scan[function, Permutations @ list];

ScanPermutations[function_, list_] :=
  Do[
    Scan[
      function @ Permute[list, #] &,
      GroupElements[
        SymmetricGroup @ Length @ list,
        Range[8! * (i - 1) + 1, 8! * i]
      ]
    ],
    {i, Length[list]! / 8!}
  ];

The code splits up the permutations in blocks of $8! = 40320$ whenever the given list has more than 8 elements. This is about ten times faster than calling GroupElements to generate just one permutation at a time.

Printing e.g. the permutations of {1,2,3} can now be done as follows:

ScanPermutations[Print, {1,2,3}]
{1,2,3}
{1,3,2}
{2,1,3}
{3,1,2}
{2,3,1}
{3,2,1}

Going one step further, we may use Reap and Sow to select permutations that match a given criterion:

SelectPermutations[list_, crit_] :=
  First @ Last @ Reap @ ScanPermutations[
    If[crit @ #, Sow @ #] &,
    list
  ];

And finally, here's a small example that puts the above in action. It selects the even permutations of {1,2,3}:

SelectPermutations[{1, 2, 3}, Signature[#] === 1 &]
{{1, 2, 3}, {2, 3, 1}, {3, 1, 2}}
$\endgroup$
1
  • 2
    $\begingroup$ I appreciate this late answer. This technique looks promising and has several attractive properties when compared with some of the answers above. Neat stuff. Thanks! $\endgroup$
    – Harold
    Commented Oct 1, 2014 at 22:41
12
$\begingroup$

Chunks of permutations

Here is a LibraryFunction implementation of "CoolMulti" algorithm generating permutations of multisets. The algorithm is described in:

  • A. Williams: Loopless generation of multiset permutations by prefix shifts
/* permutations.c */
#include "WolframLibrary.h"

DLLEXPORT mint WolframLibrary_getVersion() {
    return WolframLibraryVersion;
}
DLLEXPORT int WolframLibrary_initialize(WolframLibraryData libData) {
    return LIBRARY_NO_ERROR;
}
DLLEXPORT void WolframLibrary_uninitialize(WolframLibraryData libData) {}

DLLEXPORT int nextPermutationsChunk(
        WolframLibraryData libData, mint Argc, MArgument *Args, MArgument Res
) {
    /* Values tensor: multiset of integers in non-increasing order. */
    MTensor valuesT = MArgument_getMTensor(Args[0]);
    /* Actual data of values tensor. */
    mint* values = libData->MTensor_getIntegerData(valuesT);
    /* Number of elements in multiset. */
    mint n = libData->MTensor_getDimensions(valuesT)[0];

    /* `stateT` tensor: `{next1, next2, ..., head, ref}`. */
    MTensor stateT = MArgument_getMTensor(Args[1]);
    /*
     * First `n` elements of `next` array contain indices of next nodes
     * in emulated linked list. Other elements of `stateT` tensor are used
     * only through direct pointers.
     */
    mint* next = libData->MTensor_getIntegerData(stateT);
    /* Pointer to index of head node. */
    mint* head = next + n;
    /* Pointer to index of reference node. */
    mint* ref = head + 1;

    /* Number of permutations in returned chunk. */
    mint chunkSize = MArgument_getInteger(Args[2]);
    /* Dimensions of returned `chunk` tensor. */
    mint chunkDims[2] = {chunkSize, n};
    /* 2 dimentional tensor with chunk of permutations to be returned. */
    MTensor chunkT;
    libData->MTensor_new(MType_Integer, 2, chunkDims, &chunkT);
    /* Actual data of the chunk tensor. */
    mint* chunk = libData->MTensor_getIntegerData(chunkT);

    mint i;
    for (i = 0; i < chunkSize; i++) {
        /*
         * Based on:
         * Aaron Williams. 2009. Loopless generation of multiset permutations
         * using a constant number of variables by prefix shifts.
         * http://webhome.csc.uvic.ca/~haron/CoolMulti.pdf
         */
        mint afterRef = next[*ref];
        mint localRef;
        if (next[afterRef] < n && values[*ref] >= values[next[afterRef]]) {
            localRef = afterRef;
        } else {
            localRef = *ref;
        }
        mint newHead = next[localRef];

        next[localRef] = next[newHead];
        next[newHead] = *head;

        if (values[newHead] < values[*head]) {
            *ref = newHead;
        }
        *head = newHead;

        /* Populate i-th permutation in chunk. */
        mint j, index;
        for (j = 0, index = *head; j < n; j++) {
            chunk[i*n + j] = values[index];
            index = next[index];
        }
    }

    /* Return control over state tensor back to Wolfram Language. */
    libData->MTensor_disown(stateT);

    /* Set chunk tensor as returned value of LibraryFunction. */
    MArgument_setMTensor(Res, chunkT);

    return LIBRARY_NO_ERROR;
}

Save above code in permutations.c file in same directory as current notebook, or paste it as a string, instead of {"permutations.c"}, as first argument of CreateLibrary in code below. Pass, in "CompileOptions", appropriate optimization flags for your compiler, the ones below are for GCC.

Needs@"CCompilerDriver`"
SetDirectory@NotebookDirectory[];
CreateLibrary[{"permutations.c"}, "permutations"(*,
    "CompileOptions" -> "-Wall -march=native -O3"*)
]
nextPermutationsChunk = LibraryFunctionLoad[%, "nextPermutationsChunk",
    {{Integer, 1, "Constant"}, {Integer, 1, "Shared"}, Integer},
    {Integer, 2}
]

Original "CoolMulti" algorithm works by manipulating singly linked list containing elements of multiset. To generate next permutation algorithm requires said linked list and pointers to two of list nodes from previous iteration. In this implementation we emulate singly linked list using two arrays.

Constant array of values is passed as first argument of nextPermutationsChunk. It should be non-empty list of integers in non-increasing order.

Array representing "state" that needs to be passed between iterations, containing indices of "next" nodes, of emulated linked list, and indices of two "special" nodes, needed between iterations, is passed as second argument of nextPermutationsChunk. It should contain n+2 (where n is number of elements in multiset) non-negative integers.

Third argument of nextPermutationsChunk is number of permutations in returned chunk.

nextPermutationsChunk function returns list containing chunk of permutations of given multiset. As a side effect it modifies list given in second argument, so that it correctly represents "state of generator" after generating returned chunk.

Initial "generator state" is list containing consecutive integers from 1 to n, index of "head node": 0, and index of penultimate node: n-2. For example let's generate all 12 permutations of {4, 2, 1, 1} multiset in two chunks of 5 elements and third chunk containing 2 elements.

values = {4, 2, 1, 1};
state = Join[Range@4, {0, 2}];
(* {1, 2, 3, 4, 0, 2} *)

nextPermutationsChunk[values, state, 5]
(* {{1, 4, 2, 1}, {4, 1, 2, 1}, {1, 4, 1, 2}, {1, 1, 4, 2}, {4, 1, 1, 2}} *)
state
(* {3, 4, 1, 2, 0, 2} *)

nextPermutationsChunk[values, state, 5]
(* {{2, 4, 1, 1}, {1, 2, 4, 1}, {2, 1, 4, 1}, {1, 2, 1, 4}, {1, 1, 2, 4}} *)
state
(* {4, 0, 1, 2, 3, 2} *)

nextPermutationsChunk[values, state, 2]
(* {{2, 1, 1, 4}, {4, 2, 1, 1}} *)
state
(* {1, 3, 4, 2, 0, 2} *)

With simple example from Leonid's answer we get:

values = Reverse@Range@14;
state = Join[Range@14, {0, 12}];
ctr = 0;
While[(nxt = nextPermutationsChunk[values, state, 10000]) =!= {} && ctr < 5*10^6,
    ctr += Length[nxt]
]; // AbsoluteTiming
ctr
(* {0.098408, Null} *)
(* 5000000 *)

which is ten times faster than JLink` based solution, that on my computer gives {1.00089, Null}, but of course, as is pointed out in said answer, processing time is anyway likely to dominate over generation time.


All permutations

Using our LibraryFunction we can define a function that returns all permutations.

permutations // ClearAll
permutations[expr : _@Repeated[_, {0, 1}]] := {expr}
permutations[list_List /; VectorQ[Unevaluated@list, IntegerQ]] :=
    With[{n = Length@list},
        nextPermutationsChunk[
            Reverse@Sort@list,
            Join[Range@n, {0, n - 2}],
            Multinomial @@ Tally[list][[All, 2]]
        ]
    ]
permutations[expr_ /; Not@AtomQ@Unevaluated@expr] :=
    With[{n = Length@expr, tallied = Sort@Tally@(List @@ expr)},
        With[{unique = Head@expr @@ tallied[[All, 1]]},
            unique[[#]] & /@ nextPermutationsChunk[
                Flatten@Reverse@
                    MapIndexed[ConstantArray[First@#2, Last@#1]&, tallied],
                Join[Range@n, {0, n - 2}],
                Multinomial @@ tallied[[All, 2]]
            ]
        ]
    ]

Check that it generates same permutations as built-in Permutations for integer lists:

And @@ (Sort@permutations@# === Sort@Permutations@# & /@ Join @@ (Tuples[Range@#, #] & /@ Range@6))
(* True *)

and arbitrary symbolic expressions:

ClearAll[f, g]
And @@ (Sort@permutations@# === Sort@Permutations@# & /@ Join @@ (g @@@ Tuples[f /@ Range@#, #] & /@ Range@6))
(* True *)

For lists of unique integers speed and memory usage of our permutations is almost the same as of built-in Permutations:

tmp = Join[ConstantArray[0, 1], Range@9];
Permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.095855, 290304520} *)
permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.09785,  290305440} *)

For larger lists with more duplicates our new permutations becomes much faster than built-in:

tmp = Join[ConstantArray[0, 21], Range@5];
Permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {1.02257,  1641870472} *)
permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.563034, 1641870848} *)

tmp = Join[ConstantArray[0, 117], Range@3];
Permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {5.73228,  1617644568} *)
permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.575785, 1617641984} *)

tmp = Join[ConstantArray[0, 15000], Range@1];
Permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {2446.53,  1800960704} *)
permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.766165, 1800481744} *)

For symbolic expressions with unique arguments built-in Permutations are few times faster and use less memory.

tmp = {a, b, c, d, e, f, g, h, i};
Permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.122717,  43546144} *)
permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {0.435274, 127736464} *)

Since we don't have direct access to symbolic expressions from library link permutations use top level mapping and unpacks one level of integer array returned by LibraryFunction. I've tried using WSTP in library link, but all solutions I tried were slower than top level mapping used in permutations.

Nonetheless for sufficiently large array with many duplicates permutations catches up:

tmp = Join[ConstantArray[a, 90], {b, c, d}];
Permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {2.89692,  654026416} *)
permutations@tmp // MaxMemoryUsed // AbsoluteTiming (* {2.74743, 1376560792} *)

Lazy permutations

Using Streaming` package we can implement lazy permutations:

Needs@"Streaming`"

lazyPermutationsInternal = Function[{values, n, length, chunkSize, opts},
    Module[{active = False, left = length, state = Join[Range@n, {0, n - 2}]},
        LazyListCreate[
            IteratorCreate[
                ListIterator,
                (active = True)&,
                With[{realChunkSize = Min[chunkSize, left]},
                    If[realChunkSize === 0,
                        {}
                    (* else *),
                        left -= realChunkSize;
                        nextPermutationsChunk[values, state, realChunkSize]
                    ]
                ]&,
                TrueQ@active&,
                Remove[active, left, state]&
            ],
            chunkSize,
            opts,
            "Length" -> length,
            "FiniteList" -> True
        ]
    ]
];

lazyPermutations // ClearAll
lazyPermutations[
    expr : _@Repeated[_, {0, 1}],
    chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]
] :=
    LazyListCreate[{expr}, chunkSize, opts]
lazyPermutations[
    list_List /; VectorQ[Unevaluated@list, IntegerQ],
    chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]
] :=
    lazyPermutationsInternal[
        Reverse@Sort@list,
        Length@list,
        Multinomial @@ Tally[list][[All, 2]],
        chunkSize,
        {opts}
    ]
lazyPermutations[
    expr_ /; Not@AtomQ@Unevaluated@expr,
    chunkSize : _Integer?Positive : 10^5, opts : OptionsPattern[]
] :=
    With[{tallied = Sort@Tally@(List @@ expr)},
        With[{unique = Head@expr @@ tallied[[All, 1]]},
            unique[[#]] & /@ lazyPermutationsInternal[
                Flatten@Reverse@
                    MapIndexed[ConstantArray[First@#2, Last@#1]&, tallied],
                Length@expr,
                Multinomial @@ tallied[[All, 2]],
                chunkSize,
                {opts}
            ]
        ]
    ]

Check that lazyPermutations generates same permutations as built-in Permutations for integer lists:

And @@ (LazyListBlock[Sort@Normal@lazyPermutations[#, 50] === Sort@Permutations@#] & /@ Join @@ (Tuples[Range@#, #] & /@ Range@5))
(* True *)

and arbitrary symbolic expressions:

ClearAll[f, g]
And @@ (LazyListBlock[Sort@Normal@lazyPermutations[#, 50] === Sort@Permutations@#] & /@ Join @@ (g @@@ Tuples[f /@ Range@#, #] & /@ Range@5))
(* True *)

lazyPermutations generates LazyList of permutations that can be iterated in chunks of given length. Here we generate LazyList of permutations of f[a, d, a, b] with 3-element chunks:

lazyPermutations[f[a, d, a, b], 3]
(* « LazyList[f[a,d,b,a],f[d,a,b,a],f[a,d,a,b],...] » *)
Scan[Print, %]
(* f[a,d,b,a]
   f[d,a,b,a]
   f[a,d,a,b]
   f[a,a,d,b]
   f[d,a,a,b]
   f[b,d,a,a]
   f[a,b,d,a]
   f[b,a,d,a]
   f[a,b,a,d]
   f[a,a,b,d]
   f[b,a,a,d]
   f[d,b,a,a] *)

We can control memory needed to iterate over LazyList by setting appropriate chunk size:

Scan[Identity, lazyPermutations[Range@11, 10^4]] // MaxMemoryUsed // AbsoluteTiming // LazyListBlock
(* {55.3121,   9360664} *)
Scan[Identity, lazyPermutations[Range@11, 10^5]] // MaxMemoryUsed // AbsoluteTiming // LazyListBlock
(* {40.8806,  38275680} *)
Scan[Identity, lazyPermutations[Range@11, 10^6]] // MaxMemoryUsed // AbsoluteTiming // LazyListBlock
(* {38.9187, 376096448} *)
$\endgroup$
10
$\begingroup$

Consider than the permutations of {1, 2, 3, 4, 5} are each of the permutations of {1, 2, 3, 4} with 5 inserted at each possible place. One can therefore examine the permutations of {1, 2, 3, 4, 5} in blocks like this:

p4 = Permutations@Range@4;

Table[
  ReplaceList[x, {h___, t___} :> {h, 5, t}],
  {x, p4}
]

For example, making a certain selection:

Table[
  Select[
    ReplaceList[x, {h___, t___} :> {h, 5, t}],
    # + #2 - #3*#4/#5 > 7 & @@ # &
  ],
  {x, p4}
]

The same can be applied to the permutations of Range@12.

$\endgroup$
9
$\begingroup$

Migrated from here

Mathematica is quite likely to be slow at generating permutations. Here is a Java-based approach. It is based on Java reloader and this nice Java code, which I slightly extended with a getNextMultiple method.

So, load the Java reloader first (run the code from that post). Then, execute the following:

JCompileLoad @ "
 import java.util.Arrays;

 public class PermutationGenerator {
  private int[] array;
  private int firstNum;
  private boolean firstReady = false;

  public PermutationGenerator(int n, int firstNum_) {
     if (n < 1) {
         throw new IllegalArgumentException(\"The n must be min. 1\");
     }
     firstNum = firstNum_;
     array = new int[n];
     reset();
  }

  public void reset() {
     for (int i = 0; i < array.length; i++) {
         array[i] = i + firstNum;
     }
     firstReady = false;
  }

  public boolean hasMore() {
     boolean end = firstReady;
     for (int i = 1; i < array.length; i++) {
         end = end && array[i] < array[i-1];
     }
     return !end;
  }

  public int[] getNext() {

     if (!firstReady) {
         firstReady = true;
         return array;
     }

     int temp;
     int j = array.length - 2;
     int k = array.length - 1;

     // Find largest index j with a[j] < a[j+1]

     for (;array[j] > array[j+1]; j--);

     // Find index k such that a[k] is smallest integer
     // greater than a[j] to the right of a[j]

     for (;array[j] > array[k]; k--);

     // Interchange a[j] and a[k]

     temp = array[k];
     array[k] = array[j];
     array[j] = temp;

     // Put tail end of permutation after jth position in increasing order

     int r = array.length - 1;
     int s = j + 1;

     while (r > s) {
         temp = array[s];
         array[s++] = array[r];
         array[r--] = temp;
     }

     return array;
  } // getNext()

public int[][] getNextMultiple(int n) {
    int[][] result = new int[n][];
    int i = 0;
    int len = array.length;
    while(hasMore()){
        if(i>=n) break;
        result[i++] = Arrays.copyOf(getNext(), len);
    }
    if(i == n){ return result;}
    int[][] truncated = new int[i][];
    for(int j=0;j<i;j++){
        truncated[j] = result[j];
    }
    return truncated;
}

  // For testing of the PermutationGenerator class
  public static void main(String[] args) {
     PermutationGenerator pg = new PermutationGenerator(3, 1);

     while (pg.hasMore()) {
         int[] temp =  pg.getNext();
         for (int i = 0; i < temp.length; i++) {
             System.out.print(temp[i] + \" \");
         }
         System.out.println();
     }
  }

}";

Now, we can try things. Create the permutation generator:

gen = JavaNew["PermutationGenerator", 14, 1]

Generate first 5 permutations:

gen@getNextMultiple[10]

(*
  {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14}, 
  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 13}, 
  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 12, 14}, 
  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 13, 14, 12}, 
  {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 14, 12, 13}}
*)

Let's get more serious. On my machine, I do generate 5 million permutations per second (including the data transfer time):

ctr = 0;
gen@reset[]
While[(nxt = gen@getNextMultiple[10000]) =!= {} && ctr < 5*10^6, 
  ctr += Length[nxt]
]; // AbsoluteTiming
ctr

(* {1.040131, Null} *)

(* 5000000 *)

This means that you'll need a few hours to process them all - I only count the time to generate them:

14!/(5*10^6*3600) // N

(* 4.84324 *)

In fact, following halirutan's argument, if we just increment a counter on the top-level, it will take longer than generating permutations:

testctr = 0;
Do[testctr++, {5*10^6}] // AbsoluteTiming

(* {2.361853, Null} *)

But because we generate permutations in large batches, there is some chance that you could write a compiled (and possibly parallelized) function to process them also in batches, and that one may be much faster.

Still, hours, not years ). But of course, I agree with halirutan, that if your function which you plan to apply to permutations is any intensive computationally, this won't help you. So, the goal of my post was to show how to reduce the overhead of permutation generation, only. Whether or not this will be useful will entirely depend on what sort of function you will apply to them.

$\endgroup$
6
$\begingroup$

Well from computational point of view, if you wanted the whole list or some part of it then the size of output would be the main problem.

Assuming plaintext output is used ... by my rough estimation it would take over 225 GB (gigabytes) to store the whole list on a disk. Furthermore it would take about 4 days to compute them all on this laptop.

You need data like:

filename = "f://out.txt";
seed = Range[12];
size = 10;

Recursive method call like:

seed = MPermutations[seed, size, filename]

Other useful lines:

FilePrint[filename]
DeleteFile[filename]

This is an example of what MPermutations could look like if flat file output is used. It computes next count number of permutations of list, prints them in filename and returns the last element.

MPermutations[list_, count_Integer, filename_String] := Module[
  {n = 1, current = list},
  OpenAppend[filename];
  While[
   n < count, current = NextPermutation[current];
   n++;
   Write[filename, current]];
  Close[filename];
  current]

Return result is needed to make recursive method call possible. Recomputing that line will add next size amount of results to file (assuming that data is in a different cell).

$\endgroup$
6
$\begingroup$

A lot of times, we don't necessarily need all the permutations, the real purpose is to select a part of interest according to condition. I wrote a function selectPermutations, it's quite efficient and take up very little memory, sometimes would be useful.

ClearAll[selectPermutations];
Options[selectPermutations]={CompilationTarget->"WVM"};
selectPermutations[nn_Integer,cond_,A_List,OptionsPattern[]]:=
Module[{F,var,mm,iter,cTarget,cTime,cf,ans},
  F=Symbol@*FromLetterNumber;
  var=Variables[A];
  mm=Length[var];
  iter=Table[{F[ii+1],If[Select[cond&&Unequal@@var//Not//LogicalExpand,
    F[ii]==Last@Sort@Cases[#,_Symbol,-1]&],0,Evaluate@If[ii<mm,nn,1]]},{ii,mm}];
  cTarget=OptionValue[CompilationTarget];
  cTime=First@AbsoluteTiming[
  cf=Compile[{{a,_Integer}},
    Module[{B=Internal`Bag[Rest@{0}]},
      Do[Internal`StuffBag[B,#,1],##2];Internal`BagPart[B,All]],
        RuntimeAttributes->{Listable},RuntimeOptions->"Speed",CompilationTarget->cTarget
  ]&[A,Sequence@@iter]];
  If[cTarget=="C",Print["Compilation time: ",cTime]];
  Print["Executing time: ",First@AbsoluteTiming[
  ans=ArrayReshape[#,{Length[Flatten@#]/#2,#2}]&[cf[Range[nn]],Length[A]]]];
  ans
]

Performance testing

I am using Mathematica 12 on a laptop running Windows 10 (64-bit), CPU: Intel i7-9750H, RAM: 32GB.

Case 1

Want to compute the permutations of {1, 2, ..., 11} with only 3 GB of memory

cond = And @@ (Abs[#-#2]!=1& @@@ {{a,b},{a,c},{b,c},{b,f},{b,e},{c,e},{c,f},{c,g},{d,f},
  {d,g},{e,b},{e,a},{e,f},{e,i},{e,h},{f,g},{f,j},{f,i},{f,h},{g,i},{g,j},{h,i},{h,k},{i,j},
  {i,k},{j,k}});

selectPermutations[11,cond,{a,b,c,d,e,f,g,h,i,j,k}]//Length

Output:

Executing time: 0.0915339
23632

Using Select and Permutations

(*cond=Replace[cond, x_Symbol :> Indexed[#, LetterNumber@ToString@x], -1]*)

cf = With[{cond=cond},
Compile[{},
  Select[Permutations[Range[11]],
    Block[{a=#[[1]],b=#[[2]],c=#[[3]],d=#[[4]],e=#[[5]],f=#[[6]],g=#[[7]],h=#[[8]],i=#[[9]],
      j=#[[10]],k=#[[11]]}, cond]&],
       CompilationTarget->"C",RuntimeOptions->"Speed"
  ]
];
cf[]//Length//AbsoluteTiming

Output:

{4.29387, 23632}


Case 2

https://brg.a2hosted.com/?page_id=5489
Just a casual example. Actually, this program can be optimized, you can solve a diophantine equation first, and then traverse the 7 variables, this is not the point.

cond = (a+c+i==g+j+l&&a+h+l==g+i+k&&b+i+j==c+f+g&&b+d+i==e+g+h&&d+f+l==e+g+k)&&
  (a+c+i<a+h+l<b+i+j<b+d+i<d+f+l);

selectPermutations[12,cond,{a,b,c,d,e,f,g,h,i,j,k,l},CompilationTarget->"C"]

Output:

Compilation time: 0.36566
Executing time: 0.340409
{{4, 8, 5, 9, 12, 11, 1, 6, 2, 7, 10, 3}}

Using Select and Permutations

cf = With[{cond=cond},
Compile[{},
  Select[Permutations[Range[12]],
    Block[{a=#[[1]],b=#[[2]],c=#[[3]],d=#[[4]],e=#[[5]],f=#[[6]],g=#[[7]],h=#[[8]],i=#[[9]],
      j=#[[10]],k=#[[11]],l=#[[12]]},cond]&],
        CompilationTarget->"C",RuntimeOptions->"Speed"
  ]
];
cf[]//AbsoluteTiming

Ouput:

{412.059, {{4, 8, 5, 9, 12, 11, 1, 6, 2, 7, 10, 3}}}

1211 times faster.


Case 3

Magic hexagon
The task is to place the integers 1,...,19 in the following grid so that the sum of all numbers in a straight line (there are lines of length 3, 4 and 5) is equal.
enter image description here
Easy to prove, the magic constant is 38. First we can solve following equations, then traverse the 7 variables.

Solve[{x1+x2+x3==x4+x5+x6+x7==x8+x9+x10+x11+x12==x13+x14+x15+x16==x17+x18+x19==x1+x4+x8==
  x2+x5+x9+x13==x3+x6+x10+x14+x17==x7+x11+x15+x18==x12+x16+x19==x3+x7+x12==x2+x6+x11+x16==
  x1+x5+x10+x15+x19==x4+x9+x14+x18==x8+x13+x17==38},Array[Symbol["x"<>ToString@#] &, 19],
    Integers,GeneratedParameters->Symbol@*FromLetterNumber][[1]]//Normal//Column

Output

x1->a
x2->b
x3->38-a-b
x4->c
x5->d
x6->e
x7->38-c-d-e
x8->38-a-c
x9->f
x10->g
x11->38-b-d-e-f-g
x12->-38+a+b+c+d+e
x13->38-b-d-f
x14->38-c-d-e-f-g
x15->-38+b+c+d+e+f
x16->d+f+g
x17->-38+a+b+c+d+f
x18->d+e+g
x19->76-a-b-c-2 d-e-f-g

A={a,b,38-a-b,c,d,e,38-c-d-e,38-a-c,f,g,38-b-d-e-f-g,-38+a+b+c+d+e,38-b-d-f,
   38-c-d-e-f-g,-38+b+c+d+e+f,d+f+g,-38+a+b+c+d+f,d+e+g,76-a-b-c-2 d-e-f-g};
cond=LogicalExpand[Unequal@@A]&&And@@Thread[1<=Complement[A,Variables[A]]<=19];
selectPermutations[19,cond,A,CompilationTarget->"C"]

Output:

Compilation time: 0.656829
Executing time: 0.0106396
{{3, 17, 18, 19, 7, 1, 11, 16, 2, 5, 6, 9, 12, 4, 8, 14, 10,13, 15},
{3, 19, 16, 17, 7, 2, 12, 18, 1, 5, 4, 10, 11, 6, 8, 13, 9, 14, 15},
{9, 11, 18, 14, 6, 1, 17, 15, 8, 5, 7, 3, 13, 4, 2, 19, 10, 12, 16},
{9, 14, 15, 11, 6, 8, 13, 18, 1, 5, 4, 10, 17, 7, 2, 12, 3, 19, 16},
{10, 12, 16, 13, 4, 2, 19, 15, 8, 5, 7, 3, 14, 6, 1, 17, 9, 11, 18},
{10, 13, 15, 12, 4, 8, 14, 16, 2, 5, 6, 9, 19, 7, 1, 11, 3, 17, 18},
{15, 13, 10, 14, 8, 4, 12, 9, 6, 5, 2, 16, 11, 1, 7, 19, 18, 17, 3},
{15, 14, 9, 13, 8, 6, 11, 10, 4, 5, 1, 18, 12, 2, 7, 17, 16, 19, 3},
{16, 12, 10, 19, 2, 4, 13, 3, 7, 5, 8, 15, 17, 1, 6, 14, 18, 11, 9},
{16, 19, 3, 12, 2, 7, 17, 10, 4, 5, 1, 18, 13, 8, 6, 11, 15, 14, 9},
{18, 11, 9, 17, 1, 6, 14, 3, 7, 5, 8, 15, 19, 2, 4, 13, 16, 12, 10},
{18, 17, 3, 11, 1, 7, 19, 9, 6, 5, 2, 16, 14, 8, 4, 12, 15, 13, 10}}

Using Select and Permutations

cf = With[{cond=cond},
Compile[{{v,_Integer,2}},
  Select[v,Block[{a=#[[1]],b=#[[2]],c=#[[3]],d=#[[4]],e=#[[5]],f=#[[6]],g=#[[7]]},cond]&],
  CompilationTarget->"C",RuntimeOptions->"Speed"
  ]
];
Function[{a,b,c,d,e,f,g},Evaluate@A] @@@ cf[Permutations[Range[19],{7}]]//AbsoluteTiming

Output:

{27.3219, {{3, 17, 18, 19, 7, 1, 11, 16, 2, 5, 6, 9, 12, 4, 8, 14, 10,13, 15},
{3, 19, 16, 17, 7, 2, 12, 18, 1, 5, 4, 10, 11, 6, 8, 13, 9, 14, 15},
{9, 11, 18, 14, 6, 1, 17, 15, 8, 5, 7, 3, 13, 4, 2, 19, 10, 12, 16},
{9, 14, 15, 11, 6, 8, 13, 18, 1, 5, 4, 10, 17, 7, 2, 12, 3, 19, 16},
{10, 12, 16, 13, 4, 2, 19, 15, 8, 5, 7, 3, 14, 6, 1, 17, 9, 11, 18},
{10, 13, 15, 12, 4, 8, 14, 16, 2, 5, 6, 9, 19, 7, 1, 11, 3, 17, 18},
{15, 13, 10, 14, 8, 4, 12, 9, 6, 5, 2, 16, 11, 1, 7, 19, 18, 17, 3},
{15, 14, 9, 13, 8, 6, 11, 10, 4, 5, 1, 18, 12, 2, 7, 17, 16, 19, 3},
{16, 12, 10, 19, 2, 4, 13, 3, 7, 5, 8, 15, 17, 1, 6, 14, 18, 11, 9},
{16, 19, 3, 12, 2, 7, 17, 10, 4, 5, 1, 18, 13, 8, 6, 11, 15, 14, 9},
{18, 11, 9, 17, 1, 6, 14, 3, 7, 5, 8, 15, 19, 2, 4, 13, 16, 12, 10},
{18, 17, 3, 11, 1, 7, 19, 9, 6, 5, 2, 16, 14, 8, 4, 12, 15, 13, 10}}}

2570 times faster.

Related link: Using Meta programming + Compile to speed up Backtracking algorithm

$\endgroup$
1
  • $\begingroup$ The link in Case 2 is dead; what's the original problem in that link? $\endgroup$
    – user688486
    Commented Dec 26, 2023 at 5:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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