function c = solver(a)
global randfact
randfact=0;
c = subsolver(a);
fp = sum(a(c,3));
ft = sum(a(:,3));
da=a(:,4)-a(:,3);
if fp/ft<.85 & length(c)>10,
atmp = a;
atmp(:,3)=atmp(:,3).^3.0;
freightMean = median(a(find(a(:,3)>0),3));
highFreight = find(a(:,3)>freightMean);
atmp(highFreight,3) = 1.6*atmp(highFreight,3);
randfact=7;
c2 = subsolver(atmp);
score2 = sum(da(c2(2:end-1)));
score1 = sum(da(c(2:end-1)));
if score2<score1,
c = c2;
end
end
function c = subsolver(a)
iworth=find(a(:,3));
mfreight=sum(a(iworth,3))/length(iworth);
ia=find( a(:,4)>a(1,4)*.0714 | a(:,3)>=(mfreight-1.9)*.074);
a=a(ia,:);
xy = (a(:,1)-a(1)) + 1i*(a(:,2)-a(1,2));
if ~(abs(xy(2:end)) <= a(1,4))
% Not enough fuel to go anywhere, stay home (Ave)
c=[1 1];
return
end
n=size(a,1);
da=a(:,4)-a(:,3);
X = xy(:, ones(n,1));
dist = abs(X-X.');
dist0=dist;
sfr=sum(a(:,3))*0.9;
[c,score1]=routefinder(a,dist,xy);
[c,score1]=addgroup(a,c,dist,da,xy);
sfrc=sum(a(c,3));
if sfrc<sfr
[c,score1]=findenum(a,dist,da,sfr);
[c,score1]=addgroup(a,c,dist,da,xy);
end
sfrc=sum(a(c,3));
if sfrc<sfr
if n<19
c2=turbodiesel(a,dist0);
else
c2=diesel2(a,dist0);
end
[c2,score2]=addgroup(a,c2,dist,da,xy);
if(score2<score1)
c=c2;
score1=score2;
end
end
sfrc=sum(a(c,3));
if sfrc<sfr
c3 = petrol3(a,dist0);
[c2,score2]=addgroup(a,c2,dist,da,xy);
if(score2<score1)
c=c2;
score1=score2;
end
end
if da(1)<score1
% dummy solution is the best(!!)
c=[1 1];
end
c=ia(c);
function [c,score]=addgroup(a,c,dist,da,xy)
score=sum(da(c(2:end-1)));
[c,score]=optim(a,c,score,dist,da);
da1=da;
da1(c)=0;
if any(da1<0)
if numel(c)>4
c=savegas(a,c,score,xy,da);
end
score=sum(da(c(2:end-1)));
[c,score]=optim(a,c,score,dist,da);
c = addGas(a(:,4),c,dist);
c = addFreight(a,c,dist);
score=sum(da(c(2:end-1)));
[c,score]=optim(a,c,score,dist,da);
end
function c = addGas(a,c,dist)
cont=1;
c=c(:);
while cont
gasv=a;
gasv(c)=0;
fIdx = find(gasv);
if isempty(fIdx), return; end % NO GAS LEFT TO ADD
Nc = length(c);
N = length(fIdx);
% CALCULATE EXCESS GAS AT EACH POINT IN PATH
d=ones(Nc-1,1);
for i=1:Nc-1 % faster in Matlab6.5 than sub2ind.. ?
d(i)=dist(c(i),c(i+1));
end
eg = cumsum(a(c(1:Nc-1))-d);
dispgas=zeros(Nc-1,1);
dispgas(Nc-1)=eg(Nc-1);
for counter=Nc-2:-1:1
dispgas(counter)=min(dispgas(counter+1),eg(counter));
end
% SORT FREIGHT BY VALUE
% TRY TO FIT BIGGEST FIRST
ad = dist(c,fIdx);
h = dispgas(:,ones(1,N)) - (ad(1:Nc-1,:) + ad(2:Nc,:) - d(:,ones(1,N)));
proximity = max(h);
[ignore,idx] = sort(-a(fIdx).*proximity(:));
cont=0;
fIdx = fIdx(idx);
for n=1:N,
i=fIdx(n);
% SET OF EXCESS DISTANCES ADDED BY VISITING THIS POINT
ad = dist(c,i);
ed = ad(1:Nc-1) + ad(2:Nc) - d;
if(any(dispgas>=ed))
h=dispgas-ed;
[ignore,closest] = max(h);
if ed(closest)<a(fIdx(n))
if h(closest)>=0
cont=1;
d1=dist(c(closest),i);
d2=dist(i,c(closest+1));
d0=dist(c(closest),c(closest+1));
c = [c(1:closest); i; c(closest+1:Nc)];
d=[d(1:closest-1);d1;d2;d(closest+1:Nc-1)];
eg = cumsum(a(c(1:Nc))-d);
dispgas=zeros(Nc,1);
dispgas(Nc)=eg(Nc);
for counter=Nc-1:-1:1
dispgas(counter)=min(dispgas(counter+1),eg(counter));
end
Nc=Nc+1;
end
end
end % if any dispgas
end % for n
end % while cont
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% addFreight
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = addFreight(a,c,dist)
cont=1;
c=c(:);
a(c,3)=0;
while cont
Nc = length(c);
fIdx = find(a(:,3));
if isempty(fIdx), return; end % NO FREIGHT LEFT TO ADD
N = length(fIdx);
% CALCULATE EXCESS GAS AT EACH POINT IN PATH
d=ones(Nc-1,1);
for i=1:Nc-1 % faster in Matlab6.5 than sub2ind.. ?
d(i)=dist(c(i),c(i+1));
end
eg = cumsum(a(c(1:Nc-1),4)-d);
dispgas=zeros(Nc-1,1);
dispgas(Nc-1)=eg(Nc-1);
for counter=Nc-2:-1:1
dispgas(counter)=min(dispgas(counter+1),eg(counter));
end
% SORT FREIGHT BY VALUE
% TRY TO FIT BIGGEST FIRST
ad = dist(c,fIdx);
h = dispgas(:,ones(1,N)) - (ad(1:Nc-1,:) + ad(2:Nc,:) - d(:,ones(1,N)));
proximity = max(h);
[ignore,idx] = sort(-a(fIdx,3).*proximity(:));
%[ignore,idx] = sort(-a(fIdx,3));
fIdx = fIdx(idx);
cont=0;
for n=1:N
i=fIdx(n);
% SET OF EXCESS DISTANCES ADDED BY VISITING THIS POINT
ad = dist(c,i);
ed = ad(1:Nc-1) + ad(2:Nc) - d;
%%%%%%%%%%%%%%%%%%% This seems to be quicker than the max...
if any(dispgas>=ed)
%%%%%%%%%%%%%%%%%%%
h=dispgas-ed;
[ignore,closest] = max(h);
if h(closest)>=0
cont=1;
d1=dist(c(closest),i);
d2=dist(i,c(closest+1));
d0=dist(c(closest),c(closest+1));
a(i,3)=0;
c = [c(1:closest); i; c(closest+1:Nc)];
d=[d(1:closest-1);d1;d2;d(closest+1:Nc-1)];
eg=[eg(1:closest-1);eg(closest)+d0-d1;eg(closest:Nc-1)-(d1+d2-d0)];
dispgas=zeros(Nc,1);
dispgas(Nc)=eg(Nc);
for counter=Nc-1:-1:1
dispgas(counter)=min(dispgas(counter+1),eg(counter));
end;
Nc=Nc+1;
end
end
end
end
function c = diesel2(a,Dmat)
global randfact;
n=size(a,1);
rand('state',1+randfact);
freightv=a(:,3);
petrolv=a(:,4);
Dmat1=Dmat(:,1);
pos=1;
jdest=1;
totalpetrol=sum(petrolv);
unitpetrol = totalpetrol/(n-1);
petroln=petrolv/totalpetrol;
totalfreight=sum(freightv);
freightn=freightv/totalfreight;
petrol = petrolv(1);
petrolv(1)=0; freightv(1)=-1;
imax=[0;0];
notbeenthere=logical(ones(n,1));
notbeenthere(1)=0;
c=ones(n+1,1);
notnearlythere=notbeenthere;
while 1
jdest=jdest+1;
% is it possible to go anywhere?
dpos=Dmat(:,pos);
bpossible=(petrol >= dpos & notbeenthere);
if ~(bpossible)
break; % fuel low, go home if possible
end
% is it safe to go there?
safe = find(bpossible & (petrol + petrolv >= dpos + Dmat1));
if (isempty(safe))
possible=find(bpossible);
% is the next station after this one safe? (Ave)
safe=zeros(n,1);
for li=1:length(possible)
la=possible(li);
lanotbeenthere=notbeenthere;
lanotbeenthere(la)=0;
ladpos=Dmat(la,pos)+Dmat(:,la);
lapetrol=petrol+petrolv(la);
lapossible=(lapetrol>=ladpos & lanotbeenthere);
if any(lapossible & lapetrol+petrolv>=ladpos+Dmat1)
safe(la)=1;
break;
end
end
safe=find(safe);
end
if isempty(safe)
break;
end
possible=safe;
lm=length(possible);
if (lm>1)
Dpossible=Dmat(possible,pos);
metric = ( (petroln(possible) - Dpossible/totalpetrol) * 2*(n-jdest+1)/n*(1+unitpetrol/petrol)^5 ...
+ freightn(possible)) ...
./ Dpossible;
sm = [-realmax; -realmax];
for imetric=1:lm
mi=metric(imetric);
if(mi>sm(2))
sm(1)=sm(2); imax(1)=imax(2);
sm(2)=mi; imax(2)=imetric;
elseif mi>sm(1)
sm(1)=mi; imax(1)=imetric;
end
end
rn=0;
if sm(1)/sm(2)>.995
rn=fix(rand+.5);
end
next=possible(imax(end-rn));
else
next=possible(1);
end
notnearlythere(next)=0;
r=Dmat(pos,next);
% find places that are on the way
ontheway = find( (Dmat(:,next) + Dmat(:, pos) < 1.088*r & notnearlythere & bpossible) );
if(~isempty(ontheway))
dp = Dmat(ontheway,pos); dn=Dmat(ontheway,next);
ontheway = ontheway( find ( ...
dp+dn+Dmat(1,next) <= petrol+petrolv(ontheway)+petrolv(next) & ... % can get home from
((dp+dn-r)<petrolv(ontheway) | freightv(ontheway)>0 ) ) ... % net freight or petrol gain
);
[maxf,im]=max(freightv(ontheway));
if (~isempty(ontheway))
next=ontheway(im(1));
end
end
c(jdest)=next;
petrol=petrol + petrolv(next) - Dmat(pos,next);
pos=next;
notbeenthere(pos)=0;
notnearlythere=notbeenthere;
end
c(jdest)=1;
c=c(1:jdest);
function cb=turbodiesel(a,dist)
freightv=a(:,3);
petrolv=a(:,4);
da=freightv-petrolv;
n=size(a,1);
Nstck=1000;
stck=zeros(Nstck,1);
nstck=1;
c=ones(n+1,1);
cb=[1 1];
avail=logical(ones(n,1));
GAS=ones(n,1);
stck(2)=1;
nstck=2;
totalpetrol=sum(petrolv);
unitpetrol = totalpetrol/(n-1);
petroln=petrolv/totalpetrol;
totalfreight=sum(freightv);
freightn=freightv/totalfreight;
gas=petrolv(1);
petrolv(1)=0; freightv(1)=-1;
dist1=dist(:,1);
if n<14
IMAX=3;
else
IMAX=2;
end
fgbest=-inf;
loop=1;
while 1
current=stck(nstck);
if current
c(loop)=current;
avail(current)=0;
if loop>1
gas=gas-dist(c(loop-1),current)+petrolv(current);
%test
if gas>=dist(current)
fg=sum(da(c(2:loop)));
if fgbest<fg % best
cb=c(1:loop+1);
cb(loop+1)=1;
fgbest=fg;
end
end
if (da(avail)<=0)
break;
end
end
dpos=dist(:,current);
bpossible=(gas >= dpos & avail);
nok=1;
if any(bpossible)
% is it safe to go there?
safe = find(bpossible & (gas + petrolv >= dpos + dist1));
if isempty(safe)
reachable=find(bpossible);
% is the next station after this one safe? (Ave)
bsafe=logical(zeros(n,1));
for li=1:length(reachable)
la=reachable(li);
lanotbeenthere=avail;
lanotbeenthere(la)=0;
ladpos=dist(la,current)+dist(:,la);
lapetrol=gas+petrolv(la);
lapossible=(lapetrol>=ladpos & lanotbeenthere);
if any(lapossible & lapetrol+petrolv>=ladpos+dist1)
bsafe(la)=1;
break;
end
end
safe=find(bsafe);
if isempty(safe)
safe=find(bpossible);
Dpossible=dist(safe,current);
metric = ( (petroln(safe) - Dpossible/totalpetrol) * 2*(n-loop+1)/n*(1+unitpetrol/gas)^5 ...
+ freightn(safe)) ...
./ Dpossible;
[sm,imax]=sort(-metric);
if nstck>16
IMAX2=1;
else
IMAX2=IMAX;
end
safe=safe(imax(1:min(end,IMAX2)));
end
end
if ~isempty(safe)
reachable=safe;
Dpossible=dist(reachable,current);
metric = ( (petroln(reachable) - Dpossible/totalpetrol) * 2*(n-loop+1)/n*(1+unitpetrol/gas)^5 ...
+ freightn(reachable)) ...
./ Dpossible;
[sm,imax]=sort(-metric);
%imax=imax(sm<sm(1)*0.9); this doesn't work because of possible negative and positive values
if nstck>30
IMAX2=1;
else
IMAX2=IMAX;
end
reachable=reachable(imax(1:min(end,IMAX2)));
notnearlythere=avail;
for i=1:length(reachable)
next=reachable(i);
r=dist(current,next);
% find places that are on the way
ontheway = find( (dist(:,next) + dist(:, current) < 1.088*r & notnearlythere & bpossible) );
if(~isempty(ontheway))
dp = dist(ontheway,current); dn=dist(ontheway,next);
ontheway = ontheway( find ( ...
dp+dn+dist(1,next) <= gas+petrolv(ontheway)+petrolv(next) & ... % can get home from
((dp+dn-r)<petrolv(ontheway) | freightv(ontheway)>0 ) ) ... % net freight or petrol gain
);
[maxf,im]=max(freightv(ontheway));
if (~isempty(ontheway))
reachable(i)=ontheway(im(1));
end
end
notnearlythere(next)=0;
end % for all reachables
if ~isempty(reachable)
stck(nstck)=0;
stck(nstck+1:nstck+length(reachable))=reachable;
nstck=nstck+length(reachable);
GAS(loop)=gas;
loop=loop+1;
nok=0;
end % not empty reachable
end % no safe (or almost safe) found
end % not totally impossible
if nok
if dist(current)<=gas % normally always?
fg=sum(da(c(2:loop)));
if fgbest<fg % best
cb=c(1:loop+1);
cb(loop+1)=1;
fgbest=fg;
end
end
nstck=nstck-1;
avail(current)=1;
end
else
loop=loop-1;
if loop<=1
break;
end
avail(c(loop))=1;
nstck=nstck-1;
gas=GAS(loop-1);
end % current
end % while
function c = petrol3(a,Dmat)
n=size(a,1);
freightv=a(:,3);
petrolv=a(:,4);
petrol=0;
freight=0;
pos=1;
jdest=1;
totalpetrol=sum(petrolv);
unitpetrol = totalpetrol/(n-1);
petroln=petrolv/totalpetrol;
totalfreight=sum(freightv);
freightn=freightv/totalfreight;
petrol = petrol + petrolv(1);
freight = freight + freightv(1);
petrolv(1)=0; freightv(1)=-1;
notbeenthere=logical(ones(n,1));
notbeenthere(1)=0;
c=ones(1,n+1);
Dp = max(Dmat(:,1) - petrolv,0);
while 1
jdest=jdest+1;
pmD = petrol - Dmat(:,pos);
possible = find( ...
(pmD >= Dp) & ...
notbeenthere ...
);
if isempty(possible)
break
end
Dpossible=Dmat(possible,pos);
metric = ( (petroln(possible) - Dpossible/totalpetrol) * 2*(n-jdest+1)/n*(1+unitpetrol/petrol)^5 ...
+ freightn(possible)) ...
./ Dpossible;
imax= 1;maxm = -inf;
for count = 1:numel(metric)
if metric(count)>maxm
imax = count;
maxm = metric(count);
end
end
% [maxm,imax]=max(metric);
next=possible(imax);
c(jdest)=next;
petrol=petrol + petrolv(next) - Dmat(pos,next);
freight=freight+freightv(next);
pos=next;
notbeenthere(pos)=0;
end
c=c(1:jdest);
function [cbs,fgsbest]=findenum(a,dist,da,sfr)
n=size(a,1);
weightfreight=[ones(1,7)*50 zeros(1,n)];
Nstck=n*9;
stck=zeros(Nstck,1);
nstck=1;
c=ones(n+1,1);
cbs=[1 1]; % short-lines
fgsbest=da(1);
avail=logical(ones(n,1));
GAS=c;
gas=a(1,4);
GAS(1)=gas;
stck(2)=1;
nstck=2;
loop=1;
cnt=0;
opt=1; % not yet used
while 1
current=stck(nstck);
if current
c(loop)=current;
avail(current)=0;
if loop>1
gas=gas-dist(c(loop-1),current)+a(current,4);
end
if dist(current)<=gas % can return home
fg=sum(da(c(2:loop)));
if fgsbest>fg
cbs=c(1:loop+1);
cbs(loop+1)=1;
fgsbest=fg;
end
end
i=find(avail);
reachable=i(dist(i,current)<=gas);
if isempty(reachable)
nstck=nstck-1;
avail(current)=1;
cnt=cnt+1;
if cnt>38
opt=0;
break;
end
else
metric=dist(reachable,current)+a(reachable,3)*weightfreight(loop);
[dsts,i]=sort(metric);
stck(nstck)=0;
nn=min(length(reachable),fix(rand+4.5));
stck(nstck+1:nstck+nn)=reachable(i(nn:-1:1));
nstck=nstck+nn;
GAS(loop)=gas;
loop=loop+1;
end
else
loop=loop-1;
avail(c(loop))=1;
if loop==1
break;
end
nstck=nstck-1;
gas=GAS(loop-1);
end
if loop==n
if abs(sfr-sum(a(c,3)))<1e-10
break
end
end
end
function [c,s]=optim(a,c,s,dist,da)
if length(c)<4
return
end
opt=1;
g0=ones(length(c)-1,1);
g0(1)=a(1,4)-dist(c(2));
f=a(c,4);
n=length(c);
for i=2:n-1
g0(i)=g0(i-1)+f(i)-dist(c(i),c(i+1));
end
dg=g0;
while opt
opt=0;
da1=da(c);
i=find(da1(2:end-1)>0)+1;
if isempty(i)
return
end
for k=1:length(i)
m=i(k);
dg(k)=f(m)-dist(c(m-1),c(m))-dist(c(m),c(m+1))+dist(c(m-1),c(m+1));
end
j=find(dg(1:length(i))<=g0(end));
if isempty(j)
return
end
i=i(j);
dg=dg(j);
[da11,j]=sort(-da1(i));
i=i(j);
dg=dg(j);
for j=1:length(i)
if (g0(i(j):end)>=dg(j))
s=s+da11(j);
g0=[g0(1:i(j)-2);g0(i(j):end)-dg(j)];
c(i(j))=[];
f(i(j))=[];
opt=1;
break;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% RouteFinder Version 1.0
% by Yi Cao, 08/04/2003
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c,s]=routefinder(a,mD,xY)
aN = size(a,1); % Constant section
onesN = ones(aN,1);
aF = a(:,3);
aG = a(:,4);
tF = sum(aF);
da = aG - aF;
vD = mD(:,1);
tL = 5;
C=ones(tL,aN+1);
K=ones(tL,1);
S=zeros(tL,1);
for L=1:tL,
I=1;
c=C(L,:);
g=0;
f=aF;
rf=tF;
for k=2:aN
g=g+aG(I);
d=g-mD(:,I);
d(c(1:k-1))=-inf;
rf=rf-f(I);
f(I)=0;
if rf<10*eps & g >= mD(I,1),
k=k-1;
break;
end
if L==2,
[Y,I]=max(d);
elseif L==5
[Y,I]=max(d+sign(d).*aG);
elseif L==4
if k > max(0.04*aN,5),
[Y,I]=max(d);
else
[Y,I]=max(0.4*d+sign(d).*aG);
end
elseif L==3
xy=abs(xY-sum(xY.*f)/(aN-k+1)/(rf+eps))+eps;
[Y,I]=max(d+0.15*sign(d).*aG./xy);
else
p=d>0;
af=f.*p;
cf=sum(af);
if cf > 10*eps,
xy=abs(xY-sum(af.*xY)/sum(p)/(cf+eps))+eps;
[Y,I]=max(d+sign(d).*af./xy);
else
xy=abs(xY-sum(xY.*f)/(aN-k+1)/(rf+eps))+eps;
[Y,I]=max(d);
end
end
if Y<0,
k=k-1;
break;
end
g=d(I);
c(k) = I;
end
C(L,:)=c;
K(L)=k;
if k/aN > 0.9,
break;
end
end
for L=1:tL,
k=K(L);
c=C(L,1:k+1);
if k<2,
c=[1 1];
i2=2;
g1=aG(1);
else,
vG=aG(c);
d=diag(mD(c,c),1);
g1 = [aG(1);cumsum(vG(1:k-1))-cumsum(d(1:k-1))+vG(2:k)-mD(c(2:k),1)];
i1 = find(g1>=0);
i2 = i1(end) + 1;
c = c(1:i2);
c(i2) = 1;
end
K(L)=i2;
C(L,1:i2)=c;
S(L)=sum(da(c(2:i2)));
end
[s,L]=min(S);
c=C(L,1:K(L));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by Yi Cao
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% by Mr.Bond
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c0=savegas(a,c,score,xy,da)
% see if we can save gas by going A->C->B instead of A->B->C
c0 =c;
n = numel(c);
m=n-1;
g0 = cumsum(a(c(1:m),4)-abs(diff(xy(c))));
% more iterations takes too much time but there may be a smarter way
if n<150, OUTER_MAX=4; else OUTER_MAX=2; end
for outer = 2:OUTER_MAX
for point = 2:2:n-outer
c = c0;
c([point:point+outer-1]) =c([point+outer-1:-1:point]);
g = cumsum(a(c(1:m),4)-abs(diff(xy(c))));
if g(m)>g0(m)
if any(g<0)
else
c0=c;
g0=g;
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
