Winner Paulo Uribe (dm1)

Finish 2004-04-28 09:00:00 UTC

relaxp1

by Paulo Uribe

Status: Failed
Results:

Based on: relaxp (diff)

Comments
Please login or create a profile.
Code
function d = solver(map,n)

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

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

if s>801 && s<1999
    d1 = crystallise6(map,n);
    d1 = trade(d1,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 ~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

% 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
                d(:,(j-1)*k+1:j*k) = d1 + (j-1)*n1;
            end
            return
        end
    end
end

% is this an all-integer map?
maps=map*goal;
allint=false;
for i=1:min(s1,10);
    for j=1:min(s2,10);
        if abs(maps(i,j)-round(maps(i,j)))>1e-12
            allint = true;
            break
        end
    end
    if allint
        break
    end
end
if ~allint
    minscore = n*abs(goal-round(goal))/goal + 1e-12;
else
    minscore = 0;
end

len = s1*s2;
limiet=n/570;
limiet = max(limiet,minscore);
np13 = 9;
np4 = 8;
np5 = 86;
np = np13 + np4 + np5;
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*2);

stop = 0;

pat3x3a = [1  8  9; ...
        2  7  6; ... 
        3  4  5]; 
pat3x3b = [1  2  9; ...
        4  3  8; ... 
        5  6  7]; 
pat4x4(:,:,1) = [ 1  8  9 16
         2  7 10 15
         3  6 11 14
         4  5 12 13];
pat4x4(:,:,2) = [ 1 14 15 16
         2 13 12 11
         3  6  7 10
         4  5  8  9];
pat4x4(:,:,3) = [ 1 14 15 16
         2 13 10  9
         3 12 11  8
         4  5  6  7];
pat4x4(:,:,4) = [ 1 12 13 16
         2 11 14 15
         3 10  9  8
         4  5  6  7];
pat4x4(:,:,5) = [ 1  4  5 16
         2  3  6 15
         9  8  7 14
        10 11 12 13];
pat4x4(:,:,6) = [ 1  2 15 16
         4  3 14 13
         5  8  9 12
         6  7 10 11];
pat4x4(:,:,7) = [ 1  2  3 16
         8  7  4 15
         9  6  5 14
        10 11 12 13];
pat4x4(:,:,8) = [ 1  2  3 16
         6  5  4 15
         7 10 11 14
         8  9 12 13];
pat5x5(:,:, 1) = [ 1 10 11 24 25
         2  9 12 23 22
         3  8 13 20 21
         4  7 14 19 18
         5  6 15 16 17];
pat5x5(:,:, 2) = [ 1 10 11 24 25
         2  9 12 23 22
         3  8 13 14 21
         4  7 16 15 20
         5  6 17 18 19];
pat5x5(:,:, 3) = [ 1 10 11 12 25
         2  9 14 13 24
         3  8 15 22 23
         4  7 16 21 20
         5  6 17 18 19];
pat5x5(:,:, 4) = [ 1 10 11 12 25
         2  9 14 13 24
         3  8 15 16 23
         4  7 18 17 22
         5  6 19 20 21];
pat5x5(:,:, 5) = [ 1 12 13 14 25
         2 11 10 15 24
         3  8  9 16 23
         4  7 18 17 22
         5  6 19 20 21];
pat5x5(:,:, 6) = [ 1 22 23 24 25
         2 21 20 19 18
         3  8  9 16 17
         4  7 10 15 14
         5  6 11 12 13];
pat5x5(:,:, 7) = [ 1 22 23 24 25
         2 21 20 19 18
         3  8  9 10 17
         4  7 12 11 16
         5  6 13 14 15];
pat5x5(:,:, 8) = [ 1 22 23 24 25
         2 21 18 17 16
         3 20 19 14 15
         4  7  8 13 12
         5  6  9 10 11];
pat5x5(:,:, 9) = [ 1 12 13 24 25
         2 11 14 23 22
         3 10 15 16 21
         4  9  8 17 20
         5  6  7 18 19];
pat5x5(:,:,10) = [ 1 12 13 14 25
         2 11 16 15 24
         3 10 17 18 23
         4  9  8 19 22
         5  6  7 20 21];
pat5x5(:,:,11) = [ 1 14 15 16 25
         2 13 12 17 24
         3 10 11 18 23
         4  9  8 19 22
         5  6  7 20 21];
pat5x5(:,:,12) = [ 1 22 23 24 25
         2 21 20 19 18
         3 10 11 12 17
         4  9  8 13 16
         5  6  7 14 15];
pat5x5(:,:,13) = [ 1 22 23 24 25
         2 21 20 13 12
         3 18 19 14 11
         4 17 16 15 10
         5  6  7  8  9];
pat5x5(:,:,14) = [ 1 22 23 24 25
         2 21 14 13 12
         3 20 15 16 11
         4 19 18 17 10
         5  6  7  8  9];
pat5x5(:,:,15) = [ 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];
pat5x5(:,:,16) = [ 1 22 23 24 25
         2 21 16 15 14
         3 20 17 12 13
         4 19 18 11 10
         5  6  7  8  9];
pat5x5(:,:,17) = [ 1 16 17 24 25
         2 15 18 23 22
         3 14 19 20 21
         4 13 12 11 10
         5  6  7  8  9];
pat5x5(:,:,18) = [ 1 16 17 18 25
         2 15 20 19 24
         3 14 21 22 23
         4 13 12 11 10
         5  6  7  8  9];
pat5x5(:,:,19) = [ 1 18 19 20 25
         2 17 16 21 24
         3 14 15 22 23
         4 13 12 11 10
         5  6  7  8  9];
pat5x5(:,:,20) = [ 1 22 23 24 25
         2 21 20 19 18
         3 14 15 16 17
         4 13 12 11 10
         5  6  7  8  9];
pat5x5(:,:,21) = [ 1  6  7 24 25
         2  5  8 23 22
         3  4  9 20 21
        12 11 10 19 18
        13 14 15 16 17];
pat5x5(:,:,22) = [ 1  6  7 24 25
         2  5  8 23 22
         3  4  9 10 21
        14 13 12 11 20
        15 16 17 18 19];
pat5x5(:,:,23) = [ 1  6  7  8 25
         2  5 10  9 24
         3  4 11 22 23
        14 13 12 21 20
        15 16 17 18 19];
pat5x5(:,:,24) = [ 1  6  7  8 25
         2  5 10  9 24
         3  4 11 12 23
        16 15 14 13 22
        17 18 19 20 21];
pat5x5(:,:,25) = [ 1 14 15 16 25
         2 13 12 17 24
         3  4 11 18 23
         6  5 10 19 22
         7  8  9 20 21];
pat5x5(:,:,26) = [ 1 22 23 24 25
         2 21 20 19 18
         3  4 11 12 17
         6  5 10 13 16
         7  8  9 14 15];
pat5x5(:,:,27) = [ 1 22 23 24 25
         2 21 20 17 16
         3  4 19 18 15
         6  5 10 11 14
         7  8  9 12 13];
pat5x5(:,:,28) = [ 1 20 21 22 25
         2 19 18 23 24
         3  4 17 16 15
         6  5 10 11 14
         7  8  9 12 13];
pat5x5(:,:,29) = [ 1 22 23 24 25
         2 21 20 15 14
         3  4 19 16 13
         6  5 18 17 12
         7  8  9 10 11];
pat5x5(:,:,30) = [ 1 20 21 22 25
         2 19 18 23 24
         3  4 17 14 13
         6  5 16 15 12
         7  8  9 10 11];
pat5x5(:,:,31) = [ 1 18 19 20 25
         2 17 16 21 24
         3  4 15 22 23
         6  5 14 13 12
         7  8  9 10 11];
pat5x5(:,:,32) = [ 1 22 23 24 25
         2 21 20 19 18
         3  4 15 16 17
         6  5 14 13 12
         7  8  9 10 11];
pat5x5(:,:,33) = [ 1  8  9 10 25
         2  7  6 11 24
         3  4  5 12 23
        16 15 14 13 22
        17 18 19 20 21];
pat5x5(:,:,34) = [ 1 22 23 24 25
         2 21 20 19 18
         3  4  5 16 17
         8  7  6 15 14
         9 10 11 12 13];
pat5x5(:,:,35) = [ 1 22 23 24 25
         2 21 20 19 18
         3  4  5  6 17
        10  9  8  7 16
        11 12 13 14 15];
pat5x5(:,:,36) = [ 1  4  5 24 25
         2  3  6 23 22
        11 10  7 20 21
        12  9  8 19 18
        13 14 15 16 17];
pat5x5(:,:,37) = [ 1  4  5 24 25
         2  3  6 23 22
         9  8  7 20 21
        10 13 14 19 18
        11 12 15 16 17];
pat5x5(:,:,38) = [ 1  4  5 24 25
         2  3  6 23 22
        13 12  7  8 21
        14 11 10  9 20
        15 16 17 18 19];
pat5x5(:,:,39) = [ 1  4  5  6 25
         2  3  8  7 24
        13 12  9 22 23
        14 11 10 21 20
        15 16 17 18 19];
pat5x5(:,:,40) = [ 1  4  5  6 25
         2  3  8  7 24
        11 10  9 22 23
        12 15 16 21 20
        13 14 17 18 19];
pat5x5(:,:,41) = [ 1  4  5  6 25
         2  3  8  7 24
        15 14  9 10 23
        16 13 12 11 22
        17 18 19 20 21];
pat5x5(:,:,42) = [ 1  2 13 14 25
         4  3 12 15 24
         5 10 11 16 23
         6  9 18 17 22
         7  8 19 20 21];
pat5x5(:,:,43) = [ 1  2 23 24 25
         4  3 22 21 20
         5 10 11 18 19
         6  9 12 17 16
         7  8 13 14 15];
pat5x5(:,:,44) = [ 1  2 23 24 25
         4  3 22 21 20
         5 10 11 12 19
         6  9 14 13 18
         7  8 15 16 17];
pat5x5(:,:,45) = [ 1  2 15 16 25
         4  3 14 17 24
         5 12 13 18 23
         6 11 10 19 22
         7  8  9 20 21];
pat5x5(:,:,46) = [ 1  2 23 24 25
         4  3 22 21 20
         5 12 13 14 19
         6 11 10 15 18
         7  8  9 16 17];
pat5x5(:,:,47) = [ 1  2 23 24 25
         4  3 22 15 14
         5 20 21 16 13
         6 19 18 17 12
         7  8  9 10 11];
pat5x5(:,:,48) = [ 1  2 21 22 25
         4  3 20 23 24
         5 18 19 14 13
         6 17 16 15 12
         7  8  9 10 11];
pat5x5(:,:,49) = [ 1  2 19 20 25
         4  3 18 21 24
         5 16 17 22 23
         6 15 14 13 12
         7  8  9 10 11];
pat5x5(:,:,50) = [ 1  2 23 24 25
         4  3 22 21 20
         5 16 17 18 19
         6 15 14 13 12
         7  8  9 10 11];
pat5x5(:,:,51) = [ 1  2 15 16 25
         4  3 14 17 24
         5  6 13 18 23
         8  7 12 19 22
         9 10 11 20 21];
pat5x5(:,:,52) = [ 1  2 23 24 25
         4  3 22 21 20
         5  6 13 14 19
         8  7 12 15 18
         9 10 11 16 17];
pat5x5(:,:,53) = [ 1  2 23 24 25
         4  3 22 19 18
         5  6 21 20 17
         8  7 12 13 16
         9 10 11 14 15];
pat5x5(:,:,54) = [ 1  2 21 22 25
         4  3 20 23 24
         5  6 19 18 17
         8  7 12 13 16
         9 10 11 14 15];
pat5x5(:,:,55) = [ 1  2 23 24 25
         4  3 22 17 16
         5  6 21 18 15
         8  7 20 19 14
         9 10 11 12 13];
pat5x5(:,:,56) = [ 1  2 21 22 25
         4  3 20 23 24
         5  6 19 16 15
         8  7 18 17 14
         9 10 11 12 13];
pat5x5(:,:,57) = [ 1  2 19 20 25
         4  3 18 21 24
         5  6 17 22 23
         8  7 16 15 14
         9 10 11 12 13];
pat5x5(:,:,58) = [ 1  2 23 24 25
         4  3 22 21 20
         5  6 17 18 19
         8  7 16 15 14
         9 10 11 12 13];
pat5x5(:,:,59) = [ 1  2  9 10 25
         4  3  8 11 24
         5  6  7 12 23
        16 15 14 13 22
        17 18 19 20 21];
pat5x5(:,:,60) = [ 1  2 23 24 25
         4  3 22 21 20
         5  6  7 18 19
        10  9  8 17 16
        11 12 13 14 15];
pat5x5(:,:,61) = [ 1  2 23 24 25
         4  3 22 21 20
         5  6  7  8 19
        12 11 10  9 18
        13 14 15 16 17];
pat5x5(:,:,62) = [ 1  2  3 24 25
        10  9  4 23 22
        11  8  5 20 21
        12  7  6 19 18
        13 14 15 16 17];
pat5x5(:,:,63) = [ 1  2  3 24 25
         8  7  4 23 22
         9  6  5 20 21
        10 13 14 19 18
        11 12 15 16 17];
pat5x5(:,:,64) = [ 1  2  3 24 25
        12 11  4 23 22
        13 10  5  6 21
        14  9  8  7 20
        15 16 17 18 19];
pat5x5(:,:,65) = [ 1  2  3 24 25
         6  5  4 23 22
         7 12 13 20 21
         8 11 14 19 18
         9 10 15 16 17];
pat5x5(:,:,66) = [ 1  2  3 24 25
         6  5  4 23 22
         7 12 13 14 21
         8 11 16 15 20
         9 10 17 18 19];
pat5x5(:,:,67) = [ 1  2  3 24 25
         6  5  4 23 22
         7 14 15 16 21
         8 13 12 17 20
         9 10 11 18 19];
pat5x5(:,:,68) = [ 1  2  3 24 25
         6  5  4 23 22
         7 18 19 20 21
         8 17 16 15 14
         9 10 11 12 13];
pat5x5(:,:,69) = [ 1  2  3 24 25
         6  5  4 23 22
         7  8 15 16 21
        10  9 14 17 20
        11 12 13 18 19];
pat5x5(:,:,70) = [ 1  2  3 24 25
         6  5  4 23 22
         7  8 19 20 21
        10  9 18 17 16
        11 12 13 14 15];
pat5x5(:,:,71) = [ 1  2  3 24 25
         6  5  4 23 22
         7  8  9 20 21
        12 11 10 19 18
        13 14 15 16 17];
pat5x5(:,:,72) = [ 1  2  3 24 25
         6  5  4 23 22
         7  8  9 10 21
        14 13 12 11 20
        15 16 17 18 19];
pat5x5(:,:,73) = [ 1  2  3  4 25
        12 11 10  5 24
        13 14  9  6 23
        16 15  8  7 22
        17 18 19 20 21];
pat5x5(:,:,74) = [ 1  2  3  4 25
        14 13 12  5 24
        15 10 11  6 23
        16  9  8  7 22
        17 18 19 20 21];
pat5x5(:,:,75) = [ 1  2  3  4 25
        10  9  8  5 24
        11 12  7  6 23
        14 13 18 19 22
        15 16 17 20 21];
pat5x5(:,:,76) = [ 1  2  3  4 25
        12 11  6  5 24
        13 10  7 22 23
        14  9  8 21 20
        15 16 17 18 19];
pat5x5(:,:,77) = [ 1  2  3  4 25
        10  9  6  5 24
        11  8  7 22 23
        12 15 16 21 20
        13 14 17 18 19];
pat5x5(:,:,78) = [ 1  2  3  4 25
        14 13  6  5 24
        15 12  7  8 23
        16 11 10  9 22
        17 18 19 20 21];
pat5x5(:,:,79) = [ 1  2  3  4 25
         8  7  6  5 24
         9 14 15 22 23
        10 13 16 21 20
        11 12 17 18 19];
pat5x5(:,:,80) = [ 1  2  3  4 25
         8  7  6  5 24
         9 14 15 16 23
        10 13 18 17 22
        11 12 19 20 21];
pat5x5(:,:,81) = [ 1  2  3  4 25
         8  7  6  5 24
         9 16 17 18 23
        10 15 14 19 22
        11 12 13 20 21];
pat5x5(:,:,82) = [ 1  2  3  4 25
         8  7  6  5 24
         9 20 21 22 23
        10 19 18 17 16
        11 12 13 14 15];
pat5x5(:,:,83) = [ 1  2  3  4 25
         8  7  6  5 24
         9 10 17 18 23
        12 11 16 19 22
        13 14 15 20 21];
pat5x5(:,:,84) = [ 1  2  3  4 25
         8  7  6  5 24
         9 10 21 22 23
        12 11 20 19 18
        13 14 15 16 17];
pat5x5(:,:,85) = [ 1  2  3  4 25
         8  7  6  5 24
         9 10 11 22 23
        14 13 12 21 20
        15 16 17 18 19];
pat5x5(:,:,86) = [ 1  2  3  4 25
         8  7  6  5 24
         9 10 11 12 23
        16 15 14 13 22
        17 18 19 20 21];
nmap0 = { D{1}*0+n, D{2}*0+n };
mj = 1;
for i = 1:8
    mj = ~mj;    
    nmap = nmap0{real(mj)+1};
    for j=0:np-1
        idx=j+1;
        if ~mj
            sz = siz;
        else
            sz = siz2;
        end
        idx2 = idx + mj*np;
        if ( ~mj | s1 ~= s2 )
            if idx==1
                ind{idx2} = refind1(sz);
            elseif idx==2
                ind{idx2} = refind2(sz);
            elseif idx==3
                ind{idx2} = refind3(sz,2);
            elseif idx==4
                ind{idx2} = refind3(sz,3);
            elseif idx==5
                ind{idx2} = refind3(sz,4);
            elseif idx==6
                ind{idx2} = refind3(sz,5);         
            elseif idx==7
                ind{idx2}= refind4(sz);
            elseif idx==8
                ind{idx2} = refind5(sz,pat3x3a);
            elseif idx==9
                ind{idx2} = refind5(sz,pat3x3b);
            elseif idx > (np13) & idx <= (np13 + np4)
                ind{idx2} = refind5(sz,pat4x4(:,:,idx-np13));
            elseif idx > (np13 + np4)
                ind{idx2} = refind5(sz,pat5x5(:,:,idx-np13-np4));
            end            
        else
            ind{np+j+1} = ind{j+1};
        end
        
        d1 = pattern(D{i},n,len,ind{idx2},nmap);
        
        Score1 = myscore(D{i},d1,n);
        if Score1<Score
            Score=Score1;
            d=d1;
            dec=i;

        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 < 1e-2)
    d = mysolver(map, n, d, siz);
    s = myscore(map,d,n);
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)

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;
        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;

iadd=0;

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
    
    iadded=[];
    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( apop + map(halo(ih)) - idtarg);
                if testval < best
                    iadd = ih;
                    best = testval;
                end
            end
            i1=halo(iadd);
            if abs(dpop(idis)+map(i1)-targ) < abs(dpop(idis)-targ)   % will this point make things better?
                iadded = [iadded find(oldhalo==i1)];
                d(i1)=idis;     % add the pixel
                if iadd==1
                    halo=halo(2:nhalo);
                elseif iadd==nhalo
                    halo=halo(1:nhalo-1);
                else
                    halo=[halo(1:iadd-1) halo(iadd+1:nhalo)];
                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)
                iadd=iadded(ia);
                i1 = oldhalo(iadd);
                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
            iadded=[];
            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)
            iadd=iadded(ia);
            i1 = oldhalo(iadd);
            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
                        iadd=ih;
                        mmax=map(halo(ih));
                    end
                end
                i1=halo(iadd);
                d(i1)=idis;
                halo(iadd)=[];
                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
                        iadd=ih;
                        mmax=map(halo(ih));
                    end
                end
                i1=halo(iadd);
                d(i1)=idis;
                halo(iadd)=[];
                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

%unpad
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
traded = zeros(r,c);

% 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 is 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
    tradevalue = map(ndx1);
    tradeben = min(tradevalue,min(err1,abs(err2)));
    ndx = [tradevalue < max(err1,abs(err2))]; % No min benefit
    ndx1  = ndx1(ndx);
    if isempty(ndx1)
        notdone = 0;
        continue
    end
    dist1 = dist1(ndx);
    dist2 = dist2(ndx);
    err1  = err1(ndx);
    err2  = err2(ndx);
    tradevalue = tradevalue(ndx);
    tradeben = tradeben(ndx);
    
    % Make as many trades as possible
    [v,ndx] = sort(tradeben);
    k = 0;
    for i = ndx(end:-1:1) % Max benefit first
        
        if traded(ndx1(i)) % can only trade each square once ever
            k = k + 1; % This square has been traded too many times
            continue
        end
        
        % Does this trade hurt score and is trading still valid for dist2?
        if (tradeben(i) > min( abs( err([dist1(i),dist2(i)]) ) ))
            continue
        end
        
        % Consider this square traded
        traded(ndx1(i)) = 1;
        
        % Make the trade
        d(ndx1(i)) = dist2(i);
        
        % Check if districts are still contiguous
        if isNCont(d==dist1(i),[r c])
            d(ndx1(i)) = dist1(i); % Undo trade
            continue
        end
        
        % Update vars if needed
        if tradevalue(i) > 0
            err(dist1(i)) = err(dist1(i)) - tradevalue(i); % the giver
            err(dist2(i)) = err(dist2(i)) + tradevalue(i); % the taker
            if (err(dist1(i)) < 0) || (err(dist2(i)) > 0)
                break % These districts shouldn't be traded any more, so recompute stats
            end
        end
        
    end
    if k == length(ndx)
        notdone = 0; % All squares have been traded too many times so stop trading
    end
end