Winner GUO (asdf1)

2004-11-10 09:00:00 UTC

# lean94g

Status: Failed
Results:

Based on: lean94f (diff)
Basis for: lean94h (diff)

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_1
moves2 = solver_1(Af,Ai,w);

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

if res-optcost < 30
return
end

moves2 = reversemov(solver_1(Af,Ai,w));

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

if res-optcost < 30
return
end

% try solver_3
moves2 = solver_3(Ai,Af,w));

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

if res-optcost < 30
return
end

moves2 = solver_SFAL(Af,Ai,w);

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

if res-optcost < 100
return
end

% not solved yet? then get desperate
while ~isfinite(res)
moves = solver_1(Ai,Af,w);
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));

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 = [];

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 = [ceil(2*log(n)):(n-1)]

choosew = sum(blockers).*invw;
for i = 1:leaveout
end

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

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

% solve the subproblem
[mov2,ok] = solver_linspac(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

Ap(ipos(i)) = i;
end

if ok
partialmov = trymov;
ok = 1;
mov = partialmov;
return
end
else
dropw = blockers.*invw;
partialmov(find(partialmov(:,choose)),:) = [];
end
end
ok = 0;
end

mov(find(mov(:,c)),:) = [];
nmov = size(mov,1);
[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);

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

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

for k=1:(nmov+1)
R = zeros(size(A));
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);

```