Re: Special permutation pattern ascending groups sum

*To*: mathgroup at smc.vnet.net*Subject*: [mg113265] Re: Special permutation pattern ascending groups sum*From*: Daniel Lichtblau <danl at wolfram.com>*Date*: Thu, 21 Oct 2010 07:02:16 -0400 (EDT)

Daniel Lichtblau wrote: > ----- Original Message ----- >> From: "me you" <coconut_dj at yahoo.com> >> To: mathgroup at smc.vnet.net >> Sent: Sunday, October 17, 2010 5:06:25 AM >> Subject: [mg113187] Special permutation pattern ascending groups sum of subset >> hi, >> >> i'm very new to Mathematica and i would like to ask for some help, i >> need some code which can handle with the following problem >> >> i have a vector of 9 elements {0, 0, 1, 2, 3, 4, 5, 7, 8} >> i have k = 3, k is the number of subsets >> each subset should contain elements in ascending order >> i need to get all the permutation that fit the following conditions >> a1 <= a2 <= a3 >> a2 <= a5 <= a6 >> a7 <= a8 <= a9 >> a1 <= a4 <= a7 >> a1 + a2 + a3 = 10 >> a4 + a5 + a6 = 10 >> a7 + a8 + a9 = 10 >> >> possible solution would be >> {0, 2, 8, 0, 3, 7, 1, 4, 5} >> and >> {0, 3, 7, 0,2, 8, 1, 4, 5} >> >> the script should be very dynamic because i would like to use it on a >> largeset of numbers (40) >> somehow before we can check for the sum of elements we need to >> generate just those permutation where elements in subset sare >> ascending order (a1<=a2<=a3) (a4<=a6<=a6) (a7<=a8<=a9) >> >> Thanks > > If I follow correctly, your aj's are meant as for a permutation where set element j goes to aj. If so, I think you want a4 rather than a2 for the middle inequality constraint. > > But it can be handled without finding explicit permutations that require ordering constraints. Can instead treat as a 0-1 assignment problem. Set up a 9x3 array of variables, where a[j,k] is one if element j is in subset k. We can then set up various constraints and solve as an integer linear constraint satisfaction problem. > > In[174]:= set = {0, 0, 1, 2, 3, 4, 5, 7, 8}; > n = Length[set]; > k = 3; > ssize = n/3; > target = 10; > > In[184]:= vars = Array[a, {k, n}]; > fvars = Flatten[vars]; > > Now we require that each variable be 0 or 1 (using inequalities and a domain restriction later). Also need rows to sum to the number of items in each subset (3, in this example). And we need that certain linear equalities be met (subsets sum to 10). > > Finally, in order to cut down on essentially equivalent solutions, we force the last element into the first subset, and also decree that the second subset receive either the next to last or third to last unless the last subset receives neither. This tends to work for the example in question (we get a unique solution). Maybe there is a better way of forcing solutions to be essentially unique, I'm not sure. (By essentially unique, I mean we do not get solutions that are simply subset reorderings of one another. If the meaning of this is not clear, try the setup but without c5 below.) > > In[201]:= c1 = Map[0 <= # <= 1 &, fvars]; > c2 = Map[Total[#] == ssize &, vars]; > c3 = Map[Total[#] == 1 &, Transpose[vars]]; > c4 = Map[set.# == target &, vars]; > c5 = {vars[[1, -1]] == 1, vars[[2, -2]] >= vars[[3, -2]], > Total[vars[[2, -3 ;; -2]]] >= Total[vars[[3, -3 ;; -2]]]}; > constraints = Flatten[Join[c1, c2, c3, c4, c5]]; > > In[207]:= Timing[ > solns = Map[Pick[set, #, 1] &, > vars /. {ToRules[Reduce[constraints, fvars, Integers]]}, {2}]] > > Out[207]= {0.016, {{{0, 2, 8}, {0, 3, 7}, {1, 4, 5}}, {{0, 2, 8}, {0, > 3, 7}, {1, 4, 5}}}} > > Daniel Lichtblau > Wolfram Research In private email I was asked whether this approach would handle somewhat larger examples. It gets complicated, so I thought I'd post a partial result. With an explicit branch-and-prune loop for ILP solving we can get individual solutions, Finding ALL solutions might be a challenge. It is easy to modify an instance solver for this, but if there are many solutions, and for some examples this seems to be the case, then it is not obvious to me whether or when it will terminate. There are a few reasons that an explicit loop (rather than invocation to FindInstance or Reduce) is required. One is that the built in functions use exact simplex algorithm to solve relaxed linear programming problems. We can instead use machine precision and hope we do not mess up. In practice, this is robust for most 0-1 ILPs I run across. It certainly seems to find solutions for this class of problem. Another reason is that these problems have a "smart" variable ordering, and FindInstance et al will not recognize that. Specifically, for relaxed problem solutions that are not integral, we want to branch on those that correspond to larger values in the set. This tends to tighten the candidate sums quickly. To give an idea of the effect, my first naive attempt that did not use this variable ordering approach handled a problem with Range[25] split into 5 sets in around three hours (just to find one solution). With smart ordering it came down to a fraction of a second. That was enough improvement to make me think, at first, that I'd introduced some bugs. Another issue is that we want to prevent essential duplicates of solutions, that is, separation into subsets that are permutations of one another. A method that works reasonably well is to try to force last elements of each subset to be in descending order. I do not know how to guarantee this, but if we simply force that the largest element be in subset#1, and that the second largest not be in subsets 3 and after, the third largest not be in subsets 4 and after, etc. This is what triangleConstraints below is enforcing. Below is the code. I direct LinearProgramming to use some fast library interior point code for solving the relaxed LPs. triangleConstraints[vars_, k_, n_] := Table[vars[[j, i]] == 0, {j, 2, k}, {i, n - j + 2, n}] SetOptions[LinearProgramming, Method -> "CLP"]; subsetSums[set_, k_, target_, nsols_: 1] := Module[{n = Length[set], a, ssize, c1, c2, c3, c4, c5, constraints, vars, fvars, jj = 0, stack, iter = 0, minval, sign, var, badvar, newcon, vals, solns = {}, soln, eps = 10^(-7)}, ssize = n/k; vars = Array[a, {k, n}]; fvars = Reverse[Flatten[Transpose[vars]]]; c1 = Map[0 <= # <= 1 &, fvars]; c2 = Map[Total[#] == ssize &, vars]; c3 = Map[Total[#] == 1 &, Transpose[vars]]; c4 = Map[set.# == target &, vars]; c5 = Join[{a[1, n] == 1}, triangleConstraints[vars, k, n]]; constraints = Flatten[Join[c1, c2, c3, c4, c5]]; stack = {constraints, {}}; While[stack =!= {}, iter++; If[Mod[iter, 200] == 0, Print[{iter}]]; constraints = stack[[1]]; stack = stack[[2]]; sign = RandomChoice[{-1, 1}]; var = sign*fvars[[Mod[iter, Length[fvars], 1]]]; Quiet[vals = NMinimize[{var, constraints}, fvars]]; If[Head[vals] == NMinimize || ! FreeQ[vals, Indeterminate], Continue[]]; minval = First[vals]; vals = Chop[vals[[2]], eps]; badvar = Position[fvars /. vals, a_ /; Chop[a - Round[a], eps] != 0, {1}, 1, Heads -> False]; If[badvar == {}, jj++; soln = vals; soln = soln /. (a_ -> b_) :> a -> Round[b]; solns = {solns, soln}; Print["new solution ", {jj, Map[Pick[set, #, 1] &, vars /. soln]}]; If[jj == nsols, Break[], Continue[]]; ]; badvar = badvar[[1, 1]]; badvar = fvars[[badvar]]; newcon = Sequence[]; If[sign == 1, If[Ceiling[minval - eps] == 1, newcon = var == 1]; ,(*else*) If[Ceiling[minval - eps] == 0, newcon = var == 0]; ]; stack = {{Join[constraints, {badvar == 0, newcon}]}, stack}; stack = {{Join[constraints, {badvar == 1, newcon}]}, stack}; ]; solns = Union[Partition[Flatten[solns], n*k]]; {iter, jj, Map[Pick[set, #, 1] &, vars /. solns, {2}]}] Example: We'll split Range[36] into 6 subsets of equal totals, and request five solutions. In[124]:= set = Range[36]; k = 6; target = Total[set]/k; In[127]:= Timing[subsetSums[set, k, target, 5]] During evaluation of In[127]:= new solution {1,{{13,14,15,16,17,36},{1,4,5,32,34,35},{8,9,10,22,29,33},{3,7,21,25,27,28},{11,12,19,20,23,26},{2,6,18,24,30,31}}} During evaluation of In[127]:= new solution {2,{{13,14,15,16,17,36},{2,3,5,32,34,35},{8,9,10,22,29,33},{4,6,21,25,27,28},{11,12,19,20,23,26},{1,7,18,24,30,31}}} During evaluation of In[127]:= new solution {3,{{13,14,15,16,17,36},{1,2,7,32,34,35},{8,9,10,22,29,33},{4,6,21,25,27,28},{11,12,19,20,23,26},{3,5,18,24,30,31}}} During evaluation of In[127]:= new solution {4,{{12,13,15,16,19,36},{2,3,5,32,34,35},{8,9,10,22,29,33},{4,6,21,25,27,28},{11,14,17,20,23,26},{1,7,18,24,30,31}}} During evaluation of In[127]:= new solution {5,{{11,14,15,16,19,36},{2,3,5,32,34,35},{8,9,10,22,29,33},{4,6,21,25,27,28},{12,13,17,20,23,26},{1,7,18,24,30,31}}} Out[127]= {1.86072, {63, 5, {{{11, 14, 15, 16, 19, 36}, {2, 3, 5, 32, 34, 35}, {8, 9, 10, 22, 29, 33}, {4, 6, 21, 25, 27, 28}, {12, 13, 17, 20, 23, 26}, {1, 7, 18, 24, 30, 31}}, {{12, 13, 15, 16, 19, 36}, {2, 3, 5, 32, 34, 35}, {8, 9, 10, 22, 29, 33}, {4, 6, 21, 25, 27, 28}, {11, 14, 17, 20, 23, 26}, {1, 7, 18, 24, 30, 31}}, {{13, 14, 15, 16, 17, 36}, {1, 2, 7, 32, 34, 35}, {8, 9, 10, 22, 29, 33}, {4, 6, 21, 25, 27, 28}, {11, 12, 19, 20, 23, 26}, {3, 5, 18, 24, 30, 31}}, {{13, 14, 15, 16, 17, 36}, {1, 4, 5, 32, 34, 35}, {8, 9, 10, 22, 29, 33}, {3, 7, 21, 25, 27, 28}, {11, 12, 19, 20, 23, 26}, {2, 6, 18, 24, 30, 31}}, {{13, 14, 15, 16, 17, 36}, {2, 3, 5, 32, 34, 35}, {8, 9, 10, 22, 29, 33}, {4, 6, 21, 25, 27, 28}, {11, 12, 19, 20, 23, 26}, {1, 7, 18, 24, 30, 31}}}}} How long might it take to find all solutions? I don't know. Hours, at least. But what will be the eventual appropriate unit of measure? Days? Months? Lifetimes of universe? Daniel Lichtblau Wolfram Research