Services & Resources / Wolfram Forums / MathGroup Archive

MathGroup Archive 2004

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Re: Re: Re: optimally picking one element from each list

I had a hard-wired debug print in that version; here's the code without it.

Clear[optimalSplit]
Options[optimalSplit] = {Verbosity -> Brief};
optimalSplit::usage =
   "optimalSplit[listOfSets]
       finds an optimal choice from each set,\n minimizing \
cost=Length@Split@choices and returning {cost,choices}.\nOptions:\n \
Verbosity->Brief (default),\n Verbosity->Tableau gives an L.P. tableau,\n \
Verbosity->Equations gives the constraints,\n Verbosity->Both gives both."
optimalSplit[test_List, opts___Rule] := Block[
     {verbosity, elements, intersections, x, y, z, nodes, arcs, n, a, rules, \
m, b, c},
     verbosity = Verbosity /. Join[{opts}, Options[optimalSplit]];
     elements = Flatten[MapIndexed[Thread[x[First@#2, #1]] &, test], 1];
     intersections = z @@@ Flatten[MapIndexed[Thread[{First@#2, #1}] \
&, Intersection @@@ Partition[test, 2, 1]], 1];
     nodes = Join[{"
     obj"}, x /@ Range@Length@test, elements, {dummy}]; Print[nodes];
     arcs = Join[elements, y @@@ elements, intersections, {"RHS"}];
     h[i_, j_] := {nodes[[i]], arcs[[j]]};
     rules = {{i_, i_} -> 1, {x[i_], x[i_,
     j_]} -> -1, {
         i_x, i_} -> 1, {x[i__], z[
           i__]} -> -1, {x[i_, j_], z[ii_, j_]} /; ii == i - 1 ->
         1, {x[1], "RHS"} -> -1, {"obj", y[i_, _]} /; i < Length[
             test] -> 1, {x[i__], y[i__]} -> -1, {x[i_], y[ii_, j_]} /; ii ==
           i - 1 -> 1, {
         dummy, y[Length@test, _]} -> 1, {dummy, "RHS"} -> 1, {_, _} -> 0};
     {n, a} = Length /@ {nodes, arcs};
     m = Array[{nodes[[#1]], arcs[[#2]]} &, {n, a}] /. rules;
     If[MemberQ[{Tableau, Both}, verbosity],
            Print@TableForm[Transpose[m], TableHeadings -> Reverse@ {nodes,
           arcs}]];
     b = m[[Range[2, n], -1]];
     c = m[[1, Range[a - 1]]];
     m = m[[Range[2, n], Range[a - 1]]];
     If[MemberQ[{Equations, Both}, verbosity], Print[
     Thread[m.Most@arcs == b[[All]]] // ColumnForm];
       Print["objective = ", b.Rest[nodes]]
       ];
     b = Thread[{b, 0}];
     soln = Transpose@{Most@arcs, LinearProgramming[c, m, b]};
     {Length@
       Split@#, #} &[Sort[soln /. {{_, 0} -> Sequence[], {y[__],
           1} -> Sequence[], {x[i_, j_], 1} -> {
             i, j}, {z[i_, j_], 1} -> {i + 1, j}}][[All, -1]]]
     ]

Bobby

On Wed, 26 May 2004 14:41:19 -0500, DrBob <drbob at bigfoot.com> wrote:

> The "weights" you mention don't appear in my solution, I think, and the right-hand-sides from mine don't seem to appear in yours... unless I'm misunderstanding you, and the two are actually the same. Your weights seem to be 1 when my RHS is 0, and vice-versa, which seems odd. Calling them weights makes me think perhaps you're setting up the dual problem, which would be perfectly legitimate. Perhaps that's the difference.
>
> Let me describe my solution in words, and maybe you'll be able to set me straight. (At this point, two monitors are very useful!)
>
> 1) optimalSplit's Tableau display shows a Transpose of the actual tableau solved. Here's another version of the code (cleaned up just a little from last time). Run it on a small example, and I'll walk you through the Tableau.
>
> Clear[optimalSplit]
> Options[optimalSplit] = {Verbosity -> Brief};
> optimalSplit::usage =
>    "optimalSplit[listOfSets]
>        finds an optimal choice from each set,\n minimizing \
> cost=Length@Split@choices and returning {cost,choices}.\nOptions:\n \
> Verbosity->Brief (default),\n Verbosity->Tableau gives an L.P. tableau,\n \
> Verbosity->Equations gives the constraints,\n Verbosity->Both gives both."
> optimalSplit[test_List, opts___Rule] := Block[
>      {verbosity, elements, intersections, x, y, z, nodes, arcs, n, a, rules, \
> m, b, c},
>      verbosity = Verbosity /. Join[{opts}, Options[optimalSplit]];
>      elements = Flatten[MapIndexed[Thread[x[First@#2, #1]] &, test], 1];
>      intersections = z @@@ Flatten[MapIndexed[Thread[{First@#2, #1}] \
> &, Intersection @@@ Partition[test, 2, 1]], 1];
>      nodes = Join[{"
>      obj"}, x /@ Range@Length@test, elements, {dummy}]; Print[nodes];
>      arcs = Join[elements, y @@@ elements, intersections, {"RHS"}];
>      h[i_, j_] := {nodes[[i]], arcs[[j]]};
>      rules = {{i_, i_} -> 1, {x[i_], x[i_,
>      j_]} -> -1, {
>          i_x, i_} -> 1, {x[i__], z[
>            i__]} -> -1, {x[i_, j_], z[ii_, j_]} /; ii == i - 1 ->
>          1, {x[1], "RHS"} -> -1, {"obj", y[i_, _]} /; i < Length[
>              test] -> 1, {x[i__], y[i__]} -> -1, {x[i_], y[ii_, j_]} /; ii ==
>            i - 1 -> 1, {
>          dummy, y[Length@test, _]} -> 1, {dummy, "RHS"} -> 1, {_, _} -> 0};
>      {n, a} = Length /@ {nodes, arcs};
>      m = Array[{nodes[[#1]], arcs[[#2]]} &, {n, a}] /. rules;
>      If[MemberQ[{Tableau, Both}, verbosity],
>             Print@TableForm[Transpose[m], TableHeadings -> Reverse@ {nodes,
>            arcs}]];
>      Print@Transpose@m;
>      b = m[[Range[2, n], -1]];
>      c = m[[1, Range[a - 1]]];
>      m = m[[Range[2, n], Range[a - 1]]];
>      If[MemberQ[{Equations,
>      Both}, verbosity], Print[Thread[m.Most@arcs == b[[All]]] // ColumnForm];
>        Print["objective = ", b.Rest[nodes]]
>        ];
>      b = Thread[{b, 0}];
>      soln = Transpose@{Most@arcs, LinearProgramming[c, m, b]};
>      {Length@Split@#, #} &[Sort[soln /. {{_, 0} ->
>             Sequence[], {y[__], 1} ->
>              Sequence[], {x[i_, j_], 1} -> {i, j}, {z[i_, j_], 1} -> {i + 1,
>      j}}][[All, -1]]]
>      ]
>
> test = Array[Union@Array[Random[Integer, 10] &, 4] &, 3]
> optimalSplit[test, Verbosity -> Both]
>
> 2) I'll speak in terms of the actual Tableau, even though you're looking at its Transpose, because that way we stay closer to the theory. It has a row for each node (constraint) of the network, a column for each arc (variable), a column for the supply/demand (RHS) at each node, and a row for the cost (objective) function. Each variable/arc/column has -1 in the row corresponding to the node the arc comes from, +1 in the row corresponding to the node it goes to, and (in the "obj" row) a unit-cost (distance) for using the arc.
>
> 3) The node (constraint) names are straightforward; there's a node x[i] for each set, and a node x[i,j] for each element j in set i. The RHS at each node is the sum of flow on arcs entering the node minus arcs leaving the node. The "dummy" (demand) node is where flow ends up in the end; it has RHS = +1, because one unit of flow will end up there. The first set has RHS = -1; this is where flow starts (supply).
>
> 4) I reused x[i,j] as an arc name, to mean the arc from x[i] to node x[i,j]. These arcs have zero cost, since a choice is mandatory.
>
> 5) Each node x[i,j] (i<n) is the start of an arc y[i,j], leading to node x[i+1] (the next set). These arcs have cost = 1, because they represent new terms in the Split of a solution.
>
> 6) Whenever j is in the intersection of set i and set i+1, I've added an arc z[i,j] from node x[i,j] to node x[i+1,j]. These arcs represent selecting node j from both sets, so the cost is 0.
>
> 7) A different network problem might have a second "dummy" node at the beginning, but in this case, the first set serves the purpose. The "dummy" at the end was necessary to collect flow from all the x[n,i].
>
> 8) Flow starts at the supply node (the first set) and finds a path to the demand node (dummy). At set 1, a choice of arcs x[1,j] is required, so flow goes to one of the elements. From there, it either goes over y[1,j] to the next set (causing a Split) or goes over z[1,j] to element x[2,j], avoiding the choice at set 2. Things go on in that way; some of the sets aren't visited by any flow.
>
> 9) In the optimal solution, for each i, a unique x[i,j] or z[i,j] is 1. In either case, that means element j was selected.
>
> I hope that's of use to someone!
>
> Bobby
>
> On Wed, 26 May 2004 08:54:21 -0500, Daniel Lichtblau <danl at wolfram.com> wrote:
>
>> DrBob wrote:
>>> My second solution may or may not be equivalent to what you're proposing
>>> -- you lost me in the details -- but the thought process behind it seems
>>> more straightforward, at least. I drew the problem as a network, named
>>> the nodes and arcs, and set up the LP tableau from that, using rules
>>> that paid attention to the names I'd chosen. Because it comes from a
>>> shortest-path network problem, it's guaranteed to be unimodular, hence
>>> basic solutions are integer. It doesn't begin as a quadratic problem, so
>>> things are a little less complicated. In the optimalSplit routine, I
>>> left all that transparent, rather than optimizing the setup process,
>>> which might obscure the thought process.
>>>
>>> I'm thinking of writing a threshold shortest-path code like the one I
>>> wrote in Fortran for my dissertation, years ago. That would simplify the
>>> setup, as I'd eliminate the tableau and input the arcs directly. It
>>> might speed up the solution, as well, although LinearProgramming is
>>> certainly fast already.
>>>
>>> Bobby
>>
>>
>> I know but little about network optimization but I see a way to view my
>> setup that may show it to be the same as yours.
>>
>> For each subset we form a column, with a "vertex" in the column for each
>> element in the subset (so each vertex corresponds to a value from the
>> universal set, and vertices in different columns can share such values).
>>
>> We form directed weighted arcs between all vertices of adjacent columns,
>> such that:
>>
>> (1) All weights are nonnegative.
>> (2) The sum of weights out of each column (except the last) is 1.
>> (3) The sum of weights into each column (except the first) is 1.
>> (4) For each vertex not in the first or last column, the sum of weights
>> of incoming arcs equals the sum of weights of outgoing arcs.
>>
>> This corresponds to the constraints I imposed.
>>
>> The function to minimize is the sum of arc weights between vertices that
>> do not correspond to the same value. This was the gist of my
>> KroneckerDelta usage.
>>
>> Anyway, this is a simpler way to set up the LP problem I formulated.
>> It's nicer than the roundabout way I first did it, but, as I said, I
>> don't tend to think in terms of network problems. If, as I suspect, this
>> reformulation corresponds to yours then the two are in fact the same
>> approach.
>>
>> Regardless, my feeling is that for largish problems the setup speed will
>> be small compared to the solving step. Hence it is best to have a
>> formulation that is most natural from the user's point of view, with
>> translation to an LP done, if need be, by the program. That at least is
>> my quick opinion on solver design; feel free to ignore it.
>>
>>
>> Daniel
>>
>>
>
>
>



-- 
Using M2, Opera's revolutionary e-mail client: http://www.opera.com/m2/
  • Prev by Date: Re: Uniform design
  • Next by Date: 3D fitting of data points
  • Previous by thread: Re: Re: Re: Re: optimally picking one element from each list
  • Next by thread: Re: Re: Re: optimally picking one element from each list