Diffing "just checking" and "checking just checking"

Title: just checking checking just checking
Author: Alfonso Nieto-Castanon Amitabh Verma
Submitted: 2011-11-05 04:10:11 UTC 2011-11-05 05:06:45 UTC
Status: Passed Passed
Score: 658097.0 659388.0
Result: 65778395 (cyc: 33, node: 7406) 65908024 (cyc: 33, node: 7406)
CPU Time: 63.708 57.541
Code:
function [moves, vine] = solver(board, limit)
[moves,vine,scr] = alfi(board,limit);
if sum(board(:))-scr > 1e4
if sum(board(:))-scr > 1e5
    [moves_a,vine_a,scr_a] = nick(board,limit);    
    if scr_a>scr
        moves = moves_a;
        vine = vine_a;
    end
end

end

%% nick
function [moves, vine,scr] = nick(board, limit)
% Added filter to crosseyed
%
% N. Howe

doclusters=any(any(board(:,2:end)==board(:,1:end-1)))||any(any(board(1:end-1,:)==board(2:end,:))); %%%%
if doclusters, [vine1,scr1] = kurt(board);
else vine1=[];scr1=-inf; 
end

%moves1 = [];
[moves2,vine2,scr2] = pusher(board,limit);
[nrow,ncol] = size(board);
ntot = nrow*ncol;
[moves3,vine3,scr3] = pusher(board',limit);
[moves4,vine4,scr4] = pusher(rot90(board,2),limit);
[moves5,vine5,scr5] = pusher(rot90(board,3),limit);
[moves6,vine6,scr6] = dofilter(board,limit);
[scr,bestpusher]=max([scr1 scr2 scr3 scr4 scr5 scr6]);
switch bestpusher
    case 1
        moves = [];
        vine = vine1;
    case 2
        moves = moves2;
        vine = vine2;
    case 3
        id = reshape(1:ntot,nrow,ncol);
        id3 = id';
        moves = id3(moves3);
        vine = id3(vine3);
    case 4
        id = reshape(1:ntot,nrow,ncol);
        id4 = rot90(id,2);
        moves = id4(moves4);
        vine = id4(vine4);
    case 5
        id = reshape(1:ntot,nrow,ncol);
        id5 = rot90(id,3);
        moves = id5(moves5);
        vine = id5(vine5);
    case 6
        moves = moves6;
        vine = vine6;
end;
end

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

moves = zeros(limit,2);
[nrow,ncol] = size(board);
if (limit < nrow)
    scr = -inf;
    moves = [];
    vine = [];
    return;
end;
ntot = nrow*ncol;
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));
        while (cx~=tx)||(cy~=ty)
            dx = sign(tx-cx);
            dy = sign(ty-cy);
            if (dx~=0)&&(dy~=0)
                if cboard(cy,cx+dx)<=cboard(cy+dy,cx)
                    dy = 0;
                else
                    dx = 0;
                end;
            end;
            if (rank(cy+dy,cx+dx)>0)
                vx = cx+dx;
                vy = cy+dy;
                if (dx==0)
                    if (vx>1)&&(cboard(vy,vx)>cboard(vy,vx-1))&&((vx==ncol)||(cboard(vy,vx+1)>cboard(vy,vx-1)))
                        % move left
                        steps = steps+1;
                        if (steps > limit)
                            break;
                        end;
                        moves(steps,1) = pos(vy,vx);
                        moves(steps,2) = pos(vy,vx-1);
                        cboard(vy,vx-1) = cboard(vy,vx);
                        cboard(vy,vx) = 0;
                        if (rank(vy,vx-1)>0)
                            remain(rank(vy,vx-1)) = 0;
                        end;
                        rank(vy,vx-1) = rank(vy,vx);
                        rank(vy,vx) = 0;
                        %imagesc(cboard);
                        %axis image;
                        %drawnow;
                    elseif (vx < ncol)&(cboard(vy,vx)>cboard(vy,vx+1))
                        % move right
                        steps = steps+1;
                        if (steps > limit)
                            break;
                        end;
                        moves(steps,1) = pos(vy,vx);
                        moves(steps,2) = pos(vy,vx+1);
                        cboard(vy,vx+1) = cboard(vy,vx);
                        cboard(vy,vx) = 0;
                        if (rank(vy,vx+1)>0)
                            remain(rank(vy,vx+1)) = 0;
                        end;
                        rank(vy,vx+1) = rank(vy,vx);
                        rank(vy,vx) = 0;
                        %imagesc(cboard);
                        %axis image;
                        %drawnow;
                    else
                        remain(rank(cy+dy,cx+dx)) = 0;
                    end;
                else
                    if (vy>1)&(cboard(vy,vx)>cboard(vy-1,vx))&((vy==nrow)|(cboard(vy+1,vx)>cboard(vy-1,vx)))
                        % move north
                        steps = steps+1;
                        if (steps > limit)
                            break;
                        end;
                        moves(steps,1) = pos(vy,vx);
                        moves(steps,2) = pos(vy-1,vx);
                        cboard(vy-1,vx) = cboard(vy,vx);
                        cboard(vy,vx) = 0;
                        if (rank(vy-1,vx)>0)
                            remain(rank(vy-1,vx)) = 0;
                        end;
                        rank(vy-1,vx) = rank(vy,vx);
                        rank(vy,vx) = 0;
                        %imagesc(cboard);
                        %axis image;
                        %drawnow;
                    elseif (vy < nrow)&(cboard(vy,vx)>cboard(vy+1,vx))
                        % move right
                        steps = steps+1;
                        if (steps > limit)
                            break;
                        end;
                        moves(steps,1) = pos(vy,vx);
                        moves(steps,2) = pos(vy+1,vx);
                        cboard(vy+1,vx) = cboard(vy,vx);
                        cboard(vy,vx) = 0;
                        if (rank(vy+1,vx)>0)
                            remain(rank(vy+1,vx)) = 0;
                        end;
                        rank(vy+1,vx) = rank(vy,vx);
                        rank(vy,vx) = 0;
                        %imagesc(cboard);
                        %axis image;
                        %drawnow;
                    else
                        remain(rank(cy+dy,cx+dx)) = 0;
                    end;
                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+dy,cx+dx) = i; %cboard(cy,cx);
            rank(cy,cx) = 0;
            cx = cx+dx;
            cy = cy+dy;
            %imagesc(cboard);
            %axis image;
            %drawnow;
        end;
        if (steps > limit)
            break;
        end;
        rank(ty,tx) = i;
        nextfit = nextfit+1;
    end;
end;
moves = moves(1:min(steps,limit),:);
[vine,scr] = kurt(cboard);
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
        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;
        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;
        %imagesc(board);
        %axis image;
        %drawnow;
    end;
    [vine,scr] = sebastian(board);
else
    scr = -inf;
    moves = [];
    vine = [];
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);
end
% post-processing moves
if limit>0
    [board,iC,moves,limit] = postprocess(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

function [board,tiC,moves,limit] = postprocess(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

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;
% D=Dijkstra(tP,tiC(1));
% maskout=tP&~isinf(D)&board>board(tiC(1));
% nmaskout=sum(D(maskout));
for ndir=1:4
    K=size(tiC,2);
    switch(ndir)
        case 1, 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)));
        case 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)));
        case 3, 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)));
        case 4, 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)));
    end
    %for n2=1:numel(idx1),
    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');
        %idxZ=find(~tP(tempA)|maskout(tempA),1,'first');
        switch(ndir)
            case {1,2},
                %idxA=find(board(tiC(2*idx1(n2)-1))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2))),1,'first');
                [nill,idxA]=max((1./(abs(board(tiC(2*idx1(n2)-1))-board(tempA))+1)).*(board(tiC(2*idx1(n2)-1))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)))));
                if ~nill,idxA=[]; end
            case {3,4},
                %idxA=find(board(tiC(2*idx1(n2)))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)+1)),1,'first');
                [nill,idxA]=max((1./(abs(board(tiC(2*idx1(n2)))-board(tempA))+1)).*(board(tiC(2*idx1(n2)))>=board(tempA)&board(tempA)>=board(tiC(2*idx1(n2)+1))));
                if ~nill,idxA=[]; end
        end
        if ~isempty(idxA)&&idxA<idxZ
            idxZ=find(~tP(tempB),1,'first');
            %idxZ=find(~tP(tempB)|maskout(tempA),1,'first');
            switch(ndir)
                case {1,2},
                    idxB=find(board(tempA(idxA))>=board(tempB)&board(tempB)>=board(tiC(2*idx1(n2))),1,'first');
                case {3,4},
                    idxB=find(board(tempA(idxA))>=board(tempB)&board(tempB)>=board(tiC(2*idx1(n2)+1)),1,'first');
            end
            if ~isempty(idxB)&&idxB<idxZ&&idxA+idxB-2<=limit,%-nmaskout,%.0*boardab(tiC(1))*(n+m)/2, %%%
                ok=1;
                switch(ndir)
                    case 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)];
                    case 2, 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)];
                    case 3, 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)];
                    case 4, 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)];
                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);
                %break;
            end
        end
    end
end
end

function [board,moves,tiC,tP,limit,ok]=postprocess_endpoint(board,moves,tiC,tP,limit,type)
[m,n]=size(board);
ok=1;
switch(type)
    case 1
        mask=tP&board<2*board(tiC(1));
        D=Dijkstra(mask,tiC(1));
        idx1=find(mask&D-1<=limit&board>=board(tiC(1)));
        if isempty(idx1),
            mask=tP;
            D=Dijkstra(mask,tiC(1));
            idx1=find(mask&D-1<=limit&board>=board(tiC(1)));
            if isempty(idx1),ok=0;end
        end
        if ok
            [nill,idx2]=min(board(idx1)+1e-3*D(idx1)); %%%t3
            idx0=idx1(idx2);
            limit=limit-(D(idx0)-1);
            for n2=D(idx0)-1:-1:1
                if      D(idx0+1)==n2, tidx0=idx0+1;
                elseif  D(idx0-1)==n2, tidx0=idx0-1;
                elseif  D(idx0+m)==n2, tidx0=idx0+m;
                elseif  D(idx0-m)==n2, tidx0=idx0-m;
                end
                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));
        idx1=find(mask&D-1<=limit&board<=board(tiC(end))&board>0);
        if isempty(idx1),ok=0;end
        if ok
            [nill,idx2]=min(-board(idx1)+1e-3*D(idx1)); %%%t3
            idx0=idx1(idx2);
            limit=limit-(D(idx0)-1);
            for n2=D(idx0)-1:-1:1
                if      D(idx0+1)==n2, tidx0=idx0+1;
                elseif  D(idx0-1)==n2, tidx0=idx0-1;
                elseif  D(idx0+m)==n2, tidx0=idx0+m;
                elseif  D(idx0-m)==n2, tidx0=idx0-m;
                end
                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_pixel(board,limit,boardab)
[m,n]=size(board);
cD=board;
IDX=zeros(m,n);
touched=board>0;
valid=touched;
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(:));
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

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);
        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]=size(A);
mn=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;
    % 	if A(i1+1); i1=i1+1;
    %     elseif A(i1-1); i1=i1-1;
    %     elseif A(i1+m); i1=i1+m;
    %     elseif A(i1-m); i1=i1-m;
    %     end
end
D(i1)=1;
for n1=2:mn,
    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

function [vine,scr] = kurt(board)
doclusters=any(any(board(:,2:end)==board(:,1:end-1)))||any(any(board(1:end-1,:)==board(2:end,:))); 
if doclusters, 
    [vine,scr] = real_kurt(board);
else
    [m,n]=size(board);
    board=[zeros(1,n+2);[zeros(m,1),board,zeros(m,1)];zeros(1,n+2)];
    vine=bellman_pixel(board);
    [iC1,iC2]=ind2sub([m+2,n+2],vine);
    vine=sub2ind([m,n],iC1-1,iC2-1);
    board=board(2:end-1,2:end-1);
    scr=sum(board(vine));
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));

inhoud = unique(A)';
t=sort(A(:));
tt = cumsum(t);
ttt = find(diff([t; 0]));
tttt = tt(ttt);
t2 = sparse(0);
t2(inhoud+1) = tttt;
H = full(t2(A+1));

if any(size(H) ~= size(A))
    H=H';
end


B=A;


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

Hmod = zeros(m,n);


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)
        verplaatsindex = IND(i-1,j);
        if ~any(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];
                %                 Hmod(verplaatsindex) = sum(A(C{verplaatsindex})==A(verplaatsindex));
            end
        end
    end
    
    if i<size(A,1) && A(i+1,j)<=A(i,j)
        verplaatsindex = IND(i+1,j);
        if ~any(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];
                %                 Hmod(verplaatsindex) = sum(A(C{verplaatsindex})==A(verplaatsindex));
            end
        end
    end
    if j<size(A,2) && A(i,j+1)<=A(i,j)
        verplaatsindex = IND(i,j+1);
        if ~any(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];
                %                 Hmod(verplaatsindex) = sum(A(C{verplaatsindex})==A(verplaatsindex));
            end
        end
    end
    if 1<j && A(i,j-1)<=A(i,j)
        verplaatsindex = IND(i,j-1);
        if ~any(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];
                %                 Hmod(verplaatsindex) = sum(A(C{verplaatsindex})==A(verplaatsindex));
            end
        end
    end
    
    updated(i,j)=false;
    
    
    
end

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

end

function [vine, score] = sebastian(board)

[M,N] = size(board);
OFF = M*N+1;

iboard = reshape(1:M*N,[M,N]);
 
n = [OFF*ones(1,N); iboard(1:end-1,:)];
s = [iboard(2:end,:); OFF*ones(1,N)];

board = [board(:);0];
w = [OFF*ones(M,1);(1:M*N-M)';OFF];
e = [(M+1:M*N)';OFF*ones(M,1);OFF];
n = [n(:);OFF];
s = [s(:);OFF];

w(board<board(w)) = OFF;
e(board<board(e)) = OFF;
n(board<board(n)) = OFF;
s(board<board(s)) = OFF;
wens = [w(:) e(:) n(:) s(:)];

[B,IX] = sort(board(:));

lastMove = nan(size(board));
sumMoves = zeros(size(board));

for i=IX'
    
    [C,I] = max(sumMoves(wens(i,:)));
    lastMove(i) = wens(i,I);
    sumMoves(i) = board(i) +  C;
    
end

[score,endPoint] = max(sumMoves);

vine = endPoint;
while sumMoves(vine(end))
    vine(end+1) = lastMove(vine(end));
end
vine = vine(end-1:-1:1);

end