Winner GUO (asdf1)

Finish 2004-11-10 09:00:00 UTC

Lean8/2

by Niilo Sirola

Status: Passed
Results: 57309
CPU Time: 36.0634
Score: 57.355
Submitted at: 2004-11-09 16:28:54 UTC
Scored at: 2004-11-09 16:34:33 UTC

Current Rank: 841st
Based on: Lean8++ (diff)

Comments
Please login or create a profile.
Code
function moves = solver(Ai,Af,w)

global Dperfect

w = 1+abs(w);
n = length(w);
ipos=zeros(n,2);
fpos=zeros(n,2);

for i=1:n
    [ipos(i,1) ipos(i,2)] = find(Ai==i);
    [fpos(i,1) fpos(i,2)] = find(Af==i);
end

% optimum number of moves
optmove = sum(sum(abs(fpos-ipos)));

% cost of the optimum move
optcost = sum(w.*sum(abs(fpos-ipos),2));

Dperfect = optmove;

% try different solvers with fastest first

moves = zeros(0,2);
res = inf;

if all(Ai(:)==Af(:))
    return
end

% try solver_linspac
moves2 = solver_linspac(Ai,Af,w);
res2 = grade2(Ai,Af,w,moves2);
    
if res2<res
    moves = moves2;
    res = res2;
end

if res-optcost < 100
    return
end

% try backwards
moves2 = reversemov(solver_linspac(Af,Ai,w));
res2 = grade2(Ai,Af,w,moves2);

if res2<res
    moves = moves2;
    res = res2;
end

if res-optcost < 100
    return
end

% try solver_3
moves2 = solver_3(Ai,Af,w);
res2 = grade2(Ai,Af,w,moves2);
    
if res2<res
    moves = moves2;
    res = res2;
end

if res-optcost < 100
    return
end

% try backwards
moves2 = reversemov(solver_3(Af,Ai,w));
res2 = grade2(Ai,Af,w,moves2);

if res2<res
    moves = moves2;
    res = res2;
end

if res-optcost < 100
    return
end

% try solver_SFAL
moves2 = solver_SFAL(Ai,Af,w);
res2 = grade2(Ai,Af,w,moves2);

if res2<res
    moves = moves2;
    res = res2;
end

if res-optcost < 100
    return
end

% try backwards
moves2 = reversemov(solver_SFAL(Af,Ai,w));
res2 = grade2(Ai,Af,w,moves2);

if res2<res
    moves = moves2;
    res = res2;
end

% not solved yet? then get desperate
while ~isfinite(res)
    moves = solver_1(Ai,Af,w);
    res = grade2(Ai,Af,w,moves);
end

%--------------------------------
function mov = reversemov(mov)

o = [3 4 1 2]'; % for turning the moves backwards
if isempty(mov)
    return
end

mov = mov(end:-1:1,:);
mov(:,2)=o(mov(:,2));


function score = grade2(ai,af,w,mv)


if isempty(mv) 
    score = inf;
    return
end

% N=1, E=2, S=3, W=4
I = [0  1  0 -1];
J = [1  0 -1  0];

% First column of moves should be integer values from 1 to N
if any(round(mv(:,1)) ~= mv(:,1)) || any(mv(:,1) < 1) || any(mv(:,1) > length(w))
    score = inf;
    return
end

% Second column of moves should be integer values from 1 to 4
if any(round(mv(:,2)) ~= mv(:,2)) || any(mv(:,2) < 1) || any(mv(:,2) > 4)
    score = inf;
    return
end

% Step through the moves and make sure there are no illegal moves and
% verify that the final location is the specified location
a = ai;
for k = 1:size(mv,1)
    [row, col] = find(a==mv(k,1));
    a(row,col) = 0;
    % new position of the block
    row = row + I(mv(k,2));
    col = col + J(mv(k,2));
    % check edges
    if (row < 1) || (col < 1) || (row > size(a,1)) || (col > size(a,2))
        score = inf;
        return
    end
    % check collisions
    if (a(row,col) ~= 0)
        score = inf;
        return
    end
    % place the block in the new position
    a(row, col) = mv(k,1);
end

% Make sure we end up in the correct location
if any(a(:) ~= af(:))
        score = inf;
        return
end

% compute the score
score = sum(w(mv(:,1)));


% -----------------------------------
% Niilo Sirola's solver
function [moves,ok] = solver_SFAL(A,B,w0)

ym = size(A,1);
n = length(w0);

M = [ym 1 -ym -1];

w = abs(w0(:)')+.1;

[mov, ok] = mainsolver(A,B,w);

if ok
    mov = 1*(mov==M(1)) + 2*(mov==M(2)) + 3*(mov==M(3)) + 4*(mov==M(4));
    moves = [(mov>0)*(1:n)' sum(mov,2)];
else
    moves = [];
end


function [mov, ok] = mainsolver(A,B,w)

global Dperfect

ym = size(A,1);
M = [ym 1 -ym -1];
n = length(w);
mov = [];

mask = true(1,n);

ipos = zeros(1,n);
fpos = zeros(1,n);
for i=1:n
    ipos(i) = find(A==i);
    fpos(i) = find(B==i);
end

perfmoves = abs(rem(fpos-ipos,ym)) + abs(floor(fpos-ipos)/ym);

invw = max(w)+1-w;

blockers = zeros(size(w));

for leaveout = [0 ceil(2*log(n)):(n-1)]

    mmask = mask;
    choosew = sum(blockers).*invw;
    for i = 1:leaveout
        randi = rand*sum(mmask.*choosew);
        choose = find(cumsum(mmask.*choosew)>=randi,1,'first');
        mmask(choose) = false;
    end

    Ap = zeros(size(A));
    Bp = zeros(size(B));
    for i = find(mmask)
        Ap(ipos(i)) = i;
        Bp(B==i) = i;
    end

    % form a subproblem
    ind2 = find(mmask);
    A2 = zeros(size(A));
    B2 = zeros(size(B));
    w2 = zeros(sum(mmask),1);
    for i=1:length(ind2)
       A2(ipos(ind2(i))) = i;
       B2(B==ind2(i)) = i;
       w2(i) = w(ind2(i));
    end

    Dperfect = sum(perfmoves(mmask));
    % solve the subproblem
    [mov2,ok] = solver_3(A2,B2,w2);

    % convert moves to my format
    nmov = size(mov2,1);
    partialmov = zeros(nmov,n);
    for i=1:nmov
        partialmov(i,ind2(mov2(i,1))) = M(mov2(i,2));
    end
    
    blockers = blockers + sum(findoverlaps(cumsum([ipos;partialmov])));

    if ~ok
        continue
    end

    for i = find(mask)
        Ap(ipos(i)) = i;
    end

    for i = 1:3*sum(mask)
        randi = rand*sum(w.*(~mmask & mask));
        choose = find(cumsum(w.*(~mmask & mask))>=randi,1,'first');
        [trymov, ok] = imoves(choose,partialmov,Ap,B,w,mmask);
        if ok
            partialmov = trymov;
            mmask(choose) = true;
            if all(mmask==mask)
                ok = 1;
                mov = partialmov;
                return
            end
        else
            dropw = blockers.*invw;
            randi = rand*sum(dropw.*(mmask&mask));
            choose = find(cumsum(dropw.*(mmask&mask))>=randi,1,'first');
            mmask(choose) = false;
            partialmov(find(partialmov(:,choose)),:) = [];
        end
    end
    ok = 0;
end

function [mov, ok] = imoves(c,mov,A,B,w,mask)

mov(find(mov(:,c)),:) = [];
nmov = size(mov,1);
n = length(mask);
[crow,ccol] = find(A==c);
[trow,tcol] = find(B==c);

if nmov == 0 && crow==trow && ccol==tcol
    ok = 1;
    return
end

[rows,cols] = size(A);
slices = false(rows, cols, nmov+1);
ipos = zeros(1,n);

for i=find(mask)
    ipos(i) = find(A(:)==i);
end

pos = cumsum([ipos; mov],1);

for k=1:(nmov+1)
    R = zeros(size(A));
    R(pos(k,mask)) = -1;
    slices(:,:,k) = (R==-1);
end

target = sub2ind(size(slices),trow,tcol,nmov+1);
optlen = abs(crow-trow) + abs(ccol-tcol);

[path] = dijkstra(numel(slices),find(A==c), target, slices, optlen);

if isempty(path)
    ok = 0;
    return
end

for i=length(path):-1:2
    if path(i)-path(i-1) < rows*cols
        k = floor((path(i-1)-1)/(rows*cols))+1;
        mov = [mov(1:(k-1),:); zeros(1,n); mov(k:end,:)];
        mov(k,c) = path(i)-path(i-1);
    end
end

ok = 1;


function overlaps = findoverlaps(pos)

[npos,n] = size(pos);
cols = find(all(pos));
pos = pos(:,cols);

A = zeros(npos,length(cols));

[sortpos, ind] = sort(pos,2);
ind1 = ind(:,1:end-1);
ind2 = ind(:,2:end);
dp = diff(sortpos,1,2)==0;

for i=1:npos
    A(i,ind1(i,dp(i,:))) = true;
    A(i,ind2(i,dp(i,:))) = true;
end

overlaps = zeros(npos,n);
overlaps(:,cols) = A;

function [path] = dijkstra(n, s, d, R, optlen)

RELAX = 1;

[rows,cols,depth] = size(R);
R = R(:);

visited = false(1,n);

distance = inf*ones(1,n);
parent = zeros(1,n);

distance(s) = 0;

stack = zeros(1,n);
stack(1)=s;
next = 1;
last = 1;

for i = 1:(n-1),
    if next>last
        break
    else
        u = stack(next);
        next = next + 1;
    end

    visited(u) = true;
    ndx = u - 1;
    k = floor(ndx/(rows*cols)) + 1;
    ndx = rem(ndx, rows*cols);
    col = floor(ndx/rows) + 1;
    ndx = rem(ndx, rows);
    row = floor(ndx) + 1;
    v = u + 1;
    if v>0 && v<=n && ~R(v) && row<rows
        if distance(u) + 1 < distance(v)
            distance(v) = distance(u) + 1;
            parent(v) = u;
            if ~visited(v)
                if (v==d) && (distance(d) <= RELAX*optlen)
                    break
                end
                last = last + 1;
                stack(last) = v;
            end
        end
    end
    v = u - 1;
    if v>0 && v<=n && ~R(v) && row>1
        if distance(u) + 1 < distance(v)
            distance(v) = distance(u) + 1;
            parent(v) = u;

            if ~visited(v)
                if (v==d) && (distance(d) <= RELAX*optlen)
                    break
                end
                last = last + 1;
                stack(last) = v;
            end
        end
    end
    v = u + rows;
    if v>0 && v<=n && ~R(v) && col<cols
        if distance(u) + 1 < distance(v)
            distance(v) = distance(u) + 1;
            parent(v) = u;
            if ~visited(v)
                if (v==d) && (distance(d) <= RELAX*optlen)
                    break
                end
                last = last + 1;
                stack(last) = v;
            end
        end
    end
    v = u - rows;
    if v>0 && v<=n && ~R(v) && col>1
        if distance(u) + 1 < distance(v)
            distance(v) = distance(u) + 1;
            parent(v) = u;
            if ~visited(v)
                if (v==d) && (distance(d) <= RELAX*optlen)
                    break
                end
                last = last + 1;
                stack(last) = v;
            end
        end
    end

    v = u + rows*cols;
    if v>0 && v<=n && ~R(v) && k < depth
        if distance(u) < distance(v)
            distance(v) = distance(u);
            parent(v) = u;
            if ~visited(v)
                if (v==d) && (distance(d) <= RELAX*optlen)

                    break
                end
                last = last + 1;
                stack(last) = v;
            end
        end
    end
end

path = [];
if parent(d) ~= 0
    t = d;
    path = d;
    while t ~= s
        p = parent(t);
        path = [p path];
        t = p;
    end
end



% ----------------------------------------------
% Stijn Helsen's solver
function [mv,perfectMV] = solver_linspac(ai,af,w,states)
% (argument 'states' not used)

global Ac Ar m2 Dperfect

if nargin<4
    states=1111;
end

nBlocks = length(w);
[m,n] = size(ai);

m2=m+2;
n2=n+2;

A = -ones(m2,n2);
Af = A;
A( 2:m+1, 2:n+1 ) = ai;
Af( 2:m+1, 2:n+1 ) = af;

Ac=1:n2;
Ac=Ac(ones(m2,1),:);
Ar=(1:m2)';
Ar=Ar(:,ones(n2,1));

Pi=w;
Pf=w;
for i=m2+2:numel(A)-m2-1
    if A(i)>0, Pi(A(i)) = i; end
    if Af(i)>0, Pf(Af(i)) = i; end
end

P=Pi;

nmv=1;
mv = zeros(500,2);
nNOK=sum(P~=Pf);

Paths=zeros(m+n,nBlocks);
lPaths=w;
fPaths=w;
Pend=w;
bOK=w;
obs=zeros(nBlocks,2);

nmv0=0;
while nmv0<nmv&&nNOK
    obs(:)=0;
    nmv0=nmv;
    for i=1:nBlocks
        if P(i)==Pf(i)
            lPaths(i)=0;
            fPaths(i)=0;
            bOK(i)=1;
        else
            [P1,f1]=SearchPath(A,P(i),Pf(i));
            if isempty(P1)
                lPaths(i)=0;
                fPaths(i)=0;
                obs(i,1:length(f1))=f1;
            else
                if isempty(f1)
                    fPaths(i)=1;
                else
                    fPaths(i)=0;
                    obs(i,1:length(f1))=f1;
                end
                lPaths(i)=length(P1);
                Paths(1:lPaths(i),i)=P1;
                Pend(i)=P1(end);
            end
            bOK(i)=0;
        end
    end

    iP=find(~bOK&lPaths);
    PCol=zeros(length(iP));
    L=lPaths(iP);
    for i=1:length(iP)

        Pe=Pend(iP(i));
        for j=1:length(iP)
            if i~=j
                lj=L(j);
                PCol(i,j)=any(Paths(1:lj,iP(j))==Pe);
            end
        end
    end
    sPCol=sum(PCol,2);
    pOK=find(sPCol==0);
    pNOK=find(sPCol~=0);
    if isempty(pOK)
        if length(pNOK)==1
            pOK(end+1)=pNOK;
        elseif ~isempty(pNOK)
            if length(pNOK)>1&&any(fPaths(iP(pNOK)))
                pNOK(~fPaths(iP(pNOK)))=[];
            end
            iNOK1=pNOK(find(sPCol(pNOK)==1));
            for i=iNOK1'
                j=find(PCol(i,:));
                jj=iP(j);
                p=iP(i);
                pe=Pend(p);
                A(P(p))=0;
                A(pe)=p;
                [P1,f1]=SearchPath(A,P(jj),Pf(jj));
                A(P(p))=p;
                A(pe)=0;
                if isempty(f1)
                    pOK(end+1)=i;
                    break
                end
            end
        end
    end
    if length(pOK)>1
        obs1=abs(obs(:));
        temp = zeros(1,length(obs1));
        i=1;
        q=0;
        pOK1=pOK;
        while i<=length(obs1)
            temp(obs1==obs1(i))=1;
            if sum(temp)>1
                j=find(obs1==obs1(i));
                obs1(j(2:end))=[];
            end
            if any(obs1(i)==pOK)
                q=1;
                j=find(obs1(i)==pOK);
                pOK1(j)=0;
            end
            i=i+1;
        end
        if q
            pOK(pOK1~=0)=[];
        end
        if length(pOK)>1&&any(fPaths(iP(pOK)))
            pOK(~fPaths(iP(pOK)))=[];
        end
    end
    j=1;
    while j<=length(pOK)
        i=pOK(j);
        b=iP(i);
        k=nmv+1:nmv+L(i);
        mv(k,1)=b;
        mv(k,2)=Paths(1:L(i),b);
        nmv=nmv+L(i);
        A(P(b))=0;
        A(Pend(b))=b;
        P(b)=Pend(b);
        if fPaths(b)
            nNOK=nNOK-1;
        end
        j=j+1;
    end
end

P=Pi;
for i=2:nmv
    b=mv(i);
    p=mv(i,2);
    dp=p-P(b);
    if dp==1
        mv(i,2)=2;
    elseif dp==-1
        mv(i,2)=4;
    elseif dp==m2
        mv(i,2)=1;
    else
        mv(i,2)=3;
    end
    P(b)=p;
end
mv=mv(2:nmv,:);

if nNOK
    perfectMV = 0;
else
    perfectMV = 1;
end


function [P,stopped]=SearchPath(A, p1, p2)
global Ac Ar m2
c1=Ac(p1);c2=Ac(p2);
r1=Ar(p1);r2=Ar(p2);

stopped=[];
P=zeros(sum(size(A)),1);
if r1>r2
    d_r=-1;
elseif r1<r2
    d_r=1;
else
    d_r=0;
end
if c1>c2
    d_c=-1;
elseif c1<c2
    d_c=1;
else
    d_c=0;
end
n=0;
p=p1;
if d_r==0||d_c==0
    while p~=p2
        np=p+d_c*m2+d_r;
        if A(np)
            stopped=A(np);
            break
        end
        p=np;
        n=n+1;
        P(n)=p;
    end
else



    Ah=A;
    c=c2;
    i1=p2-d_r;
    r_1=r2-d_r;
    while c~=c1
        r=r_1;
        r_1=0;
        i2=i1;
        while r~=r1
            if Ah(i2)
                Afil=-Ah(i2);
                if r==r2

                    r_1=r2;
                    i1=i2-d_c*m2;
                else
                    r_1=r+d_r;
                    i1=i2-d_c*m2+d_r;
                end
                r=r-d_r;
                i2=i2-d_r;
                while r~=r1
                    if Ah(i2)==0
                        Ah(i2)=Afil;
                    end
                    r=r-d_r;
                    i2=i2-d_r;
                end
                if Ah(i2)==0
                    Ah(i2)=Afil;
                end
                break;
            end
            r=r-d_r;
            i2=i2-d_r;
        end
        if r_1==0&&Ah(i2)
            i1=i2-d_c*m2+d_r;
            r_1=r+d_r;
        end
        c=c-d_c;
        if r_1==0
            break;
        end
    end
    r=r2;

    i1=p2-d_c*m2;
    c_1=c2-d_c;
    while r~=r1
        c=c_1;
        c_1=0;
        i2=i1;
        while c~=c1
            if Ah(i2)
                if Ah(i2)<0
                    Afil=Ah(i2);
                else
                    Afil=-Ah(i2);
                end
                if c==c2
                    c_1=c2;
                    i1=i2-d_r;
                else
                    c_1=c+d_c;
                    i1=i2-d_r+d_c*m2;
                end
                c=c-d_c;
                i2=i2-d_c*m2;
                while c~=c1
                    if Ah(i2)==0
                        Ah(i2)=Afil;
                    end
                    c=c-d_c;
                    i2=i2-d_c*m2;
                end
                if Ah(i2)==0
                    Ah(i2)=Afil;
                end

                break;
            end
            c=c-d_c;
            i2=i2-d_c*m2;
        end
        if c_1==0&&Ah(i2)
            i1=i2-d_r+d_c*m2;
            c_1=c+d_c;
        end
        r=r-d_r;
        if c_1==0
            break;
        end
    end

    while p~=p2
        if c1~=c2&&Ah(p+m2*d_c)==0
            di=m2*d_c;
            dc=d_c;
            dr=0;
        elseif r1~=r2&&Ah(p+d_r)==0
            di=d_r;
            dc=0;
            dr=d_r;
        else
            if c1==c2
                stopped=Ah(p+d_r);
            elseif r1==r2
                stopped=Ah(p+m2*d_c);
            else
                stopped=[Ah(p+d_r) Ah(p+m2*d_c)];
            end
            break;
        end
        p=p+di;
        n=n+1;
        P(n)=p;
        r1=r1+dr;
        c1=c1+dc;
    end
end
P=P(1:n);





% solver by wu zhili
% ---------------------------------------
function [mv,ok] = solver_3(ai,af,w,states)
global mv a airow aicol afrow afcol rowDis colDis I J b m n wn Dperfect

if nargin<4
    states=1111;
end

mv = [];
aiMap = ai>0;
a = ai;
[m,n] = size(aiMap);
aiBlock = ai(aiMap);
nBlocks = length(aiBlock);

[aiOrder, aiPos] = sort(aiBlock);
[airow,aicol] = find(aiMap);
airow = airow(aiPos)';
aicol = aicol(aiPos)';
w0 = w;
w = w(aiPos);
wn = w;

afMap = af>0;
[afOrder, afPos ] = sort( af(afMap) );
[afrow,afcol] = find(afMap);
afrow = afrow(afPos)';
afcol = afcol(afPos)';
rowDis = airow-afrow;
colDis = aicol-afcol;
rowDis0 = rowDis;
colDis0 = colDis;
[nw,nd]=sort(-w);
ds = sum(abs(w.*rowDis') + abs(w.*colDis'));
I = [0  1  0 -1];
J = [1  0 -1  0];
recurP = [1:4; 3 4 1 2; 1:4];
iter = 0;
while ( sum(abs(rowDis) + abs(colDis))>0 ) > 0 && iter<2
    iter = iter + 1;
    for irr = 1:nBlocks
        ir = nd(irr);
        priority = 0;
        len=0;
        iterr = 0;
        while abs(rowDis(ir))+abs(colDis(ir)) > 0 && iterr < 10
            iterr = iterr+1;
            b=ai*0;
            x = airow(ir);
            y = aicol(ir);
            b(x,y) = ir;

            if ~chainMV( ir,priority )
                break;
            end;
            len = len+1;

            if len>3
                if all( mv(end:-1:end-2,1) == ir ) && ...
                        any( all( repmat( mv(end:-1:end-2,2), 1, 4 ) == recurP ) )
                    priority = Inf;
                    mv(end-1:end,:) = [];
                end;
            end
        end
    end
end
res = sum(abs(rowDis) + abs(colDis));
if res > 0
    ok = false;
else
    ok = true;
end


function isMovable = chainMV( ir , priority)
global mv a airow aicol afrow afcol rowDis colDis I J b m n wn
isMovable = 0;
x = airow(ir);
y = aicol(ir);
tx = afrow(ir);
ty = afcol(ir);
disr = rowDis(ir);
disc = colDis(ir);
if ( abs( disr ) + abs( disc ) > 0 ) || ( priority > wn(ir) )

    mvDirects = ( disr * I + disc * J );
    [goodMv, mvInd] = sort(mvDirects);
    for k = 1:2

        i = I(mvInd(k));
        j = J(mvInd(k));
        nx = x + i;
        ny = y + j;
        succMv = 0;
        if (nx>0) && (nx<=m) && (ny>0) && (ny<=n)
            if a(nx,ny) == 0
                succMv = 1;
            else
                if b(nx,ny) == 0
                    b(nx,ny) = 1;
                    succMv = chainMV( a(nx,ny) , priority );
                else
                    b(nx,ny) = 0;
                end;
            end;
        end;
        if succMv == 1

            airow(ir) = nx;
            aicol(ir) = ny;
            rowDis(ir) = rowDis(ir) + I(mvInd(k));
            colDis(ir) = colDis(ir) + J(mvInd(k));
            a(x,y) = 0;
            a(nx,ny) = ir;
            mv(end+1,:) = [ir,mvInd(k)];
            isMovable = 1;
            break;
        end;
    end;
end;

% it would be nice to give credit of this solver to the original author
%-------------------------------------------------
function move = solver_1(aInit, aFinal, wt, bscore)

if nargin<4
    bscore = inf;
end

allmoves  = [];
wtinitori = wt;
aInitori  = aInit;

sortbydist   = 0;
sortbyweight = 1;

recur     = 24;
minw = 10000;

if sortbyweight
    [tmp inds] = sort(-wt);
    wtinit = wtinitori;
end;

lenwt = length(wt);
absx = zeros(lenwt,1);
absy = absx;
for index = 1:lenwt
    [hvix hviy] = find(aInit  == index);
    [hvfx hvfy] = find(aFinal == index);

    absx(index) = abs(hvix-hvfx);
    absy(index) = abs(hviy-hvfy);
end;

if sortbydist
    dist = abs(absx) + abs(absy);

    [tmp inds] = sort(-dist);
    wtinit = wtinitori;
end;

for index = 1:length(inds)
    hv = inds(index);
    [hvix hviy] = find(aInit  == hv);
    [hvfx hvfy] = find(aFinal == hv);

    dist = abs(hvix-hvfx)+abs(hviy-hvfy);

    wt(hv) = -Inf;

    wtinit(hv) = wtinit(hv)+minw+10000;
    [move cost aInit] = movefrompos(hv, [hvix hviy], [hvfx hvfy], aInit, aFinal, wtinit, 1, dist+5);

    if isinf(cost) || isempty(move)
    else
        move = [move(:,1) 3*abs(move(:,4))-move(:,4)+2*abs(move(:,5))-move(:,5) ];
        allmoves = [allmoves; move];
        if ( sum(wtinitori(allmoves(:,1)))>bscore )
            move=[];
            return;
        end
    end;
end;

count   = 1;
oldinds = [];
while ~all(aFinal(:) == aInit(:))

    sortbyweight = sortbydist;
    sortbydist   = ~sortbyweight;

    if sortbyweight
        inds = find(aFinal ~= aInit);
        inds = aFinal(inds);
        inds = inds(inds > 0);
        [tmp indices] = sort(-wt(inds));

        inds = inds(indices);
    end;

    if sortbydist
        for index = 1:length(wt)
            [hvix hviy] = find(aInit  == index);
            [hvfx hvfy] = find(aFinal == index);
            dist(index) = abs(hvix-hvfx)+abs(hviy-hvfy);
        end;
        [tmp inds] = sort(-dist);
        firstzeros = find(tmp == 0);
        inds(firstzeros:end) = [];
    end;

    wtinit = wtinitori;
    notmove = setdiff(oldinds, inds);
    for hv = 1:length(notmove)
        wtinit(notmove(hv)) = wtinit(notmove(hv))+minw+10000;
    end;


    for hind = 1:length(inds)
        hv = inds(hind);
        [hvix hviy] = find(aInit  == hv);
        [hvfx hvfy] = find(aFinal == hv);

        dist = abs(hvix-hvfx)+abs(hviy-hvfy);

        wtinit(hv) = wtinit(hv)+minw+10000;
        [move cost aInit] = movefrompos(hv, [hvix hviy], [hvfx hvfy], aInit, aFinal, wtinit, 1, dist+6);
        if hind > 1
            wtinit(inds(hind-1)) = wtinit(inds(hind-1))-minw-10000;
        end;
        if isinf(cost) || isempty(move)
        else
            move = [move(:,1) 3*abs(move(:,4))-move(:,4)+2*abs(move(:,5))-move(:,5) ];
            allmoves = [allmoves; move];
            if ( sum(wtinitori(allmoves(:,1)))>bscore )
                move=[];
                return;
            end
            move_ok=1;
        end;
    end;
    oldinds = inds;
    count = count+1;
end;

move = allmoves;

function [move, cost, curposnew] = movefrompos(item, pinit, pfin, curpos, finpos, wt, strategy, recur)

cost = 0;
curposnew = curpos;
if all(pinit == pfin), move = []; return; end;
if recur == 0, move = []; cost = 0; return; end;

diffx = pfin(1)-pinit(1);
diffy = pfin(2)-pinit(2);
dirx = sign(diffx);
diry = sign(diffy);
cost = Inf;

if (strategy < 10 && abs(diffx) > abs(diffy)) || (strategy == 11 && abs(diffx) < abs(diffy))
    if dirx && curpos(pinit(1)+dirx, pinit(2)) <= 0
        [move cost curposnew] = onemove(item, pinit, [dirx 0], pfin, curpos, finpos, wt, strategy, recur);
        if length(move) > 0 && cost/size(move,1) == mod(wt(curpos(pinit(1),pinit(2))),10000), return; end;
    end;
    if diry && curpos(pinit(1), pinit(2)+diry) <= 0
        [move1 cost1 curposnew1] = onemove(item, pinit, [0 diry], pfin, curpos, finpos, wt, strategy, recur);
        if length(move1) > 0 && cost1/size(move1,1) == mod(wt(curpos(pinit(1),pinit(2))),10000)
            move=move1;
            cost=cost1;
            curposnew=curposnew1;
            return;
        end;
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
else


    if diry && curpos(pinit(1), pinit(2)+diry) <= 0
        [move cost curposnew] = onemove(item, pinit, [0 diry], pfin, curpos, finpos, wt, strategy, recur);
        if length(move) > 0 && cost/size(move,1) == mod(wt(curpos(pinit(1),pinit(2))),10000), return; end;
    end;
    if dirx && curpos(pinit(1)+dirx, pinit(2)) <= 0
        [move1 cost1 curposnew1] = onemove(item, pinit, [dirx 0], pfin, curpos, finpos, wt, strategy, recur);
        if length(move1) > 0 && cost1/size(move1,1) == mod(wt(curpos(pinit(1),pinit(2))),10000)
            move=move1;
            cost=cost1;
            curposnew=curposnew1;
            return;
        end;
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
end;


if isinf(cost) || cost == 0
    if dirx && diry && curpos(pinit(1), pinit(2)+diry) > 0 && curpos(pinit(1)+dirx, pinit(2)) > 0
        [move1 cost1 curposnew1] = onemove(item, pinit, [dirx 0], pfin, curpos, finpos, wt, strategy, recur);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
        [move1 cost1 curposnew1] = onemove(item, pinit, [0 diry], pfin, curpos, finpos, wt, strategy, recur);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    else
        if dirx && curpos(pinit(1)+dirx, pinit(2)) > 0
            [move1 cost1 curposnew1] = onemove(item, pinit, [dirx 0], pfin, curpos, finpos, wt, strategy, recur);
            if cost1<cost
                move = move1;
                cost = cost1;
                curposnew = curposnew1;
            end
        end;
        if diry && curpos(pinit(1), pinit(2)+diry) > 0
            [move1 cost1 curposnew1] = onemove(item, pinit, [0 diry], pfin, curpos, finpos, wt, strategy, recur);
            if cost1<cost
                move = move1;
                cost = cost1;
                curposnew = curposnew1;
            end
        end;
    end;

end;


if isinf(cost) && recur > 1 && abs(strategy) < 5


    mv  = [dirx diry];
    mv1 = abs(mv)-1;
    if any(mv1)
        mv2 = abs(mv1);
        mv3 = -mv;
        ok3 = 1;
    else
        mv1 = [-mv(1) 0];
        mv2 = [0 -mv(2)];
        ok3 = 0;
    end;
    sz  = size(curpos);
    strategy = strategy + sign(strategy)*1;

    if all(pinit+mv1 > 0) && all(pinit+mv1 <= sz) && ~curpos(pinit(1)+mv1(1),pinit(2)+mv1(2))
        [move1 cost1 curposnew1] = onemove(item, pinit, mv1, pfin, curpos, finpos, wt, strategy, recur-1);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
    if all(pinit+mv2 > 0) && all(pinit+mv2 <= sz) && ~curpos(pinit(1)+mv2(1),pinit(2)+mv2(2))
        [move1 cost1 curposnew1] = onemove(item, pinit, mv2, pfin, curpos, finpos, wt, strategy, recur-1);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
    if ok3 && all(pinit+mv3 > 0) && all(pinit+mv3 <= sz) && ~curpos(pinit(1)+mv3(1),pinit(2)+mv3(2))
        [move1 cost1 curposnew1] = onemove(item, pinit, mv3, pfin, curpos, finpos, wt, strategy, recur-1);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
    if all(pinit+mv1 > 0) && all(pinit+mv1 <= sz) && curpos(pinit(1)+mv1(1),pinit(2)+mv1(2)) > 0
        [move1 cost1 curposnew1] = onemove(item, pinit, mv1, pfin, curpos, finpos, wt, strategy, recur-1);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
    if all(pinit+mv2 > 0) && all(pinit+mv2 <= sz) && curpos(pinit(1)+mv2(1),pinit(2)+mv2(2)) > 0
        [move1 cost1 curposnew1] = onemove(item, pinit, mv2, pfin, curpos, finpos, wt, strategy, recur-1);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;
    if ok3 && all(pinit+mv3 > 0) && all(pinit+mv3 <= sz) && curpos(pinit(1)+mv3(1),pinit(2)+mv3(2)) > 0
        [move1 cost1 curposnew1] = onemove(item, pinit, mv3, pfin, curpos, finpos, wt, strategy, recur-1);
        if cost1<cost
            move = move1;
            cost = cost1;
            curposnew = curposnew1;
        end
    end;

end;

if isinf(cost)
    move=[];
    curposnew = curpos;
end;


function [move, cost, curpos] = onemove(item, pinit, movedir, pfin, curpos, finpos, wt, strategy, recur)



curcost  = mod(wt(curpos(pinit(1), pinit(2))), 10000);
move3    = [];
costnew3 = 0;
newpos   = pinit+movedir;
if curpos(newpos(1), newpos(2)) > 0

    if wt(curpos(newpos(1), newpos(2))) > 8000
        cost = Inf;
        move = [];
        return;
    end;

    if wt(curpos(newpos(1), newpos(2))) > 2*curcost && recur > 0 && strategy > 0
        tmppos1 =  abs(movedir)-1   + pinit;
        tmppos2 =  abs(movedir)-1   + pinit;
        tmppos3 = movedir + tmppos1;
        tmppos4 = movedir + tmppos2;
        sz = size(curpos);
        wt2       = [0; wt];
        curpos2 = curpos;
        curpos2(curpos2 < 0) = 0;
        curpos2 = curpos2+1;
        tmpcost = [Inf Inf];
        if all(tmppos1 > 0) && all(tmppos1 <= sz)
            if all(tmppos3 > 0) && all(tmppos3 <= sz)
                tmpcost(1) = wt2(curpos2(tmppos1(1), tmppos1(2))) + wt2(curpos2(tmppos1(1), tmppos1(2))) + 2*curcost;
            end;
        end;
        tmpcost(2)=tmpcost(1);

        if tmpcost(2) < wt(curpos(newpos(1), newpos(2)))
            [move3 costnew3 curpostmp] = onemove(item, pinit, tmppos2-pinit, pfin, curpos, finpos, wt, strategy, -1);
            if ~isinf(costnew3)
                [move4 costnew4 curpostmp] = onemove(item, tmppos2, tmppos4-tmppos2, pfin, curpostmp, finpos, wt, strategy, -1);
                if ~isinf(costnew4)
                    [move5 costnew5 curpostmp] = movefrompos(item, tmppos4, pfin, curpostmp, finpos, wt, -strategy, recur-1);
                    cost = costnew3 + costnew4 + costnew5;
                    move = [move3; move4; move5];
                    if ~isinf(cost), curpos = curpostmp; end;
                    return;
                end;
            end;
        end;
    end;

    newitem = curpos(newpos(1), newpos(2));
    [pfin2(1) pfin2(2)] = find(finpos == newitem);


    dirx2 = sign(pfin2(1)-newpos(1));

    diry2 = sign(pfin2(2)-newpos(2));
    wt(newitem) = wt(newitem) + 20000;
    if all( [dirx2 diry2] == -movedir ) || all(newpos == pfin2)
        costnew3 = Inf;
    else

        [move3 costnew3 curposnew] = movefrompos(newitem, newpos, pfin2, curpos, finpos, wt, strategy, -1);
        if costnew3 == 0, costnew3 = Inf; end;
        if ~isinf(costnew3), curpos = curposnew; end;
    end;

    if isinf(costnew3)
        wtdir = [ Inf Inf Inf Inf ];
        diffpos = pfin - newpos;
        if (diffpos(1) && movedir(1)) || (diffpos(2) && movedir(2))
            mv1 = abs(movedir)-1;
        else
            if any(sign(diffpos) > 0)
                mv1 = abs(movedir)-1;
            else
                mv1 = abs(abs(movedir)-1);
            end;
        end;
        mv2 = -mv1;
        mv3 = movedir;
        sz  = size(curpos);
        if all(newpos+mv1 > 0) && all(newpos+mv1 <= sz)

            if curpos(newpos(1)+mv1(1),newpos(2)+mv1(2)) <= 0
                [move3 costnew3 curpos] = movefrompos(newitem, newpos, newpos+mv1, curpos, finpos, wt, strategy, -1);
            else wtdir(1) = wt(curpos(newpos(1)+mv1(1),newpos(2)+mv1(2)));
            end;
        end;
        if isinf(costnew3) && all(newpos+mv2 > 0) && all(newpos+mv2 <= sz)
            if curpos(newpos(1)+mv2(1),newpos(2)+mv2(2)) <= 0
                [move3 costnew3 curpos] = movefrompos(newitem, newpos, newpos+mv2, curpos, finpos, wt, strategy, -1);
            else wtdir(2) = wt(curpos(newpos(1)+mv2(1),newpos(2)+mv2(2)));
            end;
        end;
        if isinf(costnew3) && all(newpos+mv3 > 0) && all(newpos+mv3 <= sz)
            if curpos(newpos(1)+mv3(1),newpos(2)+mv3(2)) <= 0
                [move3 costnew3 curpos] = movefrompos(newitem, newpos, newpos+mv3, curpos, finpos, wt, strategy, -1);
            else wtdir(3) = wt(curpos(newpos(1)+mv3(1),newpos(2)+mv3(2)));
            end;
        end;


        if isinf(costnew3)
            [wtdir si] = sort(wtdir);
            for index = 1:length(wtdir)
                if ~isinf(wtdir(index)) && isinf(costnew3)
                    switch si(index),
                        case 1, [move3 costnew3 curpos] = movefrompos(newitem, newpos, newpos+mv1, curpos, finpos, wt, strategy, -1);
                        case 2, [move3 costnew3 curpos] = movefrompos(newitem, newpos, newpos+mv2, curpos, finpos, wt, strategy, -1);
                        case 3, [move3 costnew3 curpos] = movefrompos(newitem, newpos, newpos+mv3, curpos, finpos, wt, strategy, -1);
                    end;
                end;
            end;
        end;
    end;
    wt(newitem) = wt(newitem) - 20000;


    if isinf(costnew3) || costnew3 == 0
        cost = Inf;
        move = [];
        return;
    end;
end;


if curpos(newpos(1), newpos(2)) == -item
    cost = Inf;
    move = [];
    return;
end;

curpos(newpos(1), newpos(2)) =  curpos(pinit(1), pinit(2));
curpos(pinit(1) , pinit(2) ) = -curpos(pinit(1), pinit(2));

if recur > 0
    [move2 costnew curposnew] = movefrompos(item, pinit+movedir, pfin, curpos, finpos, wt, strategy, recur-1);
else
    [move2 costnew curposnew] = movefrompos(item, pinit+movedir, pfin, curpos, finpos, wt, strategy, 0);
end;
curpos = curposnew;

if curpos(pinit(1), pinit(2))<0
    curpos(pinit(1), pinit(2)) = 0;
end;


move = [ move3; item pinit movedir; move2];
cost = costnew3+costnew+curcost;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function ndx = sub2ind(siz,varargin)
siz = [siz ones(1,nargin-length(siz)-1)];
mt = cellfun('isempty',varargin);
if any(mt)
    ndx = zeros(~mt);
    return;
end
k = [1 cumprod(siz(1:end-1))];
ndx = 1;
for i = 1:length(siz),
    v = varargin{i};
    ndx = ndx + (v-1)*k(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function P = perms(V)
V = V(:).';
n = length(V);
if n <= 1, P = V;
    return;
end
q = perms(1:n-1);
m = size(q,1);
P = zeros(n*m,n);
P(1:m,:) = [n * ones(m,1) q];
for i = n-1:-1:1
    t = q;
    t(t == i) = n;
    P((n-i)*m+1:(n-i+1)*m,:) = [i*ones(m,1) t];
end
P = V(P);