Solver V5

function [xy,hist,hist2] = solver(kd,bx)

% set parameters
alpha  = 0.95;
beta   = 0.9;
N      = size(kd,1);
M      = size(find(kd>0),1)/2;
P      = 1;
hist2  = [];
rand('state',sum(100*clock));

% modifications
pd        = sqrt((bx(2)-bx(1))^2+(bx(4)-bx(3))^2);
kd(kd>pd) = pd;
meankd    = sum(sum(kd.*(kd~=-1)))/M;
eta       = meankd^1.8;

% initialize
if N == 1
   xy = [bx(1),bx(3)];
   return
end
z     = bx(1)+rand(N,1)*(bx(2)-bx(1))+j*(bx(3)+rand(N,1)*(bx(4)-bx(3)));
Z     = z(:,ones(N,1));
dZ    = Z-transpose(Z);
err   = (abs(dZ)-kd).*(kd~=-1);
ztmp  = z;
etmp  = sum(abs(err(:)));
hist  = etmp;
oldgr = zeros(N,1);
e = etmp;

% simple search
for m = 1:120
   
   % compute gradient
   gr   = beta*sum(err.*dZ,2)+(1-beta)*oldgr;
   norm = sqrt(sum(abs(gr).^2));
   if norm == 0
      xy = [real(z),imag(z)];
      return
   end
   gr   = gr/norm;
   var  = norm/M/M;
   
   hist2 = [hist2,norm/M/M];
   
   % update state 
   zz = z-rndrayl(eta,N,1).*(rndnorm(real(gr),var,N,1)+j*rndnorm(imag(gr),var,N,1));
   zr = min(max(real(zz),bx(1)),bx(2));
   zi = min(max(imag(zz),bx(3)),bx(4));
   zz = zr+j*zi;
   
   % compute error
   Z   = zz(:,ones(N,1));
   dZ  = Z-transpose(Z);
   err = (abs(dZ)-kd).*(kd~=-1);
   e   = sum(abs(err(:)));
   
   % remember best state
   if (e > hist(end)) & (hist(end) < etmp)
      ztmp = z;
      etmp = hist(end);
   end
   z = zz;
   
   if m < 80
      eta = alpha*eta;
   end
   
   oldgr = gr;
   hist  = [hist;e];
end

if e > etmp
   z = ztmp;
   e = etmp;
else
   etmp = e;
end

% start local search
%kd1  = kd(:);
%kd2  = kd>0;
%kd2  = kd2(:);
%Z    = z(:,ones(N,1));
%dZ   = Z-transpose(Z);
%gr   = sum((abs(dZ)-kd).*(kd~=-1).*dZ./(abs(dZ)+1e-12),2);
%%gr   = real(gr).*~(real(gr) > 0 & real(z) == bx(1)).*~(real(gr) < 0 & real(z) == bx(2))+...	% gradient correction
%%   j*imag(gr).*~(imag(gr) > 0 & imag(z) == bx(3)).*~(imag(gr) < 0 & imag(z) == bx(4));
%norm = sqrt(sum(abs(gr).^2));
%if norm == 0
%   xy = [real(z),imag(z)];
%   return
%end
%gr   = gr/norm;
%var  = norm/M/M*sqrt(sum(abs(err),2));
%var  = var(:,ones(P-1,1));

% process local search
%z     = [z,z(:,ones(P-1,1))-rndrayl(pd/7.5,N,P-1).*(rndnorm(real(gr(:,ones(P-1,1))),var,N,P-1)+...
%      j*rndnorm(imag(gr(:,ones(P-1,1))),var,N,P-1))];
%zr    = min(max(real(z),bx(1)),bx(2));
%zi    = min(max(imag(z),bx(3)),bx(4));
%z     = zr+j*zi;
%zz1   = kron(ones(N,1),z)-kron(z,ones(N,1));
%kron(eye(N),ones(1,N));
%zz2   = sum((abs(abs(zz1)-kd1(:,ones(P,1)))).*kd2(:,ones(P,1)));
%[m,n] = min(zz2);
%err   = reshape((abs(zz1(:,n))-kd1).*kd2,N,N);
%z     = z(:,n);

xy    = [real(z),imag(z)];

% * * * Functions * * * %

function x = rndnorm(mu,sigma,m,n)

x = randn(m,n).*sigma+mu;

function x = rndrayl(b,m,n)

x = (rndnorm(0,b,m,n).^2+rndnorm(0,b,m,n).^2);