I do not exactly understand why ArgMin
fails, but I suspect a difficulty in optimizing over the integers only.
An alternative is to use free numerical optimization over the reals (i.e. NArgMin
), then rounding this result both up and down to integer values, constructing all possible tuples of those integers, which is a much smaller search space, and selecting the tuple that best represents the ratios in the original population. This is what I do below. Note that I also rewrote the conditions in what seems to me a more natural way, and localized the variables inside a Module
:
ClearAll[dean2]
dean2[populations_, seats_] := Module[
{nonzeropops, vars, objective, constraintOnTotal, constraintOnSign, min,
viableoptions, best},
nonzeropops = DeleteCases[populations, 0];
vars = Array[x, Length[nonzeropops]];
objective = Sum[
Log[(x[i]!)^3 2^(x[i] - 1)/(x[i]^2 nonzeropops[[i]]^(x[i] - 1) ((2 x[i] - 1)!))],
{i, Length[nonzeropops]}];
constraintOnTotal = Total[vars] == seats;
constraintOnSign = And @@ Thread[vars >= 0];
min = NArgMin[{objective, constraintOnTotal && constraintOnSign}, vars];
viableoptions = Select[
Tuples@Transpose@Through[{Floor, Ceiling}[min]],
Total[#] == seats &];
best = First@
Nearest[(Ratios /@ viableoptions) -> viableoptions, Ratios[nonzeropops]];
Fold[Insert[#1, 0, #2] &, best, Position[populations, 0]]
]
This works fine even for very large seat numbers:
dean2[{0, 40, 20, 12, 0}, 21]
(* Out: {0, 11, 6, 4, 0} in 0.073 s *)
dean2[{0, 40, 20, 12, 0}, 210]
(* Out: {0, 117, 58, 35, 0} in 0.1s *)
dean2[{0, 40, 20, 12, 0}, 21000]
(* Out: {0, 11667, 5833, 3500, 0} in 0.14 s *)
Perhaps more directly, one can implement the spirit of Dean's apportionment method by minimizing the differences in the ratios of the numbers of representatives apportioned to each population, and the ratios of the populations themselves:
ClearAll[deanAlt]
deanAlt[populations_?(VectorQ[#, (# >= 0 &)] &), seats_?(# > 0 &)] :=
Module[{nonzeropops, vars, min, viableoptions, best},
nonzeropops = DeleteCases[populations, 0];
vars = Array[s, Length[nonzeropops]];
min = ArgMin[{Norm[First@*Ratios /@ Partition[vars, 2, 1] - Ratios[nonzeropops]], Total[vars] == seats}, vars];
viableoptions = Select[Tuples@Transpose@Through[{Floor, Ceiling}[min]], Total[#] == seats &];
best = First@Nearest[(Ratios /@ viableoptions) -> viableoptions, Ratios[nonzeropops]];
Fold[Insert[#1, 0, #2] &, best, Position[populations, 0]]
]
This reproduces the results from the original dean
, and from dean2
above, and it is much, much faster.