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