Winner Ave (More tricks)

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

another heuristic?

by Jack Snoeyink

Status: Passed
Results: 1687.271K raw score
CPU Time: 101.536
Score: 1710.17
Submitted at: 2003-04-09 16:46:46 UTC
Scored at: 2003-04-09 17:20:41 UTC

Current Rank: 66th
Based on: Add more gas and freight correct (diff)
Basis for: tuning again b (diff)

Comments
Please login or create a profile.
Code
function c = solverq13(a)


global randfact
randfact=0;
c = [1 1];
xy = (a(:,1)-a(1)) + 1i*(a(:,2)-a(1,2));
d = abs(xy);
gi = find((d <= a(1,4)) & (a(:,4)>0)); % reachable fuel
if length(gi) <= 2
    if length(gi) <= 1
        ia = find(d<=a(1,4)/2);
        a = a(ia,:);
        if (size(a,1)>1)
            X = xy(ia,ones(size(ia)));
            X = abs(X.'-X); %dist
            c = addFreight(a,c,X);
            c = ia(c);
        end
        return 
    else % only 2 fuel spots
        d2 = abs(xy-xy(gi(2)));
        g = a(1,4)-d2(1)+a(gi(2),4);
        ia = find((d <= a(1,4)) | (d2 <= g)); % reachable
        if sum(a(ia,4)>0) <=2 
            a = a(ia,:);
            if sum(a(2:end,3))>1e6
                c = subsolver(a);
                c = ia(c);
            end 
            return
        end
    end
else % some fuel, is it enough?
    g = sum(a(gi,4));
    g2 = sum(a((d <=g) & (a(:,4)>0),4));
    if g2-g < 0.1 
        ia = find(d<=g);
        if (.7* size(a,1) > length(ia))
            a = a(ia,:);
            if sum(a(2:end,3))>1e-2
                c = subsolver(a);
                c = ia(c);
            end
            return
        end
    end
end

c = subsolver(a);
ft = sum(a(:,3));    

   da=a(:,4)-a(:,3);
score1 = sum(da(c(1:end-1)));
randfact=9;
fp = ft - sum(a(c,3));
if (fp  > 1e3) && (fp < 8e3) && (length(c)>10),
    gi = a(:,4)>0;
        atmp = a;
        atmp(gi,3)=800*atmp(gi,4);
        c2 = subsolver(atmp);
        score2 = sum(da(c2(1:end-1)));
        if score2<score1,
            c = c2;
            score1 = score2;
            opt = 1;
        end
    end
fp = sum(a(c,3));
if fp/ft<.9275 && numel(c)>2
  X = xy(:, ones(size(a,1),1));
  dist = abs(X-X.');
  c = addGas(a(:,4),c,dist);
  c = addFreight(a,c,dist);
  fp = sum(a(c,3));

  if fp/ft<.9275 && numel(c)>10
 
    atmp = a;
    pp=find(a(:,3)>0);
    freightMean = sum(a(pp,3))/numel(d)*2.56;
    highFreight = find(a(:,3)>freightMean);
    randfact=7;
    for k=1:4
      atmp(pp,3)=a(pp,3).^(k+1);
      atmp(highFreight,3) = 1.57*atmp(highFreight,3);
      c2 = subsolver(atmp);
      c2 = addGas(a(:,4),c2,dist);
      c2 = addFreight(a,c2,dist);
      score2 = sum(da(c2(2:end-1)));
      score1 = sum(da(c(2:end-1)));
      if score2<score1,
        c = c2;
        if k<3
          fp = sum(a(c,3));
        end
      end
      if fp/ft>=.925-.1*k
        break
      end
    end
  end
end

function c = subsolver(a)

iworth=find(a(:,3));

mfreight=sum(a(iworth,3))/length(iworth);

ia=find( a(:,4)>a(1,4)*.0712 | a(:,3)>=(mfreight-1.9)*.075);

a=a(ia,:);
xy = (a(:,1)-a(1)) + 1i*(a(:,2)-a(1,2));

n=size(a,1);
da=a(:,4)-a(:,3);
X = xy(:, ones(n,1));
dist = abs(X-X.');
sfr=sum(a(:,3))*.9;
[c,score1]=routefinder(a,dist,xy);
[c,score1]=optim(a,c,score1,dist,da);
da1=da;
da1(c)=0;
if any(da1<0)
  if numel(c)>6
    c=movepoints(a,c,dist);
    c=savegas(a,c,score1,dist,da);
  end
  c = addGas(a(:,4),c,dist);
  c = addFreight(a,c,dist);
  [c,score1]=optim(a,c,score1,dist,da);
end
sfrc=sum(a(c,3));
if sfrc<sfr
  [c2,score2]=findenum(a,dist,da);
  %[c2,score2]=optim(a,c2,score2,dist,da);
  da1=da;
  da1(c2)=0;
  if score2<score1 && any(da1<0)
    if numel(c2)>6
      c2=movepoints(a,c2,dist);
      c2=savegas(a,c2,score2,dist,da);
    end
    c2 = addGas(a(:,4),c2,dist);
    c2 = addFreight(a,c2,dist);
    [c2,score2]=optim(a,c2,score2,dist,da);
  end
  if score2 < score1,
    c = c2;
    score1 = score2;
  end
sfrc=sum(a(c,3));
end
if sfrc<sfr
  c2=diesel2(a,dist);
  score2 = sum(da(c2(2:end-1)));
  c3 = petrol3(a,dist);
  score3 = sum(da(c3(2:end-1)));
  if(score2<score1)
    [c,score1]=optim(a,c2,score2,dist,da);
  end
  if(score3<score1)
    [c,score1]=optim(a,c3,score3,dist,da);
  end
  da1=da;
  da1(c)=0;
  if any(da1<0)
    if numel(c)>6
      c=movepoints(a,c,dist);
      c=savegas(a,c,score1,dist,da);
    end
    c = addGas(a(:,4),c,dist);
    c = addFreight(a,c,dist);
    [c,score1]=optim(a,c,score1,dist,da);
  end
end
da1=da;
da1(c)=0;
if any(da1<0)
  if numel(c)>6
    c=movepoints(a,c,dist);
    c=savegas(a,c,score1,dist,da);
  end
  c = addGas(a(:,4),c,dist);

  c = addFreight(a,c,dist);
  [c,score1]=optim(a,c,score1,dist,da);
end
if da(1)<score1
  % dummy solution is the best(!!)
  c=[1 1];
end
c=ia(c);

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  % must be true
        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  % must be true from any(dispgas>=ed)
      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 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)
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;
STOP=false;
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;
        STOP=true;
      end
    else
      metric=dist(reachable,current)+a(reachable,3)*weightfreight(loop);
      [dsts,i]=sort(metric);
      stck(nstck)=0;
      if STOP
        stck(nstck+1)=reachable(i(1));
        nstck=nstck+1;
      else
        nn=min(length(reachable),fix(rand+4.5));
        stck(nstck+1:nstck+nn)=reachable(i(nn:-1:1));
        nstck=nstck+nn;
      end
      GAS(loop)=gas;
      loop=loop+1;
    end
  else
    loop=loop-1;
    avail(c(loop))=1;
    if loop==1
      break;
    end
    nstck=nstck-1;
    if STOP
      while loop>5
        while stck(nstck)
          nstck=nstck-1;
        end
        loop=loop-1;
        nstck=nstck-1;
      end
    end
    gas=GAS(loop-1);
  end
end

function [c,s]=optim(a,c,s,dist,da)
n = length(c);
if n<4
  return
end
mi=1; 
d = zeros(n-1,1);
eg = d;
while mi~=0 % while removed a point 
  n = length(c);
  n1=n-1;
  d(1) = dist(c(2));
  eg(1) = a(1,4)-d(1);
  for i=2:n1 
    d(i)=dist(c(i),c(i+1));
    eg(i) = eg(i-1) + a(c(i),4)-d(i);
  end
  %   xy = a(:,1)+1i*a(:,2);
  %e = cumsum(a(c(1:end-1),4)-abs(diff(xy(c))));
  %all(abs(e-eg(1:n1))<1e-9)
  mg = eg(n1); % min spare gas
  mx = 0; % max improvement
  mi = 0; 
  for i = n1:-1:2
    if(mg>eg(i))
      mg = eg(i);
    end
    g = da(c(i));
    if g >=0
      delta = -d(i-1)-d(i)+dist(c(i-1),c(i+1))+g;
      if delta<=0
        c(i) = [];
        d(i-1) = dist(c(i-1),c(i));
        s = s - g;
        mi = -1;
        mg = 0; %mg-delta; % increases spare gas
      elseif (delta<=mg) && (g>mx) % possible max adjustment
        mx = g;
        mi = i;
      end
    end
  end
  if mi>0
    c(mi) = [];
    s = s - mx;
  end
end

function c=savegas(a,c,score,dist,da)

% more iterations takes too much time but there may be a smarter way
n = numel(c);
m=n-1;
d=zeros(m,1);
eg = d;
for outer = [14:-2:2 3:2:9 8:-2:2 2 2]
  d(1) = dist(c(2));
  eg(1) = a(1,4)-d(1);
  for i=2:m 
    d(i)=dist(c(i),c(i+1));
    eg(i) = eg(i-1) + a(c(i),4)-d(i);
  end
  
  for point = n-outer:-1:2
    ia = c(point-1);
    ib = c(point);
    ic = c(point+outer-1);
    id = c(point+outer);
    % The JIT is a strange beast
    d1 = dist(ia,ib)+dist(ic,id);
    d2 = dist(ia,ic)+dist(ib,id);
    if d1>d2
      
      g = eg(point-1)+dist(ia,ib)-dist(ia,ic);
      
      for i = point+outer-1:-1:point+1
        if g<0 
          break;
        end
        g = g + a(c(i),4)-dist(c(i),c(i-1));
      end
      if g<0 
        continue;
      end
      g = g + a(ib,4)-dist(ib,id);
      c(point:point+outer-1) = c(point+outer-1:-1:point);
      g = eg(point-1)+dist(ia,ib)-dist(ia,ic);
      for ii=point-1:point+outer-1
        eg(ii) = g; 
        g = g + a(c(ii),4) - dist(c(ii),c(ii+1));
      end
    end
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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);
eps10=10*eps;
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<eps10 && 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.17*sign(d).*aG./xy);
    else
      p=d>0;
      af=f.*p;
      cf=sum(af);
      if cf > eps10,
        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;
    S(L)=sum(da(c(2:i2-1)));
  end
  K(L)=i2;
  C(L,1:i2)=c;
end
[s,L]=min(S);
c=C(L,1:K(L));

function c=movepoints(a,c,dist)
if length(c)<4
  return
end
c=c(:);
n=length(c);
opt=1;
nopt=1;
while opt && nopt<4
  opt=0;
  g0=ones(n-1,1);
  g0(1)=a(1,4)-dist(c(2));
  f=a(c,4);
  dxy=zeros(n-1,1);
  d2xy=zeros(n-2,1);
  dxy(1)=dist(c(2));
  for i=2:n-1
    dxy(i)=dist(c(i),c(i+1));
    d2xy(i-1)=dist(c(i-1),c(i+1));
    g0(i)=g0(i-1)+f(i)-dxy(i);
  end
  lim=sum(dxy)/(n-1)*2;
  ddxy=dxy(1:end-1)+dxy(2:end)-d2xy;
  icorn=find(ddxy>lim);
  if isempty(icorn)
    break;
  end
  [dum,i]=sort(-ddxy(icorn));
  icorn=icorn(i);
  for i=1:length(icorn)
    k=icorn(i)+1;
    ddxy2=dist(c(1:end-1),c(k))+dist(c(2:end),c(k))-dxy;
    ddxy2(k)=NaN;
    ddxy2(k-1)=NaN;
    [mn,l]=min(ddxy2);
    if mn<ddxy(k-1)
      opt=0;
      if k<l
        if ddxy(k-1)>=a(c(k),4)
          opt=1;    % you gain more then what you loose, so no problem
        elseif g0(k+1:l)>=a(c(k),4)-ddxy(k-1)
          opt=1;
        end
        if opt
          c0=c;
          c=[c(1:k-1);c(k+1:l);c(k);c(l+1:n)];
        end
      else
        if l==1
          g01=a(1,4)-dist(c(k));
        else
          g01=g0(l-1)+a(c(l),4)-dist(c(l),c(k));
        end
        if g01>=0
          dg0=dxy(l)-dist(c(k),c(l))-dist(c(k),c(l+1))+a(c(k),4);
          if all(g0(l:k-1)+dg0>=0)
            opt=1;
          end
        end
        if opt
          c0=c;
          c=[c(1:l);c(k);c(l+1:k-1);c(k+1:n)];
        end
      end
      if opt
        break;    % currently just one correction/loop
        nopt=nopt+1;
      end
    end
  end
end