Finish 2012-04-11 12:00:00 UTC

The Last Monkey

by Fel

Status: Passed
Results: 108982 (cyc: 10, node: 5749)
CPU Time: 65.49
Score: 11017.0
Submitted at: 2012-04-11 15:58:01 UTC
Scored at: 2012-04-11 18:27:34 UTC

Current Rank: 343rd (Highest: 337th )
Based on: Defense (diff)

Comments
Please login or create a profile.
Code
function [board, orientation] = solver(tiles, boardSize)
% First initialize the board
[board orientation] = curr_best(tiles, boardSize);

% Now optimize until convergence
[board orientation] = checkerboard_outer(tiles, board, orientation);
end

% Quickly compute the score of a board
function score = compute_score(tiles, board, vertical, horizontal)
% Compute the board score
score = sum(abs(diff(vertical(:,:)))) + sum(abs(diff(horizontal(:,:))));

% Add on the score of the remaining tiles
ind = 1:size(tiles, 1);
ind(board(board~=0)) = [];
score = score + sum(reshape(tiles(ind,:), [], 1));
end

% Compute the edges of the tiles in the board
function [vertical horizontal] = compute_edges(tiles, board, orient)
% Put all the tiles in the correct orientation and location
sz = size(board);
M = board(:) == 0;
if any(M)
    orient(end+1) = 1;
    tiles(end+1,:) = 0;
    board(M) = numel(orient);
end
tiles = tiles(board(:),:);
orient = orient(board(:));
tiles = reshape(rotate_tiles(tiles, orient), [sz 4]);

% Compute the vertical array
vertical = zeros(2, sz(1), sz(2)+1); 
vertical(2,:,1:end-1) = tiles(:,:,4);
vertical(1,:,2:end) = tiles(:,:,2);

% Compute the horizontal array
horizontal = zeros(2, sz(1)+1, sz(2));
horizontal(2,1:end-1,:) = tiles(:,:,1);
horizontal(1,2:end,:) = tiles(:,:,3);
end

% Rotate tiles so that the first value is the north value
function tiles = rotate_tiles(tiles, orient)
shift = [0 1 2 3; 1 2 3 0; 2 3 0 1; 3 0 1 2] * size(tiles, 1);
tiles = tiles(bsxfun(@plus, shift(orient,:), (1:size(tiles, 1))'));
end

% Optimize until convergence using the checkerboard update approach
function [B O score] = checkerboard_outer(tiles, B, O)
max_iters = floor(1e5 ./ (numel(B) * size(tiles, 1)));
if ~max_iters
    return
end
topleft = true;
last_score = 1e100;
iters = 0;
% Compute the edges of the tiles in the board
[vertical horizontal] = compute_edges(tiles, B, O);
% Compute the score
score = compute_score(tiles, B, vertical, horizontal);
% Loop until convergence
while score && score < last_score && iters < max_iters
    iters = iters + 1;
    last_score = score;
    % Do the checkerboard optimization
    [B O] = checkerboard_inner(tiles, B, O, topleft, vertical, horizontal);
    % Switch which squares are black
    topleft = ~topleft;
    
    % Compute the edges of the tiles in the board
    [vertical horizontal] = compute_edges(tiles, B, O);
    % Compute the score
    score = compute_score(tiles, B, vertical, horizontal);
    %assert(score <= last_score);
end
end

% Compute a globally optimal placement of tiles in all black squares of a
% checkerboard, given the tiles in all the white squares
function [B O] = checkerboard_inner(tiles, B, O, topleft, vertical, horizontal)
% Compute which tiles are fixed (i.e. white). topleft indicates whether the
% top left square is fixed.
sz = size(B);
fixed = true(sz);
fixed(topleft+1:2:end,1:2:end) = false;
fixed(~topleft+1:2:end,2:2:end) = false;
free = ~fixed;

% Extract the fixed board
fb = [reshape(horizontal(1,[free; false(1, sz(2))]), [], 1) ...
      reshape(vertical(2,[false(sz(1), 1) free]), [], 1) ...
      reshape(horizontal(2,[false(1, sz(2)); free]), [], 1) ...
      reshape(vertical(1,[free false(sz(1), 1)]), [], 1)];
fb = shiftdim(fb, -1);

% Extract the remaining tiles
fixed = B(fixed);
ind = 1:size(tiles, 1);
ind(fixed(fixed~=0)) = [];
rt = reshape(tiles(ind,:), [], 1, 4);

% Compute the scores of all remaining tiles in all locations and orientations
costMat = sum(abs(bsxfun(@minus, rt, fb)), 3);
costMat(:,:,2:4) = 0;
for a = 2:4
    rt = circshift(rt, [0 0 -1]);
    costMat(:,:,a) = sum(abs(bsxfun(@minus, rt, fb)), 3);
end

% Compute the best orientation and score for each tile in each location
[costMat orient] = min(costMat, [], 3);
costMat = bsxfun(@minus, bsxfun(@minus, costMat, sum(rt, 3)), sum(abs(fb), 3));
costMat(costMat >= 0) = Inf;

% Find the best placement of tiles
A = munkres(costMat);

% Compute the outputs
M = find(A ~= 0);
costMat = costMat(M,:)';
costMat = costMat(A(M)+(0:size(costMat, 1):numel(costMat)-1));
M = M(costMat < Inf);
A = A(M);
ind = ind(M);
free = find(free);
B(free) = 0;
B(free(A)) = ind;
orient = orient(M,:)';
O(ind) = orient(A+(0:size(orient, 1):numel(orient)-1));
end

function [assignment,cost] = munkres(costMat)
% MUNKRES   Munkres (Hungarian) Algorithm for Linear Assignment Problem. 
%
% [ASSIGN,COST] = munkres(COSTMAT) returns the optimal column indices,
% ASSIGN assigned to each row and the minimum COST based on the assignment
% problem represented by the COSTMAT, where the (i,j)th element represents the cost to assign the jth
% job to the ith worker.
%
% Partial assignment: This code can identify a partial assignment is a full
% assignment is not feasible. For a partial assignment, there are some
% zero elements in the returning assignment vector, which indicate
% un-assigned tasks. The cost returned only contains the cost of partially
% assigned tasks.

% This is vectorized implementation of the algorithm. It is the fastest
% among all Matlab implementations of the algorithm.

% Examples
% Example 1: a 5 x 5 example
%{
[assignment,cost] = munkres(magic(5));
disp(assignment); % 3 2 1 5 4
disp(cost); %15
%}
% Example 2: 400 x 400 random data
%{
n=400;
A=rand(n);
tic
[a,b]=munkres(A);
toc                 % about 2 seconds 
%}
% Example 3: rectangular assignment with inf costs
%{
A=rand(10,7);
A(A>0.7)=Inf;
[a,b]=munkres(A);
%}
% Example 4: an example of partial assignment
%{
A = [1 3 Inf; Inf Inf 5; Inf Inf 0.5]; 
[a,b]=munkres(A)
%}
% a = [1 0 3]
% b = 1.5
% Reference:
% "Munkres' Assignment Algorithm, Modified for Rectangular Matrices", 
% http://csclab.murraystate.edu/bob.pilgrim/445/munkres.html

% version 2.3 by Yi Cao at Cranfield University on 11th September 2011

assignment = zeros(1,size(costMat,1));
cost = 0;

validMat = costMat == costMat & costMat < Inf;
bigM = 10^(ceil(log10(sum(costMat(validMat))))+1);
costMat(~validMat) = bigM;

% costMat(costMat~=costMat)=Inf;
% validMat = costMat<Inf;
validCol = any(validMat,1);
validRow = any(validMat,2);

nRows = sum(validRow);
nCols = sum(validCol);
n     = max(nRows,nCols);
% if ~n
%     return
% end

maxv=10*max(costMat(validMat));

dMat = zeros(n) + maxv;
dMat(1:nRows,1:nCols) = costMat(validRow,validCol);

%*************************************************
% Munkres' Assignment Algorithm starts here
%*************************************************

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   STEP 1: Subtract the row minimum from each row.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
minR = min(dMat,[],2);
minC = min(bsxfun(@minus, dMat, minR));

%**************************************************************************  
%   STEP 2: Find a zero of dMat. If there are no starred zeros in its
%           column or row start the zero. Repeat for each zero
%**************************************************************************
zP = dMat == bsxfun(@plus, minC, minR);

starZ = zeros(n,1);
while any(zP(:))
    [r,c]=find(zP,1);
    starZ(r)=c;
    zP(r,:)=false;
    zP(:,c)=false;
end

while 1
%**************************************************************************
%   STEP 3: Cover each column with a starred zero. If all the columns are
%           covered then the matching is maximum
%**************************************************************************
    if all(starZ>0)
        break
    end
    uncoverColumn = true(n,1);
    uncoverColumn(starZ(starZ>0))=false;
    uncoverRow = true(n,1);
    primeZ = zeros(n,1);
    Mask = false(size(dMat));
    Mask(uncoverRow,uncoverColumn) = dMat(uncoverRow,uncoverColumn)==bsxfun(@plus,minR(uncoverRow),minC(uncoverColumn));
    while 1
        %**************************************************************************
        %   STEP 4: Find a noncovered zero and prime it.  If there is no starred
        %           zero in the row containing this primed zero, Go to Step 5.  
        %           Otherwise, cover this row and uncover the column containing 
        %           the starred zero. Continue in this manner until there are no 
        %           uncovered zeros left. Save the smallest uncovered value and 
        %           Go to Step 6.
        %**************************************************************************
        Step = 6;
        while 1
            [uZr uZc] = find(Mask, 1);
            if isempty(uZr)
                break;
            end
            primeZ(uZr) = uZc;
            stz = starZ(uZr);
            if ~stz
                Step = 5;
                break;
            end
            uncoverRow(uZr) = false;
            Mask(uZr,:) = false;
            Mask(uncoverRow,stz) = dMat(uncoverRow,stz) == minR(uncoverRow) + minC(stz);
            uncoverColumn(stz) = true;
        end
        if Step == 6
            % *************************************************************************
            % STEP 6: Add the minimum uncovered value to every element of each covered
            %         row, and subtract it from every element of each uncovered column.
            %         Return to Step 4 without altering any stars, primes, or covered lines.
            %**************************************************************************
            T = bsxfun(@minus, bsxfun(@minus, dMat(uncoverRow,uncoverColumn), minR(uncoverRow)), minC(uncoverColumn));
            minval = min(T(:));
            Mask(uncoverRow,uncoverColumn) = T==minval;
            minC(uncoverColumn) = minC(uncoverColumn) + minval;
            minR(~uncoverRow) = minR(~uncoverRow) - minval;
        else
            break
        end
    end
    %**************************************************************************
    % STEP 5:
    %  Construct a series of alternating primed and starred zeros as
    %  follows:
    %  Let Z0 represent the uncovered primed zero found in Step 4.
    %  Let Z1 denote the starred zero in the column of Z0 (if any).
    %  Let Z2 denote the primed zero in the row of Z1 (there will always
    %  be one).  Continue until the series terminates at a primed zero
    %  that has no starred zero in its column.  Unstar each starred
    %  zero of the series, star each primed zero of the series, erase
    %  all primes and uncover every line in the matrix.  Return to Step 3.
    %**************************************************************************
    rowZ1 = find(starZ==uZc);
    starZ(uZr)=uZc;
    while rowZ1>0
        starZ(rowZ1)=0;
        uZc = primeZ(rowZ1);
        uZr = rowZ1;
        rowZ1 = find(starZ==uZc);
        starZ(uZr)=uZc;
    end
end

% Cost of assignment
rowIdx = find(validRow);
colIdx = find(validCol);
starZ = starZ(1:nRows);
vIdx = starZ <= nCols;
assignment(rowIdx(vIdx)) = colIdx(starZ(vIdx));
pass = assignment(assignment>0);
pass(~diag(validMat(assignment>0,pass))) = 0;
assignment(assignment>0) = pass;
cost = trace(costMat(assignment>0,assignment(assignment>0)));
end

function [board, orientation] = curr_best(tiles, boardSize)

[board,orientation,solved,nrow,ncol,ntiles,nElem] = board_perfect_snake(boardSize,tiles);
if solved;return;end

cs=hankel(1:4,0:3);
played = tiles;
numUnplayed = ntiles-nElem;

sums = sum(tiles,2);
if numUnplayed > 0
    [sorted,si] = sort(sums);
	numUnplayedp1=numUnplayed+1;
    equals = find(sorted(numUnplayedp1)==sorted);
    if nnz(sorted(numUnplayedp1) == sorted(numUnplayedp1:end)) ~= length(equals)
        temp_si = si(equals);
        [~,mi] = sort(max(tiles(temp_si,:),[],2));
        si(equals(end:-1:1))= temp_si(mi);
    end
    played(si(1:numUnplayed),:) = Inf(numUnplayed,4);
    ntiles=nElem;
    height=nrow;
else
%    height=max(min(ceil(sqrt(ntiles)), nrow),ceil(ntiles/ncol));
    height=ceil(max(min(sqrt(ntiles), nrow),ntiles/ncol));
end

[y,x]=ind2sub([height,ncol],1:ntiles);
ind = sub2ind(boardSize,y,x);
board(ind) = -1;
last = x(end);

b =repmat(board,[1,1,18]);
o = ones(size(tiles,1), 18);

NZ0 = prod(played,2);

[b(:,:,1),o(:,1),list] = place([1 height],1:last,b(:,:,1),played,played,o(:,1),ntiles, cs,nrow,ncol,false,NZ0);
[b(:,:,1),o(:,1),list] = place(height:-1:1,[last 1],b(:,:,1),played,list,o(:,1),ntiles, cs,nrow,ncol,false,NZ0);
newList(:,:,1) = list;

[b(:,:,2),o(:,2),list] = place(height:-1:1,last,b(:,:,2),played,played,o(:,2),ntiles, cs,nrow,ncol,false,NZ0);
[b(:,:,2),o(:,2),list] = place(height,1:last,b(:,:,2),played,list,o(:,2),ntiles, cs,nrow,ncol,false,NZ0);
[b(:,:,2),o(:,2),list] = place(1:height-1,1,b(:,:,2),played,list,o(:,2),ntiles, cs,nrow,ncol,false,NZ0);
[b(:,:,2),o(:,2),list] = place(1,2:last,b(:,:,2),played,list,o(:,2),ntiles, cs,nrow,ncol,false,NZ0);
newList(:,:,2) = list;

for i=[2 5]
    c = floor(i/5) + 1;
    [b(:,:,i+1),o(:,i+1)] = place(height-1:-1:2,ncol-1:-1:2,b(:,:,c),played,newList(:,:,c),o(:,c),ntiles , cs,nrow,ncol,true,NZ0);
    [b(:,:,i+2),o(:,i+2)] = place(2:height-1,2:ncol-1,b(:,:,c),played,newList(:,:,c),o(:,c),ntiles, cs,nrow,ncol,true,NZ0);
    [b(:,:,i+3),o(:,i+3)] = place(2:height-1,ncol-1:-1:2,b(:,:,c),played,newList(:,:,c),o(:,c),ntiles, cs,nrow,ncol,true,NZ0);
end

[b(:,:,9),o(:,9)] = place(height-1:-1:2,2:ncol-1,b(:,:,2),played,newList(:,:,2),o(:,2),ntiles, cs,nrow,ncol,true,NZ0);

[b(:,:,1),o(:,1)] = place(height:-1:1,1:ncol,b(:,:,12),played,played,o(:,12),ntiles,cs,nrow,ncol,-1,NZ0);
[b(:,:,2 ), o(:,2 )] = calibrating_pure_werner(played, boardSize,nrow,ncol,true);
S = 9;

if nElem < 1000
    for ii = 1:6
	Sii=S+ii;
        [b(:,:,Sii),o(:,Sii)] = BiswasMichaelC(height,last, played, played, o(:,Sii), ntiles, cs,nrow,ncol);
    end
   
    S=Sii;
    if nElem < 200
        S = S +1;
        [b(:,:,S ), o(:,S )] = calibrating_pure_werner(tiles, boardSize,nrow,ncol,false);
     
    end
end

overall = zeros(S,1);
parfor i=1:S
    overall(i) = overallScore(b(:,:,i),nrow,ncol,o(:,i),tiles, cs);
end
[~,mo] = min(overall);
[board,orientation]=localopt(b(:,:,mo), o(:,mo),tiles,cs,nrow,ncol);
end

function [board,orientation,gain]=localopt(board, orientation,tiles,cs,r,c)
gain=0;
n = size(tiles,1);
didx = [3;4;1;2];
uidx = [1;2;3;4];
lidx = [4;1;2;3];
ridx = [2;3;4;1];
bp = board>0;
bt = board(bp);
te = [zeros(1,4);tiles];
oe = [1;orientation];
ub = [zeros(1,c);board(1:end-1,:)];
db = [board(2:end,:);zeros(1,c)];
lb = [zeros(r,1) board(:,1:end-1)];
rb = [board(:,2:end) zeros(r,1)];
ut = ub(bp)+1;
dt = db(bp)+1;
lt = lb(bp)+1;
rt = rb(bp)+1;
np1=n+1;
nv = [te(ut(:)+(didx(oe(ut))-1)*(np1)),...
      te(rt(:)+(lidx(oe(rt))-1)*(np1)),...
      te(dt(:)+(uidx(oe(dt))-1)*(np1)),...
      te(lt(:)+(ridx(oe(lt))-1)*(np1))];
bt=bt(:);
bv = tiles(bt(:,ones(1,4))+(cs(orientation(bt),:)-1)*n);
for i=1:2
pv = sum(abs(bv-nv),2);
[mpv,id]=max(pv);
while mpv
    bid=bt(id);
    bid1=bid+1;
    pv(id)=0;
    obid = orientation(bid);
    bx = tiles(bid,:);
    p0 = sum(abs(bx(cs)-nv(id+zeros(4,1),:)),2);
    [pmin,ido]=min(p0);
    gain=gain+(p0(obid)-pmin);
    if p0(obid)>pmin
        bv(id,:)=bx(cs(ido,:));
        orientation(bid)=ido;
        u1=ut==bid1;
        d1=dt==bid1;
        l1=lt==bid1;
        r1=rt==bid1;
        nv(u1,1) = tiles(bid,didx(ido));
        nv(d1,3) = tiles(bid,uidx(ido));
        nv(l1,4) = tiles(bid,ridx(ido));
        nv(r1,2) = tiles(bid,lidx(ido));
        L=u1|d1|l1|r1;
        pv(L) = sum(abs(bv(L,:)-nv(L,:)),2);
    end
    [mpv,id]=max(pv);
end
pv = sum(abs(bv-nv),2);
[mpv,id]=max(pv);
N=numel(pv);
tv = tiles(bt,:);
k=0;
while mpv>0
    bid = bt(id);
    bid1=bid+1;
    nvid = zeros(1,16);
    nvid(:) = nv(id+zeros(4,1),:);
    u1=ut==bid1;
    d1=dt==bid1;
    l1=lt==bid1;
    r1=rt==bid1;
    L=~(u1|d1|l1|r1);
    L(id)=false;
    qv = Inf(N,4);
    tmp=abs(tv(L,cs)-nvid(ones(sum(L),1),:));
    qv(L,:) = [tmp(:,1)+tmp(:,5)+tmp(:,9)+tmp(:,13),tmp(:,2)+tmp(:,6)+tmp(:,10)+tmp(:,14),tmp(:,3)+tmp(:,7)+tmp(:,11)+tmp(:,15),tmp(:,4)+tmp(:,8)+tmp(:,12)+tmp(:,16)];
    [iv,to] = min(qv(:));
    jd = mod(to-1,N)+1;
    bjd = bt(jd);
    bjd1=bjd+1;
    % objd = orientation(bjd);
    jo = ceil(to/N);
    nvjd = nv(jd+zeros(4,1),:);
    bvid = zeros(4);
    bvid(:) = tiles(bid,cs);
    [jv,io]=min(sum(abs(bvid - nvjd),2));
    jpv = sum(abs(bv(jd,:) - nv(jd,:)));
    if iv+jv < mpv+jpv
        bt(jd) = bid;
        bt(id) = bjd;
        orientation(bid) = io;
        orientation(bjd) = jo;
        board(bp)=bt;
        bv(id,:) = tv(jd,cs(jo,:));
        bv(jd,:) = tv(id,cs(io,:));
        bx = tv(id,:);
        tv(id,:) = tv(jd,:);
        tv(jd,:) = bx;
        nv(u1,1) = tiles(bjd,didx(jo));
        nv(d1,3) = tiles(bjd,uidx(jo));
        nv(l1,4) = tiles(bjd,ridx(jo));
        nv(r1,2) = tiles(bjd,lidx(jo));
        u2=ut==bjd1;
        d2=dt==bjd1;
        l2=lt==bjd1;
        r2=rt==bjd1;
        nv(u2,1) = tiles(bid,didx(io));
        nv(d2,3) = tiles(bid,uidx(io));
        nv(l2,4) = tiles(bid,ridx(io));
        nv(r2,2) = tiles(bid,lidx(io));
        L=u1|d1|l1|r1|u2|d2|l2|r2;
        pv(L) = sum(abs(bv(L,:)-nv(L,:)),2);
        uid = ut==bid1;
        ujd = ut==bjd1;
        ut(uid) = bjd1;
        ut(ujd) = bid1;
        did = dt==bid1;
        djd = dt==bjd1;
        dt(did) = bjd1;
        dt(djd) = bid1;
        lid = lt==bid1;
        ljd = lt==bjd1;
        lt(lid) = bjd1;
        lt(ljd) = bid1;
        rid = rt==bid1;
        rjd = rt==bjd1;
        rt(rid) = bjd1;
        rt(rjd) = bid1;
        k=k+1;
    end
    pv(id)=0;
    pv(jd)=0;
    [mpv,id]=max(pv);
end
end
end

function score = overallScore(board,nrow,ncol, orientation, tiles, cs)

tile = zeros(nrow*4,ncol);
tmp = tiles.';
for y = 1:nrow
    for x = 1:ncol
        if board(y,x)
            tile((y-1)*4+1:y*4,x)=tmp(cs(orientation(board(y,x)),:),board(y,x));
        end
    end
end

out_tiles = 1:size(tiles,1);
in_tiles = board(board>0);
out_tiles(in_tiles) = [];

score = sum(sum(abs(tile(5:4:end,1:end)-tile(3:4:end-4,1:end))))+ ...
    sum(sum(abs(tile(2:4:end,1:end-1)-tile(4:4:end,2:end))))+ ...
    sum(tile(4:4:end,1))+sum(tile(2:4:end,end))+sum(tile(end-1,1:end))+ ...
    sum(tile(1,1:end))+sum(sum(tiles(out_tiles,:)));
end

function [board, ort, lst] = BiswasMichaelC(height, width, tiles, lst, ort, ntiles, cs,nrow,ncol)
[cols,rows]=meshgrid(1:width,1:height);
center_x = width/2 + 0.5;
center_y = height/2 + 0.5;
dist_center = ((cols - center_x).^2 + (rows - center_y).^2).^(1/2.3);
dist_center = dist_center + rand(size(dist_center))*0.03; %param
[~, sort_index] = sort(dist_center(:));
loop_order = flipud(sort_index)';

board = zeros(nrow,ncol);
brd = zeros(height,width)-1;

if ntiles < height*width
    loop_order(~ismember(loop_order,1:ntiles)) = [];
    brd(ntiles+1:end) = 0;
end
score = inf(size(lst,1),1,4);
mst   = ~isinf(lst(:,1));
for ii = loop_order
    y = rows(ii);
    x = cols(ii);
    if (brd(y,x)==-1)
        [n,e,s,w] = findBoundaries(brd, tiles, ort, y, x, cs,height,width);
        
        nesw = cat(3, [n e s w], [w n e s], [s w n e],[e s w n]);
        score=score+inf;
        score(mst,:,:) = sum(abs(int32(bsxfun(@minus,lst(mst,:),nesw))),2);
        
        [ms, o] = find(score==min(score(:)));
        which = round(rand*(length(ms)-1)+1);
        t = ms(which);
        
        ort(t) = o(which);
        brd(y,x) = t;
%        lst(t,:) = inf(1,4);
        mst(t)   = false;
    end
    lst(~mst,:) = inf;
end
board(1:height,1:width) = brd;
end

function [brd, ort, lst] = place(rowRange, colRange, brd, tiles, lst, ort, ntiles, cs,nrow,ncol,flag,NZ)
k     = 1;
score = inf(size(lst,1),1,4);
mst   = ~isinf(lst(:,1));
for y=rowRange
    for x=colRange
        if (brd(y,x)==-1)
            [n,e,s,w] = findBoundaries(brd, tiles, ort, y, x, cs,nrow,ncol);
            nesw  = cat(3, [n e s w], [w n e s], [s w n e],[e s w n]);
            score = score+inf;
            score(mst,:,:) = sum(abs(int32(bsxfun(@minus,lst(mst,:),nesw))),2);
            
            [TT, O] = find(score==min(score(:)));
            %
            if flag==-1
                idx_t = ceil(length(TT)/1.5);%round(rand*(length(ms)-1)+1);
            elseif flag
                [~,idx_t] = max(NZ(TT));
            else
                [~,idx_t] = min(NZ(TT));
            end
            t = TT(idx_t);
            ort(t) = O(idx_t);
            brd(y,x)=t;
%            lst(t,:)= inf(1,4);
            mst(t)  = false;     
            k = k + 1;
            if (k > ntiles)
                return;
            end
        end
    end
end
lst(~mst,:)=inf;
end

function [north, east, south, west] = findBoundaries(board, tiles, orientation, row, col, cs,nrow,ncol)
north = boundary(board, tiles, orientation, row-1,col ,3, cs,nrow,ncol );
south = boundary(board, tiles, orientation, row+1,col ,1, cs,nrow,ncol );
west = boundary(board, tiles, orientation, row ,col-1,2, cs,nrow,ncol );
east = boundary(board, tiles, orientation, row ,col+1,4, cs,nrow,ncol );
end

function value = boundary(board, tiles, orientation, row,col,id, cs,nrow,ncol )

if row<1 || col<1 || row>nrow || col>ncol || board(row,col) == 0
    value = 0;
elseif board(row,col) == -1
    value = nan;
else
    value = tiles(board(row,col),cs(orientation(board(row,col)),id));
end

end


function [board, orientation] = calibrating_pure_werner(tiles, boardSize,nrows,ncols,flag)
board = zeros(boardSize);
ntiles = size(tiles, 1);
orientation = ones(ntiles, 1);
boarde = board;

% Kudos to Martin F. for quadrupling tiles
tiles = [tiles; tiles(:,[2:4,1]); tiles(:,[3:4,1:2]); tiles(:,[4,1:3])];

tilestaken = false(4*ntiles, 1);

if flag
boarde(nrows,1) = 4;
else
boarde(1,1) = 4;
end

[nextrow, nextcol,v]  = nextfield();
while v && ~all(tilestaken);     
    tix = findbesttilefor(nextrow, nextcol,tiles,board,nrows,ncols,tilestaken,flag);
    placetile(tix, nextrow, nextcol);
    [nextrow, nextcol,v]  = nextfield();
end

    function [nextrowv, nextcolv,v] = nextfield()
        [v, ix] = max(boarde(:));
          nextrowv=mod(ix-1,nrows)+1; nextcolv= ceil(ix/nrows);
    end
    function placetile(tileix, row, col)
        board(row, col) = tileix;
        boarde(row, col) = 0;
        rowm1=row > 1;
        colm1= col > 1;
        rowmnrows=row < nrows;
        colmncols=col < ncols;
        rowp1=row+1;
        colp1=col+1;
        coln1=col-1;
        rown1=row-1;
        
        method4a(tileix, row, col,rowm1,colm1,rowmnrows ,colmncols, colp1,coln1,rowp1,rown1);
        method4b( rowm1,colm1,colmncols, colp1,coln1, rown1);
        method4c( colm1,rowmnrows,colmncols, colp1,coln1,rowp1 );
          
        poptile(tileix);
                
        function method4a(tileix, row, col,rowm1,colm1,rowmnrows,colmncols, colp1,coln1,rowp1,rown1)
                     
            if rowm1 && ~board(rown1,col ); boarde(rown1, col ) = boarde(rown1, col ) + tiles(tileix, 1); end;
            if rowmnrows && ~board(rowp1,col ); boarde(rowp1, col ) = boarde(rowp1, col ) + tiles(tileix, 3); end;
            if colm1 && ~board(row ,coln1); boarde(row , coln1) = boarde(row , coln1) + tiles(tileix, 4); end;
            if colmncols && ~board(row ,colp1); boarde(row , colp1) = boarde(row , colp1) + tiles(tileix, 2); end;
        end
        function method4b( rowm1, colm1,colmncols, colp1,coln1, rown1)
            if rowm1 && colm1 && ~board(rown1,coln1); boarde(rown1, coln1) = boarde(rown1, coln1) + 0.1; end;
            if rowm1 && colmncols && ~board(rown1,colp1); boarde(rown1, colp1) = boarde(rown1, colp1) + 0.1; end;
        end
        function method4c( colm1,rowmnrows,colmncols, colp1,coln1,rowp1 )
            if rowmnrows && colm1 && ~board(rowp1,coln1); boarde(rowp1, coln1) = boarde(rowp1, coln1) + 0.1; end;
            if rowmnrows && colmncols && ~board(rowp1,colp1); boarde(rowp1, colp1) = boarde(rowp1, colp1) + 0.1; end;
        end
        
        function poptile(tileix )
            ix = mod(tileix-1, ntiles)+1;
            orientation(ix) = ceil(tileix/ntiles);
            tilestaken(ix+ntiles*(0:3), :) = true;
        end
    end

boardzeros = board == 0;
board = mod(board-1, ntiles)+1;
board(boardzeros) = 0;

end

function ret = findbesttilefor(r, c,tiles,board,nrows,ncols,tilestaken,flag)
  function myret=spycat(x,y,z,tiles,board)
        if board(x,y) > 0; myret = tiles(board(x,y), z);
        else myret = nan;
        end
    end 

v = zeros(1,4);
if r > 1;
    v(1) = spycat(r-1,c , 3,tiles,board);
end
if c < ncols;
    v(2)= spycat(r,c+1,4,tiles,board)  ;
end
if r < nrows;
    v(3) = spycat(r+1,c,1,tiles,board);
end
if c > 1;
    v(4)=spycat(r,c-1,2,tiles,board) ;
end

 mask = ~isnan(v);
    tilecosts = sum(abs(bsxfun(@minus, tiles(:,mask), v(mask))), 2);
    tilecosts(tilestaken) = inf;
if flag
    [~, ret] = min(tilecosts - 1e-10*sum(tiles(:, isnan(v)), 2));
else
    [~, ret] = min(tilecosts-sum(tiles,2).*1e-6);
end

end

% Richard routines
function [board,orientation,solved,nr,nc,numtiles,nrnc]=board_perfect_snake(boardSize,tiles)
% check for perfect unique board
solved=false;
board = zeros(boardSize,'single');
numtiles=size(tiles,1);

orientation = ones(numtiles, 1);

[nr nc]=size(board);
nrnc=nr*nc;
if numtiles~=nrnc,return;end

 

% minimum required #tiles with a zero
if size(find( tiles==0),1)<2*(nr+nc-2),return;end % No solution

otiles=tiles;

[tv,orientation,otiles,valid]= ...
    find_corner(tiles,otiles,orientation);
if ~valid,return;end
board(1)=tv; % 2 pts

for i=1:nr
    
    if i>1
        tile=board(i-1);
        
        [board,orientation,otiles,valid]= ...
            find_tiles(board,orientation,tiles,otiles,i,1,nrnc,tile,3);
        % [board,orientation,otiles,valid]= ...
        % first_column(board,orientation,otiles,tiles,i,j,nelem,tile,edge);
        if ~valid,return;end
        
    end % i>1
    
    for j=2:nc
        tile=board(i,j-1);
        
        [board,orientation,otiles,valid]= ...
            find_tiles(board,orientation,tiles,otiles,i,j,nrnc,tile,2);
        if ~valid,return;end
        
    end % j 2:nc
    
end % i 1:nr

solved=true;

end


function [tv,orientation,otiles,valid]= ...
    find_corner(tiles,otiles,orientation)

valid=false;
cs=[4 1 2 3]; % 4 pts define outside

for i=1:4
    adj= tiles(:,i) + tiles(:,cs(i) );
    tv=find(adj==0,1); % Corner found
    if ~isempty(tv)
        orientation(tv)=i;
        otiles(tv,:)=rot_tile(i,tiles(tv,:)); % or to i 2 pts
        valid=true;
        return; % corner found and rotated to position
    end
end

% No corner found
end % find corner

function [one_tile]=rot_tile(rot,tile)
one_tile=circshift(tile,[0 1-rot]);
end

function [board,orientation,otiles,valid]= ...
    find_tiles(board,orientation,tiles,otiles,i,j,nrnc,tile,ul)
% [board,orientation,otiles,valid]= ...
% first_column(board,orientation,otiles,tiles,i,j,nrnc,tile,UL);

valid=false; % could remove valid=false from returns

tvec=find(otiles==otiles(tile,ul));
if size(tvec,1)~=2,return;end % Non Unique solutions
tvecn1=tvec-1;
tnum=mod(tvecn1,nrnc)+1;
trot=floor((tvecn1)/nrnc)+1;


if tnum(1)~=tile % use
    tv=tnum(1);
    rot=trot(1);
else % use 2nd
    tv=tnum(2);
    rot=trot(2);
end

% for found pushed to L(4) 4:1,3:4,2:3,1:2
if j>1 % doing rows
    rot = mod(rot,4)+1;
end
board(i,j)=tv;
otiles(tv,:)=rot_tile(rot,tiles(tv,:));
orientation(tv)=rot;

valid=true;
end