Make a map of squares

I'm trying to use Mathematica to make a map similar to the one here

i.e. a US map where each state is represented by a square of a given size, located as close as possible to its true geographical position while not overlapping with other squares. I get the data for the squares as follows:

usa = Import[
"Data"];
transform[s_] := StringTrim[s, Whitespace ~~ "(" ~~ ___ ~~ ")"]
polygons =
"PlacemarkNames" /. usa[[1]]] -> ("Geometry" /. usa[[1]])];
(* Remove Alaska and Hawai *)
polygons =
Join[{polygons[[1]]}, polygons[[3 ;; 10]], polygons[[12 ;;]]];
centers = First /@ polygons[[All, 2, 1]];
area[poly_] :=
Apply[Plus,
Flatten[First@poly Map[({1, -1} Reverse[#] &),
RotateLeft[First@poly]]]]/2
areas = Table[Total[area /@ polygons[[i, 2, 2 ;;]]], {i, 1, 48}];
sizes = Sqrt /@ areas;


which uses the true area of the state as square area, which is enough for now. I can plot them:

Graphics[Table[{RGBColor[RandomReal[], RandomReal[], RandomReal[]],
Rectangle[{centers[[i, 1]] - sizes[[i]]/2,
centers[[i, 2]] - sizes[[i]]/2},
{centers[[i, 1]] + sizes[[i]]/2,
centers[[i, 2]] + sizes[[i]]/2}]}, {i, 1, 48}]]


Starting now, I demonstrate how far I've managed to go on my own. I want to move the squares to non-overlapping positions. I try to minimize the following score function:

dist[c1_, c2_, s1_, s2_] := Module[{dx, dy},
dx = Abs[c1[1] - c2[1]]; dy = Abs[c1[2] - c2[2]];
Max[0, dx - s1 - s2] + Max[0, dy - s1 - s2]
];
res = FindMinimum[
Sum[Sum[dist[c[i], c[j], sizes[[i]], sizes[[j]]], {j, i + 1, 48}], {i, 1, 48}],
Flatten[Table[{c[i][j], centers[[i, j]]}, {i, 1, 48}, {j, 1, 2}], 1]]


but plotting this gives terrible results:

Graphics[Table[{RGBColor[RandomReal[], RandomReal[], RandomReal[]],
Rectangle[{c[i][1] - sizes[[i]]/2, c[i][2] - sizes[[i]]/2},
{c[i][1] + sizes[[i]]/2, c[i][2] + sizes[[i]]/2}]}, {i, 1, 48}] /. res[[2]]]


I have not been able to see what is the issue with my code, because the score function clearly shouldn't be zero for this particular solution, but I am at loss to understand why. Any hint about that, or advice on tackling the problem in a different way, is welcome!

-
Looks like a job for belisarius' answer to spread out rectangles to make them non-overlapping. Perhaps also subject it to a shape constraint (like in Heike's answer to the word cloud question) – R. M. Sep 10 '12 at 21:00
@R.M yeah, I saw the word-cloud question (and answers). I'd prefer something that clearly uses Mathematica's optimizers rather than a specially crafted algorithm, as it is also intended as a demonstration of Mathematica's power to solve real-life tasks… – F'x Sep 10 '12 at 21:03
Did you see the link to belisarius' answer? There's a link to the code in the comments. Looks like it's a pretty straightforward application of that... I'll leave this for him to answer it for after all, he did all the hard work :) – R. M. Sep 10 '12 at 21:06
I had the same buttonbar problem tonight and I saw this complaint elsewhere too. – Sjoerd C. de Vries Sep 10 '12 at 21:09
@R.M The code in that answer is crap, because it was done (almost) procedurally on purpose. I'll see if I can remember what I did and recast it functionally – Dr. belisarius Sep 11 '12 at 5:35

First set up the lower 48 lists, along with a set of weighted offset vectors. I'll use these for initial points in a constrained optimization. They help insofar as they can be adjusted to satisfy, or almost satisfy, the inequality constraints (which, inexplicably, I call "penalty" below).

sl48 = Take[sizes, 48];
cenl48 = Take[centers, 48];
middle = Total[cenl48]/48;
wmean = sl48.cenl48/Total[sl48]
offsets = Map[# - wmean &, cenl48];


Our function will measure how far the new centers are from the initial ones. This is needed because otherwise we tend to get a result that is a tight fit of weighted squares, but bears little resemblance to the USA map.

f[centers : {{_?NumberQ, _?NumberQ} ..},
orig : {{_?NumberQ, _?NumberQ} ..}] := With[
{diffs = centers - orig},
Total@Map[#.# &, diffs]]

overlap[c1_, c2_, s1_, s2_] := Min[(s1 + s2)/2 - Abs[c1 - c2]] <= 0

penalty[centers_, wts_] :=

Flatten[Table[
overlap[centers[[i]], centers[[j]], wts[[i]], wts[[j]]], {i,
47}, {j, i + 1, 48}]]

cvars = Array[c, {48, 2}];

Timing[
res = FindMinimum[{f[cvars, cenl48], penalty[cvars, sizes]},
Transpose[{Flatten[cvars], Flatten[cenl48 + offsets]}],
MaxIterations -> 50, Method -> "InteriorPoint"];]

(* Out[344]= {136.37, Null} *)


We can plot this as follows.

colors = Table[
RGBColor[RandomReal[], RandomReal[], RandomReal[]], {48}];
rex2 = Transpose[{(cvars /. res[[2]]) - sl48/2, (cvars /. res[[2]]) +
sl48/2}];
Graphics[Transpose[{colors, Rectangle @@@ rex2}]]


-
I guess I should mention that this can probably be done as well or better, and much faster, if recast as a linear programming problem using Manhattan distances rather than Euclidean in the objective function (and slightly changing the constraints an explicitly linear formulation). – Daniel Lichtblau Sep 12 '12 at 15:20