# Find permutations with constraints without using Permutations

EDIT: To clarify, the bottleneck right now is available RAM, so any answer should keep that in mind (I cannot store all T! lists of length T and filter out those that satisfy the condition a posteriori.)

I want to find all permutations of the elements of Range[0,T-1] that satisfy a condition, but where T may be too large for Permutations to be useable: first generating and storing all permutations simply consumes too much RAM. The condition is always such that cond = {c,c,...,c[T]} means that the first element of the permutation must be larger than or equal to c, the second element must be larger than or equal to c etc. The condition is sorted in increasing order, and we can assume that the condition is not so strict that no permutations survive.

I have managed to implement what I want, but in a very procedural way using a recursive function (the details here are not that important):

recuPerm[level_] :=
If[level == 0,
res[] = Total[avail];
Sow[res],
((res[[level + 1]] = #; avail[[First[#] + 1]] = 0;
recuPerm[level - 1];
avail[[
First[#] + 1]] = #) & /@ ({(allow[[level + 1]].avail)} /.
Plus -> Sequence));
]


and I call it from the wrapper function:

listPerm[T_, cond_] :=
Block[
{a, avail, allow, res = ConstantArray[1, T], rip},
avail = a /@ Range[0, T - 1];
allow =
Table[PadLeft[ConstantArray[1, T - cond[[i]]], T], {i,
T}];
rip = Reap[recuPerm[T - 1]][];
If[rip == {}, {}, rip[]]
]


(The dummy head a is simply there so I can use Total and Dot in order to pick out allowed elements.)

Do you know of an approach that is more functional in nature and/or can better take advantage of the strengths of Mathematica? If it's more memory efficient (or faster) than my (unelegant) attempt then that's of course a bonus!

• Could you give us an order of magnitude for T ? – A.G. Feb 24 '15 at 22:17
• Well, the number of permutations grows so fast with T right, so I have only dared to try T about 11 or 12. But I am also mainly interested in constrains that are such that, if I subtract the constraint from an allowed permutation, I would get something whose Total has the same order of magnitude as T. Take for instance T = 8 and cond = {0,1,2,2,3,3,4,5}. Contrast this to no constraints, cond = ConstantArray[0,T], then any perm has Total[perm - cond] = T(T-1)/2. – Marius Ladegård Meyer Feb 25 '15 at 6:44
• Thank you djp and rasher for some nice answers! I started from rashers solution and made a compiled versjon which was naturally even a bit faster. – Marius Ladegård Meyer Mar 10 '15 at 8:23

This is pretty functional:

f = Module[{comps, r = Reverse@Range[#2, #1 - 1]},
comps[l1_, l2_] := Join @@ Map[Thread[{Sequence @@ #, Complement[l2, #]}] &, l1];
Reverse /@ Fold[comps, Transpose@{First@r}, Rest@r]] &;


This is about 10-15% faster, very slightly higher memory use (but still far below your current solution):

fz = With[{r = Reverse@Range[#2, #1 - 1]},
{#1, Outer[Complement, {#2}, #1, 1][]}]) &,
Transpose@{r[]}, Rest@r][[All, -1 ;; 1 ;; -1]]] &;


Comparing and including djp's interesting solution:

t = 11;
c = cond = {0, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6};

lp = listPerm[t, c]; // Timing // First

Timing[
possibleElements = Range[cond, Length@cond - 1];
fn[{{last_}}] := {{last}};
fn[{most__List, last_List}] := Table[{fn[{most} /. i -> Sequence[]], i}, {i, last}];
intermediate = fn[possibleElements];
result = intermediate //. {x : {{__Integer} ..}, i_Integer} :> Sequence @@ (Append[#, i] & /@ x);
] // First

fr = f[t, c]; // Timing // First

fzr = fz[t, c]; // Timing // First

(lp /. a[x_] :> x) == result == fr == fzr

(*

8.704856

6.661243

2.839218

2.464816

True

*)


Timings on an old netbook, but ~3X faster. Memory utilization s/b close to optimal: it never grows the intermediate results list to more than the ultimate results list length. Fails gracefully - if restrictions have no results, it returns no permutations (your current code, I'm sure you're aware, goes bonkers ;-} ). This s/b easy to adapt to a staggered restriction range, that is, differing lower and upper bounds for each position, should you so desire.

Neat puzzle, BTW...

Side note: You'll beat these by doing things iteratively... bodging up a function generator that builds such a solution based on parameters was 60% faster than my own fastest on some quick tests...

• This would probably deserve a question in itself, but how would you go about commenting/documenting the code so as to make it better understandable? – A.G. Mar 3 '15 at 16:12
• @a.g.. Something like "(* Given cap and list of restrictions, produces all permutations with restructions *)". In all seriouslness, though, were it my code for my use, I'd do nothing - it's self-evident. If I knew others would be using/reading/trying to understand it, something between a version of the terse comment and a step-by-step explanation. Terse Mathematica is nothing like terse APL... – ciao Mar 3 '15 at 20:45
• @A.G. please ask it. I have struggled with this dozens of times --- my collaborators probably have sought counselling over some of my more opaque code. – djp Mar 4 '15 at 2:09

This is a question of finding a suitable algorithm rather than use of Mathematica. The challenge is to generate only the permutations that will be used, rather than generate all permutations and filter those that satisfy a criterion. Fun problem. Here's my solution:

cond = {0, 0, 1};
possibleElements = Range[cond, Length@cond-1];
(* {{0, 1, 2}, {0, 1, 2}, {1, 2}} *)


So we need to generate all series {a,b,c,...} where each element appears precisely once.

I wrote a recursive function fn to take the last element off the end of the list, and remove it from the preceding elements. That is:

fn[{a,b,c}, {a,b,c}, {b,c}] ->
{
{fn[{a, c}, {a, c}], b},
{fn[{a, b}, {a, b}], c}
}

fn[{most__List, last_List}] :=
Table[{fn[{most} /. i -> Sequence[]], i}, {i, last}]


We also want to terminate the recursion:

fn[{{a}}] -> {{a}}
fn[{{last_}}] := {{last}};


The output of fn looks like this:

intermediate = fn[possibleElements]

(* {{{{{{2}}, 0}, {{{0}}, 2}}, 1}, {{{{{1}}, 0}, {{{0}}, 1}}, 2}} *)


To reduce this, we need to make the transformation, and do so repeatedly.

{ x:{Lists of Integers}, i_Integer } -> {{x1, i}, {x2, i}, ...}
result = intermediate //. {x : {{__Integer} ..}, i_Integer} :>
Sequence @@ (Append[#, i] & /@ x)


All together:

cond = {0, 1, 2, 3, 4, 5, 6, 6, 8, 8, 10, 11, 12, 13, 14};
possibleElements = Range[cond, Length@cond-1]
fn[{{last_}}] := {{last}};
fn[{most__List, last_List}] :=
Table[{fn[{most} /. i -> Sequence[]], i}, {i, last}]

intermediate = fn[possibleElements]
result = intermediate //. {x : {{__Integer} ..}, i_Integer} :>
Sequence @@ (Append[#, i] & /@ x)


It isn't a function and doesn't fail gracefully (thanks @rasher), but you can add those if you like. I think @rasher's algorithm is the same as mine.

• Good use of MMA patterns, but as written returns trash if conditions result in no valid permutations... – ciao Mar 3 '15 at 8:42
• @rasher Thanks. I think your algorithm is much the same? – djp Mar 3 '15 at 9:45
• Yes, there's no magic bullet (a la factorial numbering systems for "normal" perms), so any algo. is going to have a similar tack. How it's done is where memory/performance bennies might be found. – ciao Mar 3 '15 at 20:39