function [f t X] = pc(dist,X0)

% Given a distance matrix (not necessarily Euclidean) D and a initial set of
% seeding point set X0, this function finds a set of points X such that
% Euclidean distance among the new points X is as close as to D.

% Problem: min_X sum_{i,j} ||D_ij - d(x_i,x_j)|| where d(x_i,x_j) computes the
% Euclidean distance between points x_i and x_j.

% Test usage:
% tic; n=10;d=5;
% Z = rand(d,n);
% X0 =rand(d,n);
% D = distAllPairsL2(Z);
% [f,t,X] = pc(D,X0);
% plot(t,f);

% Note that in the case of l_1 norm, mean step(line 43) can be replace by
% computing the median, Weisfield's algorithm and "computeIntersections"
% (line 42) function by the funciton that computes the intersections of the sphere
% with the line using the input space's geometry

% Author: Arvind Agarwal (http://www.cs.utah.edu/~arvind)
% Date Created: Feb 06, 2010


[d ,n] = size(X0);

% Initialize the palceCenter Algorithm 
X=X0;

% compute the residue
Dn = distAllPairsL2(X);
c = sum(sum((dist-Dn).^2));
fprintf(1,'Initial cost : %g\n',c);

for iter = 1:100 %increase the numiter if you want
    f(iter) = c;
    t(iter) = toc;
    for ii = 1:n
        for it=1:floor(sqrt(iter))
            Z = computeIntersections(X,ii,dist);
            p = mean(Z,2);
            X = [X(:,1:ii-1) p X(:,ii+1:n)];
        end
    end
    Dn = distAllPairsL2(X);
    c = sum(sum((dist-Dn).^2));
    fprintf(1,'Cost at iter %d : %g \t %g\n',iter,c,toc);
end

%----- compute the intersecting points ---------------
function Z = computeIntersections(X,ii,dist)
[d ,n] = size(X);
V = X-repmat(X(:,ii),1,n);
D = sqrt(sum(V.*V,1))';
Q = dist(ii,:)'./D;
Z = X - V.*repmat(Q',d,1);
Z(:,ii)=[];

%--------------- some useful function not required for the algorithm------------%
function K=distAllPairsL2(X)
ps = X'*X; 
n = size(ps,1);
normx = repmat(sum(X'.^2,2),1,n);
K = real(sqrt(-2*ps + normx + normx'));
% making diagonal distance zero
K = K - diag(diag(K));