Winner Paulo Uribe (dm1)

2004-04-28 09:00:00 UTC

# Imagine that

by cyclist

Status: Passed
Results: 0.219%
CPU Time: 9.894
Score: 4.5487
Submitted at: 2004-04-28 15:49:35 UTC
Scored at: 2004-04-28 15:50:12 UTC

Current Rank: 24th
Based on: chkIntMap (diff)

cyclist
28 Apr 2004
Inspired by vague memories that "i" and "j" are not the best variable names to use, because Matlab needs to check that the user does not mean sqrt(-1).
Code
```function d = solver(map,n)

goal = sum(map(:))/n;
map = map/goal;

d = solver0(map,n,goal);
s = myscore(map,d,n)*goal;

if s>801 && s<1999
d1=crystallise6(map,n);
s1=myscore(map,d1,n)*goal;
if s1<s
d=d1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = solver0(map,n,goal)
[s1,s2] = size(map);

d = solver1(map,n,s1,s2,goal);
s=myscore(map,d,n)*goal;
if s<21
return
end

if ~any(map(1,:)) || ~any(map(s1,:))
notemptyrows = find(any(map,2));
row1=notemptyrows(1);
rowe=notemptyrows(end);
trim=1;
else
trim=0;
row1=1;
rowe=s1;
end

if ~any(map(:,1)) || ~any(map(:,s2))
notemptycols = find(any(map));
col1=notemptycols(1);
cole=notemptycols(end);
trim=1;
else
trim=0;
col1=1;
cole=s2;
end
if trim
d1=d;
d1(row1:rowe,col1:cole)=solver1(map(row1:rowe,col1:cole),n,rowe-row1+1,cole-col1+1,goal);
sc1=myscore(map,d1,n)*goal;
if sc1<s
s=s1;
d1(1:row1-1,:)=d1(row1,col1);
d1(rowe+1:s1,:)=d1(rowe,cole);
d1(:,1:col1-1)=d1(row1,col1);
d1(:,cole+1:s2)=d1(rowe,cole);
d=d1;
end
end

if s<80|| (s>185&&s<6000)||s>8000
return
end

% Use prince of Darkness
if s1*s2 == 20
d_2  =  ozzy(map,n);
if myscore(map,d_2,n)*goal<s
d = d_2;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = solver1(map,n,s1,s2,goal)
is1=s1:-1:1;
is2=s2:-1:1;

for ndivs = [23,19,17,13,11,7,5,3,2]
if (rem(n,ndivs) == 0) && (rem(s1,ndivs) == 0)
k = s1/ndivs; mm = map(1:k,:);
idivs = (1:k)'; idivs = idivs(:,ones(1,ndivs));
MM = mm(idivs,:);
if isequal(map,MM)
n1 = n/ndivs;
d1 = solver1(mm,n1,k,s2,goal);
d = map;
for ccc = 1:ndivs
d((ccc-1)*k+1:ccc*k,:) = d1 + (ccc-1)*n1;
end
return
end
end
if (rem(n,ndivs) == 0) && (rem(s2,ndivs) == 0)
k = s2/ndivs; mm = map(:,1:k);
idivs = (1:k)'; idivs = idivs(:,ones(1,ndivs));
MM = mm(:,idivs);
if isequal(map,MM)
n1 = n/ndivs;
d1 = solver1(mm,n1,s1,k,goal);
d = map;
for ccc = 1:ndivs
jk=ccc*k;
d(:,jk-k+1:jk) = d1 + (ccc-1)*n1;
end
return
end
end
end

% is this an all-integer map?
if any( any( rem(map*goal,1)>0 ) )
minscore=0;
else
minscore = n*abs(goal-round(goal))/goal + 1e-13;
end

len = s1*s2;
limiet=n/570;
limiet = max(limiet,minscore);
np = 11;
D = cell(1,8);
siz = [s1 s2];
siz2= [s2 s1];
D{1} = map;
D{2} = map';
D{3} = map(is1,:);
D{4} = D{2}(is2,:);
D{5} = map(:,is2);
D{6} = D{2}(:,is1);
D{7} = D{3}(:,is2);
D{8} = D{4}(:,is1);
Score = inf;
dec =0;
ind = cell(1,np*8);

stop = 0;

pat3x3a = [1  8  9; ...
2  7  6; ...
3  4  5];
pat3x3b = [1  2  9; ...
4  3  8; ...
5  6  7];
pat4x4 = [1  2 15 16; ...
4  3 14 13; ...
5  8  9 12; ...
6  7 10 11];
pat5x5 = [1 20 21 22 25; ...
2 19 18 23 24; ...
3 16 17 12 11; ...
4 15 14 13 10; ...
5  6  7  8  9];

%mj = 1;
init_idx=zeros(22);

for bbb = 1:8
nmap = D{bbb}*0+n;
%    mj=~mj;
mj = rem(bbb-1,2);
for ccc=[0:5 7:np-2]

idx=ccc+mj*np+1;
if ~init_idx(idx)
if idx==9
ind{idx} = refind1(siz);
elseif idx==11
ind{idx} = refind2(siz);
elseif idx==4
ind{idx} = refind3(siz,2);
elseif idx==3
ind{idx} = refind3(siz,3);
elseif idx==8
ind{idx}= refind4(siz);
elseif idx==1
ind{idx} = refind5(siz,pat3x3a);
elseif idx==2
ind{idx} = refind5(siz,pat3x3b);
elseif idx==7
ind{idx} = refind5(siz,pat5x5);
elseif idx==5
ind{idx} = refind5(siz,pat4x4);
elseif idx==6
ind{idx} = refind3(siz,4);
elseif idx==10
ind{idx} = refind3(siz,5);
end

if s1 ~= s2
if idx==20
ind{idx} = refind1(siz2);
elseif idx==22
ind{idx} = refind2(siz2);
elseif idx==15
ind{idx} = refind3(siz2,2);
elseif idx==14
ind{idx} = refind3(siz2,3);
elseif idx==19
ind{idx} = refind4(siz2);
elseif idx==12
ind{idx} = refind5(siz2,pat3x3a);
elseif idx==13
ind{idx} = refind5(siz2,pat3x3b);
elseif idx==18
ind{idx} = refind5(siz2,pat5x5);
elseif idx==16
ind{idx} = refind5(siz2,pat4x4);
elseif idx==17
ind{idx} = refind3(siz2,4);
elseif idx==21
ind{idx} = refind3(siz2,5);
end
else
for  jj = 1:np
ind{np+jj} = ind{jj};
init_idx(np+jj)=1;
end
end
init_idx(idx)=1;
end

d1 = pattern(D{bbb},n,len,ind{idx},nmap);

Score1 = myscore(D{bbb},d1,n);
if Score1<Score
Score=Score1;
d=d1;
dec=bbb;
if Score1 < limiet
dec = bbb;
stop = 1;
break
end
end
end
if stop
break;
end
end

switch dec;
case 2
d = d';
case 3
d = d(is1,:);
case 4
d = d(is2,:)';
case 5
d = d(:,is2);
case 6
d = d(:,is1)';
case 7
d = d(is1,is2);
case 8
d = d(is2,is1)';
end
if Score*n < 5e-2, return; end
d = mysolver(map, n, d, siz,goal);
%s=myscore(map,d,n);

%if (s>80&&s<185)
%    if s1*s2==20
%        d  =  ozzy(map,n);
%    end
%    return
%end
%if s>6000&&s<8000
%    if s1*s2==20
%        d  =  ozzy(map,n);
%    end
%
%end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = refind1(s)
nel = s(1)*s(2);
X = reshape(1:nel,s(1),s(2));
X(:,2:2:s(2)) = X(s(1):-1:1,2:2:s(2));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = refind2(siz)
a=siz(1);b=siz(2);
d = zeros(a,b);
ss = a*b;
upper = 1;
lower = a;
jj = 1;
kk = b;
ii = 1;

while ii <= ss,
dd = lower-upper;
d(upper:lower,jj) = (ii:ii+dd)';
ii = ii + dd + 1;
jj = jj + 1;
if ii > ss, break; end

dd = kk-jj;
d(lower,jj:kk) = (ii:ii+dd);
ii = ii + dd + 1;
lower = lower - 1;
if ii > ss, break; end

dd = lower-upper;
d(lower:-1:upper,kk) = (ii:ii+dd)';
ii = ii + dd + 1;
kk = kk - 1;
if ii > ss, break; end

dd = kk-jj;
d(upper,kk:-1:jj) = (ii:ii+dd);
ii = ii + dd + 1;
upper = upper + 1;
if ii > ss, break; end
end

X = zeros(1,ss);
for dec = 1:ss,
X(d(dec)) = dec;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function X = refind3(siz,nr)
a=siz(1);b=siz(2);
d = zeros(a,b);
len = a*b;
bbb = 1;

dir = 1;
row = 1;
col = 1;
c = 1;

if mod(b,2) == 0,
d(:,1) = (a:-1:1)';
bbb = bbb + a;
c = 2;
col = 2;
row = 1;
end

q = mod(a,nr);
while q
if dir > 0
d(row,c:b) = bbb:bbb+b-c;
else
d(row,b:-1:c) = bbb:bbb+b-c;
end
bbb = bbb+b-c+1;
col = col + (b-c)*dir;
dir = -dir;
row = row + 1;
q = q - 1;
end

while bbb <= len,
while bbb <= len
for g = 0:nr-1
d(row+g,col) = bbb+g;
end
bbb = bbb + nr;
if (col == c && dir == -1) || (col == b && dir == 1)
break;
end
col = col + dir;
for g = 0:nr-1
d(row+nr-1-g,col) = bbb+g;
end
bbb = bbb + nr;
col = col + dir;
end
row = row + nr;
dir = -dir;
end

X = zeros(1,len);
for dec = 1:len,
X(d(dec)) = dec;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Indices = refind4(siz)
s1 = siz(1);
s2 = siz(2);
bbb = 1;
d = zeros(s1,s2);
r = s1;
c = s2;

x = 1;
y = 1;
while r >= 3 && c >= 2
d(y,x:end) = bbb:bbb+c-1;
bbb = bbb + c;
d(y+1,end:-1:x) = bbb:bbb+c-1;
bbb = bbb + c;
y = y + 2;
r = r - 2;
d(y:end,x) = (bbb:bbb+r-1)';
bbb = bbb + r;
d(end:-1:y,x+1) = (bbb:bbb+r-1)';
bbb = bbb + r;
x = x + 2;
c = c - 2;
end

d(y:end,x:end) = reshape((1:r*c)+bbb-1,c,r)';
d(y+1:2:end,end:-1:x) = d(y+1:2:end,x:end);

l = s1*s2;
Indices = zeros(1,l);
for bbb = 1:l
Indices(d(bbb)) = bbb;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Indices = refind5(siz,p1)
s1 = siz(1);
s2 = siz(2);
d = zeros(s1,s2);

[gx,gy]=size(p1);
p2 = p1(:,gy:-1:1);

ny = floor((s1-1)/(2*gy))*2;
nx = floor((s2-1)/gx);
bbb = 1;
r = 1;
c = 1;
dir = gx;
for y = 1:ny
for x = 1:nx
if dir > 0
d(r:r+gy-1,c:c+gx-1) = p1+(bbb-1);
else
d(r:r+gy-1,c:c+gx-1) = p2+(bbb-1);
end
bbb = bbb + gx*gy;
c = c + dir;
end
c = c - dir*(gx-1)/gx;
if dir < 0
d(r:r+gy-1,c) = (bbb:bbb+gy-1)';
else
d(r:r+gy-1,c+gx-1) = (bbb:bbb+gy-1)';
end
bbb = bbb + gy;
dir = -dir;
r = r + gy;
end

dr = s1 - r + 1;
dc = gx*nx+1;
d(r:end,1:dc) = reshape(bbb:bbb+dr*dc-1,dr,dc);
d(r:end,2:2:dc) = d(end:-1:r,2:2:dc);
if dc < s2
c = dc + 1;
bbb = bbb + dr*dc;
dc = s2 - c + 1;
d(1:end,c:s2) = reshape(bbb:bbb+s1*dc-1,s1,dc);
d(1:end,c:2:end) = d(end:-1:1,c:2:end);
end

l = s1*s2;
Indices = zeros(1,l);
for bbb = 1:l
Indices(d(bbb)) = bbb;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = pattern(map,n,l,ind,d)
% function d = pattern(map,n,l,ind)
% d = map*0+n;
w = 1;
ij = 1;
ji = l;
md = max(1,floor((n-1)/2));
k = l+1;
for bbb = 1:l,
ccc = ind(bbb);
m = map(ccc);
if ji<=n-ij || w+w < m ,
if ij == md,
ij = ij + 1;
k = bbb;
break;
end
ij = ij + 1;
w = 1-m;
else
w = w - m;
end
ji = ji-1;
d(ccc) = ij;
end
w = 1;
for bbb = l:-1:k
ccc = ind(bbb);
m = map(ccc);
if ji<=n-ij || w+w < m ,
if ij == n-1,
break;
end
ij = ij + 1;
w = 1-m;
else
w = w - m;
end
ji = ji-1;
d(ccc) = ij;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = mysolver(map,n,d,siz,goal)

s1=siz(1);s2=siz(2);
rand('state',101);
w = zeros(1,n);
for bbb = 1:s1*s2
w(d(bbb)) = w(d(bbb)) + map(bbb);
end

Score = sum(abs(w-1)) / n;

nl = 3000;

if Score > 0.007241
nl = 5000;
end

pop = rand(3,nl);
pop(1,:)=pop(1,:)*(s1-2)+2;
pop(2,:)=pop(2,:)*(s2-2)+2;
pop(3,:)=pop(3,:)*4;
pop=double(uint16(pop));

p=uint16(d);

lastcheck = 0;
bestscore = Score;

for kk=1:nl

y1 = pop(1,kk);
x1 = pop(2,kk);
y2 = y1;
x2 = x1;

switch pop(3,kk)
case 0
y2 = y2 + 1;
case 1
y2 = y2 - 1;
case 2
x2 = x2 + 1;
case 3
x2 = x2 - 1;
end

row1 = p(y1,x1);
rowe = p(y2,x2);
if row1 == rowe
continue;
end

a = w(row1);
b = w(rowe);
da = b-a;
if da>0
p(y2,x2) = row1;
m = map(y2,x2);

if da <= abs(da-2*m) || isNCont(p==rowe,siz)
p(y2,x2) = rowe;
else
w(rowe) = b - m;
w(row1) = a + m;
end
else
p(y1,x1) = rowe;
m = map(y1,x1);

if da >= -abs(da+2*m) || isNCont(p==row1,siz)
p(y1,x1) = row1;
else
w(rowe) = b + m;
w(row1) = a - m;
end
end

if kk-lastcheck > 500
s = sum(abs(w-1)) / n;
if s - bestscore > -0.000001
break;
else

bestscore = s;
lastcheck = kk;
end
end
end
d=double(p);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function tf = isNCont(a,siz)

s1=siz(1);s2=siz(2);

[ii,jj] = find(a);
tf = length(ii);
if tf==0;return;end
dec = 1;
a(ii(1),jj(1)) = 0;
tf = tf - 1;
while dec
r = ii(dec);
c = jj(dec);
dec = dec-1;
if (r > 1) && a(r-1,c)
dec = dec+1;
ii(dec) = r-1;
jj(dec) = c;
a(r-1,c) = 0;
tf = tf - 1;
end
if (r < s1) && a(r+1,c)
dec = dec+1;
ii(dec) = r+1;
jj(dec) = c;
a(r+1,c) = 0;
tf = tf - 1;
end
if (c > 1) && a(r,c-1)
dec = dec+1;
ii(dec) = r;
jj(dec) = c-1;
a(r,c-1) = 0;
tf = tf - 1;
end
if (c < s2) && a(r,c+1)
dec = dec+1;
ii(dec) = r;
jj(dec) = c+1;
a(r,c+1) = 0;
tf = tf - 1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Score = myscore(ma,in,n)
Total = zeros(1,n);
for bbb = 1:numel(ma),
Total(in(bbb)) = Total(in(bbb)) + ma(bbb);
end
Score = abs(Total(1)-1);
for k=2:n
Score = Score + abs(Total(k)-1);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = ozzy(map,n)

[xn,yn] = size(map);
num = xn*yn;

d=zeros(xn,yn);
[e,w]=sort(-map(:));
d(w(1:n))=1:n;

liberties = zeros(n,num);
libnz=zeros(n,1);
groupsize(1:n) = map(w(1:n));

ym=1:yn;
ymesh=ym(ones(xn,1),:);
xm=(1:xn)';
xmesh=xm(:,ones(1,yn));

for bbb=1:n
wi=w(bbb);
nzi=libnz(bbb);
x = xmesh(wi); y = ymesh(wi);
if x>1 && ~d(x-1,y),
nzi=nzi+1;
liberties(bbb,nzi)=wi-1;
end
if x<xn && ~d(x+1,y)
nzi=nzi+1;
liberties(bbb,nzi)=wi+1;
end
if y>1 && ~d(x,y-1)
nzi=nzi+1;
liberties(bbb,nzi)=wi-xn;
end
if y<yn && ~d(x,y+1)
nzi=nzi+1;
liberties(bbb,nzi)=wi+xn;
end
libnz(bbb)=nzi;
end

num=num-n;
while num
if all(groupsize==inf)
break
end
[smallest,bbb]=min(groupsize);
nzi = libnz(bbb);
notzero = (liberties(bbb,1:nzi));
test=~d(notzero);
if ~any(test) || ~nzi
groupsize(bbb)=inf;
continue;
end
if ~all(test),
nzi=sum(test);
notzero = notzero(test);
end

notzero = sort(notzero);
choice = map(notzero);
[tmp,move]=max(choice);
nzm = notzero(move);
liberties(bbb,1:nzi)=notzero;
if ~d(nzm)
x = xmesh(nzm); y = ymesh(nzm);
if x>1 && ~d(x-1,y)
nzi=nzi+1;
liberties(bbb,nzi)=nzm-1;
end
if x<xn && ~d(x+1,y)
nzi=nzi+1;
liberties(bbb,nzi)=nzm+1;
end
if y>1 && ~d(x,y-1)
nzi=nzi+1;
liberties(bbb,nzi)=nzm-xn;
end
if y<yn && ~d(x,y+1)
nzi=nzi+1;
liberties(bbb,nzi)=nzm+xn;
end
num = num - 1;
d(nzm) = bbb;
groupsize(bbb) = groupsize(bbb)+ tmp; %map(nzm);
liberties(bbb,move:nzi-1)=liberties(bbb,move+1:nzi);
libnz(bbb)=nzi-1;
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = crystallise6(map,n)

[m1,m2] = size(map);
m=m1*m2;
m11=m1+2;
m22=m2+2;
targ = 1;

dpop=zeros(n,1);    % population of district
dsq = zeros(n,1);   % number of squares in district

ihalos = zeros(n,m);
nhalos = zeros(n,1);

% pad d matrix with ones
e=ones(m11,m22);
e(2:m1+1,2:m2+1)=0;
d=e;
e(:)=10000;
e(2:m1+1,2:m2+1)=map;
map=e;
hmap=e;

apop=0;  %assigned population to date

for idis=1:n
% choose starting pixel - highest unassigned pixel

[mmax,imax]=max(~d(:).*map(:));
i1=imax;
d(i1)=idis;
dpop(idis)=map(i1);
apop = apop+dpop(idis);

dsq(idis)=1;
% make the halo list
nhalo=0;
halo=zeros(1,m);
notstuck=1;

hmap(:)=d(:);

% make the halo for the first point
if ~hmap(i1-1)
nhalo=nhalo+1;
halo(nhalo)=i1-1;
hmap(i1-1) = 1;
end
if ~hmap(i1+1)
nhalo=nhalo+1;
halo(nhalo)=i1+1;
hmap(i1+1) = 1;
end
if ~hmap(i1-m11)
nhalo=nhalo+1;
halo(nhalo)=i1-m11;
hmap(i1-m11)=1;
end
if ~hmap(i1+m11)
nhalo=nhalo+1;
halo(nhalo)=i1+m11;
hmap(i1+m11)=1;
end
if nhalo>0
halo=halo(1:nhalo);
end

oldhalo=halo;

idtarg=idis*targ;

updatehalo=0; halo_up=1;

while notstuck   &&   apop < idtarg
if nhalo>0                  % if there is the possibility of adding a point
halo=halo(1:nhalo);
notstuck=1;
% choose the best halo point to add, and add it
best = Inf;
for ih=1:nhalo
%testval = abs(dpop(idis) + map(halo(ih)) - targ);
testval = abs( apop + map(halo(ih)) - idtarg);
if testval < best
best = testval;
end
end
if abs(dpop(idis)+map(i1)-targ) < abs(dpop(idis)-targ)   % will this point make things better?
d(i1)=idis;     % add the pixel
halo=halo(2:nhalo);
halo=halo(1:nhalo-1);
else
end
halo_up=0;
nhalo=nhalo-1;
dpop(idis)  =  dpop(idis) + map(i1);
apop  =  apop + map(i1);
dsq(idis)=dsq(idis)+1;
updatehalo=0;
elseif ~halo_up             % if the halo is not up to date, update it and try again
updatehalo=1;
else
notstuck=0;   % go no further
end
elseif ~halo_up               % if the halo is not up to date, update it and try again
updatehalo=1;
else
notstuck=0;
end   % if nhalo>0

% constuct the halo of the last added point, [i1]
if updatehalo
for ia = 1:length(iadded)
if ~hmap(i1-1)
nhalo=nhalo+1;
halo(nhalo)=i1-1;
hmap(i1-1) = 1;
end
if ~hmap(i1+1)
nhalo=nhalo+1;
halo(nhalo)=i1+1;
hmap(i1+1) = 1;
end
if ~hmap(i1-m11)
nhalo=nhalo+1;
halo(nhalo)=i1-m11;
hmap(i1-m11)=1;
end
if ~hmap(i1+m11)
nhalo=nhalo+1;
halo(nhalo)=i1+m11;
hmap(i1+m11)=1;
end
if nhalo>0
halo=halo(1:nhalo);
end
end
oldhalo=halo;
halo_up=1;
end

end     %while - working through this district and its halo

% Finish with this district:
% constuct the halo of the last added point, [i1,j2] - it won't
% have been updated inside the while loop

% update teh halo before leaving this district
if ~halo_up
for ia = 1:length(iadded)
if ~hmap(i1-1)
nhalo=nhalo+1;
halo(nhalo)=i1-1;
hmap(i1-1) = 1;
end
if ~hmap(i1+1)
nhalo=nhalo+1;
halo(nhalo)=i1+1;
hmap(i1+1) = 1;
end
if ~hmap(i1-m11)
nhalo=nhalo+1;
halo(nhalo)=i1-m11;
hmap(i1-m11)=1;
end
if ~hmap(i1+m11)
nhalo=nhalo+1;
halo(nhalo)=i1+m11;
hmap(i1+m11)=1;
end
if nhalo>0
halo=halo(1:nhalo);
end
end
end

nhalos(idis)=nhalo;
if nhalo>0
ihalos(idis,1:nhalo)=halo;
end
end % loop over districts

% UNASSIGNED SQUARES - GROW DISTRICTS THAT ARE TOO SMALL

% unassigned squares? grow all districts as far as possible, starting with
% the smallest
if sum(dsq)<m
[sss,idsort] = sort(dpop-targ);
idsort=idsort(sss<0);
for k=1:length(idsort);
hmap=d;
idis=idsort(k);
% Retrieve the halo for district idis, and delete entries that have
% become occupied since halo was constructed.
halo = ihalos(idis,:);
nhalo = nhalos(idis);
for ih = 1:nhalo
if d(halo(ih))
halo(ih)=0;
end
end
halo(halo==0)=[];
nhalo = length(halo);
nhalos(idis)=nhalo;
if nhalos(idis)>0
ihalos(idis,1:nhalos(idis))=halo;
end
for ih = 1:nhalo
hmap(halo(ih))=1;
end
% halo is now updated. Start growing the district as before, this time
% until it gets stuck.
notstuck=1; updatehalo=0;
if nhalo==0
notstuck=0;
end

while dpop(idis)<targ&&notstuck
if updatehalo     % construct the halo of the last added point, [i1,j2]
if ~hmap(i1-1)
nhalo=nhalo+1;
halo(nhalo)=i1-1;
hmap(i1-1) = 1;
end
if ~hmap(i1+1)
nhalo=nhalo+1;
halo(nhalo)=i1+1;
hmap(i1+1) = 1;
end
if ~hmap(i1-m11)
nhalo=nhalo+1;
halo(nhalo)=i1-m11;
hmap(i1-m11)=1;
end
if ~hmap(i1+m11)
nhalo=nhalo+1;
halo(nhalo)=i1+m11;
hmap(i1+m11)=1;
end
end
if nhalo>0
notstuck=1;
nhalo=length(halo);
% choose the best halo point to add, and add it
mmax = -1;
for ih=1:nhalo
if map(halo(ih))>mmax
mmax=map(halo(ih));
end
end
d(i1)=idis;
nhalo=nhalo-1;
apop  =  apop + map(i1);
dpop(idis)=dpop(idis)+map(i1);
dsq(idis)=dsq(idis)+1;
updatehalo=1;        % a square has been added - next iteration will need halo to be reconstructed
else
notstuck=0;
end
end     % while growing district in second pass

% Finish with this district:
% construct the halo of the last added point, [i1,j2] - it won't
% have been updated inside the while loop
if updatehalo
if ~hmap(i1-1)
nhalo=nhalo+1;
halo(nhalo)=i1-1;
hmap(i1-1) = 1;
end
if ~hmap(i1+1)
nhalo=nhalo+1;
halo(nhalo)=i1+1;
hmap(i1+1) = 1;
end
if ~hmap(i1-m11)
nhalo=nhalo+1;
halo(nhalo)=i1-m11;
hmap(i1-m11)=1;
end
if ~hmap(i1+m11)
nhalo=nhalo+1;
halo(nhalo)=i1+m11;
hmap(i1+m11)=1;
end
if nhalo>0
halo=halo(1:nhalo);
end
end

% store the halo
nhalos(idis)=nhalo;
if nhalo>0
ihalos(idis,1:nhalo)=halo;
end

end      % phase 2 loop over districts
end

% UNASSIGNED SQUARES - GROW ANY DISTRICTS

% unassigned squares? grow all districts as far as possible, starting with
% the smallest
if sum(dsq)<m
[sss,idsort] = sort(dpop-targ);
for k=1:n
hmap=d;
idis=idsort(k);
% Retrieve the halo for district idis, and delete entries that have
% become occupied since halo was constructed.
nhalo = nhalos(idis);
halo = ihalos(idis,1:nhalo);
for ih = 1:nhalo
if d(halo(ih))
halo(ih)=0;
end
end
halo(halo==0)=[];
nhalo = length(halo);
nhalos(idis)=nhalo;
if nhalos(idis)>0
ihalos(idis,1:nhalo)=halo;
end
for ih = 1:nhalo
hmap(halo(ih))=1;
end
% halo is now updated. Start growing the district as before, this time
% until it gets stuck.
notstuck=1; updatehalo=0;
if nhalo==0
notstuck=0;
end

while notstuck
if updatehalo     % construct the halo of the last added point, [i1,j2]
if ~hmap(i1-1)
nhalo=nhalo+1;
halo(nhalo)=i1-1;
hmap(i1-1) = 1;
end
if ~hmap(i1+1)
nhalo=nhalo+1;
halo(nhalo)=i1+1;
hmap(i1+1) = 1;
end
if ~hmap(i1-m11)
nhalo=nhalo+1;
halo(nhalo)=i1-m11;
hmap(i1-m11)=1;
end
if ~hmap(i1+m11)
nhalo=nhalo+1;
halo(nhalo)=i1+m11;
hmap(i1+m11)=1;
end
end
if nhalo>0
notstuck=1;
nhalo=length(halo);	%?????????
% choose the best halo point to add, and add it
mmax = -1;
for ih=1:nhalo
if map(halo(ih))>mmax
mmax=map(halo(ih));
end
end
d(i1)=idis;
nhalo=nhalo-1;
dpop(idis)=dpop(idis)+map(i1);
dsq(idis)=dsq(idis)+1;
updatehalo=1;        % a square has been added - next iteration will need halo to be reconstructed
else
notstuck=0;
end
end     % while growing district in second pass

end      % loop over districts
end

d = d(2:m1+1,2:m2+1);

% Empty districts: try to create single-square districts that don't cause
% contiguity problems
emptydis = find(~dsq);
for k = 1:length(emptydis)
idis = emptydis(k);
for i1=1:m1
for j1=1:m2
if ~critcon(d,i1,j1,m1,m2)
id = d(i1,j1);
d(i1,j1) = idis;
dpop(id) = dpop(id)-map(i1,j1);
dpop(idis) = dpop(idis)+map(i1,j1);
idis=idis+1;
if idis>n
break
end
end
end
if idis>n
break
end
end
end   % if
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d = trade(d, map, n)

% Init some vars
popcount = full(sparse(d(:),1,map(:))');
totalpop = sum(popcount);
targetpop = totalpop / n;
err = popcount - targetpop;
[r,c] = size(map);

% Trades squares between districts to improve score
notdone = 1;
while notdone

% Find district boundaries and get indices to those points
ddiff = diff(d); % Down boundary
[dy,dx] = find(ddiff);
dndx1 = (dx - 1)*r + dy; % above the boundry
dndx2 = dndx1 + 1;       % below the boundry
rdiff = diff(d,1,2); %diff(d')'; % Right boundary
[ry,rx] = find(rdiff);
rndx1 = (rx - 1)*r + ry;
rndx2 = rndx1 + r;
ndx1 = [dndx1; rndx1]'; % Combine down and right
ndx2 = [dndx2; rndx2]';
if isempty(ndx1)
notdone = 0;
continue
end

% Get district of each square

dist1 = d(ndx1);
dist2 = d(ndx2);

% Check if population errors are in opposite sense
err1 = err(dist1);
err2 = err(dist2);
ndx = sign(err1) ~= sign(err2);
ndx1  = ndx1(ndx);
if isempty(ndx1)
notdone = 0;
continue
end
ndx2  = ndx2(ndx);
dist1 = dist1(ndx);
dist2 = dist2(ndx);
err1  = err1(ndx);
err2  = err2(ndx);

% Swap vars so dist1 (the giver) always gives a square to dist2 (the taker)
ndx = err1 < 0;
ndx1(ndx)  = ndx2(ndx);
%ndx2(ndx)  = dummy; % never used after this, so dropped
dummy      = dist1(ndx);
dist1(ndx) = dist2(ndx);
dist2(ndx) = dummy;
dummy      = err1(ndx);
err1(ndx)  = err2(ndx);
err2(ndx)  = dummy;

% Check if score will be improved
ndx = tradevalue < max(err1,abs(err2)); % do examples
ndx1  = ndx1(ndx);
if isempty(ndx1)
notdone = 0;
continue
end
dist1 = dist1(ndx);
dist2 = dist2(ndx);
err1  = err1(ndx);
err2  = err2(ndx);

% Check if trade benefit is large enough to be worthwhile
%ndx = [tradeben > 0.00001*totalpop]; % Limit amount of trading
ndx = tradeben > 0.0000011*totalpop; % Limit amount of trading
%ndx = [tradeben > 0]; % Make all trades - way too slow!
ndx1  = ndx1(ndx);
if isempty(ndx1)
notdone = 0;
continue

end
dist1 = dist1(ndx);
dist2 = dist2(ndx);

% Make as many trades as possible but don't affect a single district more
% than once during this loop
%ndx = flplr(ndx); % Max benefit first
k = 0;
for bbb = ndx(end:-1:1)

continue % One or both districts already traded in this iteration
end
if traded(ndx1(bbb)) % can only trade each square once ever
k = k + 1; % This square has been traded too many times
continue
end

% Consider this square traded

% Make the trade
d(ndx1(bbb)) = dist2(bbb);

% Check if districts are still contiguous
%		if ~isconnected(d,dist1(bbb))
if isNCont(d==dist1(bbb),[r c])
d(ndx1(bbb)) = dist1(bbb); % Undo trade
continue
end

% Update vars
err(dist1(bbb)) = err(dist1(bbb)) - tradevalue(bbb); % the giver
err(dist2(bbb)) = err(dist2(bbb)) + tradevalue(bbb); % the taker