Optimizing the performance of an algorithm

Question

First - a bit of an introduction. If you're only interested in the code, you can skip this section.

The following question is drawn from Dennis E. Shasha's Puzzling Adventures, and is listed under the question Strategic Bullying.

We have a set of entities, each of them with a power represented by a number. A coalition's power will be the sum of its members' powers. Coalition A will destroy Coalition B if A's power is greater than that of B, but A will back down if B can obtain enough new allies such that B's power is greater than or equal to A. When a coalition is beaten, it is destroyed and its wealth (proportional to its power) is distributed among the winners. The winners do not lose any of their power. Each of the entities plays selfishly, striving above all to avoid certain destruction but also to acquire wealth.

Based on the description of the question above, a set of entities is stable if there is a strong enough coalition to prevent any war and every member of that coalition has an interest in participating in the coalition.

---

Mathematically speaking, a set of entities with distinct integer powers is:

• Unstable if the most powerful entity has more than half the total power of all entities, since it would destroy all other entities without being challenged;
• Stable if the most powerful entity has exactly half of the total power of all entities, since it would not be able to destroy the other entities (they have an incentive to form a coalition to stop the most powerful entity;
• and otherwise, Stable if and only if ALL subset of entities formed by removing any coalition of entities with less than half the power of all entities is Unstable.

(Thanks to Xerxes for pointing out the error in my original formulation).

What is the most efficient code we can write to check the stability of a given set of entities?

Now, I've written a function to check the stability of any given set of entities:

fstable = Which[
   2*Max[#1] > Total[#1], 0, 
   2*Max[#1] == Total[#1], 1,
   True, BitXor[Max[fstable /@ Subsets[#1, {Length[#1] - 1}]], 1]
   ] &;

Edit:The above code doesn't cover all the cases in the question as it assumed that one only needed to consider the stability of any subset formed by removing one element - look at Xerxes' answer below for a correct and faster implementation.

Examples of implementation:

fstable[{3, 4, 5, 6}]
0

fstable[{3, 4, 5}]
1

However, this function does not run very fast, and in fact even checking the stability of a set of 10 entities takes some time:

Timing[fstable[Range[10]]]
{30.966199, 0}

Is there any way we can improve this, or will we require some additional mathematical insight to allow us to eliminate specific cases?

---

Note: The question the author actually asks is for us to find the largest set of entities with distinct integer powers and with the most powerful entity having a power of 21. To solve the problem by brute force, I would define the function using the idiom f[x_]:=f[x] to make use of Mathematica's ability to store function evaluations. However, the total number of function evaluations to be solved is quite large, so I'm not sure whether doing so will cause the kernel to run out of memory.

Answer 1

I think the key here is to abort evaluation as soon as we have enough information to judge the stability of the group. There's no need to evaluate the stability of all removed-member subsets if any one of them turns out to be stable. We implement the short-circuit using a Throw-Catch pair. We also want to avoid constructing the full Subsets if we can avoid it, so we use the third parameter inside a Do loop to construct them one by one. I also memoized it, since even for 21, this will only require a couple million entries.

fstable[{_}] := True
fstable[set_] := (fstable[set] = 
   Module[{m = Max[set], t = Total[set]}, 
    Catch[Which[m > t - m, False, 2 m == t, True, True, 
      Do[If[fstable[s], Throw[False]], {s, 
        Subsets[set, {Length[set] - 1}]}]; True]]])

This gives a few orders of magnitude speed-up.

Timing[oldfstable[Range[11]] // Positive]
(* {140.236764, True} *)
Timing[fstable[Range[11]]]
(* {0.024002, True} *)

For 21, we get

Timing[fstable[Range[21]]]
(* {1.664104, True} *)

Of course, as far as I can tell, fstable[Range[n]] is equivalent to OddQ[n], which is even faster, but this solution will work for more arbitrary cases.

EDIT: Update to address the concerns in my comment above (i.e. that we have failed to consider multi-element coalitions being knocked out simultaneously in a conflict).

This case is much more difficult than the other, so I'm employing a few more tricks. Firstly, we memoize onto a different head fs instead of fstable. This allows us to quickly scan new sets to see if they have been previously resolved without doing any hard work in case they haven't. This means that we can short-circuit very quickly if we already know a set from previous calculations. Secondly, we sort the knocked-out subsets according to how much power remains, and we make sure to only consider cases where the knocked-out coalition is actually weaker than the remaining one.

fstable[{_}] := True
fstable[set_] := 
 If[MatchQ[fs[set], True | False], 
  fs[set], (fs[set] = 
    Module[{m = Max[set], t = Total[set]}, 
     Catch[Which[m > t - m, False, m == t - m, True, True, 
       With[{ss = 
          SortBy[Select[Subsets[set, {2, Length[set] - 1}], 
            Total[#] > t/2 &], -Total[#] &]}, 
        Do[If[fs[s] === True, Throw[False]], {s, ss}]; 
        Do[If[fstable[s], Throw[False]], {s, ss}]; True]]]])]

Then, we can use this to try to address the given problem. What is the largest stable set of non-repeated integer elements less than and including 21?

Timing[Block[{max = 21}, best = {};
 Do[Block[{S = Subsets[Range[max - 1], {0, max - 1}, {si}][[1]]}, 
   If[Length[S] >= Length[best] && fstable[Append[S, max]], 
    best = Append[S, max]]], {si, 2^(max - 1)}]; best]]
(* 34096.250882, {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 14, 20, 21} *)

Yikes! That's a long time. We're at the very limits of what brute force can do. Probably there is a better technique that I'm not seeing.

SECOND EDIT: When thinking about why certain configurations are stable, it's sometimes unsatisfying to know that a configuration is unstable without knowing what it destabilizes into. We know that a configuration can be no more than one iteration from stable, since any entity that would be eliminated in subsequent steps would want to oppose the instability from the start. Thus, I replace False with the stable configuration each set devolves into. According to the rules, we choose the one in which each participant in the winning coalition receives the largest amount of loot. This means we can't use one of our short circuits, so it makes the code slower.

fstable[{_}] := True
fstable[set_] := 
 If[MatchQ[fs[set], True | False], 
  fs[set], (fs[set] = 
    Module[{m = Max[set], t = Total[set]}, 
     Catch[Which[m > t - m, {m}, m == t - m, True, True, 
       With[{ss = 
          SortBy[Select[Subsets[set, {2, Length[set] - 1}], 
            Total[#] > t/2 &], N[(Total[#] - t)/Length[#]] &]}, 
        Do[If[fstable[s], Throw[s]], {s, ss}]; True]]]])]

Here's an interesting case, which we wrongly treated in the comments:

fstable[{1, 2, 3, 4, 5}]
(* {1, 3, 4} *)

3 and 4 won't form a coalition with 5 to knock out puny {1,2} when they can ally with 1 to knock out the more valuable {2,5} group.