Finish 2011-11-09 12:00:00 UTC

Twiddling my phone

by Michael Bindschadler

Status: Passed
Results: 69719198 (cyc: 14, node: 11515)
CPU Time: 76.867
Score: 697403.0
Submitted at: 2011-11-08 23:02:49 UTC
Scored at: 2011-11-08 23:04:22 UTC

Current Rank: 897th (Highest: 620th )
Based on: twiddling twiddling (diff)
Basis for: SOC FTW (hopefully) (diff)

Comments
Please login or create a profile.
Code
function [moves, vine] = solver(board, limit)

% based on "m"
% (http://www.mathworks.com/matlabcentral/contest/contests/34/submissions/64931)

%Allocate result cell
resultcell(:,3) = num2cell([-Inf*ones(10,1);0]);

%Get an index matrix for all those flipud/fliplr/rot90 cases --> not
%elegant at all, maybe incorporate this into the functions later
[nrow,ncol] = size(board);
ntot = nrow*ncol;
id = reshape(1:ntot,nrow,ncol);

%% Mike Add
limit = max(limit-100,0);
limit4soc = min(100,limit);

if limit4soc==0
    [moves,vine] = alfi_soc(board,limit);
    %disp('Winner is alfi_soc')
    return
end
%% END Mike add

limitlimit = 15100;
if limit < 15100
    [resultcell{1,1} resultcell{1,2} resultcell{1,3}] = alfi(board, limit);
    resultcell{1,4} = 'alfi';
end

largestValue = max(board(:));

if largestValue > 200 %|| limit>10000
    if (limit > 0)
        [resultcell{3,1} resultcell{3,2} resultcell{3,3}] = robert(board, limit);
        resultcell{3,4} = 'robert';
        [movestemp vinetemp resultcell{2,3}] = robert(fliplr(board),limit); % monotonous should work good on flipping lr/ud as well
        resultcell{2,4} = 'robert_fliplr';
        idtemp = fliplr(id);
        resultcell{2,1} = idtemp(movestemp);
        resultcell{2,2} = idtemp(vinetemp);
        [movestemp vinetemp resultcell{10,3}] = robert(flipud(board),limit);
        resultcell{10,4} = 'robert_flipud';
        idtemp = flipud(id);
        resultcell{10,1} = idtemp(movestemp);
        resultcell{10,2} = idtemp(vinetemp);
        [resultcell{9,1} resultcell{9,2} resultcell{9,3}] = dofilter(board,limit);
        resultcell{9,4} = 'dofilter';
        
    end
    
    doclusters=any(any(board(:,2:end)==board(:,1:end-1)))||any(any(board(1:end-1,:)==board(2:end,:))); %%%%
    if doclusters && limit<numel(board)
        [resultcell{4,2} resultcell{4,3}] = real_kurt(board);
        resultcell{4,4} = 'real_kurt';        
    end
    
    [resultcell{5,1} resultcell{5,2} resultcell{5,3}] = nick(board,limit);
    resultcell{5,4} = 'nick';
    [movestemp vinetemp resultcell{7,3}] = nick(rot90(board,2),limit);
    resultcell{7,4} = 'nick_rot90_2';
    idtemp = rot90(id,2);
    resultcell{7,1} = idtemp(movestemp);
    resultcell{7,2} = idtemp(vinetemp);
    if max(cell2mat(resultcell([5 7],3))) < 2 * max(cell2mat(resultcell([1 2 3 4 9 10],3)))
        [movestemp vinetemp resultcell{6,3}] = nick(board',limit);
        resultcell{6,4} = 'nick_transpose';
        idtemp = id';
        resultcell{6,1} = idtemp(movestemp);
        resultcell{6,2} = idtemp(vinetemp);
        [movestemp vinetemp resultcell{8,3}] = nick(rot90(board,3),limit);
        resultcell{8,4} = 'nick_rot90_3';
        idtemp = rot90(id,3);
        resultcell{8,1} = idtemp(movestemp);
        resultcell{8,2} = idtemp(vinetemp);
    end
end

%% Mike Add
siz = size(board);
ind_board = id;%reshape(1:prod(siz),siz);
north_neighbor = shiftnD(ind_board,1,0);
south_neighbor = shiftnD(ind_board,-1,0);
east_neighbor = shiftnD(ind_board,-2,0);
west_neighbor = shiftnD(ind_board,2,0);
neighbors = cat(3,north_neighbor,south_neighbor,east_neighbor,west_neighbor); %these are indices, not values
n_lut = reshape(neighbors,prod(siz),4);

    function board = do_moves(board,moves)
        % Perform moves on the board:
        for abc = 1:size(moves,1)
            board(moves(abc,:)) = [0 board(moves(abc,1))];
        end
    end
oldresultcell = resultcell;

alfi_soc_run_flag = false;    
for i = 1:size(resultcell,1)
    if ~isempty(resultcell{i,2})
        [new_moves,resultcell{i,2},success_flag] = convert2soc(resultcell{i,2},limit4soc,board,siz,n_lut);
        if success_flag
            resultcell{i,1} = [resultcell{i,1};new_moves];
            mod_board = do_moves(board,resultcell{i,1});
            resultcell{i,3} = sum(mod_board(resultcell{i,2}));
        else
            if ~alfi_soc_run_flag
                [alfi_soc_moves,alfi_soc_vine,alfi_soc_score] = alfi_soc(board,limit);
                alfi_soc_run_flag = true;
            end
            resultcell{i,1} = alfi_soc_moves;
            resultcell{i,2} = alfi_soc_vine;
            resultcell{i,3} = alfi_soc_score;
            resultcell{i,4} = 'alfi_soc';
        end
    end
end
%% End Mike Add

% [old_best,old_best_index] = max(cell2mat(oldresultcell(:,3)));
% disp(['Old winner was: ', oldresultcell{old_best_index,4}]);

[soc_best,best_index]=max(cell2mat(resultcell(:,3)));
% disp(['New Winner is: ', resultcell{best_index,4}])

moves = resultcell{best_index,1};
vine = resultcell{best_index,2};

if vine(1)~=1
    vine(1)=1; %SOC compliance insurance!
end

% if size(moves,1) < limit
%     [umm_moves, umm_vine] = use_more_moves(board, moves, vine, limit);
% 
%     if umm_vine(1)==1
%         moves = umm_moves;
%         vine = umm_vine;
%     end
% end
% 

end


%% nick
function [moves, vine, scr] = nick(board, limit)

moves = zeros(limit,2);
[nrow,ncol] = size(board);
ntot = nrow*ncol;

if (limit < nrow-1)
    scr = -inf;
    moves = [];
    vine = [];
    return;
end;

xg = repmat(1:ncol,nrow,1);
yg = repmat((1:nrow)',1,ncol);
pos = reshape(1:ntot,nrow,ncol);
fitloc = reshape(1:ntot,nrow,ncol);
fitloc(:,2:2:end) = flipud(fitloc(:,2:2:end));
%[s,r] = sort(board(:),'descend');
[s,r] = sort(-board(:));
s=-s;
rank = zeros(nrow,ncol);
rank(r) = 1:ntot;
remain = true(1,ntot);
nextfit = 1;
steps = 0;
cboard = board;
for i = 1:ntot
    if remain(i)
        [cy,cx] = find(rank==i);
        tx = xg(fitloc(nextfit));
        ty = yg(fitloc(nextfit));
        [cboard,moves,remain,rank,steps]=oftl2(cx,tx,cy,ty,cboard,rank,steps,moves,pos,remain,ncol,nrow,limit,s,i);
        if (steps > limit)
            break;
        end;
        rank(ty,tx) = i;
        nextfit = nextfit+1;
    end;
end;
moves = moves(1:min(steps,limit),:);
[vine,scr] = real_kurt(cboard);
end

% %% michael
% % http://www.mathworks.com/matlabcentral/contest/contests/34/submissions/64578
% function [moves, vine] = use_more_moves(board, moves, vine, limit)
% 
% m = size(board,1);
% 
% for i = 1:size(moves,1)
%     board(moves(i,:)) = [0 board(moves(i,1))];
% end
% 
% extra_moves = limit - size(moves,1);
% % disp('Moves left over')
% larger_i = find(board >= max(board(vine)));
% larger_i = setdiff(larger_i,vine);
% if ~isempty(larger_i)
%     end_row = mod(vine(end),m);
%     if end_row == 0
%         end_row = m;
%     end
%     end_col = ceil(vine(end)/m);
%     
%     larger_rows = mod(larger_i,m);
%     larger_rows(larger_rows == 0) = m;
%     larger_cols = ceil(larger_i/m);
%     
%     moves_req = abs(end_row - larger_rows) + abs(end_col - larger_cols) - 1;
%     possible_moves = larger_i(moves_req < extra_moves);
%     moves_req = moves_req(moves_req < extra_moves);
%     
%     while ~isempty(possible_moves) & extra_moves > 0
%         [junk, sort_i] = sort(moves_req);
%         
%         possible_moves = possible_moves(sort_i);
%         
%         test_move = possible_moves(1);
%         test_row = mod(test_move,m);
%         if test_row == 0
%             test_row = m;
%         end
%         test_col = ceil(test_move/m);
%         
%         if end_row > test_row
%             move_rows = test_row:end_row;
%             move_cols = zeros(size(move_rows)) + test_col;
%         elseif end_row < test_row
%             move_rows = test_row:-1:end_row;
%             move_cols = zeros(size(move_rows)) + test_col;
%         else
%             move_rows = [];
%             move_cols = [];
%         end
%         
%         if end_col == test_col
%             move_rows(end) = [];
%             move_cols(end) = [];
%         elseif end_col > test_col
%             move_cols = [move_cols, test_col+1:end_col-1];
%             move_rows = [move_rows, zeros(size(test_col+1:end_col-1))+end_row];
%         else
%             move_cols = [move_cols, test_col-1:-1:end_col+1];
%             move_rows = [move_rows, zeros(size(test_col-1:-1:end_col+1))+end_row];
%         end
%         
%         add_moves = move_cols*m + move_rows - m;
%         add_moves = [add_moves(1:end-1)', add_moves(2:end)'];
%         
%         if any(ismember(add_moves(:),vine))
%             possible_moves(1) = [];
%             moves_req(1) = [];
%         else
%             moves = [moves; add_moves];
%             extra_moves = limit - size(moves,1);
%             vine = [vine, moves(end,2)];
%             
%             for i = 1:size(add_moves,1)
%                 board(add_moves(i,:)) = [0 board(add_moves(i,1))];
%             end
%             
%             end_row = mod(vine(end),m);
%             if end_row == 0
%                 end_row = m;
%             end
%             end_col = ceil(vine(end)/m);
%             
%             larger_i = find(board >= max(board(vine)));
%             larger_i = setdiff(larger_i,vine);
%             larger_rows = mod(larger_i,m);
%             larger_rows(larger_rows == 0) = m;
%             larger_cols = ceil(larger_i/m);
%             
%             moves_req = abs(end_row - larger_rows) + abs(end_col - larger_cols) - 1;
%             possible_moves = larger_i(moves_req < extra_moves);
%             moves_req = moves_req(moves_req < extra_moves);
%         end
%         
%     end
% end
% end

function [cboard,moves,remain,rank,steps]=oftl2(cx,tx,cy,ty,cboard,rank,steps,moves,pos,remain,ncol,nrow,limit,s,i)
while (cx~=tx)||(cy~=ty)
    dx = sign(tx-cx);
    dy = sign(ty-cy);
    if (dx~=0)&&(dy~=0)
        correction = cboard(cy,cx+dx)<=cboard(cy+dy,cx);
        dy = dy*(1-correction);
        dx = dx*correction;
    end
    if (rank(cy+dy,cx+dx)>0)
        vx = cx+dx;
        vy = cy+dy;
        if (dx==0)
            [steps,moves,cboard,rank,remain,breakmarker]=oftl5(steps,limit,moves,pos,vy,vx,cboard,rank,remain,1,0,vx,ncol,cy,dy,cx,dx);
        else
            [steps,moves,cboard,rank,remain,breakmarker]=oftl5(steps,limit,moves,pos,vy,vx,cboard,rank,remain,0,1,vy,nrow,cy,dy,cx,dx);
        end;
        if breakmarker
            break;
        end
    end;
    steps = steps+1;
    if (steps > limit)
        break;
    end;
    moves(steps,1) = pos(cy,cx);
    moves(steps,2) = pos(cy+dy,cx+dx);
    cboard(cy+dy,cx+dx) = s(i); %cboard(cy,cx);
    cboard(cy,cx) = 0;
    rank(cy,cx) = 0;
    cx = cx+dx;
    cy = cy+dy;
end;
end


function [steps,moves,cboard,rank,remain,breakmarker]=oftl5(steps,limit,moves,pos,vy,vx,cboard,rank,remain,stepx,stepy,v_main,n_main,cy,dy,cx,dx)

if (v_main>1)&&(cboard(vy,vx)>cboard(vy-(stepy~=0),vx-(stepx~=0)))&&((v_main==n_main)||(cboard(vy+(stepy~=0),vx+(stepx~=0))>cboard(vy-(stepy~=0),vx-(stepx~=0))))
    [steps,moves,cboard,rank,remain,breakmarker]=oftl3(steps,limit,moves,pos,vy,vx,cboard,rank,remain,-stepx,-stepy);
elseif (v_main < n_main)&(cboard(vy,vx)>cboard(vy+(stepy~=0),vx+(stepx~=0)))
    [steps,moves,cboard,rank,remain,breakmarker]=oftl3(steps,limit,moves,pos,vy,vx,cboard,rank,remain,stepx,stepy);
else
    remain(rank(cy+dy,cx+dx)) = 0;
    breakmarker = 0;
end;

end

function [steps,moves,cboard,rank,remain,breakmarker]=oftl3(steps,limit,moves,pos,vy,vx,cboard,rank,remain,stepx,stepy)
steps = steps+1;
if (steps > limit)
    breakmarker = 1;
else
    breakmarker = 0;
    moves(steps,1) = pos(vy,vx);
    moves(steps,2) = pos(vy+stepy,vx+stepx);
    cboard(vy+stepy,vx+stepx) = cboard(vy,vx);
    cboard(vy,vx) = 0;
    if (rank(vy+stepy,vx+stepx)>0)
        remain(rank(vy+stepy,vx+stepx)) = 0;
    end;
    rank(vy+stepy,vx+stepx) = rank(vy,vx);
    rank(vy,vx) = 0;
end
end

function [moves,vine,scr] = dofilter(board,limit)
[nrow,ncol] = size(board);
ntot = nrow*ncol;
moves = zeros(limit,2);
nmove = 0;
pos = reshape(1:ntot,nrow,ncol);
m1 = median(board,1);
m2 = median(board,2);
if (sum(sign(diff(m1)))>0.67*ncol)||(sum(sign(diff(m2)))>0.67*nrow)||(sum(sign(diff(m2)))<-0.67*nrow)||(sum(sign(diff(m2)))<-0.67*nrow)
    d = diff(board,1,2);
    mb = median(board,2);
    md = median(d,2);
    rmb = repmat(mb,1,ncol);
    rmd = repmat(md,1,ncol);
    guess = rmb+kron(md,-(ncol-1)/2:(ncol-1)/2);
    bad = (abs(board-guess)>20*abs(rmd));
    for i = nrow:-1:1
        [nmove,moves]=oftl6(i,ncol,bad,moves,nmove,limit,pos);
        if (nmove>=limit)
            break;
        end;
    end;
    moves = moves(1:min(nmove,limit),:);
    for i = 1:size(moves,1)
        board(moves(i,2)) = board(moves(i,1));
        board(moves(i,1)) = 0;
    end;
    [vine,scr] = SimpleLongestPath(board);
else
    scr = -inf;
    moves = [];
    vine = [];
end
end

function [nmove,moves]=oftl6(i,ncol,bad,moves,nmove,limit,pos)
offset = 0;
for j = ceil(ncol/2):ncol
    if bad(i,j)
        offset = offset+1;
    else
        moves(nmove+1:nmove+offset,:) = [pos(i,j:-1:j-offset+1)' pos(i,j-1:-1:j-offset)'];
        nmove = nmove+offset;
        if (nmove>=limit)
            break;
        end;
    end;
end;
offset = 0;
for j = ceil(ncol/2):-1:1
    if bad(i,j)
        offset = offset+1;
    else
        moves(nmove+1:nmove+offset,:) = [pos(i,j:1:j+offset-1)' pos(i,j+1:1:j+offset)'];
        nmove = nmove+offset;
        if (nmove>=limit)
            break;
        end;
    end;
end;

end

%% alfi
function [moves, vine, gain] = alfi(board, limit)
moves=[];
ab=accumarray(1+board(:),1);
cab=flipud(cumsum(flipud(ab)))-ab(end);
[m,n]=size(board);
board=[zeros(1,n+2);[zeros(m,1),board,zeros(m,1)];zeros(1,n+2)];
[m,n]=size(board);
boardab=cab(1+board);
if max(ab(2:end))>1,
    % compute same-valued clusters
    [BoardCluster,ClusterValue,IdxList,IdxSegments,Nclusters]=connected(board);
    ClusterSize=diff(IdxSegments);
    ClusterNeighb=neighb(BoardCluster,board,Nclusters);
    % search between-clusters
    ClustersOrder = bellman(ClusterNeighb,ClusterValue'.*sqrt(ClusterSize'));
    % search within-clusters
    iC=bellman_postprocess(ClustersOrder,IdxList,IdxSegments,BoardCluster);
else
    % search between-pixels
    %     iC = bellman_pixel(board,limit,boardab);
    iC = SimpleLongestPath( board );
    iC = iC(end:-1:1);
    iC = iC(board(iC)>0);
end
% post-processing moves
if limit>0
    [board,iC,moves,limit] = postprocess(board,iC,moves,limit,boardab);
    %    [iC1,iC2]=ind2sub([m,n],moves);
    iC1=rem(moves-1,m)+1;
    iC2=(moves-iC1)/m+1;
    %    moves=sub2ind([m-2,n-2],iC1-1,iC2-1);
    moves=iC1-1+(m-2)*(iC2-2);
end

%[iC1,iC2]=ind2sub([m,n],iC);
iC1=rem(iC-1,m)+1;
iC2=(iC-iC1)/m+1;
%vine=sub2ind([m-2,n-2],iC1-1,iC2-1);
vine=iC1-1+(m-2)*(iC2-2);
vine=vine(end:-1:1);
board=board(2:end-1,2:end-1);
gain=sum(board(vine));

end

function [board,tiC,moves,limit] = postprocess(board,tiC,moves,limit,boardab)
[m,n]=size(board);
nC=numel(tiC);
tP=board>0;
if nC>1&&limit>0, %%%
    % grow laterally
    d=abs(diff(tiC))==1;
    E=[m*~d+d; m*d+~d];
    BorderUp={repmat(1+(0:n-1)*m,[m,1]),repmat((1:m)',[1,n])};
    BorderDown={repmat(m+(0:n-1)*m,[m,1]),repmat(m*(n-1)+(1:m)',[1,n])};
    %    tP=board>0;
    tP(tiC)=false;
    for n1=1:1e2,
        [board,moves,tiC,tP,E,limit,ok]=postprocess_lateral(board,moves,tiC,tP,E,limit,BorderDown,BorderUp,boardab);
        if ~ok, break; end
    end
end
if limit>0
    % grow end-points
    %    tP=board>0;
    tP(tiC)=false;
    for n1=1:1e2,
        [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,1,m,n);
        if ~ok, break; end
    end
    for n1=1:1e3,
        [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,2,m,n);
        if ~ok, break; end
    end
end
end

function [board,moves,tiC,tP,E,limit,ok]=postprocess_lateral(board,moves,tiC,tP,E,limit,BorderDown,BorderUp,boardab)
[m,n] = size(board);
ok=0;
for ndir=1:4
    offset = ndir>2.5;
    coeff = (2*mod(ndir,2)-1);
    K=size(tiC,2);
    idx1=find(tP(tiC((1+offset):2:K-1)+coeff*E(2,(1+offset):2:K-1))&tP(tiC((2+offset):2:K)+coeff*E(2,(1+offset):2:K-1)));
    for n2=numel(idx1):-1:1, %%%
        switch(ndir)
            case 1,
                tempA=tiC(2*idx1(n2)-1)+E(2,2*idx1(n2)-1):E(2,2*idx1(n2)-1):BorderDown{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)-1));
                tempB=tiC(2*idx1(n2))+E(2,2*idx1(n2)-1):E(2,2*idx1(n2)-1):BorderDown{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)));
            case 2,
                tempA=tiC(2*idx1(n2)-1)-E(2,2*idx1(n2)-1):-E(2,2*idx1(n2)-1):BorderUp{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)-1));
                tempB=tiC(2*idx1(n2))-E(2,2*idx1(n2)-1):-E(2,2*idx1(n2)-1):BorderUp{1+(E(2,2*idx1(n2)-1)==m)}(tiC(2*idx1(n2)));
            case 3,
                tempA=tiC(2*idx1(n2))+E(2,2*idx1(n2)):E(2,2*idx1(n2)):BorderDown{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)));
                tempB=tiC(2*idx1(n2)+1)+E(2,2*idx1(n2)):E(2,2*idx1(n2)):BorderDown{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)+1));
            case 4,
                tempA=tiC(2*idx1(n2))-E(2,2*idx1(n2)):-E(2,2*idx1(n2)):BorderUp{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)));
                tempB=tiC(2*idx1(n2)+1)-E(2,2*idx1(n2)):-E(2,2*idx1(n2)):BorderUp{1+(E(2,2*idx1(n2))==m)}(tiC(2*idx1(n2)+1));
        end
        idxZ=find(~tP(tempA),1,'first');
        [nill,idxA]=max((1./(abs(board(tiC(2*idx1(n2)-1+offset))-board(tempA))+1)).*(board(tiC(2*idx1(n2)-1+offset))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)+offset))));
        if ~nill,idxA=[]; end
        if ~isempty(idxA)&&idxA<idxZ
            idxZ=find(~tP(tempB),1,'first');
            idxB=find(board(tempA(idxA))>=board(tempB)&board(tempB)>=board(tiC(2*idx1(n2)+offset)),1,'first');
            if ~isempty(idxB)&&idxB<idxZ&&idxA+idxB-2<=limit
                ok=1;
                tiC=[tiC(1:2*idx1(n2)-1+offset),tiC(2*idx1(n2)-1+offset)+coeff*E(2,2*idx1(n2)-1+offset), tiC(2*idx1(n2)+offset)+coeff*E(2,2*idx1(n2)-1+offset), tiC((2*idx1(n2)+offset):end)];
                E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
                tP(tiC)=false;
                newmoveA=[tempA(2:idxA)',tempA(1:idxA-1)'];
                newmoveB=[tempB(2:idxB)',tempB(1:idxB-1)'];
                board(newmoveA(:,2))=board(tempA(idxA));
                board(newmoveB(:,2))=board(tempB(idxB));
                board(newmoveA(:,1))=0;
                board(newmoveB(:,1))=0;
                moves=[moves;flipud(newmoveA);flipud(newmoveB)];
                limit=limit-size(newmoveA,1)-size(newmoveB,1);
            end
        end
    end
end
end

function [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,type,m,n)
[m,n]=size(board);
ok=1;
switch(type)
    case 1
        mask=tP&board<2*board(tiC(1));
        D=Dijkstra(mask,tiC(1),m,n);
        idx1=find(mask&D-1<=limit&board>=board(tiC(1)));
        if isempty(idx1),
            mask=tP;
            D=Dijkstra(mask,tiC(1),m,n);
            idx1=find(mask&D-1<=limit&board>=board(tiC(1)));
            if isempty(idx1),ok=0;end
        end
        if ok
            if limit>150
                [nill,idx2]=min(board(idx1)+1e-3*D(idx1)); %%%t3
            else
                [nill,idx2]=min(board(idx1).*(D(idx1)/limit).^0.2); %%%t3
            end
            
            idx0=idx1(idx2);
            limit=limit-(D(idx0)-1);
            for n2=D(idx0)-1:-1:1
                tidx0 = idx0 + (D(idx0+1)==n2) + (~(D(idx0+1)==n2))*(-(D(idx0-1)==n2)+(~(D(idx0-1)==n2))*(m*(D(idx0+m)==n2)+(~(D(idx0+m)==n2))*(-m*(D(idx0-m)==n2))));
                moves=[moves; idx0, tidx0];
                board(tidx0)=board(idx0);
                board(idx0)=0;
                idx0=tidx0;
            end
            tiC=[idx0,tiC];
            tP(idx0)=false;
        end
    case 2
        mask=tP;
        D=Dijkstra(mask,tiC(end),m,n);
        idx1=find(mask&D-1<=limit&board<=board(tiC(end))&board>0);
        if isempty(idx1),ok=0;end
        if ok
            if limit>300
                [nill,idx2]=min(-board(idx1)+1e-3*D(idx1)); %%%t3
            else
                [nill,idx2]=min(-board(idx1).*(limit./D(idx1)).^0.085); %%%t3
            end
            idx0=idx1(idx2);
            limit=limit-(D(idx0)-1);
            for n2=D(idx0)-1:-1:1
                tidx0 = idx0 + (D(idx0+1)==n2) + (~(D(idx0+1)==n2))*(-(D(idx0-1)==n2)+(~(D(idx0-1)==n2))*(m*(D(idx0+m)==n2)+(~(D(idx0+m)==n2))*(-m*(D(idx0-m)==n2))));
                moves=[moves; idx0, tidx0];
                board(tidx0)=board(idx0);
                board(idx0)=0;
                idx0=tidx0;
            end
            tiC=[tiC,idx0];
            tP(idx0)=false;
        end
end
end

function iC = bellman(C,D)
N=size(C,1);
c=cell(1,N); for n1=1:N,c{n1}=find(C(:,n1)); end
%C=full(C);
IDX=zeros(N,1);
cD=D;
touched=true(1,N);
while any(touched)
    for touch=find(touched),
        touched(touch)=false;
        idx=c{touch};%find(C(:,touch));
        cDnew=D(idx)+cD(touch);
        idx2=find(cD(idx)<cDnew);
        if ~isempty(idx2)
            cD(idx(idx2))=cDnew(idx2);
            touched(idx(idx2))=true;
            IDX(idx(idx2))=touch;
        end
    end
end
[nill,idx]=max(cD);
iC=zeros(N,1);
for n1=1:N,
    if idx>0
        iC(n1)=idx;
        idx=IDX(idx);
    else break;
    end
end
iC=iC(1:n1-1);
end

function  iC=bellman_postprocess(ClustersOrder,IdxList,IdxSegments,BoardCluster)
[m,n]=size(BoardCluster);
iC=[];
for n1=1:numel(ClustersOrder)
    idx=IdxList(IdxSegments(ClustersOrder(n1)):IdxSegments(ClustersOrder(n1)+1)-1);
    if numel(idx)==1, iC=[iC,idx(1)];
    else
        % longest shortest-path
        ThisCluster=false(m,n);
        ThisCluster(idx)=true;
        D=Dijkstra(ThisCluster,iC,m,n);
        if n1==numel(ClustersOrder)
            temp=D(ThisCluster);
        else
            temp=D(ThisCluster).*(BoardCluster([false(1,n);ThisCluster(1:m-1,:)])==ClustersOrder(n1+1)|BoardCluster([ThisCluster(2:m,:);false(1,n)])==ClustersOrder(n1+1)|BoardCluster([false(m,1),ThisCluster(:,1:n-1)])==ClustersOrder(n1+1)|BoardCluster([ThisCluster(:,2:n),false(m,1)])==ClustersOrder(n1+1));
        end
        [nill,idx0]=max(temp(:));
        idx0=idx(idx0);
        E=zeros(2,D(idx0)-1);
        tiC=zeros(1,D(idx0));
        tiC(end)=idx0;
        for n2=D(idx0)-1:-1:1
            if      D(idx0+1)==n2, idx0=idx0+1; E(1,n2)=1;E(2,n2)=m;
            elseif  D(idx0-1)==n2, idx0=idx0-1; E(1,n2)=1;E(2,n2)=m;
            elseif  D(idx0+m)==n2, idx0=idx0+m; E(1,n2)=m;E(2,n2)=1;
            elseif  D(idx0-m)==n2, idx0=idx0-m; E(1,n2)=m;E(2,n2)=1;
            end
            tiC(n2)=idx0;
        end
        % now grow it
        tP=ThisCluster;
        tP(tiC)=false;
        ok=1;
        while ok
            ok=0;
            K=size(tiC,2);
            idx1=find(tP(tiC(1:2:K-1)+E(2,1:2:K-1))&tP(tiC(2:2:K)+E(2,1:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)+E(2,2*idx1(n2)-1), tiC(2*idx1(n2))+E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
                E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
                tP(tiC)=false;
            end
            K=size(tiC,2);
            idx1=find(tP(tiC(1:2:K-1)-E(2,1:2:K-1))&tP(tiC(2:2:K)-E(2,1:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)-E(2,2*idx1(n2)-1), tiC(2*idx1(n2))-E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
                E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
                tP(tiC)=false;
            end
            K=size(tiC,2);
            idx1=find(tP(tiC(2:2:K-1)+E(2,2:2:K-1))&tP(tiC(3:2:K)+E(2,2:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
                E=[E(:,1:2*idx1(n2)-1),E([2,1],2*idx1(n2)), E([1,2],2*idx1(n2)), E([2,1],2*idx1(n2)), E(:,2*idx1(n2)+1:end)];
                tP(tiC)=false;
            end
            K=size(tiC,2);
            idx1=find(tP(tiC(2:2:K-1)-E(2,2:2:K-1))&tP(tiC(3:2:K)-E(2,2:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
                E=[E(:,1:2*idx1(n2)-1),E([2,1],2*idx1(n2)), E([1,2],2*idx1(n2)), E([2,1],2*idx1(n2)), E(:,2*idx1(n2)+1:end)];
                tP(tiC)=false;
            end
        end
        
        iC=[iC,tiC];
    end
end
end

function B = neighb(A,V,nA)
[m,n] = size(A);
B=sparse(A(:,1:n-1),A(:,2:n),double(V(:,1:n-1)>V(:,2:n)),nA,nA)+sparse(A(:,2:n),A(:,1:n-1),double(V(:,2:n)>V(:,1:n-1)),nA,nA)+sparse(A(1:m-1,:),A(2:m,:),double(V(1:m-1,:)>V(2:m,:)),nA,nA)+sparse(A(2:m,:),A(1:m-1,:),double(V(2:m,:)>V(1:m-1,:)),nA,nA);
B=B>0;
end

function [B,C,p,r,nc] = connected(A)
[m,n] = size(A);
N = m*n ;
K = reshape (1:N, m, n) ;
V = A(:,1:n-1) == A(:,2:n);
H = A(1:m-1,:) == A(2:m,:);
G = sparse(K([V,false(m,1)]),K([false(m,1),V]),1,N,N)+sparse(K([H; false(1,n)]),K([false(1,n); H]),1,N,N);
G = G + G' + speye(N);
[p q r] = dmperm(G);
nc = numel(r) - 1;
C = A(p(r(1:nc)));
B = ones(m, n);
for a = 2:nc
    B(p(r(a):r(a+1)-1)) = a;
end
end

function D = Dijkstra(A,i1,m,n)
D=inf(m,n);
D(~A)=0;
if isempty(i1),
    i1=find(A,1,'first');
    mode=0;
else
    i1=i1(end);
    mode=1;
end
D(i1)=1;
for n1=2:m*n,
    idx1=D(i1+1)>n1;
    idx2=D(i1-1)>n1;
    idx3=D(i1+m)>n1;
    idx4=D(i1-m)>n1;
    %i1=unique([i1(idx1)+1,i1(idx2)-1,i1(idx3)+m,i1(idx4)-m]);
    X=false(m,n);
    X(i1(idx1)+1)=true;
    X(i1(idx2)-1)=true;
    X(i1(idx3)+m)=true;
    X(i1(idx4)-m)=true;
    i1=find(X);
    if isempty(i1), break; end
    D(i1)=n1;
end
if mode, D(D>0)=D(D>0)-1; end
end

%% kurt
function [vine, bestsom] = real_kurt(A)
moves = [];
%[bestsom,bestvine] = max(A(:));

[m,n]=size(A);

% inhoud = unique(A)';
% heur = sparse(0);
% for i=inhoud;
%     heur(i+1) = sum(sum(A(A<=i)));
% end
% H=full(heur(A+1));

content = unique(A);
t=sort(A(:));
tt = cumsum(t);
%ttt = find(diff([t; 0]));
%tttt = tt(ttt);
tttt = tt(diff([t; 0])~=0);
t2 = zeros(1,size(content,1));
t2(content+1) = tttt;
H = t2(A+1);
if size(H,1) ~= size(A,1)
    H=H';
end

B=A;

IND =  reshape(1:m*n,[m,n]);
C = num2cell(IND);

updated = true(size(A));

iteration = 0;
while 1
    G=B+H;
    G=G.*updated;
    [maxg,startindex] = max(G(:));
    
    [bestsom] = max(B(:));
    
    if maxg<=bestsom
        break
    end
    
    if maxg==0
        break
    end
    
    iteration = iteration+1;
    i = mod(startindex,m);
    j = ceil(startindex/m);
    if i==0
        i=m;
    end
    
    if 1<i && A(i-1,j)<=A(i,j)
        [B,C,updated]=oftl4(IND,i,j,A,B,C,updated,-1,0);
    end
    if i<size(A,1) && A(i+1,j)<=A(i,j)
        [B,C,updated]=oftl4(IND,i,j,A,B,C,updated,1,0);
    end
    if j<size(A,2) && A(i,j+1)<=A(i,j)
        [B,C,updated]=oftl4(IND,i,j,A,B,C,updated,0,1);
    end
    if 1<j && A(i,j-1)<=A(i,j)
        [B,C,updated]=oftl4(IND,i,j,A,B,C,updated,0,-1);
    end
    
    updated(i,j)=false;
    
end

[bestsom,index]=max(B(:));
vine = fliplr(C{index});

end

function [B,C,updated]=oftl4(IND,i,j,A,B,C,updated,stepi,stepj)
verplaatsindex = IND(i+stepi,j+stepj);
if all(verplaatsindex~=C{i,j})
    [B(verplaatsindex),which] = max([B(verplaatsindex),B(i,j)+A(verplaatsindex)]);
    if which == 2
        updated(verplaatsindex)=true;
        C{verplaatsindex} = [C{i,j} verplaatsindex];
    end
end
end

function [vine bestScore] = SimpleLongestPath( board )
% Treats the board as a directed-acyclic-graph and finds the best
% possible solution. Upper-left wins to break cycles in the graph
% (when adjacent squares have the same value).
%
% This is basically a completely rewritten version of the sebastian()
% function, and is about 15 times faster than sebastian() on the contest
% machine.
%
% This function is only called 8 times presently on the sample set so it
% doesn't make much difference. I optimised it so extremely expecting to use it as
% a candidate-scoring/fitness function as part of some sort of larger algorithm.
%
% - Wesley

[rows cols] = size( board );
count = rows * cols;
sentinel = count + 1;

board = board(:);
score = [board; 0];

seq = reshape( 1:count, rows, cols );

north = [sentinel(ones( 1, cols )); seq(1:end-1,:)];
north = north(:);
north(board<score(north)) = sentinel;

south = [seq(2:end,:); sentinel(ones( 1, cols ))];
south = south(:);
south(board<score(south)) = sentinel;

west = [sentinel(ones( rows, 1 )); (1:count-rows)'];
west(board<score(west)) = sentinel;

east = [(rows+1:count)'; sentinel(ones( rows, 1 ))];
east(board<score(east)) = sentinel;

prev = zeros( sentinel, 1 );

[dum, order] = sort( board );
for i = 1:count
    curr = order(i);
    
    % Yes, this fixed size (4) sorting network really
    % is faster than using max() for some reason.
    
    northNeighbour = north(curr);
    northScore = score(northNeighbour);
    
    southNeighbour = south(curr);
    southScore = score(southNeighbour);
    
    if northScore >= southScore
        scoreNS = northScore;
        neighbourNS = northNeighbour;
    else
        scoreNS = southScore;
        neighbourNS = southNeighbour;
    end
    
    westNeighbour = west(curr);
    westScore = score(westNeighbour);
    
    eastNeighbour = east(curr);
    eastScore = score(eastNeighbour);
    
    if westScore >= eastScore
        scoreEW = westScore;
        neighbourEW = westNeighbour;
    else
        scoreEW = eastScore;
        neighbourEW = eastNeighbour;
    end
    
    if scoreEW >= scoreNS
        prev(curr) = neighbourEW;
        score(curr) = score(curr) +  scoreEW;
    else
        prev(curr) = neighbourNS;
        score(curr) = score(curr) +  scoreNS;
    end
end

[bestScore c] = max( score );
for i = 1:count
    seq(i) = c;
    c = prev(c);
    if c == sentinel
        break
    end
end
vine = seq(i:-1:1);
end

% robert
function [moves, vine, score] = robert(board, limit)

% The Fragrant Honeysuckle gives up on moving and just tries to produce
% some prettier vines 8-)

% Faster to use isconnected logic, or to embed in a zeros(m+1,n+1) and guard?


% Grow a simple vine on given board
%[vine, score, dirsimple, moves]  = robert2(board);


% Try constructing a boardwalk down left side
[movesb, boardb] = robert1(board, limit);
[vineb, vscoreb, dirsimple] = robert2(boardb);     % Rerun vine code on modified board


% Improve score on board 47 and similar that feature a steady slope and speckle

IsMonotonous = ~(any(dirsimple(vineb) == 1) && any(dirsimple(vineb) == -1));
if IsMonotonous
    % Only bother if I have a decent score already and vine goes
    % solely up or solely down.
    [movesm, boardm] = monotonous(board, vineb, dirsimple, limit);
    [vinem, vscorem] = SimpleLongestPath(boardm);     % Rerun vine code on modified board
end

% Pick best
vine = vineb;
moves = movesb;
score = vscoreb;
if IsMonotonous && vscorem > score
    vine = vinem;
    moves = movesm;
    score = vscorem;
end

end % solver function

function [moves, board] = monotonous(boardin, vine, dir, limit)

% Focus n high-scoring vines with strong vertical orientation and high limit,
% and aim to shuffle blocking values out of high-scoring rows.
% Should never damage a unidirectional vine but ineffective unless the
% board has the structure of board 47

moves = [];
board = boardin;
[rows,cols] = size(board);
boardsize   = numel(board);
lastdir = dir(vine(end));       % initialise with vine top
for i = size(vine,1)-1:-1:2     % walk down vine to polish high scores first
    if abs(dir(vine(i))) == 1   % found a vertical step
        valmin = board(vine(i)+dir(vine(i)));
        valmax = board(vine(i)); % range of non-blocking values
        %            [row,col] = ind2sub([rows,cols],vine(i));
        row=rem(vine(i)-1,rows)+1;
        col=(vine(i)-row)/rows+1;
        if lastdir > 0 && vine(i)+lastdir <= boardsize-cols
            % a right edge in mid board
            moveToCol = col+1;
            [board,moves]=oftl(board,moves,limit,row,moveToCol,valmin,valmax,cols,rows);
            row = row - 1;
            if row > 0          % Needed?
                moveToCol = col+1;
                [board,moves]=oftl(board,moves,limit,row,moveToCol,valmin,valmax,cols,rows);
            end % if second row > 0
        end % if we have a midboard right edge
        
        
        % Now repeat whole thing for left edges.  Should be a minor pickup
        % from having more moves in high-scoring areas.
        % Check lastedge logic when I make 2 upward moves in a row
        
        
        if lastdir == -cols && vine(i) > cols
            % a left edge in mid board.  Check logic.
            moveToCol = col-1;
            [board,moves]=oftl7(board,moves,limit,row,moveToCol,valmin,valmax,rows);
        end % if we have a midboard left edge
        
    end % vertical step
    lastdir = dir(vine(i));
end % for each leaf of vine

% Have now modified board, hopefully improved it in a few high-scoring cases

end % function monotonous


function [board,moves]=oftl7(board,moves,limit,row,moveToCol,valmin,valmax,rows)
while size(moves,1) < limit  % loop through all blockages (if any)
    while ( board(row,moveToCol) >= valmin && ...
            board(row,moveToCol) <= valmax    )  && moveToCol > 1
        moveToCol = moveToCol - 1;
    end
    if moveToCol > 1 % Found a blockage I'd like to replace
        moveFromCol = moveToCol-1;
        while moveFromCol > 0 && ...
                ( board(row,moveFromCol) < valmin || ...
                board(row,moveFromCol) > valmax        )
            moveFromCol = moveFromCol - 1;
        end
        if moveFromCol > 0
            % Something to replace it with exists, so do slide
            for j = moveFromCol:moveToCol-1
                if size(moves,1) < limit
                    moves = [moves; [row+rows*(j-1), row+rows*j]];
                    board(row,j+1) = board(row,j);
                    board(row,j) = 0;
                else
                    break;  % Used up all moves
                end
            end
        else
            break;  % no more nonblocks to use so row finished
        end
    else
        break;      % no blocks so row finished
    end             % (needed in case block is in other row...?)
end % First row while loop.  Break out when complete and apply
% same logic to second row

end

function [board,moves]=oftl(board,moves,limit,row,moveToCol,valmin,valmax,cols,rows)

while size(moves,1) < limit  % loop through all blockages (if any)
    while ( board(row,moveToCol) >= valmin && ...
            board(row,moveToCol) <= valmax    )  && moveToCol <= cols-1
        moveToCol = moveToCol + 1;
    end
    if moveToCol < cols % Found a blockage I'd like to replace
        moveFromCol = moveToCol+1;
        while moveFromCol <= cols && ...
                ( board(row,moveFromCol) < valmin || ...
                board(row,moveFromCol) > valmax        )
            moveFromCol = moveFromCol + 1;
        end
        if moveFromCol <= cols
            % Something to replace it with exists, so do slide
            for j = moveToCol:moveFromCol-1
                if size(moves,1) < limit
                    moves = [moves; [(row+rows*j), (row+rows*(j-1))]];
                    board(row,j) = board(row,j+1);
                    board(row,j+1) = 0;
                else
                    break;  % Used up all moves
                end
            end
        else
            break;  % no more nonblocks to use so row finished
        end
    else
        break;      % no blocks so row finished
    end             % (needed in case block is in other row...?)
end
end

function [moves, board] = robert1(boardin, limit)
% Focus n high-score boards like 8 where there is a high limit but no structure.
% Aim to construct a path from scratch. To keep the movement choices simple do
% this along an edge; a spiral in centre would require far fewer moves on
% average but would require more complex pathfinding

% Minor bug with equal values, which may upset vine constructor so add a small
% amount of noise.

moves = zeros(limit,2); mv=0;
[rows,cols] = size(boardin);
boardsize   = numel(boardin);
board = boardin - reshape(1:boardsize,size(boardin))*1e-4;
Target      = 1;                    % Pick a corner, any corner
while mv < limit         % Find another block to move
    [Biggest,BlockAt] = max(board(:));    % Move biggest
    if Biggest <= 0                 % Prob not needed
        break;
    end
    Vt=rem(Target-1,rows)+1;
    Ut=(Target-Vt)/rows+1;
    
    % First move onto a channel that will be cleared of blocks
    Vb=rem(BlockAt-1,rows)+1;
    Ub=(BlockAt-Vb)/rows+1;
    if Ub-Ut > 1     % Only if not next to target column
        Nudge = 1-mod(Vb-1,3);    % [+1,0,-1,...]
        if (Nudge ~= 0) && ~(Vb==rows && Nudge==1)
            mv=mv+1; moves(mv,1:2) = [BlockAt,BlockAt+Nudge];
            board(BlockAt+Nudge) = board(BlockAt);
            board(BlockAt) = 0;
            BlockAt = BlockAt+Nudge;
        end
    end
    
    % Now move towards target
    Vb = rem(BlockAt-1,rows)+1; U = Ut-(BlockAt-Vb)/rows-1;
    V = Vt-Vb;      % Horiz and Vert moves required
    
    while (U~=0 || V~=0) && mv<limit
        if U<-1 || V==0
            mv=mv+1; moves(mv,1:2) = [BlockAt,BlockAt-rows];
            board(BlockAt-rows) = board(BlockAt);
            board(BlockAt) = 0;
            BlockAt = BlockAt-rows;
            U = U+1;
        elseif V>0         % Just ignore erases until I get it working
            mv=mv+1; moves(mv,1:2) = [BlockAt,BlockAt+1];
            board(BlockAt+1) = board(BlockAt);
            board(BlockAt) = 0;
            BlockAt = BlockAt+1;
            V=V-1;
        else
            mv=mv+1; moves(mv,1:2) = [BlockAt,BlockAt-1];
            board(BlockAt-1) = board(BlockAt);
            board(BlockAt) = 0;
            BlockAt = BlockAt-1;
            V=V+1;
        end
    end  % Sliding this block
    board(Target) = -board(Target);  % Ignore moved blocks
    Target=Target+(mod(Ut,2))*(Vt<rows)-(~mod(Ut,2))*(Vt>1)+rows*(mod(Ut,2))*(~(Vt<rows))+rows*(~mod(Ut,2))*(~(Vt>1));
end % Finding blocks to slide
board = abs(board);
moves=moves(1:mv,:);
end % function boardwalk

function [vine, vscore, dir, moves] = robert2(board)

% Test each square on the board from low to high values to find whether vine
% can grow upwards onto an adjacent target square.

% To handle plateaux I superimpose a faint watermark on flat areas to nudge vine
% towards a space-covering pattern.  On each plateau I track 4 orientations
% separately, but whenever looking up to a higher level I can pick the best.

% Quite pretty, but as all plateaux problems in testset are low-scoring I don't
% think it makes much difference.

%

dir   = zeros(size(board));  % points way down vine, 0 at root
%   done  = false(size(board));  % have already examined this one
[rows,cols] = size(board);
boardsize   = numel(board);
watermark = reshape(1:boardsize,size(board)) * 1e-6;
watermark(:,2:2:cols) = flipud(watermark(:,2:2:cols));
board = board + watermark .* (mod(board,2)-0.5);
% superimpose a faint watermark on flat areas.  Direction should change
% between adjacent levels (but will fail if steps are even)
score = board;                % score for a vine starting and ending here
[val, ind]  = sort(board(:));

% Loop through board from low to high values to find the score associated
% with the best vine growing down from each square
% Needs more care on equal values, for now just weave up and down
%   -- Could try 4 weaves and pick best, or proper handling of plateaux

for i = 1:boardsize;
    test = ind(i);
    stepvec = [-rows,-1,1,rows];
    for step = stepvec;
        targ = test + step;
        if targ > 0 && targ <= boardsize && ...
                (abs(step)==rows || mod(min(test,targ), rows))
            % consider only adjacent squares on board
            if (board(targ) > board(test)) && (board(targ)+score(test) > score(targ))
                % Targ can grow from me (and has nothing better yet)
                score(targ) = score(test) + board(targ);
                % add test vine to targ score
                dir(targ) = -step;
                % and note where it came from
            end
        end
    end
end

% Assemble vine from top down
[dum,leaf] = max(score(:));
vine = [];
while leaf
    vine = [vine;leaf];
    leaf = (dir(leaf)~=0)*(leaf+dir(leaf));
end

vine = flipud(vine);  % return it in correct order

vscore = score(vine(end));
moves = [];
end % growvine

%% Mike Add convert2soc
function [new_moves,new_vine,success_flag] = convert2soc(vine,limit4soc,board,siz,n_lut)

% If the vine already contains 1, then just truncate the eariler part of
% the vine
if any(vine==1)
    new_vine = vine(find(vine==1):end);
    new_moves = zeros(0,2);
    success_flag = true;
    return
end

if nargin<4 || isempty(siz) || isempty(n_lut)
    siz = size(board);
    ind_board = reshape(1:prod(siz),siz);
    north_neighbor = shiftnD(ind_board,1,0);
    south_neighbor = shiftnD(ind_board,-1,0);
    east_neighbor = shiftnD(ind_board,-2,0);
    west_neighbor = shiftnD(ind_board,2,0);
    neighbors = cat(3,north_neighbor,south_neighbor,east_neighbor,west_neighbor); %these are indices, not values
    n_lut = reshape(neighbors,prod(siz),4);
end


[l2v,grown] = paths2vine(vine,siz,n_lut,limit4soc);

% We want the path which reaches the vine closest to the base without going
% over the limit of moves
attach_point = find(l2v<limit4soc-1,1,'first');
if isempty(attach_point)
    % Can't connect to vine, need emergency backup plan
%     disp('CANT REACH')
    new_moves = zeros(0,2);
    new_vine = [1]; % to conform with minimal rules
    success_flag = false;
    return
end
trunc_vine = vine(attach_point:end);
graft_length = l2v(attach_point);

graft_path = make_graft(vine(attach_point),grown,graft_length,n_lut);

graft_moves = [graft_path(end:-1:2),graft_path(end-1:-1:1)];
new_moves = [graft_moves(1,:); graft_moves];
    
new_vine = [graft_path(:); trunc_vine(:)]';
success_flag = true;

end


% idea here is to create a corridor of zeros from square one to the end of
% the already-optimized vine.  If limit is 0, then just 

% SOC ends 4pm PT TUESDAY
function [graft_path] = make_graft(attach_point_ind,grown,graft_length,n_lut)

graft_path = zeros(graft_length,1);
pathend = attach_point_ind;
level = graft_length+1;

for j = graft_length:-1:1

%while graft_path(end)~=1
    local_nodes = nb_nodes(pathend,n_lut);
    for i = 1:length(local_nodes)
        if grown(local_nodes(i))>0 && grown(local_nodes(i))<level
            pathend = local_nodes(i);
            graft_path(j) = pathend;
            level = grown(pathend);
            break
        end
    end
end



end

function [l2v,grown] = paths2vine(vine,siz,n_lut,limit)
vine_mask = false(siz);
vine_mask(vine) = true;



grown = zeros(siz);
grown(1)=1;
level = 1;
newlygrown = true(siz);
while any(newlygrown(:)) && ~(level >= limit-2)
    level = level+1;
    % grow by shifting current level
    shifts  = shiftnD(grown,-1,0) + shiftnD(grown,+1,0) + shiftnD(grown,-2,0) + shiftnD(grown,2,0);
    shifts(vine_mask) = 0; 
    newlygrown = shifts>0 & grown==0;
    grown(newlygrown) = level; 
end

grown2 = grown;
grown2(grown2==0)=Inf;
l2v = zeros(length(vine),1);
for i = 1:length(vine)
    adj_vals = nb_vals(vine(i),grown2,n_lut);
    l2v(i) = min(adj_vals);
end



end

% nb_nodes
function local_nodes = nb_nodes(node,n_lut)
local_nodes = n_lut(node,:);
local_nodes(local_nodes==0)=[] ;
end

function local_nb_vals = nb_vals(node,board,n_lut)
nb_nodes = n_lut(node,:);
nb_nodes(nb_nodes==0) = [];
local_nb_vals = board(nb_nodes);
end

function shifted_array = shiftnD(array, direction,padding)
% shifted_array = shiftnD(array, direction)
% This function should take the given array and shift it in the direction
% given.  The direction should be the number of the dimension in which to
% shift, positive if shifting in the direction of increasing indices, 
% negative if shifting in the direction of decreasing indices.

if direction ==0 %don't shift
    shifted_array = array;
    return
end

dims = size(array);

shiftdim = abs(direction); % this is the dimension to shift in
shiftdir = sign(direction); % this is the direction to shift in that dimension

if abs(direction)>length(dims)
    %error('Direction chosen is inappropriate')
end

% construct extended matrix
for i = 1:length(dims)
    c{i} = 1:dims(i);
end
c2 = c;

if shiftdir>0
    c{shiftdim} = 2:dims(shiftdim);
    c2{shiftdim} = 1:dims(shiftdim)-1;
elseif shiftdir<0
    c{shiftdim} = 1:dims(shiftdim)-1;
    c2{shiftdim} = 2:dims(shiftdim);
end

shifted_array = repmat(padding,dims);
shifted_array(c{:}) = array(c2{:});

end
%% %%%%%%% SOC FUNCTIONS %%%%%%%%%%%%%%%%%%%%%%
% I need a function for the case where a graft can't reach the vine.
% Perhaps the most likely case for this is when there are no moves allowed
% SO, we need a solver which grows the best possible vine from square onefunction [moves, vine, gain] = alfi(board, limit)
function [moves, vine, gain] = alfi_soc(board, limit)
moves=[];
ab=accumarray(1+board(:),1); % ab(i) counts the number of board locations 
                             % with value i-1 (because of the +1)
cab=flipud(cumsum(flipud(ab)))-ab(end); % computes the reverse cumulative sum of 
                             % the ab list, reverses it, and subtracts the
                             % number of times the highest value in board
                             % occurs (I wonder if it should be the number
                             % of times 0 occurs (ab(1)) instead, but I
                             % don't understand why this is here, so I'm
                             % not sure
[m,n]=size(board);
board=[zeros(1,n+2);[zeros(m,1),board,zeros(m,1)];zeros(1,n+2)]; % has the 
    % result of padding the board by a ring of zeros all the way around
[m,n]=size(board);
boardab=cab(1+board); %

if max(ab(2:end))>1,
    % If any of the ab vector is greater than 1, then there might be
    % clusters, since the same value occurs more than once in the board. We
    % don't care about ab(1) because those are zero-weight nodes
    
    % compute same-valued clusters
    [BoardCluster,ClusterValue,IdxList,IdxSegments,Nclusters]=connected(board);
    % BoardCluster is a label matrix of connected, identically-valued board
    % regions. ClusterValue is a vector of length Nclusters such that
    % ClusterValue(i) gives you the board value of one node in the region
    % described by BoardCluster==i.  IdxList and IdxSegments are p and r
    % decribed in the comments to the connected function below. Nclusters
    % is the number of clusters found. It is important to note that
    % nodes with no identically-valued neighbors are just treated as
    % clusters of size 1 throughout this section, so there are no nodes not
    % in a "cluster"
    ClusterSize=diff(IdxSegments); % vector of the number of nodes in each cluster
    ClusterNeighb=neighb_soc(BoardCluster,board,Nclusters); % Finds sparse matrix 
            % describing which clusters neighbor one another and can be moved to, i.e.
            % ClusterNeighb(i,j) = 1, iff BoardCluster==i neighbors
            % BoardCluster==j and ClusterValue(i)>ClusterValue(j)
            % (I'm not 100% sure on this, but it's something like that)
    % search between-clusters
    ClustersOrder = bellman_soc(ClusterNeighb,ClusterValue'.*sqrt(ClusterSize'));
    % search within-clusters
    iC=bellman_postprocess_soc(ClustersOrder,IdxList,IdxSegments,BoardCluster);
else
    % search between-pixels
    iC = bellman_pixel_soc(board,limit,boardab);
end
% post-processing moves
if limit>0
    [board,iC,moves,limit] = postprocess_soc(board,iC,moves,limit,boardab);
    [iC1,iC2]=ind2sub([m,n],moves);
    moves=sub2ind([m-2,n-2],iC1-1,iC2-1);
end

[iC1,iC2]=ind2sub([m,n],iC);
vine=sub2ind([m-2,n-2],iC1-1,iC2-1);
vine=vine(end:-1:1);
board=board(2:end-1,2:end-1);
gain=sum(board(vine));

end
%% postprocess_soc
function [board,tiC,moves,limit] = postprocess_soc(board,tiC,moves,limit,boardab)
[m,n]=size(board);
nC=numel(tiC);
if nC>1&&limit>0, %%%
    % grow laterally
    d=abs(diff(tiC))==1;
    E=[m*~d+d; m*d+~d];
    BorderUp={repmat(1+(0:n-1)*m,[m,1]),repmat((1:m)',[1,n])};
    BorderDown={repmat(m+(0:n-1)*m,[m,1]),repmat(m*(n-1)+(1:m)',[1,n])};
    tP=board>0;
    tP(tiC)=false;
    for n1=1:1e3,
        [board,moves,tiC,tP,E,limit,ok]=postprocess_lateral(board,moves,tiC,tP,E,limit,BorderDown,BorderUp,boardab);
        if ~ok, break; end
    end
end
if limit>0
    % grow end-points
    tP=board>0;
    tP(tiC)=false;
    for n1=1:1e3,
        [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,1);
        if ~ok, break; end
    end
%     for n1=1:1e3,
%         [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,2);
%         if ~ok, break; end
%     end
end
end
%% bellman_pixel
function iC = bellman_pixel_soc(board,limit,boardab)
[m,n]=size(board);
cD=board;
IDX=zeros(m,n);
touched = false(m,n);
touched(2,2) = true;
%touched=board>0;
%valid=touched;
valid = board>0;
while any(touched(:))
    touch=find(touched);
    touched(touch)=false;
    for dx=[1,-1,m,-m],
        cDnew=board(touch+dx)+cD(touch); %%%
        idx=find(board(touch+dx)>board(touch)&cDnew>cD(touch+dx));
        touched(touch(idx)+dx)=valid(touch(idx)+dx);
        cD(touch(idx)+dx)=cDnew(idx);
        IDX(touch(idx)+dx)=touch(idx);
    end
end
[nill,idx]=max(cD(:)-board(:));
if idx==1;
    idx=m+2; % want idx = 1 in original board
end
iC=zeros(1,m*n);
for n1=1:m*n,
    if idx>0
        iC(n1)=idx;
        idx=IDX(idx);
    else break;
    end
end
iC=iC(1:n1-1);
end
%% bellman
function iC = bellman_soc(C,D)
% C is sparse cluster neighbor matrix ClusterNeighb, and D is a weighted
% value for each cluster (single node value x sqrt(number of elements in cluster)

N=size(C,1); % Number of clusters
c=cell(1,N); 
% Construct list of neighbors for each cluster
for n1=1:N,
%    c{n1}=find(C(:,n1)); % c{i} = vector of cluster numbers which are reachable from cluster i
    c{n1}=find(C(n1,:)); % c{i} = vector of cluster numbers which are reachable from cluster i
end
%C=full(C);
IDX=zeros(N,1);
cD=D; % This is a weight to be gained by visiting a cluster
touched = false(1,N);
touched(2) = true;
%touched=true(1,N);  % start search from every cluster
while any(touched)
    for touch=find(touched),  % loop over clusters on the frontier 
        touched(touch)=false; % mark the current cluster as visited
        idx=c{touch};%find(C(:,touch)); % look up the visitable clusters for the current cluster
        cDnew=D(idx)+cD(touch); % This adds the weight of visitable clusters to the stored weight of the current cluster
        idx2=find(cD(idx)<cDnew); % This finds which scores are better than the current score
        if ~isempty(idx2)
            cD(idx(idx2))=cDnew(idx2); % For neighboring clusters where the new score is better than the old, update the stored score
            touched(idx(idx2))=true; % furthermore, you want to start working from these clusters again, since you have a better score you want to propagate
            IDX(idx(idx2))=touch; %
        end
    end
end
[nill,idx]=max(cD-D);
if idx==1
    idx=2; % can't use boarder padding, must use cluster 2
end
iC=zeros(N,1);
for n1=1:N,
    if idx>0
        iC(n1)=idx;
        idx=IDX(idx);
    else break;
    end
end
iC=iC(1:n1-1);
end
%% bellman_postprocess
function  iC=bellman_postprocess_soc(ClustersOrder,IdxList,IdxSegments,BoardCluster)
% This is the intra-cluster vine-growing algorithm
[m,n]=size(BoardCluster);
iC=[];
for n1=1:numel(ClustersOrder)
    idx=IdxList(IdxSegments(ClustersOrder(n1)):IdxSegments(ClustersOrder(n1)+1)-1);
    if numel(idx)==1, iC=[iC,idx(1)]; % if the current cluster is just one pixel, just add it to the list, no further path determinations needed
    else
        % longest shortest-path
        ThisCluster=false(m,n);
        ThisCluster(idx)=true;  % Makes ThisCluster a logical map of the current cluster
        D=Dijkstra(ThisCluster,iC,m,n); % Grows a map of number of moves to get from the 
                    % end of the current path (iC) to anywhere within the
                    % cluster
        if n1==numel(ClustersOrder)
%             temp=D(ThisCluster); % if this is the last cluster, then we don't need to find a particular endpoint, the furthest place is fine...
            temp = false(size(D(ThisCluster))); 
            temp(1) = true;
        else
            temp=D(ThisCluster).*(BoardCluster([false(1,n);ThisCluster(1:m-1,:)])==ClustersOrder(n1+1)...
                |BoardCluster([ThisCluster(2:m,:);false(1,n)])==ClustersOrder(n1+1)...
                |BoardCluster([false(m,1),ThisCluster(:,1:n-1)])==ClustersOrder(n1+1)...
                |BoardCluster([ThisCluster(:,2:n),false(m,1)])==ClustersOrder(n1+1));
            % this finds where the path across the cluster needs to exit to
            % continue the path
        end
        [nill,idx0]=max(temp(:));
        idx0=idx(idx0);
        E=zeros(2,D(idx0)-1);
        tiC=zeros(1,D(idx0));
        tiC(end)=idx0;
        % This for loop constructs a shortest path across the cluster from the
        % identified start to the identified end
        for n2=D(idx0)-1:-1:1
            if      D(idx0+1)==n2, idx0=idx0+1; E(1,n2)=1;E(2,n2)=m;
            elseif  D(idx0-1)==n2, idx0=idx0-1; E(1,n2)=1;E(2,n2)=m;
            elseif  D(idx0+m)==n2, idx0=idx0+m; E(1,n2)=m;E(2,n2)=1;
            elseif  D(idx0-m)==n2, idx0=idx0-m; E(1,n2)=m;E(2,n2)=1;
            end
            tiC(n2)=idx0;
        end
        % now grow it
        tP=ThisCluster;
        tP(tiC)=false;
        ok=1;
        while ok
            % This while loop tries to add bits to make the
            % shortest path wind across as much of the cluster as possible,
            % while maintaining the starting and end points
            ok=0;
            K=size(tiC,2);
            idx1=find(tP(tiC(1:2:K-1)+E(2,1:2:K-1))&tP(tiC(2:2:K)+E(2,1:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)+E(2,2*idx1(n2)-1), tiC(2*idx1(n2))+E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
                E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
                tP(tiC)=false;
            end
            K=size(tiC,2);
            idx1=find(tP(tiC(1:2:K-1)-E(2,1:2:K-1))&tP(tiC(2:2:K)-E(2,1:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)-1),tiC(2*idx1(n2)-1)-E(2,2*idx1(n2)-1), tiC(2*idx1(n2))-E(2,2*idx1(n2)-1), tiC(2*idx1(n2):end)];
                E=[E(:,1:2*idx1(n2)-2),E([2,1],2*idx1(n2)-1), E([1,2],2*idx1(n2)-1), E([2,1],2*idx1(n2)-1), E(:,2*idx1(n2):end)];
                tP(tiC)=false;
            end
            K=size(tiC,2);
            idx1=find(tP(tiC(2:2:K-1)+E(2,2:2:K-1))&tP(tiC(3:2:K)+E(2,2:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)+E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
                E=[E(:,1:2*idx1(n2)-1),E([2,1],2*idx1(n2)), E([1,2],2*idx1(n2)), E([2,1],2*idx1(n2)), E(:,2*idx1(n2)+1:end)];
                tP(tiC)=false;
            end
            K=size(tiC,2);
            idx1=find(tP(tiC(2:2:K-1)-E(2,2:2:K-1))&tP(tiC(3:2:K)-E(2,2:2:K-1)));
            for n2=numel(idx1):-1:1
                ok=1;
                tiC=[tiC(1:2*idx1(n2)),tiC(2*idx1(n2))-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1)-E(2,2*idx1(n2)), tiC(2*idx1(n2)+1:end)];
                E=[E(:,1:2*idx1(n2)-1),E([2,1],2*idx1(n2)), E([1,2],2*idx1(n2)), E([2,1],2*idx1(n2)), E(:,2*idx1(n2)+1:end)];
                tP(tiC)=false;
            end
        end
        
        iC=[iC,tiC];
    end
end
end
%% neighb_soc
function B = neighb_soc(A,V,nA)
% A is the label matrix BoardCluster, V is board, and nA is the number of
% clusters
[m,n] = size(A);
B=sparse(A(:,1:n-1),A(:,2:n),double(V(:,1:n-1)<V(:,2:n)),nA,nA)...
    +sparse(A(:,2:n),A(:,1:n-1),double(V(:,2:n)<V(:,1:n-1)),nA,nA)...
    +sparse(A(1:m-1,:),A(2:m,:),double(V(1:m-1,:)<V(2:m,:)),nA,nA)...
    +sparse(A(2:m,:),A(1:m-1,:),double(V(2:m,:)<V(1:m-1,:)),nA,nA);
B(1,:) = 0;
B(:,1) = 0;
B=B>0;
end