function [coord,sort_coord,permut,prox_permut,vaf,...
    alpha_multiplier,alpha_vaf,alpha_coord] = skew_symmetric_scaling(prox)

% SKEW_SYMMETRIC_SCALING performs a closed-form least squares scaling of a
% skew-symmetric proximity matrix.
%
% syntax: [coord,sort_coord,permut,prox_permut,vaf,...
%   alpha_multiplier,alpha_vaf,alpha_coord] = skew_symmetric_scaling(prox)
%
% COORD contains the least-squares coordinates which are sorted from
% smallest to largest in SORT_COORD; VAF is the variance they account for;
% PERMUT is the object permutation corresponding to SORT_COORD with
% PROX_PERMUT the reordered skew-symmetric proximity matrix.  For
% equally-spaced coordinates, ALPHA_MULTIPLIER defines the multiplicative
% constant on the integer-valued coordinates; a collection of
% equally-spaced coordinates is given by ALPHA_COORD with ALPHA_VAF the
% variance they account for.
%
% PROX is the inputted skew-symmetrix proximity matrix with 0's along the
% diagonal

%checking that imputted matrix is square

if size(prox, 1) ~= size(prox, 2)

    error('Proximity matrix needs to be square');

end

%checking that diagonal entries are all 0

if sum(diag(prox) ~= zeros(size(prox,1), 1)) ~= 0
   
    error('Diagonal entries of proximity matrix need to be 0');

end

%checking that proximity matrix is skew-symmetric

if squareform(triu(prox)) ~= -squareform(tril(prox))
   
    error('Proximity matrix must be skew-symmetric');

end


%determining the size of the proximity matrix

n = size(prox, 1);


%least-squares optimal set of coordinates

coord = (1/n) * sum(prox);


%sorting the coordinates from smallest to largest

sort_coord = sort(coord);


%optimal permuation vector (will contain optimal ordering)

permut = zeros(n, 1);


%determining the optimal ordering

for i = 1:n
    
    permut(i) = find(sort_coord(i) == coord);
    
end


%proximity matrix, now ordered optimally

prox_permut = prox(permut, permut);


%determining the variance accounted for (VAF) of the optimal solution w/o
%scaling constant

vaf = (sum(sum(prox).^2)/n) / (sum(sum(triu(prox).^2)));


%determining scaling constant alpha

suma = 0.0;
sumb = 0.0;

for i = 1:(n-1)
    
    suma = suma + sum(prox_permut(i,(i+1):n) .* (((i+1):n)-i));
    sumb = sumb + sum((((i+1):n)-i).^2);
    
end

alpha_multiplier = suma/sumb;


%determining the variance accounted for (VAF) of the optimal solution with
%scaling constant

    %denominator of VAF
denom = sum(sum(triu(prox).^2));

    %numerator of VAF
numer = 0.0;

for i = 1:(n-1)
    
    numer = numer + sum((prox_permut(i,(i+1):n) - alpha_multiplier * (((i+1):n)-i)).^2);

end

alpha_vaf = 1 - (numer / denom);
    

%determining the alpha-scaled coordinates (will be equally spaced)

alpha_coord = alpha_multiplier*((1:n) - (n*(n+1)/2) / n)';