A minimization which works but not too much

Question

Here under is a method to calculate Dean's apportionment method

dean[v_, s_] := Module[{vv = v, ss = s}, ww = DeleteCases[vv, 0];
  var = Table[x[i], {i, Length[ww]}];
  obj = Total[
    Table[Log[(x[i]!)^3 Rationalize[2^(
         x[i] - 1)]/(x[i]^2 ww[[i]]^(x[i] - 1) ((2 x[i] - 1)!))], {i, 
      Length[ww]}]];
  const = Total[Table[x[i], {i, Length[ww]}]];
  cons1 = 
   ToExpression[
    StringReplace[
     ToString[Table[x[i] >= 0, {i, Length[ww]}]], {"{" -> "", 
      "}" -> "", "," -> " &&"}]];
  int = ToExpression[
    StringReplace[
     ToString[
      Table[x[i] ∈ Integers, {i, Length[ww]}]], {"{" -> "", 
      "}" -> "", "," -> " &&"}]];
  argmin = ArgMin[{obj, const == ss && cons1 && int}, var];
  zer = Table[0, {i, 1, Length[vv] - Length[ww]}];
  fin = Join[argmin, zer];
  Fold[Insert[#1, 0, #2] &, DeleteCases[fin, 0], Position[vv, 0]]]
dean[{0, 40, 20, 12, 0}, 20]

It works nicely until s=20. But at s=21, it fails. It's the same thing if I increase the length of v. Of course, I have tried to change ArgMin to NArgMin and some methods but it doesn't ameliorate the situation. Mathematica complains of a division by 0. Is there a way to resolve this problem?

Answer 1

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.

Answer 2

If you change the constraint to "> 0" rather than ">= 0" it appears to work; however, I do not know if the results for larger values of s are reasonable since you did not provide the expected results.

dean[v_, s_] := Module[{vv = v},
  ww = DeleteCases[vv, 0];
  var = Array[x, Length[ww]];
  obj = Total[
    Table[Log[(x[i]!)^3 2^(x[i] - 
           1)/(x[i]^2 ww[[i]]^(x[i] - 1) ((2 x[i] - 1)!))], {i, 
      Length[ww]}]];
  const = Total[var];
  cons1 = And @@ Thread[var > 0];
  int = var ∈ Integers;
  argmin = ArgMin[{obj, const == s && cons1 && int}, var];
  zer = ConstantArray[0, Length[vv] - Length[ww]];
  fin = Flatten[{argmin, zer}];
  Fold[Insert[#1, 0, #2] &, DeleteCases[fin, 0], Position[vv, 0]]]

dean[{0, 40, 20, 12, 0}, #] & /@ {20, 21, 50, 100, 210}

(*  {{0, 11, 6, 3, 0}, {0, 11, 6, 4, 0}, {0, 28, 14, 8, 0}, {0, 55, 28, 
  17, 0}, {0, 117, 58, 35, 0}}  *)