Stuning10

function xy = solver(kd,bx)
npts = length(kd);
S=(kd-eye(npts))>=0;
npts0=npts;
nzS = find(sum(S));
S=S(nzS,nzS);
sI=find(S);
kd=kd(nzS,nzS);
npts = length(nzS);
if npts<2
	xy=zeros(npts,2);
elseif npts==2
	xy=[0;(bx(2)+1i*bx(4))*kd(1,2)/sqrt(bx(2)^2+bx(4)^2)];
else
nP=max(sum(S));
mx=max(kd(:))*(1+0.002*nP/npts);
dx=min(bx(2),mx);
dy=min(bx(4),mx);
D=dx+1i*dy;
pm=D/2;
xy=zeros(1,npts)+pm;
dr=abs(D);
kd(find(kd>dr))=dr;

minL = min(kd(find(triu(S,1))));
sCut = 0.1;
d=max([max(dx,dy)/10, 6*sCut/nP, minL/2]);
x=0:d:dx;
y=0:d:dy;
[xx,yy] = ndgrid(x,1i*y);
XX=xx(:)+yy(:);

nX=length(XX);
xOnes=ones(nX,1);
xZeros=zeros(nX,1);
iones=ones(1,npts);
sCut = 1;

olds = 1.0e30;
[dum,p]=sort(-max(kd));
S=S(p,p);
sI=find(S);
dX=S;
kd=kd(p,p);
[v,w]=sort(p);
for iter=1:11
  for index = 1:npts
    I = find(S(:,index)); 
    if ~isempty(I)
      st = sum(abs(abs(XX(:,ones(size(I)))-xy(xOnes,I))-kd(xZeros+index,I)),2);
      [null,minloc] = min(st);
      xy(index) = XX(minloc);
    end
  end
  [X1,X2]=meshgrid(xy); 
  dX(sI)=X1(sI)-X2(sI);
  s = sum(abs(abs(dX(sI))-kd(sI)));

  if((abs(s-olds)/(s+1)) < 0.002) | s < sCut
    break;
  end
  olds = s;
end
xy=xy.';
xy=xy(w,:);
kd=kd(w,w);
S=S(w,w);
sI=find(S);
dX=S;

alpha = 0.1;
ngrad = 60;

obj = 1.0e20*ones(ngrad,1);
xybig = zeros(npts,ngrad);

gradient = zeros(npts,1);
[X1,X2]=meshgrid(xy); 
dM=S;sM=S;
dM(sI)=X1(sI)-X2(sI);
sM(sI)=abs(dM(sI))-kd(sI);
objnew = sum(abs(sM(sI)));
gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM));
ograd=gradient;
obj(1) = objnew;
xybig(:,1) = xy;


for ij=2:ngrad  
  xy = xy - alpha  * gradient;
  xy = min(max(real(xy),0),dx) + j * min(max(imag(xy),0),dy);
  
  [X1,X2]=meshgrid(xy); 
  dM=S;sM=S;
  dM(sI)=X1(sI)-X2(sI);
  sM(sI)=abs(dM(sI))-kd(sI);
  s = sum(abs(sM(sI)));
  gradient(:) = sum(dM./(abs(dM)+eps).*sign(sM));
  
  obj(ij) = s;
  xybig(:,ij) = xy;
  if( obj(ij) > obj(ij-1) )
    xy = xybig(:,ij-1);
    alpha = alpha / 10;
    gradient=ograd;
  else
    alpha = alpha * 2;
  end
  [bestobj,index] = min(obj);
  if bestobj < sCut
    break
  end
  ograd = gradient;
end
[bestobj,index] = min(obj);
xy = xybig(:,index);
end
xy0 = zeros(npts0,1);
xy0(nzS)=xy;
xy=[real(xy0) imag(xy0)];