Finish 2008-05-07 12:00:00 UTC

tweaky bird and pushy cat 4

by srach

Status: Passed
Results: 133225.00 (cyc: 56, node: 8961)
CPU Time: 144.816
Score: 44558.1
Submitted at: 2008-05-07 11:54:55 UTC
Scored at: 2008-05-07 15:32:04 UTC

Current Rank: 3641st
Basis for: WIN (diff)

Comments
Please login or create a profile.
Code
function W = solv(B)

[W,S] = sunday(B,[1 2],[1 2],2100);


[W2,S2] = sunday(B',[1 2],[1 2],1100);
if S > S2
W = W2(:,[2 1 4 3]);
S=S2;
end

[nr,nc] = size(B);
[W2,S2] = sunday(fliplr(B),[1 2],[1 2],1100);
if S > S2
W = [W2(:,1) nc-W2(:,2)+1 W2(:,3) nc-W2(:,4)+1];
S=S2;
end

[W2,S2] = sunday(rot90(B,-1),[1 2],[1 2],1100);
if S > S2
W = [nr-W2(:,2)+1 W2(:,1) nr-W2(:,4)+1 W2(:,3)];
S=S2;
end

W=tune2(W,B);

return

W3=solver_A(B);
S3 = mygrade(B,W3);
if S3<S
    W=W3;
    S=S3;
end


end

function [W,S] = sunday(B0,x000,y000,th000)
[nR,nC]=size(B0);
Borig=nan(nR+2,nC+2);
Borig(2:end-1,2:end-1)=B0;
Bedit = Borig;

maxbridges = 4;

if size(Bedit,2) > 20
cutfirst = 4;
cutsecond = 8;
cutoff = 20;
else
cutfirst = 3;
cutsecond = 7;
cutoff = 12;
end

S = inf;
X1 = [1 2;2 1];
X2 = [1 3;3 1];
Y = [3 2 1;1 2 3];

% fprintf('\n');

for x = x000
    if x == 2
        [U BU] = phase1(Bedit,cutfirst,X1(x,:));
        [UU BUU] = phase3(BU,U,cutsecond,X2(x,:));
        [V BX] = phase3(BUU,UU,cutoff,X2(x,:));
    else
        [U BU] = phase1(Bedit,4,X1(x,:));
        [V BX] = phase3(BU,U,11,X2(x,:));
    end

    for y = y000
        %         if x == 1 && y == 2, continue, end
        %if x == 2 && y == 2 && S > th000, return, end
        W1 = phase2(Bedit,BX,V,maxbridges,cutoff,Y(y,:))-1;
        S1 = mygrade(B0,W1);
        %         fprintf('x=%d, y=%d, %d', x, y, S1);
        %         if isfinite(S), fprintf(' (%4d)\n', S1 - S), else,
        %         fprintf('\n'), end
        if S1 <= S
            S = S1;
            W = W1;
            maxbridges = maxbridges - 1;
        end
    end
end

if nR*nC > 400; return; end
 
br = sum(W(:,1)==W(:,3)&W(:,2)==W(:,4));

if br <= 4
    WSH = solverSH(B0);
    ssh = mygrade(B0,WSH);
    if ssh < S
        W = WSH;
        S = ssh;
    end
end
end

%%
function score = mygrade(B,W)
nR=size(B,1);
B(W(:,1)+(W(:,2)-1)*nR)=0;
B(W(:,3)+(W(:,4)-1)*nR)=0;
score=sum(B(:))+size(W,1)+sum(W(:,1)==W(:,3)&W(:,2)==W(:,4))*24;
end

%%
function [pincount k] = analyzeboard(B,rz)

% make a sorted list of all pins
pin = sort(B(B>0),'descend');
npins = size(pin,1);
if npins < 1
    pincount = [];
    k = 0;
    return
end

% pin, count, benefit
pincount = zeros(npins,3);
pincount(1,1) = pin(1);
k = 1;
count = 1;
for i = 2:npins
    if pin(i) == pincount(k,1)
        count = count + 1;
    else
        pincount(k,2) = count;
        k = k + 1;
        count = 1;
        pincount(k,1) = pin(i);
    end
end
pincount(k,2) = count;
pincount = pincount(1:k,:);
if rz < 3, return, end
% calculate the benefit of a path
for i = 1:k
    if pincount(i,2) >= 2
        p = pincount(i,1);
        [row col] = find(B == p);
        d = 0;
        N = size(row,1);
        for j = 2:N
            d = d + abs(row(j)-row(j-1)) + abs(col(j)-col(j-1));
        end
        pincount(i,3) = pincount(i,2)*p - 0.85 * d;
    end
end
end

%%
function [W B] = phase1(B,cutoff,rz)
% fprintf('-- phase 1 --\n');

W = [];

[pincount k] = analyzeboard(B,rz);

if k < 1
    return
end

pincount=sortrows(pincount,-rz);

for i = 1:k
    if pincount(i,2) >= 2
        p = pincount(i,1);
        [row col] = find(B == p);
        N = size(row,1);
        % find all pairwise distances
        Npairs = N*(N-1)/2;
        dist = zeros(Npairs,3);
        x = 0;
        for a = 1:N
            for b = (a+1):N
                x = x + 1;
                dist(x,1) = a;
                dist(x,2) = b;
                dist(x,3) = abs(row(a)-row(b)) + abs(col(a)-col(b));
            end
        end
        % sort by distance
        [d ix] = sort(dist(:,3));
        dist = dist(ix,:);

        % try to connect the closest pair possible
        npins = 0;
        for x = 1:Npairs
            if dist(x,3) > cutoff+1
                %                 fprintf('warning: dist = %2d\n', dist(x,3));
                break
            end
            a = dist(x,1);
            b = dist(x,2);
            %             path = simplepath([row(a); row(b)], [col(a); col(b)], -p);
            path = complexpath(B,[row(a); row(b)], [col(a); col(b)], -p, cutoff, 2*p);
            if size(path,1) > 0
                %                 fprintf('p = %d, connect a=%d, b=%d, dist = %d\n', p, a, b, dist(x,3));
                %                 fprintf('\nFound path, p = %d, r=%d c=%d, r=%d, c=%d\n', p, row(a), col(a), row(b), col(b));
                W = [W; path];

                B = addwirepath(B,path,-p);

                %                 pinlist = [row(a) col(a); row(b) col(b)];
                npins = 2;
                edit = [1:(a-1) (a+1):(b-1) (b+1):N];
                row = [row(a); row(b); row(edit)];
                col = [col(a); col(b); col(edit)];
                break
            else
                %                 fprintf('p = %d, cannot connect a=%d, b=%d, dist = %d\n', p, a, b, dist(x,3));
            end
        end

        if npins < 2
            continue
        end

        for j = 3:N
            % find all pins and wires
            [row2 col2] = find(B == -p);
            Npinwires = size(row2,1);
            %             fprintf('npins = %d, nwires = %d\n', npins, Npinwires - npins);

            % find all pairwise distances
            % a = already connected (pin or wire), b = not yet connected
            Npairs =  Npinwires * (N - npins);


            dist = zeros(Npairs,3);
            x = 0;
            for a = 1:Npinwires
                for b = (npins+1):N
                    x = x + 1;
                    dist(x,1) = a;
                    dist(x,2) = b;
                    dist(x,3) = abs(row2(a)-row(b)) + abs(col2(a)-col(b));
                end
            end
            % sort by distance
            [d ix] = sort(dist(:,3));
            dist = dist(ix,:);

            % try to connect closest pair possible
            connected = false;
            for x = 1:Npairs
                if dist(x,3) > cutoff+1
                    %                     fprintf('warning: dist = %2d\n', dist(x,3));
                    break
                end
                a = dist(x,1);
                b = dist(x,2);
                %                 path = simplepath([row2(a); row(b)], [col2(a); col(b)], -p);
                path = complexpath(B,[row2(a); row(b)], [col2(a); col(b)], -p, cutoff, p);
                if size(path,1) > 0
                    W = [W; path];

                    B = addwirepath(B,path,-p);

                    npins = npins + 1;
                    connected = true;
                    row([j b]) = row([b j]);
                    col([j b]) = col([b j]);
                    break
                end
            end

            if ~connected
                break
            end
        end

    end
end
end

%%
function [W B] = phase2(Borig,B,W,maxbridges,kappa,rz)

    function addbridgewirepath()
        for w = 1:size(path,1);
            if path(w,1) == path(w,3) % horizontal
                BH(path(w,1),path(w,2)) = false;
                BH(path(w,3),path(w,4)) = false;
                if path(w,2) == path(w,4)
                    B(path(w,1),path(w,2)) = -9999;
                end
            end
            if path(w,2) == path(w,4) % vertical
                BV(path(w,1),path(w,2)) = false;
                BV(path(w,3),path(w,4)) = false;
            end
        end
    end

% fprintf('-- phase 2 --\n');

[BH BV] = buildbridges(Borig,B,W);

[pincount k] = analyzeboard(B,rz);
if k < 1, return, end

pincount=sortrows(pincount,-rz);

for i = 1:k
    p = pincount(i,1);

    Npinwires = sum(B == -p);

    if Npinwires == 0
        %         fprintf('p = %d, no previous pinwires\n', p);

        if pincount(i,2) >= 2
            [row col] = find(B == p);
            N = size(row,1);
            % find all pairwise distances
            Npairs = N*(N-1)*.5;
            dist = zeros(Npairs,3);
            x = 0;
            for a = 1:N
                for b = (a+1):N
                    x = x + 1;
                    dist(x,1) = a;
                    dist(x,2) = b;
                    dist(x,3) = abs(row(a)-row(b)) + abs(col(a)-col(b));
                end
            end
            % sort by distance
            [d ix] = sort(dist(:,3));
            dist = dist(ix,:);

            maxstep = min((maxbridges*25)+kappa,2*p+1);

            % try to connect the closest pair possible
            connected = false;
            for x = 1:Npairs
                if dist(x,3) > maxstep+1
                    %                     fprintf('warning: dist = %2d\n', dist(x,3));
                    break
                end
                a = dist(x,1);
                b = dist(x,2);
                path = bridgepath(B,BH,BV,[row(a); row(b)], [col(a); col(b)], -p, maxbridges, kappa, ceil(1.85*p));
                if size(path,1) > 0
                    %                     fprintf('\nBRIDGE path, p = %d, r=%d c=%d, r=%d, c=%d\n', p, row(a), col(a), row(b), col(b));
                    W = [W; path];

                    B = addwirepath(B,path,-p);
                    addbridgewirepath();

                    connected = true;
                    break
                end
            end
            if ~connected
                continue
            end
        end
    end

    [row col] = find(B == p);
    Npins = size(row,1);

    maxstep = min((maxbridges*25)+kappa,p+1);

    for j = 1:Npins
        [row2 col2] = find(B == -p);
        Npinwires = size(row2,1);

        % find all pairwise distances
        % a = already connected (pin or wire), b = not yet connected
        Npairs =  Npinwires * (Npins-j+1);
        dist = zeros(Npairs,3);
        x = 0;
        for a = 1:Npinwires
            for b = j:Npins
                x = x + 1;
                dist(x,1) = a;
                dist(x,2) = b;
                dist(x,3) = abs(row2(a)-row(b)) + abs(col2(a)-col(b));
            end
        end
        % sort by distance
        [d ix] = sort(dist(:,3));
        dist = dist(ix,:);

        % try to connect closest pair possible
        connected = false;
        for x = 1:Npairs
            if dist(x,3) > maxstep+1
                %                     fprintf('warning: dist = %2d\n', dist(x,3));
                break
            end
            a = dist(x,1);
            b = dist(x,2);
            path = bridgepath(B,BH,BV,[row2(a); row(b)], [col2(a); col(b)], -p, maxbridges, kappa, p);
            if size(path,1) > 0
                %                      fprintf('\nEXTRA BRIDGE path, p = %d, r=%d c=%d, r=%d, c=%d\n', p, row2(a), col2(a), row(b), col(b));
                W = [W; path];

                B = addwirepath(B,path,-p);
                addbridgewirepath();

                connected = true;

                row([j b]) = row([b j]);
                col([j b]) = col([b j]);
                break
            end
        end

        if ~connected
            break
        end
    end

end

end

%%
function [W B] = phase3(B,W,kappa,rz)

% fprintf('-- phase 3 --\n');

[pincount k] = analyzeboard(B,rz);
if k < 1, return, end

pincount=sortrows(pincount,-rz);

for i = 1:k
    p = pincount(i,1);

    Npinwires = sum(B == -p);

    if Npinwires == 0
        %         fprintf('p = %d, no previous pinwires\n', p);

        if pincount(i,2) >= 2
            [row col] = find(B == p);
            N = size(row,1);
            % find all pairwise distances
            Npairs = N*(N-1)*.5;
            dist = zeros(Npairs,3);
            x = 0;
            for a = 1:N
                for b = (a+1):N
                    x = x + 1;
                    dist(x,1) = a;
                    dist(x,2) = b;
                    dist(x,3) = abs(row(a)-row(b)) + abs(col(a)-col(b));
                end
            end
            % sort by distance
            [d ix] = sort(dist(:,3));
            dist = dist(ix,:);

            maxstep = min(kappa,2*p+1);

            % try to connect the closest pair possible
            connected = false;
            for x = 1:Npairs
                if dist(x,3) > maxstep+1
                    %                     fprintf('warning: dist = %2d\n', dist(x,3));
                    break
                end
                a = dist(x,1);
                b = dist(x,2);
                %                 path = bridgepath(B,BH,BV,[row(a); row(b)], [col(a); col(b)], -p, maxbridges, kappa, 2*p);
                path = complexpath(B,[row(a); row(b)], [col(a); col(b)], -p, kappa, 2*p);
                if size(path,1) > 0
                    %                     fprintf('\nBRIDGE path, p = %d, r=%d c=%d, r=%d, c=%d\n', p, row(a), col(a), row(b), col(b));
                    W = [W; path];

                    B = addwirepath(B,path,-p);
                    %                     [B BH BV] = addbridgewirepath(B,BH,BV,path,-p);

                    connected = true;
                    break
                end
            end
            if ~connected
                continue
            end
        end
    end

    [row col] = find(B == p);
    Npins = size(row,1);

    maxstep = min(kappa,p+1);

    for j = 1:Npins
        [row2 col2] = find(B == -p);
        Npinwires = size(row2,1);

        % find all pairwise distances
        % a = already connected (pin or wire), b = not yet connected
        Npairs =  Npinwires * (Npins-j+1);
        dist = zeros(Npairs,3);
        x = 0;
        for a = 1:Npinwires
            for b = j:Npins
                x = x + 1;
                dist(x,1) = a;
                dist(x,2) = b;
                dist(x,3) = abs(row2(a)-row(b)) + abs(col2(a)-col(b));
            end
        end
        % sort by distance
        [d ix] = sort(dist(:,3));
        dist = dist(ix,:);

        % try to connect closest pair possible
        connected = false;
        for x = 1:Npairs
            if dist(x,3) > maxstep+1
                %                     fprintf('warning: dist = %2d\n', dist(x,3));
                break
            end
            a = dist(x,1);
            b = dist(x,2);
            %             path = bridgepath(B,BH,BV,[row2(a); row(b)], [col2(a); col(b)], -p, maxbridges, kappa, p);
            path = complexpath(B,[row2(a); row(b)], [col2(a); col(b)], -p, kappa, 2*p);
            if size(path,1) > 0
                %                      fprintf('\nEXTRA BRIDGE path, p = %d, r=%d c=%d, r=%d, c=%d\n', p, row2(a), col2(a), row(b), col(b));
                W = [W; path];

                B = addwirepath(B,path,-p);
                %                 [B BH BV] = addbridgewirepath(B,BH,BV,path,-p);

                connected = true;
                row([j b]) = row([b j]);
                col([j b]) = col([b j]);
                break
            end
        end

        if ~connected
            break
        end
    end

end

end

%%
function [BH BV] = buildbridges(Borig,B,path)
[NR NC] = size(B);
BH = Borig == 0;
BV = BH;
BH(:,[1 NC]) = false;
BV([1 NR],:) = false;
for i = 1:size(path,1)
    if path(i,1) == path(i,3) % horizontal
        BH(path(i,1),path(i,2)) = false;
        BH(path(i,3),path(i,4)) = false;
    end
    if path(i,2) == path(i,4) % vertical
        BV(path(i,1),path(i,2)) = false;
        BV(path(i,3),path(i,4)) = false;
    end
end
end

%%
function B = addwirepath(B,path,label)
B(path(1,1),path(1,2)) = label;
for i = 1:size(path,1);
    B(path(i,3),path(i,4)) = label;
end
end

%%
function path = traceback(z,PZ,NR,t)
path = zeros(t,4);
dr = mod(z,NR);
dc = ceil(z/NR);
for j = 1:t
    path(j,1:2) = [dr dc];
    z = PZ(z);
    dr = mod(z,NR);
    dc = ceil(z/NR);
    path(j,3:4) = [dr dc];
end
end


function path = complexpath(B,row,col,label,cutoff,maxpathlen)

% - complexpath -
[NR NC] = size(B);
PZ = zeros(NR,NC);
C = -ones(NR,NC);
C(row(2),col(2)) = 0; % source

% tag the targets
C(row(1),col(1)) = -2;
C( B == label ) = -2;

znext = zeros(NR*NC,1);
znext(1) = row(2) + (col(2)-1)*NR;
count = 1;
dZ = [-1 1 -NR NR];

[ign, ir] = sort(rand(1,4));
for step = 0:min(cutoff,maxpathlen)
    if count < 1, break, end
    N = count;
    z = znext;
    count = 0;
    for i = 1:N
        zi = z(i);
        for s=1:4
            Z = zi + dZ(ir(s));
            tag = C(Z);
            if tag == -2
                PZ(Z) = zi;
                path = traceback(Z,PZ,NR,step+1);
                return
            elseif tag == -1 && B(Z) == 0
                C(Z) = step+1;
                PZ(Z) = zi;
                count = count + 1; 
                znext(count) = Z;
            end
        end
    end
end
path = [];
end

%%
    function path = traceback2(Z,step,BRIDGE,zi,NR,row,col,PZ)
        PZ(Z) = zi;
        dr=mod(Z,NR);
        dc=ceil(Z/NR);
        path = zeros(step,4);
        j = 0;
        while dr ~= row || dc ~= col
            j = j + 1;
            path(j,1:2) = [dr dc];
            Z = dr + (dc-1)*NR;
            pz = PZ(Z);
            pr=mod(pz,NR);
            pc=ceil(pz/NR);
            path(j,3:4) = [pr pc];
            dr = pr;
            dc = pc;
            if BRIDGE(dr,dc)
                j = j + 1;
                path(j,:) = [dr dc dr dc];
            end
        end
        path = path(1:j,:);
    end

%%
function path = bridgepath(B,BH,BV,row,col,label,maxbridges,kappa,maxpathlen)

% - bridgepath -
% fprintf('bridgepath ...\n');
[NR NC] = size(B);
BRIDGE = false(NR,NC);
PZ = zeros(NR,NC);
C = -ones(NR,NC);
% C(row(2),col(2)) = 0; % source

% tag targets
C(row(1),col(1)) = -2;
C( B == label ) = -2;

maxstep = min((maxbridges*25)+kappa,maxpathlen+1);
nextstep = zeros(maxstep+26,1);
nextstep(1) = row(2) + (col(2)-1)*NR;
dZ = [-NR NR -1 1];

for step = 1:maxstep
    while nextstep(step)>0
        zi = nextstep(step);
        nextstep(step)=C(zi);
        for s = 1:4
            Z = zi + dZ(s);
            tag = C(Z);
            if tag==-2
                path = traceback2(Z,step,BRIDGE,zi,NR,row(2),col(2),PZ);
                return
            end
            if tag==-1
                if B(Z) == 0
                    C(Z) = nextstep(step+1);
                    nextstep(step+1) = Z;
                    PZ(Z) = zi;
                elseif (BH(Z)&&(s<3)) || (BV(Z)&&(s>2))
                    C(Z) = nextstep(step+26);
                    nextstep(step+26) = Z;
                    PZ(Z) = zi;
                    BRIDGE(Z) = true;
                end
            end
        end
    end
end
path = [];
end

%%

function w = solverSH(b)

% get pin numbers
p = unique(b);
p(1) = []; % discard zero

% count number of each pin
n = zeros(size(p));
for i = 1:length(n)
    n(i) = nnz(p(i) == b(:));
end

% ignore single pins since they can't be connected to anything
for i = 1:length(n)
    if n(i) == 1
        b(p(i) == b(:)) = -1;
    end
end

% pad board so I don't have to deal with out of bounds indexing
g = zeros(size(b)+2); % board to keep track of groups
bb = repmat(-1,size(g));
bb(2:end-1,2:end-1) = b; % full board

% loop to choose move and do it
w = [];
[rr, cc] = find(bb>0); % pins I want to connect
d = (size(bb,1)/2+0.5 - rr).^2 + (size(bb,2)/2+0.5 - cc).^2;
[d, order] = sort(d); % start in center of board and work out
order = order';
for k = 1:length(rr)-1
    bestscore = 0;
    minsteps = 32;
    for i = order
        if g(rr(i), cc(i))
            continue % this pin is already in a group, so don't try to connect
        end
        % find best route from pin to pin's group
        [score, mv, steps] = findBestMove(bb, g, rr(i), cc(i), minsteps);
        if score > bestscore
            bestscore = score;
            bestmove = mv;
            minsteps = steps;
            if minsteps < 3
                break
            end
        end
    end
    if bestscore == 0 % can't make any more connections (try to add bridges?)
        w = w - 1; % offset for padding
        return
    end
    g = doMove(g, bestmove, bb(bestmove(1,1), bestmove(1,2)));
    bb = doMove(bb, bestmove, bb(bestmove(1,1), bestmove(1,2)));
    w = [w; bestmove];
end

w = w - 1; % offset for padding
end

%%
function [bestscore, bestmove, minsteps] = findBestMove(b, g, r, c, steplimit)

bestscore = 0;
bestmove = [];

j = [1 -1 0 0];
k = [0 0 1 -1];



if ~any(g(:)==b(r,c))
    g = b; % no groups for this pin yet, so all pins are valid
end

% start at (r,c) pin and step away one unit at a time looking for this pin's group
bb = b;
bb(bb>0) = -1;
bb(r,c) = 1; % keep track of how many steps it takes to get to each position
minsteps = Inf;
p = b(r,c);
for i = 1:steplimit-1 % only allowed this many steps
    [rr, cc] = find(bb==i);
    for n = 1:length(rr)
        for m = 1:4
            thisr = rr(n) + j(m);
            thisc = cc(n) + k(m);

            if g(thisr,thisc) == p && ~(thisr == r && thisc == c)
                minsteps = i; % can link to group in this number of steps
                break
            end
            v = bb(thisr,thisc);
            if v == 0
                bb(thisr,thisc) = i+1;
            end
        end
        if minsteps < Inf
            break
        end
    end
    if minsteps < Inf
        break
    end
end
if minsteps == Inf
    return % can't reach any pins (good candidate for a bridge)
end

% score for this connection
bestscore = b(r,c) - minsteps;

% if only one step
if minsteps == 1
    bestmove = [r, c, thisr, thisc];
    return
end

% multiple steps - work backwards to define wire from group to pin
bestmove = zeros(minsteps,4);
for step = minsteps:-1:1
    for m = 1:4
        nextr = thisr + j(m);
        nextc = thisc + k(m);
        if bb(nextr, nextc) == step
            break
        end
    end
    bestmove(step,:) = [nextr, nextc, thisr, thisc];
    thisr = nextr;
    thisc = nextc;
end
end

%%
function b = doMove(b, mv, v)

% mark wires with the same number as the pins
for i = 1:size(mv,1)
    b(mv(i,1), mv(i,2)) = v;
end
b(mv(end,3), mv(end,4)) = v;

end

%***********************


function W=solver_A(B)

totEE=0;
W=zeros(0,4);
sB=size(B);
B0=zeros(sB);
idx0=find(B(:));
[b,I,J]=unique(B(idx0));

L0=zeros(size(b));
for n1=1:length(b),
	idx10{n1}=find(J==n1);
	L0(n1)=length(idx10{n1});
end
L0(L0==1)=0;
[nill,sortN]=sort(L0.*b);

for nsortN=length(sortN):-1:1+sum(L0==0),
	N=sortN(nsortN);
	L1=L0(N);
	idx1=idx10{N};

	if L1>1,
		R=cell(2*L1-1,2*L1-1);C=R;Q=R; % connecting paths
		D=zeros(2*L1-1,2*L1-1); D(:)=inf; % distances
		E=zeros(2*L1-1,2*L1-1); % gain
		IDX=zeros(2*L1-1,1); IDX(1:L1)=1; % valid clusters
		[r,c]=ind2sub(sB,idx0(idx1));
		for n1=1:L1, RR{n1}=r(n1); CC{n1}=c(n1); QQ{n1}=[]; end % clusters
		EE=zeros(2*L1-1,1); EE(1:L1)=b(N);
		WW=cell(2*L1-1,1);
		
		% first pass
		LL=2/L1/(L1-1);
		e=0; idx=[];
		for n1=1:L1, 
			for n2=n1+1:L1, 
				[R{n1,n2},C{n1,n2},Q{n1,n2},D(n1,n2)]=findpath(B,r(n1),c(n1),r(n2),c(n2)); 
				E(n1,n2)=EE(n1)+EE(n2)-D(n1,n2);
				if E(n1,n2)>e || (E(n1,n2)==e&&e>0&&rand<LL), e=E(n1,n2); idx=[n1,n2]; end
			end
		end
		K=L1; 
		if e>0, % something
			while ~isempty(idx),
				% new cluster
				K=K+1; 
				Rnow=R{idx(1),idx(2)}; Cnow=C{idx(1),idx(2)}; Qnow=Q{idx(1),idx(2)}; 
				RR{K}=cat(1,RR{idx(1)},Rnow((2:end-1)'),RR{idx(2)});
				CC{K}=cat(1,CC{idx(1)},Cnow((2:end-1)'),CC{idx(2)});
				QQ{K}=cat(1,QQ{idx(1)},Qnow,QQ{idx(2)});
				EE(K)=EE(idx(1))+EE(idx(2))-D(idx(1),idx(2));
				WW{K}=cat(1,WW{idx(1)},WW{idx(2)},[Rnow(1:end-1),Cnow(1:end-1),Rnow(2:end),Cnow(2:end)]);
				IDX(idx)=0;
				n1=K; N2=find(IDX);
				IDX(K)=1;
				if ~isempty(N2), 
					for n2=N2(:)', 
						if idx(1)<n2, mi1=idx(1);ma1=n2; else, mi1=n2;ma1=idx(1); end
						if idx(2)<n2, mi2=idx(2);ma2=n2; else, mi2=n2;ma2=idx(2); end
						if D(mi1,ma1)<D(mi2,ma2)
							R{n2,n1}=R{mi1,ma1};
							C{n2,n1}=C{mi1,ma1};
							Q{n2,n1}=Q{mi1,ma1};
							D(n2,n1)=D(mi1,ma1);
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1);
						else, 
							R{n2,n1}=R{mi2,ma2};
							C{n2,n1}=C{mi2,ma2};
							Q{n2,n1}=Q{mi2,ma2};
							D(n2,n1)=D(mi2,ma2);
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1);
						end
						if D(idx(1),idx(2))>1,
							% recompute distances
							LL=1/length(N2);
							[Rt,Ct,Qt,Dt]=findpath(B,Rnow((2:end-1)'),Cnow((2:end-1)'),RR{n2},CC{n2},1); 
							%disp([Rt,Ct]);
							if Dt<=D(n2,n1), R{n2,n1}=Rt;C{n2,n1}=Ct;Q{n2,n1}=Qt;D(n2,n1)=Dt; end
							E(n2,n1)=EE(n1)+EE(n2)-D(n2,n1);
						end
					end
					
					idx=[];
					IDX1=find(IDX);N1=length(IDX1);d=D(IDX1,IDX1);
					[sortd,idxd]=sort(d(:)); % closest path
					[r1,c1]=ind2sub([N1,N1],idxd(1:N1*(N1-1)/2));
					idxfirst=find(E(IDX1(r1)+size(E,1)*(IDX1(c1)-1))>max(EE(IDX1(r1)),EE(IDX1(c1))));  % first increasing gain
					if ~isempty(idxfirst), idx=[IDX1(r1(idxfirst(1))),IDX1(c1(idxfirst(1)))]; end
				else, idx=[]; end
			end
		end
		if K>L1,
			%for n1=1:K,EE(n1)=EE(n1)-25*sum(B0(RR{n1}+sB(1)*(CC{n1}-1))>0);end
			[maxEE,K]=max(EE(L1+1:end)); K=L1+K;
			if maxEE>0 && length(RR{K})>1,
				totEE=totEE+maxEE; 
				%disp(['On node ',num2str(b(N)),' score ',num2str(maxEE)]);
				W=cat(1,W,WW{K});
				for n3=1:length(RR{K}), 
					if ~B(RR{K}(n3),CC{K}(n3)), B(RR{K}(n3),CC{K}(n3))=-1; end; % paths
					if B0(RR{K}(n3),CC{K}(n3))>0, 
						W=cat(1,W,[RR{K}(n3),CC{K}(n3),RR{K}(n3),CC{K}(n3)]);
					end
				end
				for n3=1:size(QQ{K},1), 
					B(QQ{K}(n3,1),QQ{K}(n3,2))=1; % corners
				end
				for n3=1:length(RR{K}), 
					B0(RR{K}(n3),CC{K}(n3))=1; % paths
				end
			end
		end
	end
end
end

function [R,C,Q,L]=findpath(B,r1,c1,r2,c2,xtra);
% finds straight path from r1 c1 to r2 c2 in map B
R=[];C=[];Q=[];L=inf;
L1=length(r1); L2=length(r2);
A=abs(r1(:,ones(1,L2))-r2(:,ones(1,L1))')+abs(c1(:,ones(1,L2))-c2(:,ones(1,L1))');
[a,idx]=sort(A(:));
for n1=1:length(a),
	if L<=a(n1),return;end
	i2=ceil(idx(n1)/L1);
	i1=idx(n1)-(i2-1)*L1;
	% From B(c1(i1),r1(i1)) to B(c2(i2),r2(i2))
	if c1(i1)==c2(i2), % vertical
		C0=c1(i1);
		if r2(i2)<r1(i1), R0=(r1(i1):-1:r2(i2))';
		else, R0=(r1(i1):r2(i2))'; end
		L0=length(R0)-1;
		if L0<2||~any(B(R0(2:end-1),C0)), % clean straight
			R=R0;
			C=C0(ones(L0+1,1));
			Q=[r1(i1),c1(i1);r2(i2),c2(i2)];
			L=L0;
			return;
		end
	end
	if r1(i1)==r2(i2), % horizontal
		R0=r1(i1);
		if c2(i2)<c1(i1), C0=(c1(i1):-1:c2(i2))';
		else, C0=(c1(i1):c2(i2))'; end
		L0=length(C0)-1;
		if L0<2||~any(B(R0,C0(2:end-1))), % clean straight
			C=C0;
			R=R0(ones(L0+1,1));
			Q=[r1(i1),c1(i1);r2(i2),c2(i2)];
			L=L0;
			return;
		end
	end
	if c1(i1)<c2(i2), minc=c1(i1);maxc=c2(i2); idxc=1; else, minc=c2(i2);maxc=c1(i1); idxc=2; end
	b1=B(:,c1(i1)); b1(r1(i1))=0; 
	bb1a=cumsum(b1>0);
	bb1b=bb1a; bb1b(1)=bb1a(end);bb1b(2:end)=bb1a(end)-bb1a(1:end-1); 
	bb1=max(0,bb1a-bb1a(r1(i1)))+max(0,bb1b-bb1b(r1(i1))); 
	d1a=cumsum(b1~=0);
	d1b=d1a; d1b(1)=d1a(end);d1b(2:end)=d1a(end)-d1a(1:end-1); 
	d1=max(0,d1a-d1a(r1(i1)))+max(0,d1b-d1b(r1(i1)));
	b2=B(:,c2(i2)); b2(r2(i2))=0; 
	bb2a=cumsum(b2>0);
	bb2b=bb2a; bb2b(1)=bb2a(end);bb2b(2:end)=bb2a(end)-bb2a(1:end-1); 
	bb2=max(0,bb2a-bb2a(r2(i2)))+max(0,bb2b-bb2b(r2(i2))); 
	d2a=cumsum(b2~=0);
	d2b=d2a; d2b(1)=d2a(end);d2b(2:end)=d2a(end)-d2a(1:end-1); 
	d2=max(0,d2a-d2a(r2(i2)))+max(0,d2b-d2b(r2(i2)));
	validcorner1=find(b1==0&b2==0&bb1==0&bb2==0);
	if ~isempty(validcorner1),
		e=any(B(validcorner1,minc+1:maxc-1)>0,2);
		validcorner2=find(e==0);
		if ~isempty(validcorner2);
			f=sum(B(validcorner1(validcorner2),minc+1:maxc-1)~=0,2);
			validcorner=validcorner1(validcorner2);
			d=abs(validcorner-r1(i1))+abs(c2(i2)-c1(i1))+abs(validcorner-r2(i2))+25*(d1(validcorner)+d2(validcorner)+f);
			[md,id]=min(d);
			if md<L,
				if idxc==1, C=cat(1,c1(i1)+zeros(abs(validcorner(id)-r1(i1)),1),(c1(i1):c2(i2))',c2(i2)+zeros(abs(validcorner(id)-r2(i2)),1));
				else, 		C=cat(1,c1(i1)+zeros(abs(validcorner(id)-r1(i1)),1),(c1(i1):-1:c2(i2))',c2(i2)+zeros(abs(validcorner(id)-r2(i2)),1)); end
				if validcorner(id)>=r1(i1),
					if r2(i2)>=validcorner(id),R=cat(1,(r1(i1):validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)+1:r2(i2))');
					else, 					   R=cat(1,(r1(i1):validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)-1:-1:r2(i2))'); end
				else,
					if r2(i2)>=validcorner(id),R=cat(1,(r1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)+1:r2(i2))');
					else, 					   R=cat(1,(r1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(c2(i2)-c1(i1)),1),(validcorner(id)-1:-1:r2(i2))'); end
				end
				Q=[r1(i1),c1(i1);r2(i2),c2(i2);validcorner(id),minc; validcorner(id),maxc];
				L=md;
				if L<=a(n1),return;end
			end
		end
	end
	if r1(i1)<r2(i2), minr=r1(i1);maxr=r2(i2); idxr=1; else, minr=r2(i2);maxr=r1(i1); idxr=2; end
	b1=B(r1(i1),:)'; b1(c1(i1))=0; 
	bb1a=cumsum(b1>0);
	bb1b=bb1a; bb1b(1)=bb1a(end);bb1b(2:end)=bb1a(end)-bb1a(1:end-1); 
	bb1=max(0,bb1a-bb1a(c1(i1)))+max(0,bb1b-bb1b(c1(i1))); 
	d1a=cumsum(b1~=0);
	d1b=d1a; d1b(1)=d1a(end);d1b(2:end)=d1a(end)-d1a(1:end-1); 
	d1=max(0,d1a-d1a(c1(i1)))+max(0,d1b-d1b(c1(i1)));
	b2=B(r2(i2),:)'; b2(c2(i2))=0; 
	bb2a=cumsum(b2>0);
	bb2b=bb2a; bb2b(1)=bb2a(end);bb2b(2:end)=bb2a(end)-bb2a(1:end-1); 
	bb2=max(0,bb2a-bb2a(c2(i2)))+max(0,bb2b-bb2b(c2(i2))); 
	d2a=cumsum(b2~=0);
	d2b=d2a; d2b(1)=d2a(end);d2b(2:end)=d2a(end)-d2a(1:end-1); 
	d2=max(0,d2a-d2a(c2(i2)))+max(0,d2b-d2b(c2(i2)));
	validcorner1=find(b1==0&b2==0&bb1==0&bb2==0);
	if ~isempty(validcorner1),
		e=any(B(minr+1:maxr-1,validcorner1)>0,1)';
		validcorner2=find(e==0);
		if ~isempty(validcorner2);
			f=sum(B(minr+1:maxr-1,validcorner1(validcorner2))~=0,1)';
			validcorner=validcorner1(validcorner2);
			d=abs(validcorner-c1(i1))+abs(r2(i2)-r1(i1))+abs(validcorner-c2(i2))+25*(d1(validcorner)+d2(validcorner)+f);
			[md,id]=min(d);
			if md<L,
				if idxr==1, R=cat(1,r1(i1)+zeros(abs(validcorner(id)-c1(i1)),1),(r1(i1):r2(i2))',r2(i2)+zeros(abs(validcorner(id)-c2(i2)),1));
				else, 		R=cat(1,r1(i1)+zeros(abs(validcorner(id)-c1(i1)),1),(r1(i1):-1:r2(i2))',r2(i2)+zeros(abs(validcorner(id)-c2(i2)),1)); end
				if validcorner(id)>=c1(i1),
					if c2(i2)>=validcorner(id),C=cat(1,(c1(i1):validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)+1:c2(i2))');
					else, 					   C=cat(1,(c1(i1):validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)-1:-1:c2(i2))'); end
				else,
					if c2(i2)>=validcorner(id),C=cat(1,(c1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)+1:c2(i2))');
					else, 					   C=cat(1,(c1(i1):-1:validcorner(id))',validcorner(id)+zeros(abs(r2(i2)-r1(i1)),1),(validcorner(id)-1:-1:c2(i2))'); end
				end
				Q=[r1(i1),c1(i1);r2(i2),c2(i2);minr,validcorner(id); maxr,validcorner(id)];
				L=md;
				if L<=a(n1),return;end
			end
		end
	end
end

if isinf(L),C=[];R=[];Q=[]; end
return
end


function w=tune2(w,B);
[row,col] = size(B);
sw=size(w,1);
ww1=(w(:,2)-1)*row+w(:,1);
ww2=(w(:,4)-1)*row+w(:,3);
 
wl=find(ww1-ww2~=0);
ww=sort([ww1(wl);ww2(wl)]);
 
wcount = zeros(2*sw,2);
wcount(1,1) = ww(1);
k = 1;
count = 1;
for i = 2:size(ww),
    if ww(i) == wcount(k,1)
        count = count + 1;
    else
        wcount(k,2) = count;
        k = k + 1;
        count = 1;
        wcount(k,1) = ww(i);
    end
end
wcount(k,2) = count;
wcount = wcount(1:k,:);

mp=zeros(row,col);
mp(wcount(:,1))=wcount(:,2);

[gx,gy]=find(mp==3);
for gg=1:size(gx,1),
  r=gx(gg);
  c=gy(gg);
  
 % TR
 if c<col & r>1
    if mp(r-1,c)==2 & mp(r,c+1)==2 & mp(r-1,c+1)==1 & B(r,c+1)==0
      a=find(ww1==row*(c-1)+r & ww2==row*(c-1)+r-1);
      if size(a,1)==0, a=find(ww2==row*(c-1)+r & ww1==row*(c-1)+r-1); end
      if size(a,1)==1   
        b=[0;0];
        a=find(ww1==row*c+r-1 & ww2==row*c+r);
        if size(a,1)==0, a=find(ww2==row*c+r-1 & ww1==row*c+r); end
        if size(a,1)>0, b(1)=a; end
       
        a=find(ww1==row*(c-1)+r & ww2==row*c+r);
        if size(a,1)==0, a=find(ww2==row*(c-1)+r & ww1==row*c+r); end
        if size(a,1)>0, b(2)=a; end
       
        if all(b)
          w(b,:)=0;
          w(b(1),:)=[r-1 c r-1 c+1];
          co=find(w(:,1)~=0);
          w=w(co,:);
        end
      end  
    end
 end
 
 % TL
 if c>1 & r>1
    if mp(r-1,c)==2 & mp(r,c-1)==2 & mp(r-1,c-1)==1 & B(r,c-1)==0
      a=find(ww1==row*(c-1)+r & ww2==row*(c-1)+r-1);
      if size(a,1)==0, a=find(ww2==row*(c-1)+r & ww1==row*(c-1)+r-1); end
      if size(a,1)==1  
        b=[0;0];
        a=find(ww1==row*(c-2)+r-1 & ww2==row*(c-2)+r);
        if size(a,1)==0, a=find(ww2==row*(c-2)+r-1 & ww1==row*(c-2)+r); end
        if size(a,1)>0, b(1)=a; end
       
        a=find(ww1==row*(c-2)+r & ww2==row*(c-1)+r);
        if size(a,1)==0, a=find(ww2==row*(c-2)+r & ww1==row*(c-1)+r); end
        if size(a,1)>0, b(2)=a; end
       
        if all(b)
          w(b,:)=0;
          w(b(1),:)=[r-1 c-1 r-1 c];
          co=find(w(:,1)~=0);
          w=w(co,:);
        end
      end  
    end
 end

 % BL
  if c>1 & r<row
    if mp(r,c-1)==2 & mp(r+1,c)==2 & mp(r+1,c-1)==1 & B(r,c-1)==0
      a=find(ww1==row*(c-1)+r & ww2==row*(c-1)+r+1);
      if size(a,1)==0, a=find(ww2==row*(c-1)+r & ww1==row*(c-1)+r+1); end
      if size(a,1)==1
        b=[0;0];
        a=find(ww1==row*(c-2)+r & ww2==row*(c-2)+r+1);
        if size(a,1)==0, a=find(ww2==row*(c-2)+r & ww1==row*(c-2)+r+1); end
        if size(a,1)>0, b(1)=a; end
       
        a=find(ww1==row*(c-2)+r & ww2==row*(c-1)+r);
        if size(a,1)==0, a=find(ww2==row*(c-2)+r & ww1==row*(c-1)+r); end
        if size(a,1)>0, b(2)=a; end
       
        if all(b)
          w(b,:)=0;
          w(b(1),:)=[r+1 c-1 r+1 c];
          co=find(w(:,1)~=0);
          w=w(co,:);
        end
      end  
    end
  end
 
 % BR
  if c<col & r<row
    if mp(r,c+1)==2 & mp(r+1,c)==2 & mp(r+1,c+1)==1 & B(r,c+1)==0
      a=find(ww1==row*(c-1)+r & ww2==row*(c-1)+r+1);
      if size(a,1)==0, a=find(ww2==row*(c-1)+r & ww1==row*(c-1)+r+1); end
      if size(a,1)==1
        b=[0;0];
        a=find(ww1==row*(c-1)+r & ww2==row*c+r);
        if size(a,1)==0, a=find(ww2==row*(c-1)+r & ww1==row*c+r); end
        if size(a,1)>0, b(1)=a; end
       
        a=find(ww1==row*c+r & ww2==row*c+r+1);
        if size(a,1)==0, a=find(ww2==row*c+r & ww1==row*c+r+1); end
        if size(a,1)>0, b(2)=a; end
       
        if all(b)
          w(b,:)=0;
          w(b(1),:)=[r+1 c r+1 c+1];
          co=find(w(:,1)~=0);
          w=w(co,:);
        end
      end   
    end
  end 
end
end