Skip to main content

programming - 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[1],c[2],...,c[T]} means that the first element of the permutation must be larger than or equal to c[1], the second element must be larger than or equal to c[2] 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[[1]] = 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]][[2]];
  If[rip == {}, {}, rip[[1]]]
  ]


(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!



Answer



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]}, 
          Fold[(Join @@ MapThread[Thread[{Sequence @@ #1, #2}] &, 
               {#1, Outer[Complement, {#2}, #1, 1][[1]]}]) &, 
               Transpose@{r[[1]]}, 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...


Comments

Post a Comment

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...
Post a Comment
Read more

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...
Post a Comment
Read more

How to remap graph properties?

Graph objects support both custom properties, which do not have special meanings, and standard properties, which may be used by some functions. When importing from formats such as GraphML, we usually get a result with custom properties. What is the simplest way to remap one property to another, e.g. to remap a custom property to a standard one so it can be used with various functions? Example: Let's get Zachary's karate club network with edge weights and vertex names from here: http://nexus.igraph.org/api/dataset_info?id=1&format=html g = Import[ "http://nexus.igraph.org/api/dataset?id=1&format=GraphML", {"ZIP", "karate.GraphML"}] I can remap "name" to VertexLabels and "weights" to EdgeWeight like this: sp[prop_][g_] := SetProperty[g, prop] g2 = g // sp[EdgeWeight -> (PropertyValue[{g, #}, "weight"] & /@ EdgeList[g])] // sp[VertexLabels -> (# -> PropertyValue[{g, #}, "name"]...
Post a Comment
Read more