function c = solver(a)
iworth=find(a(:,3));
if isempty(iworth), c=[1,1]; return; end
mfreight=mean(a(iworth,3))/14.6;
ia=find( a(:,4)>a(1,4)/15 | a(:,3)>=mfreight);
a=a(ia,:);
xy = (a(:,1)-a(1)) + 1i*(a(:,2)-a(1,2));
n=size(a,1);
if ~any(a(1,4) >= abs(xy(2:n)))
% Nofreight or not enough fuel to go anywhere, stay home (Ave)
c=[1 1];
return
end
da=a(:,4)-a(:,3);
[X,Y]=meshgrid(xy);
dist = abs(X-Y);
dist0=dist;
I=1;
c=ones(1,n+1);
g=0;
for k=2:n
g=g+a(I,4);
dist(:,I)=inf;
[Y,I]=min(dist(I(1),:));
g=g-Y;
if g<0
k=k-1;
break;
end
c(k) = I;
end
c=c(1:k+1);
xy1 = xy(c);
i=sub2ind([n n],c(1:k),c(2:k+1));
g = cumsum(a(c(1:k),4)-dist0(i)');
if g(end)<0
if k<3
c=[1 1];
elseif g(end)<dist0(c(k))
% look for last point from which starting point can be reached
g=g(1:k-2)+a(c(2:k-1),4)-dist0(c(2:k-1),1);
i=find(g>=0);
if isempty(i)
c=[1 1];
else
i=i(end)+2;
c=c(1:i);
c(i)=1;
end
else
c = c(1:k+1);
c(end)=1;
end
end
score1 = sum(da(c(2:end-1)));
if sum(a(c,3))/sum(a(:,3))<0.95
c2 = diesel2(a,xy,dist0);
score2 = sum(da(c2(2:end-1)));
if(score2<score1)
c=c2;
score1=score2;
end
end
fail=0;
while ~fail %First add gas
[c,fail] = addGas(a,c,dist0);
if fail
break
end;
end
fail=0;
while ~fail %Then add freight
[c,fail] = addFreight(a,c,dist0);
if fail
break
end;
end
[c,score1]=optim(a,c,score1,xy,da);
if da(1)<score1
% dummy solution is the best(!!)
c=[1 1];
end
c=ia(c);
function [c,fail] = addGas(a,c,dist)
fail=1;
gasv=a(:,4);
gasv(c)=0;
c=c(:);
Nc = length(c);
fIdx = find(gasv>0);
if isempty(fIdx), return; end % NO GAS LEFT TO ADD
N = length(fIdx);
% CALCULATE EXCESS GAS AT EACH POINT IN PATH
i=sub2ind(size(dist),c(1:Nc-1),c(2:Nc));
d = dist(i);
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
[ignore,idx] = sort(-a(fIdx,4));
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;
h=dispgas-ed;
[ignore,closest] = max(h);
if ed(closest)<a(fIdx(n),4)
if h(closest)>=0
fail=0;
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),4)-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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% addFreight
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [c,fail] = addFreight(a,c,dist)
fail=1;
c=c(:);
Nc = length(c);
a(c,3)=0;
fIdx = find(a(:,3)>0);
if isempty(fIdx), return; end % NO FREIGHT LEFT TO ADD
N = length(fIdx);
% CALCULATE EXCESS GAS AT EACH POINT IN PATH
i=sub2ind(size(dist),c(1:Nc-1),c(2:Nc));
d = dist(i);
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
[ignore,idx] = sort(-a(fIdx,3));
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;
h=dispgas-ed;
[ignore,closest] = max(h);
if h(closest)>=0
fail=0;
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=[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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function c = diesel2(a,xy,Dmat)
n=size(a,1);
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;
notbeenthere=logical(ones(n,1));
notbeenthere(1)=0;
c=ones(n+1,1);
notnearlythere=notbeenthere;
z0=zeros(n,1);
while 1
jdest=jdest+1;
% is it possible to go anywhere?
bpossible=z0;
abpossible=0;
% is it safe to go there?
safe=z0;
safeind=1;
for i=1:n
if (petrol>=Dmat(i,pos))¬beenthere(i)
bpossible(i)=1;
abpossible=1;
if (petrol + petrolv(i)>= Dmat(i,pos) + Dmat1(i))
safe(safeind)=i;
safeind=safeind+1;
end
end
end
if abpossible==0
break; % fuel low, go home if possible
end
safe=safe(1:safeind-1);
if (isempty(safe))
possible=find(bpossible);
% is the next station after this one safe? (Ave)
safe=z0;
for li=1:length(possible)
la=possible(li);
lanotbeenthere=notbeenthere;
lanotbeenthere(la)=0;
ladpos=Dmat(la,pos)+Dmat(:,la);
lapetrol=petrol+petrolv(la);
safe(la)=0;
for i=1:n
if lanotbeenthere(i)&(lapetrol>=ladpos(i))&(lapetrol+petrolv(i)>=ladpos(i)+Dmat1(i))
safe(la)=1;
end
end
end
safe=find(safe);
end
if isempty(safe)
break;
end
maxm=-Inf;
imax=1;
for i=1:length(safe)
metric= ( (petroln(safe(i)) - Dmat(safe(i),pos)/totalpetrol) * 2*(n-jdest+1)/n*(1+unitpetrol/petrol)^5 ...
+ freightn(safe(i))) ...
./ Dmat(safe(i),pos);
if metric>maxm
maxm=metric;
imax=i;
end
end
next=safe(imax);
notnearlythere(next)=0;
r=Dmat(pos,next);
% find places that are on the way
maxf=-Inf;
im=1;
ok=0;
for i=1:n
if (((Dmat(i,next)+Dmat(i, pos)<1.088*r)¬nearlythere(i)&bpossible(i))&(Dmat(i,pos)+Dmat(i,next)+Dmat(1,next)<=petrol+petrolv(i)+petrolv(next))&(((Dmat(i,pos)+Dmat(i,next)-r)<petrolv(i))|(freightv(i)>0))&(freightv(i)>maxf))
maxf=freightv(i);
im=i;
ok=1;
end
end
if ok
next=im;
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 [c,s]=optim(a,c,s,xy,da)
opt=1;
g0 = cumsum(a(c(1:end-1),4)-abs(diff(xy(c))));
while opt
opt=0;
da1=da(c);
i=find(da1(2:end-1)>0)+1;
if isempty(i)
return
end
dg=a(c(i),4)-abs(xy(c(i))-xy(c(i-1)))-abs(xy(c(i+1))-xy(c(i)))+abs(xy(c(i+1))-xy(c(i-1)));
j=find(dg<=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 all(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))=[];
opt=1;
break;
end
end
end
|
