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);


                