Winner Paulo Uribe (dm1)

2004-04-28 09:00:00 UTC

# All Your Rectangles Are Belong To Us V2.6

Status: Passed
Results: 0.218%
CPU Time: 10.0693
Score: 4.5462
Submitted at: 2004-04-28 16:58:21 UTC
Scored at: 2004-04-28 17:05:50 UTC

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

Code
```function d = solver(map,n)

goal = sum(map(:))/n;
map = map/goal;
mapsize = prod(size(map));

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; s = s1; end

end

s_tol = 300; map_tol = 1000;

if (s > s_tol) && (mapsize < map_tol)
d2 = flowers2(map,n);
s2 = myscore(map,d2,n)*goal;
if s2 < s, d = d2; s = s2; end
end

if (s > s_tol) && (mapsize < map_tol)
d3 = seasons(d,map,n);
s3 = myscore(map,d3,n)*goal;
if s3 < s, d = d3; s = s3; end
end

% gotta beat: 0.327%, 4.92 secs

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function d = flowers2(map,N)

[m,n] = size(map); mn = m*n;

d = zeros(m,n);

% calculate index matrices
idx = d; idx(:) = [1:mn];
n1 = 1:n-1; n2 = [2:n];
m1 = 1:m-1; m2 = [2:m];

idx_up = d; idx_up(m2,:) = idx(m1,:);
idx_rt = d; idx_rt(:,n1) = idx(:,n2);
idx_dn = d; idx_dn(m1,:) = idx(m2,:);
idx_lt = d; idx_lt(:,n2) = idx(:,n1);

%%% INITIAL ALLOCATION

% divide map into roughly equal regions

rem_tol = 0.5;

Rd = m/n;
if (Rd > 1), rem_tol = 1/rem_tol; end

score_min = 1e12;

for mN = 1:min(m,N)
R = mN/ceil(N/mN);

if (abs(R - Rd) < rem_tol) && (ceil(N/mN) <= n)

nN = ceil(N/mN);

d = zeros(m,n);

%fprintf('%d  %d\n',mN,nN)

m2 = m/mN; n2 = n/nN;
M2 = ceil([0,[1:mN-1]*m2,m]); N2 = ceil([0,[1:nN-1]*n2,n]);

map_mxN = zeros(mN,nN); submap_sum = zeros(mN,nN);

% calculate the population of each region, and the maximum squares within
k = 1;
for i = 1:mN
for j = 1:nN
submap = map(M2(i)+1:M2(i+1),N2(j)+1:N2(j+1));
% choose the centre
map_mxN(i,j) = idx(round((M2(i)+1+M2(i+1))/2),round((N2(j)+1+N2(j+1))/2));
submap_sum(i,j) = sum(submap(:));
end
end

% start at the maximum square within the top N regions
[ans,submap_ssort] = sort(submap_sum(:));
map_mxN = map_mxN(submap_ssort(mN*nN:-1:mN*nN-N+1));

% allocate districts
d(map_mxN) = [1:N];

% take a guess on the max district size
D = zeros(N,ceil(mn/2));

% initialise some variables
sumD = zeros(1,N); numD = sumD; comD = sumD; jidx = zeros(1,4);

% set the starting squares in each district
for i = 1:N
D(i,1) = map_mxN(i);
sumD(i) = map(map_mxN(i));
numD(i) = 1;
end

totD = N; i = 1;

%%% GROWTH

while 1

%   fprintf('%4d / %4d    %3d / %3d\n',totD,mn,i,numD(i))

comD(i) = 1;

dn = numD(i);

for j = 1:dn % each square already occupied

jidx(1) = idx_up(D(i,j)); % the square above
jidx(2) = idx_rt(D(i,j)); % the square right
jidx(3) = idx_dn(D(i,j)); % the square below
jidx(4) = idx_lt(D(i,j)); % the square left

for k = 1:4
if jidx(k) && (d(jidx(k)) == 0) % square within bounds & not taken
comD(i) = 0;

d(jidx(k)) = i;
D(i,numD(i) + 1) = jidx(k);
sumD(i) = sumD(i) + map(jidx(k)); % increase the district population
numD(i) = numD(i) + 1; % increment the number of squares in the district
totD = totD + 1; % increment the total number of squares claimed
if (totD == mn), comD(i) = 1; break, end
end
end

if (totD == mn), break, end
end

if (totD == mn), break, end

[ans,i] = min(sumD + comD*1e12);

%   if ~rem(totD,1), plotmap(map,N,d), end

end

score = myscore(map,d,N);

if score < score_min
score_min = score;
d_min = d;
end

end
end

d = d_min;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function d = flowers3(map,N)

[m,n] = size(map); mn = m*n;

d = zeros(m,n);

% calculate index matrices
idx = d; idx(:) = [1:mn];
n1 = 1:n-1; n2 = [2:n];
m1 = 1:m-1; m2 = [2:m];

idx_up = d; idx_up(m2,:) = idx(m1,:);
idx_rt = d; idx_rt(:,n1) = idx(:,n2);
idx_dn = d; idx_dn(m1,:) = idx(m2,:);
idx_lt = d; idx_lt(:,n2) = idx(:,n1);

%%% INITIAL ALLOCATION

% divide map into roughly equal regions

rem_tol = 0.5;

Rd = m/n;
if (Rd > 1), rem_tol = 1/rem_tol; end

score_min = 1e12;

for mN = 1:min(m,N)
R = mN/ceil(N/mN);

if (abs(R - Rd) < rem_tol) && (ceil(N/mN) <= n)

nN = ceil(N/mN);

d = zeros(m,n);

%fprintf('%d  %d\n',mN,nN)

m2 = m/mN; n2 = n/nN;
M2 = ceil([0,[1:mN-1]*m2,m]); N2 = ceil([0,[1:nN-1]*n2,n]);

map_mxN = zeros(mN,nN); submap_sum = zeros(mN,nN);

% calculate the population of each region, and the maximum squares within
k = 1;
for i = 1:mN
for j = 1:nN
kidx = idx(M2(i)+1:M2(i+1),N2(j)+1:N2(j+1));
submap = map(M2(i)+1:M2(i+1),N2(j)+1:N2(j+1));
% choose the maximum
[ans,submap_mxn] = max(submap(:));
map_mxN(i,j) = kidx(submap_mxn);
submap_sum(i,j) = sum(submap(:));
end
end

% start at the maximum square within the top N regions
[ans,submap_ssort] = sort(submap_sum(:));
map_mxN = map_mxN(submap_ssort(mN*nN:-1:mN*nN-N+1));

% allocate districts
d(map_mxN) = [1:N];

% take a guess on the max district size
D = zeros(N,ceil(mn/2));

% initialise some variables
sumD = zeros(1,N); numD = sumD; comD = sumD; jidx = zeros(1,4);

% set the starting squares in each district
for i = 1:N
D(i,1) = map_mxN(i);
sumD(i) = map(map_mxN(i));
numD(i) = 1;
end

totD = N; i = 1;

%%% GROWTH

while 1

%   fprintf('%4d / %4d    %3d / %3d\n',totD,mn,i,numD(i))

comD(i) = 1; bodr_mx = -1;

dn = numD(i);

for j = 1:dn % each square already occupied

jidx(1) = idx_up(D(i,j)); % the square above
jidx(2) = idx_rt(D(i,j)); % the square right
jidx(3) = idx_dn(D(i,j)); % the square below
jidx(4) = idx_lt(D(i,j)); % the square left

for k = 1:4
if jidx(k) && (d(jidx(k)) == 0) % square within bounds & not taken
comD(i) = 0;
if (map(jidx(k)) > bodr_mx) % find the largest free bordering square
bodr_mxn = jidx(k);
bodr_mx = map(bodr_mxn);
end
%            A(i,d(jidx(k))) = 1;
end
end
end

if (~comD(i)) % take ownership of the square
d(bodr_mxn) = i;
D(i,numD(i) + 1) = bodr_mxn;
sumD(i) = sumD(i) + map(bodr_mxn); % increase the district population
numD(i) = numD(i) + 1; % increment the number of squares in the district
totD = totD + 1; % increment the total number of squares claimed
if (totD == mn), break, end
end

[ans,i] = min(sumD + comD*1e12);

%   if ~rem(totD,10), plotmap(map,N,d), end

end

score = myscore(map,d,N);

if score < score_min
score_min = score;
d_min = d;
end

end
end

d = d_min;

% -------------------------------------------------------

function d = seasons(d,map,N)

[m,n] = size(map); mn = m*n;

% calculate index matrices
idx = zeros(m,n); idx(:) = [1:mn];
n1 = 1:n-1; n2 = [2:n];
m1 = 1:m-1; m2 = [2:m];

idx_up = zeros(m,n); idx_up(m2,:) = idx(m1,:);
idx_rt = zeros(m,n); idx_rt(:,n1) = idx(:,n2);
idx_dn = zeros(m,n); idx_dn(m1,:) = idx(m2,:);
idx_lt = zeros(m,n); idx_lt(:,n2) = idx(:,n1);

idx_ur = zeros(m,n); idx_ur(m2,n1) = idx(m1,n2);
idx_rd = zeros(m,n); idx_rd(m1,n1) = idx(m2,n2);
idx_dl = zeros(m,n); idx_dl(m1,n2) = idx(m2,n1);
idx_lu = zeros(m,n); idx_lu(m2,n2) = idx(m1,n1);

D = zeros(N,ceil(mn/2));

notB = zeros(m,n);

% initialise some variables
sumD = zeros(1,N); numD = sumD; comD = sumD; jidx = zeros(1,4);

for i = 1:N
dd = find(d == i);
numD(i) = length(dd);
D(i,1:numD(i)) = dd';
sumD(i) = sum(map(dd));
end

totD = sum(map(:));
avgD = totD/N;

diffsumD = avgD - sumD;

[dsd_max,i_mxD] = max(diffsumD);

%its_tol = ceil(0.1*mn);
its_tol = fix(log(mn)*10+1);

score = zeros(1,its_tol);

score(1) = 100*sum(abs(avgD - sumD))/2/totD;
%score(1) = myscore(map,d,N);

score_tol = max(0.50*score(1),0.00001*sum(sumD)/N); its = 2;

score_min = score(1); d_min = d;

while (score(its-1) > score_tol) && (its <= its_tol) && ((its <= 5) || ((its > 5) && (any(abs(diff(score(its-1:-1:its-5))) > 1e-6))))

%   fprintf('%.4f > %.4f > %.4f    %d < %d\n',score(its-1),score_min,score_tol,its,its_tol)

num_dsd = 0;

I_mxD = zeros(N,1);
J_mxD = zeros(N,1);
I_mnD = zeros(N,1);
J_mnD = zeros(N,1);
dsmxD = zeros(N,1);
dsmnD = zeros(N,1);
mapmxD = zeros(N,1);
mapmnD = zeros(N,1);
scoreD = zeros(N,1);

d1 = d;
notB1 = notB;
D1 = D;
sumD1 = sumD;
numD1 = numD;
diffsumD1 = diffsumD;

for j_mxD = 1:numD1(i_mxD) % each square in district i_mxD

jidx(1) = idx_up(D1(i_mxD,j_mxD)); % the square above
jidx(2) = idx_rt(D1(i_mxD,j_mxD)); % the square right
jidx(3) = idx_dn(D1(i_mxD,j_mxD)); % the square below
jidx(4) = idx_lt(D1(i_mxD,j_mxD)); % the square left

for k = 1:4

if jidx(k) && map(jidx(k)) && (d(jidx(k)) ~= i_mxD) && (numD(d(jidx(k))) > 1) && (notB(jidx(k)) == 0) % square within bounds, non-zero value, bordering & not single

if k == 1  % check if expansion divides any districts
if idx_lt(jidx(k)) && idx_rt(jidx(k)) && ((~idx_up(jidx(k))) || (d(idx_up(jidx(k))) ~= d(jidx(k)))) && (d(idx_lt(jidx(k))) == d(idx_rt(jidx(k))))
continue
end
if idx_up(jidx(k)) && idx_lt(jidx(k)) && idx_lu(jidx(k)) && (d(idx_up(jidx(k))) == d(jidx(k))) && (d(idx_lt(jidx(k))) == d(jidx(k))) && (d(idx_lu(jidx(k))) ~= d(jidx(k)))
continue
end
if idx_up(jidx(k)) && idx_rt(jidx(k)) && idx_ur(jidx(k)) && (d(idx_up(jidx(k))) == d(jidx(k))) && (d(idx_rt(jidx(k))) == d(jidx(k))) && (d(idx_ur(jidx(k))) ~= d(jidx(k)))
continue
end
end
if k == 2
if idx_up(jidx(k)) && idx_dn(jidx(k)) && ((~idx_rt(jidx(k))) || (d(idx_rt(jidx(k))) ~= d(jidx(k)))) && (d(idx_up(jidx(k))) == d(idx_dn(jidx(k))))
continue
end
if idx_rt(jidx(k)) && idx_up(jidx(k)) && idx_ur(jidx(k)) && (d(idx_rt(jidx(k))) == d(jidx(k))) && (d(idx_up(jidx(k))) == d(jidx(k))) && (d(idx_ur(jidx(k))) ~= d(jidx(k)))
continue
end
if idx_rt(jidx(k)) && idx_dn(jidx(k)) && idx_rd(jidx(k)) && (d(idx_rt(jidx(k))) == d(jidx(k))) && (d(idx_dn(jidx(k))) == d(jidx(k))) && (d(idx_rd(jidx(k))) ~= d(jidx(k)))
continue
end
end
if k == 3
if idx_rt(jidx(k)) && idx_lt(jidx(k)) && ((~idx_dn(jidx(k))) || (d(idx_dn(jidx(k))) ~= d(jidx(k)))) && (d(idx_rt(jidx(k))) == d(idx_lt(jidx(k))))
continue
end
if idx_dn(jidx(k)) && idx_rt(jidx(k)) && idx_rd(jidx(k)) && (d(idx_dn(jidx(k))) == d(jidx(k))) && (d(idx_rt(jidx(k))) == d(jidx(k))) && (d(idx_rd(jidx(k))) ~= d(jidx(k)))
continue
end
if idx_dn(jidx(k)) && idx_lt(jidx(k)) && idx_dl(jidx(k)) && (d(idx_dn(jidx(k))) == d(jidx(k))) && (d(idx_lt(jidx(k))) == d(jidx(k))) && (d(idx_dl(jidx(k))) ~= d(jidx(k)))
continue
end
end
if k == 4
if idx_dn(jidx(k)) && idx_up(jidx(k)) && ((~idx_lt(jidx(k))) || (d(idx_lt(jidx(k))) ~= d(jidx(k)))) && (d(idx_dn(jidx(k))) == d(idx_up(jidx(k))))
continue
end
if idx_lt(jidx(k)) && idx_dn(jidx(k)) && idx_dl(jidx(k)) && (d(idx_lt(jidx(k))) == d(jidx(k))) && (d(idx_dn(jidx(k))) == d(jidx(k))) && (d(idx_dl(jidx(k))) ~= d(jidx(k)))
continue
end
if idx_lt(jidx(k)) && idx_up(jidx(k)) && idx_lu(jidx(k)) && (d(idx_lt(jidx(k))) == d(jidx(k))) && (d(idx_up(jidx(k))) == d(jidx(k))) && (d(idx_lu(jidx(k))) ~= d(jidx(k)))
continue
end
end

i_mnD = d(jidx(k));
j_mnD = 1;
while (D(i_mnD,j_mnD) ~= jidx(k)), j_mnD = j_mnD + 1; end

num_dsd = num_dsd + 1;

I_mxD(num_dsd) = i_mxD;
J_mxD(num_dsd) = j_mxD;
I_mnD(num_dsd) = i_mnD;
J_mnD(num_dsd) = j_mnD;
dsmxD(num_dsd) = diffsumD(i_mxD);
dsmnD(num_dsd) = diffsumD(i_mnD);
mapmxD(num_dsd) = map(D(i_mxD,j_mxD));
mapmnD(num_dsd) = map(D(i_mnD,j_mnD));

d_temp = d;
d_temp(D(i_mnD,j_mnD)) = d1(D(i_mxD,j_mxD));

sumD_tmp = sumD;
sumD_tmp(i_mxD) = sumD_tmp(i_mxD) + mapmnD(num_dsd);
sumD_tmp(i_mnD) = sumD_tmp(i_mnD) - mapmxD(num_dsd);

scoreD(num_dsd) = 100*sum(abs(avgD - sumD_tmp))/2/totD;
%            myscore(map,d_temp,N);
end
end
end

if (num_dsd > 0)
%      [ans,k] = min(abs(dsmxD(1:num_dsd) - mapmnD(1:num_dsd)) ...
%                  + abs(dsmnD(1:num_dsd) + mapmxD(1:num_dsd)));

[ans,k] = min(scoreD(1:num_dsd));

d1(D(I_mnD(k),J_mnD(k))) = d1(D(I_mxD(k),J_mxD(k)));

notB1(D(I_mnD(k),J_mnD(k))) = notB1(D(I_mnD(k),J_mnD(k))) + 1;
if (notB1(D(I_mnD(k),J_mnD(k))) > N), notB1(D(I_mnD(k),J_mnD(k))) = 0; end

D1(I_mxD(k),numD1(I_mxD(k))+1) = D1(I_mnD(k),J_mnD(k));
D1(I_mnD(k),J_mnD(k)) = D1(I_mnD(k),numD1(I_mnD(k))); D1(I_mnD(k),numD1(I_mnD(k))) = 0;

sumD1(I_mxD(k)) = sumD1(I_mxD(k)) + mapmnD(k);
sumD1(I_mnD(k)) = sumD1(I_mnD(k)) - mapmxD(k);
numD1(I_mxD(k)) = numD1(I_mxD(k)) + 1;
numD1(I_mnD(k)) = numD1(I_mnD(k)) - 1;

diffsumD1(I_mxD(k)) = diffsumD1(I_mxD(k)) - mapmnD(k);
diffsumD1(I_mnD(k)) = diffsumD1(I_mnD(k)) + mapmxD(k);
end

score(its) = 100*sum(abs(avgD - sumD1))/2/totD;
%   score(its) = myscore(map,d1,N);

if score(its) < score_min
score_min = score(its);
d_min = d1;
end

d = d1;
notB = notB1;
D = D1;
sumD = sumD1;
numD = numD1;
diffsumD = diffsumD1;

%      [dsd_max,i_mxD] = max(diffsumD);

if (num_dsd > 0)
i_mxD = I_mnD(k);
else
i_mxD = ceil(rand(1)*N);
%      [dsd_max,i_mxD] = max(diffsumD);
end

%   else
%      i_mxD = i_mxD + 1;
%      if i_mxD > N, i_mxD = 1; end

%      i_mxD = ceil(rand(1)*N);

%      d = d1;
%      notB = notB1;
%      D = D1;
%      sumD = sumD1;
%      numD = numD1;
%      diffsumD = diffsumD1;

%      [dsd_max,i_mxD] = max(diffsumD);
%   end

its = its + 1;

%   if ~rem(its,10), plotmap(map,N,d1), end

end

d = d_min;

%myscore(map,d,N)
%plotmap(map,N,d)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 j = 1:ndivs
d((j-1)*k+1:j*k,:) = d1 + (j-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 j = 1:ndivs
jk=j*k;
d(:,jk-k+1:jk) = d1 + (j-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 i = 1:8
nmap = D{i}*0+n;
%    mj=~mj;
mj = rem(i-1,2);
for j=[0:5 7:np-2]

idx=j+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{i},n,len,ind{idx},nmap);

Score1 = myscore(D{i},d1,n);
if Score1<Score
Score=Score1;
d=d1;
dec=i;
if Score1 < limiet
dec = i;
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;
i = 1;

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

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

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

while i <= len,
while i <= len
for g = 0:nr-1
d(row+g,col) = i+g;
end
i = i + 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) = i+g;
end
i = i + 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);
i = 1;
d = zeros(s1,s2);
r = s1;
c = s2;

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

d(y:end,x:end) = reshape((1:r*c)+i-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 i = 1:l
Indices(d(i)) = i;
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);
i = 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+(i-1);
else
d(r:r+gy-1,c:c+gx-1) = p2+(i-1);
end
i = i + gx*gy;
c = c + dir;
end
c = c - dir*(gx-1)/gx;
if dir < 0
d(r:r+gy-1,c) = (i:i+gy-1)';
else
d(r:r+gy-1,c+gx-1) = (i:i+gy-1)';
end
i = i + gy;
dir = -dir;
r = r + gy;
end

dr = s1 - r + 1;
dc = gx*nx+1;
d(r:end,1:dc) = reshape(i:i+dr*dc-1,dr,dc);
d(r:end,2:2:dc) = d(end:-1:r,2:2:dc);
if dc < s2
c = dc + 1;
i = i + dr*dc;
dc = s2 - c + 1;
d(1:end,c:s2) = reshape(i:i+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 i = 1:l
Indices(d(i)) = i;
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 i = 1:l,
j = ind(i);
m = map(j);
if ji<=n-ij || w+w < m ,
if ij == md,
ij = ij + 1;
k = i;
break;
end
ij = ij + 1;
w = 1-m;
else
w = w - m;
end
ji = ji-1;
d(j) = ij;
end
w = 1;
for i = l:-1:k
j = ind(i);
m = map(j);
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(j) = ij;
end

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

s1=siz(1);s2=siz(2);
rand('state',101);
w = zeros(1,n);
for i = 1:s1*s2
w(d(i)) = w(d(i)) + map(i);
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 i = 1:numel(ma),
Total(in(i)) = Total(in(i)) + ma(i);
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 i=1:n
wi=w(i);
nzi=libnz(i);
x = xmesh(wi); y = ymesh(wi);
if x>1 && ~d(x-1,y),
nzi=nzi+1;
liberties(i,nzi)=wi-1;
end
if x<xn && ~d(x+1,y)
nzi=nzi+1;
liberties(i,nzi)=wi+1;
end
if y>1 && ~d(x,y-1)
nzi=nzi+1;
liberties(i,nzi)=wi-xn;
end
if y<yn && ~d(x,y+1)
nzi=nzi+1;
liberties(i,nzi)=wi+xn;
end
libnz(i)=nzi;
end

num=num-n;
while num
if all(groupsize==inf)
break
end
[smallest,i]=min(groupsize);
nzi = libnz(i);
notzero = (liberties(i,1:nzi));
test=~d(notzero);
if ~any(test) || ~nzi
groupsize(i)=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(i,1:nzi)=notzero;
if ~d(nzm)
x = xmesh(nzm); y = ymesh(nzm);
if x>1 && ~d(x-1,y)
nzi=nzi+1;
liberties(i,nzi)=nzm-1;
end
if x<xn && ~d(x+1,y)
nzi=nzi+1;
liberties(i,nzi)=nzm+1;
end
if y>1 && ~d(x,y-1)
nzi=nzi+1;
liberties(i,nzi)=nzm-xn;
end
if y<yn && ~d(x,y+1)
nzi=nzi+1;
liberties(i,nzi)=nzm+xn;
end
num = num - 1;
d(nzm) = i;
groupsize(i) = groupsize(i)+ tmp; %map(nzm);
liberties(i,move:nzi-1)=liberties(i,move+1:nzi);
libnz(i)=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;
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?
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
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
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);
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);	%?????????
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]; % 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 i = ndx(end:-1:1)

end
k = k + 1; % This square has been traded too many times
continue
end

d(ndx1(i)) = dist2(i);

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

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