Winner Ave (More tricks)

Finish 2003-04-08 00:00:00 UTC

RouteFinder

by Yi Cao

Status: Passed
Results: 1794.707K raw score
CPU Time: 114.034
Score: 1837.47
Submitted at: 2003-04-08 04:13:18 UTC
Scored at: 2003-04-08 04:27:42 UTC

Current Rank: 410th
Basis for: RouteFinder_01 (diff)
Basis for: Getting Close (diff)

Comments
Yi Cao
08 Apr 2003
A new route finder with grouped local finer
Please login or create a profile.
Code
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<.9 & 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%