function [fit,weights,vaf] = orderpartitionfit(prox,lincon,membership)
% ORDERPARTITIONFIT provides a least-squares approximation to a proximity
% matrix based on a given collection of partitions with ordered classes.
%
% syntax: [fit,weights,vaf] = orderpartitionfit(prox,lincon,membership)
%
% PROX is the n x n input proximity matrix (with a zero main diagonal
% and a dissimilarity interpretation); LINCON is the given constraining
% linear order (a permutation of the integers from 1 to n).
% MEMBERSHIP is the m x n matrix indicating cluster membership, where
% each row corresponds to a specific ordered partition (there are
% m partitions in general);
% the columns are in the identity permutation input order used for PROX.
% FIT is an n x n matrix fitted to PROX (through least-squares) constructed
% from the nonnegative weights given in the m x 1 WEIGHTS vectors
% corresponding to each of the ordered partitions. VAF is the variance-
% accounted-for in the proximity matrix PROX by the fitted matrix FIT.
n = size(prox,1);
m = size(membership,1);
proxpermut = prox(lincon,lincon);
memberpermut = membership(1:m,lincon);
d = zeros(m,n*(n-1)/2);
for k = 1:m
index = 0;
for i = 1:(n-1)
for j = (i+1):n
index = index + 1;
if(memberpermut(k,i) == memberpermut(k,j))
d(k,index) = 0;
else
d(k,index) = 1;
end
end
end
end
index = 0;
for i = 1:(n-1)
for j = (i+1):n
index = index + 1;
proxvector(index,1) = proxpermut(i,j);
end
end
data = proxvector;
covariance = eye(n*(n-1)/2);
constraint_array = d;
constraint_constant = zeros(m,1);
equality_flag = zeros(m,1);
[solution,kuhn_tucker,interations,end_condition] = ...
dykstra(data,covariance,constraint_array,constraint_constant,equality_flag);
fitvector = data - solution;
weights = (.5)*kuhn_tucker;
fit = zeros(n,n);
meanprox = zeros(n,n);
meanproximity = (sum(sum(prox)))/(n*(n-1));
index = 0;
for i = 1:(n-1)
for j = (i+1):n
index = index +1;
meanprox(i,j) = meanproximity;
fit(i,j) = fitvector(index);
meanprox(j,i) = meanproximity;
fit(j,i) = fit(i,j);
end
end
num = sum(sum(((prox(lincon,lincon) - fit).^2)));
den = sum(sum(((prox(lincon,lincon) - meanprox).^2)));
vaf = 1 - (num/den);