Interstate 9

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

nt = size(tiles,1);
side = ceil(sqrt(nt));
if (side<min(boardSize))
    board = zeros(boardSize);
    [board(1:side,1:side),orientation] = internal(tiles,[side side]);
else
    [board,orientation] = internal(tiles,boardSize);
end;
end

function [board,orientation] = internal(tiles,boardSize)
% set up routes
nb = prod(boardSize);
nrow = boardSize(1);
ncol = boardSize(2);
route{1} = reshape(1:nb,nrow,ncol);  % down first
r = reshape(1:nb,ncol,nrow)';
[~,route{2}] = sort(r(:));  % across first
[~,route{3}] = sort(rand(1,nb));  % random
sp = spiralPath(nrow,ncol);
[~,route{4}] = sort(-sp(:));  % spiral in
[~,route{5}] = sort(sp(:));  % spiral out
dp = diagPath(nrow,ncol);
[~,route{6}] = sort(-dp(:));  % diagonal
dp2 = diagPath2(nrow,ncol);
[~,route{7}] = sort(-dp2(:));  % diagonal
sqp = squarePath(nrow,ncol);
[~,route{8}] = sort(sqp(:));  % square
sqp2 = squarePath2(nrow,ncol);
[~,route{9}] = sort(sqp2(:));  % square


nr = numel(route);
B = cell(1,nr);
O = cell(1,nr);
E = zeros(1,nr);
for i = 1:nr
    [B{i},O{i},E(i)] = routesolve(tiles,boardSize,route{i});
    if (E(i)==0)
        board = B{i};
        orientation = O{i};
        return;
    end;
end;
[~,i] = min(E);
board = B{i};
orientation = O{i};
end

function [board,orientation,err] = routesolve(tiles,boardSize,route)

nt = size(tiles,1);
board = zeros(boardSize);
nb = numel(board);
nrow = boardSize(1);
ncol = boardSize(2);
orientation = ones(nt,1);

unused = true(1,nt);
rt{4} = circshift(tiles,[0 1]);
rt{3} = circshift(tiles,[0 2]);
rt{2} = circshift(tiles,[0 3]);
rt{1} = tiles;
goalN = nan(boardSize);
goalN(1,:) = 0;
goalS = nan(boardSize);
goalS(end,:) = 0;
goalE = nan(boardSize);
goalE(:,end) = 0;
goalW = nan(boardSize);
goalW(:,1) = 0;
tsum = sum(tiles,2);
err = sum(tsum);

% greedily pick tiles
for ii = 1:min(nb,nt)
    i = route(ii);
    goal = [goalN(i),goalE(i),goalS(i),goalW(i)];
    for j = 1:4
        hd = bsxfun(@minus,rt{j}(unused,:),goal);
        hd(isnan(hd)) = 0;
        hs{j} = sum(abs(hd),2);
    end;
    hsm = cat(2,hs{:});
    hsm(hsm~=min(hsm(:))) = inf;
    hsm = bsxfun(@minus,hsm,tsum(unused));
    [~,pick] = min(hsm(:));
    [pick,ori] = ind2sub([nt-ii+1,4],pick);
    pick = find(unused,pick);
    pick = pick(end);
    board(i) = pick;
    orientation(pick) = ori;
    unused(pick) = false;
    
    % update goals
    if (i+nrow <= nb)
        goalW(i+nrow) = rt{ori}(pick,2);
    end;
    if (mod(i,nrow) ~= 0)
        goalN(i+1) = rt{ori}(pick,3);
    end;
    if (i-nrow >= 1)
        goalE(i-nrow) = rt{ori}(pick,4);
    end;
    if (mod(i-1,nrow) ~= 0)
        goalS(i-1) = rt{ori}(pick,1);
    end;
    
    % update err
    play = rt{ori}(pick,:);
    set = ~isnan(goal);
    err = err+sum(abs(play(set)-goal(set)))-sum(play(set))-sum(goal(set));
end;
end


function fitloc = spiralPath(nrow,ncol)
ntot = nrow*ncol;
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
pos = reshape(1:ntot,nrow,ncol);
nextfit = 1;
steps = 0;
radius = ceil((nrow+ncol)/12);
tops = inf;
avail = true(nrow,ncol);
sqrad = max(abs(xg-ncol/2),abs(yg-nrow/2));
%th = mod(atan2(yg-nrow/2,xg-ncol/2)+pi/2-4*pi./(sqrad+2).^2,2*pi);
th = atan2(yg-nrow/2,xg-ncol/2);
%[~,fitloc] = sort(sqrad(:)+th(:)/8);
fitn = 0;
fitloc = zeros(nrow,ncol);
for i = reshape(unique(sqrad),1,[]);
    xmp = find(sqrad==i,1);
    ring1 = find((sqrad==i)&(th>th(xmp)));
    [~,r] = sort(th(ring1));
    fitloc(ring1(r)) = fitn+1:fitn+numel(ring1);
    fitn = fitn+numel(ring1);
    ring2 = find((sqrad==i)&(th<=th(xmp)));
    [~,r] = sort(th(ring2));
    fitloc(ring2(r)) = fitn+1:fitn+numel(ring2);
    fitn = fitn+numel(ring2);
end;
%[~,sp] = sort(fitloc(:));
end


function dp = diagPath(nrow,ncol)
ntot = nrow*ncol;
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
dp = (xg+yg)+xg./ntot.^2;
end

function dp = diagPath2(nrow,ncol)
ntot = nrow*ncol;
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
dp = (xg-yg)+xg./ntot.^2;
end

function sqp = squarePath(nrow,ncol)
ntot = nrow*ncol;
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
sqp = max(xg,yg)+xg./ntot.^2+yg./ntot.^3;
end

function sqp = squarePath2(nrow,ncol)
ntot = nrow*ncol;
xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
sqp = max(ncol-xg+1,yg)+xg./ntot.^2+yg./ntot.^3;
end