# Diffing "May stinks too!" and "relax!!"

 Title: May stinks too! relax!! Author: MR Keenan Paulo Uribe Submitted: 2004-04-28 14:36:09 UTC 2004-04-28 14:56:51 UTC Status: Passed Passed Score: 7.2013 7.2166 Result: 0.197% 0.197% CPU Time: 57.6581 57.8333 Code: ```function d = solver(map,n) goal = sum(map(:))/n; map = map/goal; d = solver0(map,n,goal); d = trade(d,map,n); d = refine(map,n,d); 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 0 & length(skip) < num delta = abs(pops - mu); delta(skip) = 0; k = min(find(delta == max(delta))); [i,j] = find(B == k); [m,n] = size(B); I = [max(1,i-1); min(m,i+1); i; i]; J = [j; j; max(1,j-1); min(n,j+1)]; if pops(k) < mu f = find(pops(B(I+(J-1)*m)) > mu); else f = find(pops(B(I+(J-1)*m)) < mu); end if isempty(f) skip = [skip; k]; continue end I = I(f); J = J(f); e = mod(f-1,length(i))+1; i = i(e); j = j(e); [v,p] = sort(A(I+(J-1)*m)); I = I(p); J = J(p); p = mod(p-1,length(i))+1; i = i(p); j = j(p); if pops(k) < mu q = max(find(cumsum(v) < (mu - pops(k)))); if isempty(q) q = max(find(cumsum(v) < omega*(mu - pops(k)))); end else q = max(find(cumsum(v) < (pops(k) - mu))); if isempty(q) q = max(find(cumsum(v) < omega*(pops(k) - mu))); end end I = I(1:q); J = J(1:q); i = i(1:q); j = j(1:q); if pops(k) < mu save = B(I+(J-1)*m); B(I+(J-1)*m) = k; if ~isconnected(B) B(I+(J-1)*m) = save; skip = [skip; k]; continue end else B(i+(j-1)*m) = B(I+(J-1)*m); if ~isconnected(B) B(i+(j-1)*m) = k; skip = [skip; k]; continue end end pops = full(sparse(B(:),1,A(:))); score = scale*norm(pops-mean(pops),1); if score < S S = score; Best = B; cnt = 0; else cnt = cnt+1; end skip = []; 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 sc10 ) ) minscore=0; else minscore = n*abs(goal-round(goal))/goal + 1e-13; end len = s1*s2; limiet = n/2000; 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]; init_idx = zeros(22); for i = 1:8 nmap = D{i}*0+n; 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 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,scoreScale) 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 = 6000; 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; s1s2 = min(s1,s2); for kk=1:nl y1 = pop(1,kk); x1 = pop(2,kk); y2 = y1; x2 = x1; for i=1:s1s2 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 && i < s1s2 switch pop(3,kk) case 0 y1 = y1 + 1; if y1 == s1 y1 = 2; end case 1 y1 = y1 - 1; if y1 == 1 y1 = s1 - 1; end case 2 x1 = x1 + 1; if x1 == s2 x1 = 2; end case 3 x1 = x1 - 1; if x1 == 1 x1 = s2 - 1; end end y2 = y1; x2 = x1; else break; end end 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 > 1500 s = sum(abs(w-1)) / n; if s - bestscore == 0 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 x1 && ~d(x,y-1) nzi=nzi+1; liberties(i,nzi)=wi-xn; end if y1 && ~d(x-1,y) nzi=nzi+1; liberties(i,nzi)=nzm-1; end if x1 && ~d(x,y-1) nzi=nzi+1; liberties(i,nzi)=nzm-xn; end if y0 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)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)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)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```