turbo2

function xy = solver(kd,bx)
npts = length(kd);
if ~any(kd(:)>0),
xy = zeros(npts,2);
return
end
dx=bx(2);
dy=bx(4);
D=dx+1i*dy;
pm=D/2;
xy=zeros(1,npts)+pm;
npf=.1;
if npts<10, npf=.01; end
d=npf*max(dx,dy);
x=0:d:dx;
y=0:d:dy;
[xx,yy] = ndgrid(x,1i*y);
[nx,ny]=size(xx);
onesx=ones(1,nx);
onesy=ones(1,ny);
XX=xx(:)+yy(:);

dr=abs(D);
kd(find(kd>dr))=dr;

S=(kd-eye(npts))>=0;
nX=length(XX);
xOnes=ones(nX,1);
xZeros=zeros(nX,1);
start=1;
iones=ones(1,npts);
sCut = 2;

olds = 1.0e30;
iterbigmax = 1;

for iterbig = 1:iterbigmax
    [dum,p]=sort(-max(kd));
    S=S(p,p);
    kd=kd(p,p);
    [v,w]=sort(p);
    
    for iter=1:20 
        %	for index = 1:npts
        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
        if iter>2,
         [s,M]=findstrain(S,xy,kd);
         csum(iter)=s;
         if((abs(s-olds)/(s+1)) < 0.002) | s < sCut
            break;
         end
         olds = s;
        end
    end
    xy=xy.';
    xy=xy(w,:);
    kd=kd(w,w);
    S=S(w,w);
       
    if( iterbig < iterbigmax )
        xy=xy.';
    end
end

alpha = 0.1;
ngrad = 30;

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

gradient = zeros(npts,1);
[objnew,sM,dM]=findstrain(S,xy,kd);
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);
   
   [s,sM,dM]=findstrain(S,xy,kd);
   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);
xy=[real(xy) imag(xy)];

function [s,strainMatrix,dX]=findstrain(S,xy,kd)
sI=find(S);
strainMatrix=S;
dX=S;
[X1,X2]=meshgrid(xy); 
dX(sI)=X1(sI)-X2(sI);
strainMatrix(sI) = abs(dX(sI))-kd(sI);
s = sum(abs(strainMatrix(sI)));