Putting square pegs into square holes

Mathematica, 473 bytes

f@p_:=(s=Flatten@Round@p;v=Array[{x@#,y@#}&,n=Length@p];
  Do[w=Flatten[{g@#,h@#}&/@(b=Flatten@Position[p,x_/;Norm[x-p[[i]]]<=2,{1}])];f=Total[Norm/@(p-v)]+Total[If[#1==#2,1*^4,0]&@@@v~Subsets~{2}]/.Flatten[{x@#->g@#,y@#->h@#}&@@@w]/.Thread[Flatten@v->s];
    c=w∈Integers&&And@@MapThread[Max[#-2,0]<=#2<=Min[#+2,100]&,{Flatten@p[[b]],w}];s=Flatten@ReplacePart[s~Partition~2,Thread[b->Partition[w/.Last@Quiet@NMinimize[{f,c},w,MaxIterations->300],2]]]
    ,{i,n}]~Do~{2};s~Partition~2)

---

Before golfing:

f[p_]:=(n=Length@p;s=Flatten@Round@p;v=Array[{x[#],y[#]}&,n];
  Do[
    v2=Flatten[{x2[#],y2[#]}&/@(b=Flatten@Position[p,x_/;Norm[x-p[[i]]]<=2,{1}])];
    f2=Total[Norm/@(p-v)]+Total[If[#1==#2,1*^4,0]&@@@Subsets[v,{2}]]/.Flatten[{x[#]->x2[#],y[#]->y2[#]}&@@@v2]/.Thread[Flatten@v->s];
    c2=v2∈Integers&&And@@MapThread[Max[#-2,0]<=#2<=Min[#+2,100]&,{Flatten@p[[b]],v2}];
    s=Flatten@ReplacePart[s~Partition~2,Thread[b->Partition[v2/.Last@Quiet@NMinimize[{f2,c2},v2,MaxIterations->300],2]]];
    ,{i,n}]~Do~{2};
  s~Partition~2)

---

Explanation:

This optimization problem is not difficult to describe in Mathematica. Given a list of points p of length n,

• the variables are x[i] and y[i]: v=Array[{x[#],y[#]}&,n],
• the function to minimize is the total of displacements: f=Total[Norm/@(p-v)],
• the constraints are: c=Flatten[v]∈Integers&&And@@(Or@@Thread[#1!=#2]&@@@Subsets[v,{2}]).

And, NMinimize[{f,cons},v,MaxIterations->Infinity] will give the result. But unfortunately, such straight forward scheme seems too complicated to converge.

To work around the problem of complexity, two techniques are adopted:

• a large "interaction", If[#1==#2,1*^4,0]&, is used to avoid collision between points.
• instead of optimize all variable at the same time, we optimize on every point with their neighbors in turn.

We start from an initial guess by rounding the points. When the optimizations are done one by one, collisions are expected to be resolved, and an optimized arrangement is established.

The final solution is at least good, if not optimal. (I believe :P)

---

Result:

The result of Just for fun is shown below. Dark green points are the inputs, gray squares are the outputs, and black lines shows the displacements.

The sum of displacements is 19.4595. And the solution is

{{15,3},{5,4},{13,4},{2,5},{8,6},{21,6},{20,5},{19,5},{17,1},{17,3},{4,8},{14,6},{15,6},{13,7},{11,5},{16,5},{13,2},{22,8},{19,6},{21,7},{16,8},{12,9},{14,3},{13,5},{6,9},{10,7},{3,6},{22,7},{20,6},{8,4},{20,7},{18,4},{10,9},{17,6},{11,4},{2,8},{19,7},{22,6},{18,3},{10,8},{15,4},{10,3},{5,6},{21,8},{18,5},{2,9},{18,6},{14,8},{7,7}}

Python 3, 877 bytes

This is not a correct implementation. It fails on the second of the "further test cases", producing a solution with a total distance of 13.5325, where the solution provided needs only 13.2433. Further complicating matters is the fact that my golfed implementation doesn't match the ungolfed one I wrote first...

However, nobody else has answered, and this is too interesting a challenge to let slip past. Also, I have a picture generated from the USA data, so there's that.

The algorithm is something like this:

1. Push all points to the nearest integer coordinates (hereafter called a "square").
2. Find the square with the greatest number of points.
3. Find the lowest-cost redistribution of those points to the nine-square neighbourhood of this one, excluding any squares that have already been processed in step 2.
   • The redistribution is limited to one point per square, unless that wouldn't provide enough squares (although even then, only one point will remain on this square).
4. Repeat from step 2 until no square has more than one point.
5. Locate each of the original points, in order, and output their squares, in order.

I have absolutely no proof of optimality for any part of this algorithm, just a strong suspicion that it'll provide "pretty good" results. I think that's what we called a "heuristic algorithm" back in my uni days...!

l=len
I,G,M=-1,101,150
d=lambda x,y,X,Y:abs(x-X+1j*(y-Y))
N=(0,0),(I,0),(0,I),(1,0),(0,1),(I,I),(1,I),(1,1),(I,I)
n=lambda p,e:[(x,y)for(x,y)in(map(sum,zip(*i))for i in zip([p]*9,N))if(x,y)not in e and I<x<G and I<y<G]
def f(p):
 g={};F=[];O=[I]*l(p)
 for P in p:
  z=*map(round,P),
  if z in g:g[z]+=[P]
  else:g[z]=[P]
 while l(g)<l(p):
  L,*P=0,
  for G in g:
   if l(g[G])>l(P):L,P=G,g[G]
  o=n(L,F);h=l(o)<l(P);c=[[d(*q,*r)for r in o]for q in P];r={}
  while l(r)<l(c):
   A=B=C=M;R=S=0
   while R<l(c):
    if R not in r:
     z=min(c[R])
     if z<A:B,A=R,z;C=c[R].index(A)
    R+=1
   while S<l(c):
    if S==B:
     v=0
     while v<l(c[S]):
      if v!=C:c[S][v]=M
      v+=1
    elif C<1or not h:c[S][C]=M
    S+=1
   r[B]=C
  for q in r:
   x,y=P[q],o[r[q]]
   if y==L or y not in g:g[y]=[x]
   else:g[y]+=[x]
  F+=[L]
 for G in g:
  O[p.index(g[G][0])]=G
 return O

And the result of running it on the USA data (thanks to a utility function that turns the results into SVG):

This is slightly worse than the one the ungolfed code produced; the only visible difference is that the top-rightmost square is one further to the left in the better one.