isedm(D,tolerance_in,verbose,dim,V);

%Is real D a Euclidean Distance Matrix. -Jon Dattorro, http://convexoptimization.com
%
%[Dclosest,X,isisnot,r] = isedm(D,tolerance,verbose,dimension,V)
%
%Returns: closest EDM in Schoenberg sense (default output), 
%         a generating list X, 
%         string 'is' or 'isnot' EDM,
%         actual affine dimension r of EDM output.
%Input: candidate matrix D, 
%       optional absolute numerical tolerance for EDM determination, 
%       optional verbosity 'on' or 'off',
%       optional desired affine dim of generating list X output,
%       optional choice of 'Vn' auxiliary matrix (default) or 'V'.

function [Dclosest,X,isisnot,r] = isedm(D,tolerance_in,verbose,dim,V);

isisnot = 'is'; 
N = length(D); 

if nargin < 2 | isempty(tolerance_in)
   tolerance_in = eps;
end
tolerance = max(tolerance_in, eps*N*norm(D));
if nargin < 3 | isempty(verbose)
   verbose = 'on';
end
if nargin < 5 | isempty(V)
   use = 'Vn';
else
   use = 'V';
end

%is empty
if N < 1 
   if strcmp(verbose,'on'), disp('Input D is empty.'), end
   X = [ ];
   Dclosest = [ ];
   isisnot = 'isnot';
   r = [ ];
   return 
end 
%is square
if size(D,1) ~= size(D,2)
   if strcmp(verbose,'on'), disp('An EDM must be square.'), end
   X = [ ];
   Dclosest = [ ];
   isisnot = 'isnot';
   r = [ ];
   return
end
%is real
if ~isreal(D)
   if strcmp(verbose,'on'), disp('Because an EDM is real,'), end
   isisnot = 'isnot';
   D = real(D);
end
%is nonnegative
if sum(sum(chop(D,tolerance) < 0))
   isisnot = 'isnot';
   if strcmp(verbose,'on'), disp('Because an EDM is nonnegative,'),end
end
%is symmetric
if sum(sum(abs(chop((D - D')/2,tolerance)) > 0))
   isisnot = 'isnot';
   if strcmp(verbose,'on'), disp('Because an EDM is symmetric,'), end
   D = (D + D')/2;  %only required condition
end
%has zero diagonal
if sum(abs(diag(chop(D,tolerance))) > 0)
   isisnot = 'isnot';
   if strcmp(verbose,'on') 
      disp('Because an EDM has zero main diagonal,')
   end
end
%is EDM
if strcmp(use,'Vn')
   VDV = -Vn(N)'*D*Vn(N);
else
   VDV = -Vm(N)'*D*Vm(N);
end
[Evecs Evals] = signeig(VDV);
if ~isempty(find(chop(diag(Evals),max(tolerance_in,eps*N*normest(VDV))) < 0))
   isisnot = 'isnot';
   if strcmp(verbose,'on'), disp('Because -VDV < 0,'), end
end
if strcmp(verbose,'on')
   if strcmp(isisnot,'isnot') 
      disp('matrix input is not EDM.')
   elseif tolerance_in == eps
      disp('Matrix input is EDM to machine precision.')
   else
      disp('Matrix input is EDM to specified tolerance.')
   end
end

%find generating list 
r = max(find(chop(diag(Evals),max(tolerance_in,eps*N*normest(VDV))) > 0));
if isempty(r)
   r = 0;
end
if nargin < 4 | isempty(dim)
   dim = r;
else
   dim = round(dim);
end
t = r;
r = min(r,dim);
if r == 0 
   X = zeros(1,N);
else
   if strcmp(use,'Vn')
      X = [zeros(r,1) diag(sqrt(diag(Evals(1:r,1:r))))*Evecs(:,1:r)'];
   else
      X = [diag(sqrt(diag(Evals(1:r,1:r))))*Evecs(:,1:r)']/sqrt(2);
   end
end
if strcmp(isisnot,'isnot') | dim < t 
   Dclosest = Dx(X);
else 
   Dclosest = D; 
end
%----------------------------------

%zero entries below specified absolute tolerance threshold
function Y = chop(A,tolerance)
R = real(A);
I = imag(A);
if nargin == 1
   tolerance = max(size(A))*norm(A)*eps;
end
idR = find(abs(R) < tolerance);
idI = find(abs(I) < tolerance);
R(idR) = 0;
I(idI) = 0;
Y = R + i*I;

%-----------------------------------
function y = Vn(N)
y = [-ones(1,N-1); 
      eye(N-1)]/sqrt(2);

%-----------------------------------
%returns EDM V matrix
function V = Vm(n)
V = [eye(n)-ones(n,n)/n];

%-----------------------------------
%Make EDM from point list
function D = Dx(X)
[n,N] = size(X);
one = ones(N,1);
del = diag(X'*X);
D = del*one' + one*del' - 2*X'*X;