function W = solver(B)

%results
W=zeros(1000,4);
%actual #connections
wctr=0;

%save original matrix
% b0=B;

[nr,nc] = size(B);
% grid size
sz = nr*nc;

%pin positions
ppos = find(B);

% values of each pin
pwt = B(ppos);

%stack to save search indices
stk=zeros(nr*nc,1);
%stack ptr
tos = 0;

%search grid;
bignum=9999;
G = ones(nr,nc)*bignum;
%set all pins as obstacles
%OBSTACLE=-2 ROUTED = -1 
G(ppos) = -2;


%unique pin gp, idx(discard, membership in pu of each pwt
[pv, puind, pidx] = unique(pwt);

%group all similar pins together
[psort, ipsort]=sort(pidx);

%position of each pinval
pp = ppos(ipsort);

% count for each pin gp
pc = accumarray(pidx, ones(length(pidx),1));

% weight of each pin gp
% [wt, wtord] = sort(pv.*pc,'descend');
[wt, wtord] = sort(pv,'descend');

%start and end for each pin gp
pe = cumsum(pc);
ps = pe - pc + 1;

np = length(pv);

%row and col lookups for all indices
r = repmat((1:nr).',1,nc);
r=r(:);
c = repmat(1:nc,nr,1);
c=c(:);

%loop through all gps
for ii = 1:np
    
    ord = wtord(ii);

    %pin positions for this gp
    pind = pp(ps(ord):pe(ord));
    prow = r(pind);
    pcol = c(pind);
    
    %list of pins visited
    pdone = zeros(pc(ord),1);
    
    %source pin (index in pind)
    src = 1;
    pdone(src)=1;
    
    %number of connections in this gp = #pins -1
    for jj=2:pc(ord)        
    
    if pc(ord) >2 
        %find closest pin
        %using manhattan dist from src
        mdist = abs(prow(~pdone) - prow(src)) + abs(pcol(~pdone) - pcol(src));
    
        [sd, isd] = sort(mdist);
        
        %indices of remaining pins
        v = find(~pdone);
        tgt = v(isd(1));
        
        pdone(tgt) = 1;
    
    else
         tgt = 2;
         pdone(2) = 1;
    end
    
%     fprintf ('src(%d):%d\ttgt(%d):%d\n',B(pind(src)),pind(src),B(pind(tgt)),pind(tgt));
    
    isrc = pind(src);
    itgt = pind(tgt);        
    
    %find path from src to tgt
    % expand
    tos = 0;
    curVal =0;
    
    %all pins are obstacles except the tgt
    G(ppos)=-2;
%     G(itgt)=bignum;
    G(pind)=bignum;
    G(isrc)=0;
    
    [plen,q,G] = expand(isrc,itgt,G);
    
    %push candidates into stack
    stk(1:length(q)) = q;
    tos = length(q);
    
    if (isequal(plen,0))
    
        
        while (tos)
        
            gp = stk(tos);
            tos = tos - 1;            
            [plen,q,G] = expand(gp,itgt,G);
            
            if (plen >0)
%                 tos=0;
                break;            
            else
                if ~isempty(q)
                    %push candidates into stack
                    stk((1:length(q))+tos) = q;
                    tos = tos+length(q);
                end
            end
                        
        end
        
        
    end
    
    % backtrace

    [p,G] = traceback(isrc,itgt,plen,G);
    
%      disp(p);
    %store connections
    if (~isempty(p))
        plen = size(p,1);
        W(wctr+(1:plen),:) = p;
        wctr=wctr+plen;
    end
    % clear
    G(G>-1) = bignum;
    
    src = tgt;

    end

end

W=W(1:wctr,:);

end


%index to row,col
function [r,c] = ind2rc(ind,nr,nc)

c = ceil(ind/nr);
r = ind - (c-1)*nr;

end %ind2rc

%row col to index
function ind =rc2ind(r,c,nr,nc)


ind=0;

if (r>0 && r<=nr &&c>0 && c<=nc)
    ind = (c-1)*nr + r;
end

end

%manhattan distance
function d = md(i1, i2)

[r1,c1]=ind2rc(i1);
[r2,c2]=ind2rc(i2);

d = abs(r1-r2)+abs(c1-c2);
end %md

%check valid grid pt
function ind2 = gpat(ind, nr ,nc)

[r1,c1]=ind2rc(ind, nr, nc);

ind2=ind;

if (r1<1 || r1>nr) 
    ind2=-1; 
end

if (c1<1 || c1>nc) 
    ind2=-1; 
end

end %gridPtAt

function dind = getdirs(ind,nr,nc)
%     ind = gpat(ind,nr,nc);
    dind=zeros(4,1);

    %is a valid point
    if (ind>0)
        [r,c] = ind2rc(ind,nr,nc);
        %west
        dind(1) = rc2ind(r,c-1,nr,nc);    
        %east
        dind(2) = rc2ind(r,c+1,nr,nc);
        %north
        dind(3) = rc2ind(r-1,c,nr,nc);
        %south
        dind(4) = rc2ind(r+1,c,nr,nc);
    end

end

%expansion phase
function [plen,q,G] = expand(gp,tgt,G)

    [nr,nc]= size(G);
        
    xp = getdirs(gp, nr, nc);
    val = G(gp);
    plen = 0;
    q=zeros(4,1);
    ctr=0;
    
    for ii=1:4
        if (xp(ii) && G(xp(ii))==9999)
            G(xp(ii)) =  val+1;
            if isequal(xp(ii),tgt)
                plen=G(xp(ii));
                break;
            else
                ctr=ctr+1;
                q(ctr)=xp(ii);
            end
        end
    end
    
    if (ctr)
        q=q(1:ctr);
    else
        q=[];
    end
end

%backtrace
function [p,G] = traceback (src,tgt,plen,G)

[nr,nc]=size(G);
p=zeros(plen,4);
cur = tgt;
ctr=1;

while ~isequal(cur,src)
    
    curval = G(cur);
    %set ROUTED
    G(cur) = -1;
    nb = getdirs(cur,nr,nc);
    next = 0;
    for ii=1:4
        if (nb(ii) && G(nb(ii))>-1 && G(nb(ii)) <curval)
            next = nb(ii);
            break;
        end
    end
    
    if (~next)
%         disp('Cant traceback!');
        p=[];
        break;
    end

    [r1,c1]=ind2rc(cur,nr,nc);
    [r2,c2]=ind2rc(next,nr,nc);
    p(ctr,:)=[r1,c1,r2,c2];
    ctr=ctr+1;
    cur=next;
end

if isequal(cur,src)
    G(cur) = -1;
    p=p(1:ctr-1,:);
end


end % traceback