function [coordinates,vaf,vectors,roots] = torgerson(prox,numdim)

n = size(prox,1);

dstar = zeros(n,n);

a = (1/n)*(sum(prox.^2,2));
b = (1/n)*(sum(prox.^2,1));
c = (1/(n*n))*(sum(sum(prox.^2,2)));

for i = 1:n
    for j = 1:n
        
        dstar(i,j) = (-1/2)*((prox(i,j)^2) - a(i) - b(j) + c);
    end
end

[vectors,roots] = eig(dstar);

[vals,inds] = sort(diag(roots));

inds = flipud(inds);
vals = flipud(vals);

vectors = vectors(:,inds);
roots = diag(vals);

coordinates = vectors(:,1:numdim)*diag(sqrt(vals(1:numdim)));

distances = pdist(coordinates);

proximities = squareform(prox);

corrmatrix = corrcoef(distances',proximities');

vaf = corrmatrix(1,2)^2;

