randfsearch

function xy = solver(varargin)
if nargin==2
  kd=varargin{1};
  bx=varargin{2};
  npts = size(kd,1);
  dxy=bx([2 4])-bx([1 3]);
  if (dxy(1)<=eps & dxy(2)<=eps) | all(kd(:)<=0)
    % 0 dimensional
    xy=repmat(bx([1 3]),npts,1);
    return
  else
    [I,J,v] = find(triu(kd)>0);
    v2=kd((J-1)*npts+I).^2;
    k = 100;
    x = rand(npts,k)*(bx(2)-bx(1))+bx(1); 
    y = rand(npts,k)*(bx(4)-bx(3))+bx(3);
    [tv,id] = min(sum(abs((x(I,:)-x(J,:)).^2 ...
                +(y(I,:)-y(J,:)).^2-repmat(v2,1,k))));
    x0=x(:,id); y0=y(:,id);
    xy=[x0 y0];
  end
  [I,J,v] = find(triu(kd)>0);
  v2=kd((J-1)*npts+I).^2;
  ops.MaxIter=50;
  ops.Display='off';
  xy=fminsearch('solver',xy,ops,I,J,v2,bx,npts);
  xy=reshape(xy,[npts 2]);
else
  I=varargin{2};
  J=varargin{3};
  v2=varargin{4};
  bx=varargin{5};
  npts=varargin{6};
  xy0=varargin{1};
  if any(xy0(:,1)<bx(1) | xy0(:,1)>bx(2)) | ...
      any(xy0(:,2)<bx(3) | xy0(:,2)>bx(4))
    xy=inf;
  else
    xy = sum(abs((xy0(I,1)-xy0(J,1)).^2+(xy0(I,2)-xy0(J,2)).^2-v2));
  end
end